Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief GW using RI-RS Approximation for molecules
10 : !> \par History
11 : !> 04.2026 created [Ritaj Tyagi]
12 : ! **************************************************************************************************
13 : MODULE gw_ri_rs_large_cell_Gamma
14 : USE atomic_kind_types, ONLY: atomic_kind_type
15 : USE basis_set_types, ONLY: gto_basis_set_type
16 : USE cell_types, ONLY: cell_type,&
17 : get_cell,&
18 : pbc
19 : USE cp_dbcsr_api, ONLY: &
20 : dbcsr_add, dbcsr_binary_read, dbcsr_binary_write, dbcsr_create, dbcsr_distribution_get, &
21 : dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
22 : dbcsr_finalize, dbcsr_get_info, dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_scale, &
23 : dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
24 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
25 : copy_fm_to_dbcsr,&
26 : dbcsr_deallocate_matrix_set,&
27 : max_elements_per_block
28 : USE cp_fm_types, ONLY: cp_fm_create,&
29 : cp_fm_release,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: cp_p_file,&
34 : cp_print_key_should_output
35 : USE gw_ri_rs_compute_Z_lP_utils, ONLY: build_gram_jacobi_blas,&
36 : scale_rows_by_diag,&
37 : store_Z_lP_columns
38 : USE gw_ri_rs_grid_setup_main, ONLY: setup_ri_rs_grid
39 : USE gw_ri_rs_utils, ONLY: precompute_ri_rs_radii
40 : USE gw_tensor_large_cell_Gamma, ONLY: G_occ_vir,&
41 : compute_QP_energies,&
42 : delete_unnecessary_files,&
43 : fill_fm_Sigma_c_Gamma_time,&
44 : get_W_MIC
45 : USE gw_utils, ONLY: de_init_bs_env
46 : USE gw_utils_compute_integrals, ONLY: build_3c_integral_block_ctx,&
47 : gw_3c_ctx_create,&
48 : gw_3c_ctx_release,&
49 : gw_3c_ctx_type,&
50 : gw_3c_ws_create,&
51 : gw_3c_ws_release,&
52 : gw_3c_ws_type
53 : USE gw_utils_dbcsr, ONLY: dbcsr_contract_ABA,&
54 : hadamard_product
55 : USE gw_utils_fm, ONLY: fm_contract_ABA
56 : USE input_section_types, ONLY: section_vals_type
57 : USE kinds, ONLY: dp
58 : USE machine, ONLY: m_walltime
59 : USE message_passing, ONLY: mp_para_env_type
60 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
61 : USE orbital_pointers, ONLY: indco,&
62 : ncoset
63 : USE particle_types, ONLY: particle_type
64 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
65 : USE qs_environment_types, ONLY: get_qs_env,&
66 : qs_environment_type
67 : #include "./base/base_uses.f90"
68 :
69 : IMPLICIT NONE
70 :
71 : PRIVATE
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_ri_rs_large_cell_Gamma'
74 :
75 : PUBLIC :: gw_calc_ri_rs_large_cell_Gamma, &
76 : release_dbcsr_topology_and_matrices, &
77 : setup_square_topology
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief GW calculation using RI-RS formalism for molecules
83 : !> \param qs_env ...
84 : !> \param bs_env Band-structure environment containing GW parameters.
85 : ! **************************************************************************************************
86 :
87 0 : SUBROUTINE gw_calc_ri_rs_large_cell_Gamma(qs_env, bs_env)
88 :
89 : TYPE(qs_environment_type), POINTER :: qs_env
90 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
91 :
92 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_ri_rs_large_cell_Gamma'
93 :
94 : INTEGER :: handle
95 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_W_time
96 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
97 :
98 0 : CALL timeset(routineN, handle)
99 :
100 : !!========================================================================
101 : !! 0. Precompute AO and RI Radii
102 : !! Per-atom cutoff radii from the most diffuse Gaussian primitives in
103 : !! the AO and RI auxiliary basis sets. Stored in bs_env%ri_rs%
104 : !! radius_ao_per_atom and radius_ri_per_atom, used for sphere-cutoff
105 : !! and phi_local screening.
106 : !!========================================================================
107 0 : CALL precompute_ri_rs_radii(bs_env)
108 :
109 : !!========================================================================
110 : !! 1. Grid Generation for RI-RS
111 : !! (Modified Lebedev grids from Ivan Duchemin and Xavier Blase)
112 : !! Generate flattened 1D array of grid points for RI-RS.
113 : !! Equation: r_g(k) = R_A + r_g(A)
114 : !!========================================================================
115 0 : CALL setup_ri_rs_grid(bs_env, bs_env%ri_rs%grid_points)
116 :
117 : !!========================================================================
118 : !! 2. Atomic Basis Evaluation
119 : !! Compute values of spherical atomic basis functions at grid points.
120 : !! Expression: ϕ_μl = ϕ_μ(r_l) (mat_phi_mu_l)
121 : !!========================================================================
122 : CALL atomic_basis_at_grid_point(bs_env, bs_env%ri_rs%grid_points, &
123 0 : bs_env%ri_rs%mat_phi_mu_l)
124 :
125 : !!========================================================================
126 : !! 3. Compute RI-RS Coefficients (Z_lp)
127 : !! Solve the regularized system for each atom P, where the grid domain
128 : !! is restricted to r_l within a cutoff distance of atom P:
129 : !! a. D_ll' = [ Σ_μ ϕ_μ(r_l) ϕ_μ(r_l') ]^2 (Equation 13)
130 : !! b. D_lP = Σ_{μν} ϕ_μ(r_l) ϕ_ν(r_l) (μν|P) (Equation 15)
131 : !! c. Conditioning:
132 : !! Dvec_l = 1 / sqrt(D_ll) (Diagonal scaling vector)
133 : !! D'_ll' = Dvec_l * D_ll' * Dvec_l' + λδ_ll'
134 : !! D'_lP = Dvec_l * D_lP
135 : !! d. Solve: Σ_l' D'_ll' * Z'_l'P = D'_lP (Equation 14)
136 : !! e. Rescale: Z_lP = Z'_lP * Dvec_l (Z_lP stored in mat_Z_lP)
137 : !!========================================================================
138 : CALL compute_Z_lP(qs_env, bs_env, bs_env%ri_rs%grid_points, &
139 0 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
140 0 : bs_env%ri_rs%grid_built = .TRUE.
141 :
142 : !!========================================================================
143 : !! 4. Compute Independent-Particle Polarizability (χ)
144 : !! G^occ_µλ(i|τ|) = sum_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
145 : !! G^vir_µλ(i|τ|) = sum_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
146 : !! G^occ_ll'(i|τ|) = sum_µν ϕ_µ(r_l) G^occ_µν ϕ_ν(r_l')
147 : !! G^vir_ll'(i|τ|) = sum_µν ϕ_µ(r_l) G^vir_µν ϕ_ν(r_l')
148 : !! χ_ll'(iτ) = G^occ_ll'(i|τ|) * G^vir_ll'(i|τ|)
149 : !! χ_PQ(iτ) = sum_ll' Z_lP χ_ll'(iτ) Z_l'Q
150 : !!========================================================================
151 : CALL get_mat_chi_Gamma_tau(bs_env, bs_env%mat_chi_Gamma_tau, &
152 0 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
153 :
154 : !!========================================================================
155 : !! 5. Compute Screened Interaction (W^MIC)
156 : !! χ_PQ(iτ) -> χ_PQ(iω) -> ε_PQ(iω) -> W_PQ(iω) -> W^MIC_PQ(iτ)
157 : !!========================================================================
158 0 : CALL get_W_MIC(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_time)
159 :
160 : !!========================================================================
161 : !! 6. Compute Exact Exchange Self-Energy (Σ^x)
162 : !! D_µν = sum_n^occ C_µn C_νn
163 : !! D_ll' = sum_µν ϕ_µ(r_l) D_µν ϕ_ν(r_l')
164 : !! V^trunc_ll' = sum_PQ Z_lP V^trunc_PQ Z_l'Q
165 : !! Σ^x_ll' = D_ll' * V^trunc_ll'
166 : !! Σ^x_λσ(k=0) = -sum_ll' ϕ_λ(r_l) Σ^x_ll' ϕ_σ(r_l')
167 : !!========================================================================
168 : CALL compute_Sigma_x(bs_env, qs_env, bs_env%ri_rs%mat_phi_mu_l, &
169 0 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_x_Gamma)
170 :
171 : !!========================================================================
172 : !! 7. Compute Correlation Self-Energy (Σ^c)
173 : !! W^MIC_ll'(iτ) = sum_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
174 : !! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) * W^MIC_ll'(iτ), for τ < 0
175 : !! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) * W^MIC_ll'(iτ), for τ > 0
176 : !! Σ^c_λσ(iτ) = sum_ll' ϕ_λ(r_l) Σ^c_ll'(iτ) ϕ_σ(r_l')
177 : !!========================================================================
178 : CALL compute_Sigma_c(bs_env, fm_W_time, bs_env%ri_rs%mat_phi_mu_l, &
179 0 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_c_Gamma_time)
180 :
181 : !!========================================================================
182 : !! 8. Compute Quasiparticle Energies
183 : !! Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
184 : !! ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ) + Σ^x_nn - v^xc_nn
185 : !!========================================================================
186 0 : CALL compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
187 :
188 0 : CALL de_init_bs_env(qs_env, bs_env)
189 :
190 0 : CALL timestop(handle)
191 :
192 0 : END SUBROUTINE gw_calc_ri_rs_large_cell_Gamma
193 :
194 : ! **************************************************************************************************
195 : !> \brief Evaluates atomic basis functions on a real-space grid and builds a sparse DBCSR matrix.
196 : !> \param bs_env ...
197 : !> \param ri_rs_grid_points ...
198 : !> \param mat_phi_mu_l ...
199 : ! **************************************************************************************************
200 :
201 0 : SUBROUTINE atomic_basis_at_grid_point(bs_env, ri_rs_grid_points, mat_phi_mu_l)
202 :
203 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
204 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
205 : TYPE(dbcsr_type), INTENT(OUT) :: mat_phi_mu_l
206 :
207 : CHARACTER(LEN=*), PARAMETER :: routineN = 'atomic_basis_at_grid_point'
208 :
209 : INTEGER :: c_size, chunk_size, handle, i, i_blk, &
210 : iatom, natom, npcol, nprow, &
211 : num_grid_chunks, r_end, r_start
212 0 : INTEGER, DIMENSION(:), POINTER :: col_dist, r_blk_sizes, row_dist, sizes_AO
213 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atom_col_buffer
214 : TYPE(cell_type), POINTER :: cell
215 : TYPE(dbcsr_distribution_type) :: dbcsr_dist_ks, dist
216 : TYPE(mp_para_env_type), POINTER :: para_env
217 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
218 :
219 0 : CALL timeset(routineN, handle)
220 :
221 : ! Setup Grid Blocking
222 0 : chunk_size = max_elements_per_block
223 :
224 0 : natom = bs_env%n_atom
225 0 : cell => bs_env%ri_rs%cell
226 0 : para_env => bs_env%para_env
227 0 : particle_set => bs_env%ri_rs%particle_set
228 0 : sizes_AO => bs_env%sizes_AO
229 0 : CPASSERT(ASSOCIATED(cell))
230 0 : CPASSERT(ASSOCIATED(para_env))
231 0 : CPASSERT(ASSOCIATED(particle_set))
232 0 : CPASSERT(SIZE(sizes_AO) == natom)
233 0 : CPASSERT(SIZE(ri_rs_grid_points, 2) == bs_env%ri_rs%n_grid_points)
234 :
235 : ! =========================================================================
236 : ! 1. SETUP DBCSR MATRIX TOPOLOGY
237 : ! =========================================================================
238 :
239 : ! B. Define Row Block Sizes (Grid chunks of max size 256)
240 0 : num_grid_chunks = CEILING(REAL(bs_env%ri_rs%n_grid_points, KIND=dp)/REAL(chunk_size, KIND=dp))
241 0 : ALLOCATE (r_blk_sizes(num_grid_chunks))
242 0 : r_blk_sizes = chunk_size
243 0 : IF (MOD(bs_env%ri_rs%n_grid_points, chunk_size) /= 0) THEN
244 0 : r_blk_sizes(num_grid_chunks) = MOD(bs_env%ri_rs%n_grid_points, chunk_size)
245 : END IF
246 :
247 : ! C. Fetch CP2K's Default Process Grid Configuration
248 0 : CALL dbcsr_get_info(bs_env%mat_ao_ao%matrix, distribution=dbcsr_dist_ks)
249 0 : CALL dbcsr_distribution_get(dbcsr_dist_ks, nprows=nprow, npcols=npcol)
250 :
251 : ! D. Build Custom Mappings using Round-Robin across the 2D process grid
252 :
253 0 : ALLOCATE (row_dist(num_grid_chunks))
254 0 : DO i = 1, num_grid_chunks
255 0 : row_dist(i) = MOD(i - 1, nprow)
256 : END DO
257 :
258 0 : ALLOCATE (col_dist(natom))
259 0 : DO i = 1, natom
260 0 : col_dist(i) = MOD(i - 1, npcol)
261 : END DO
262 :
263 : ! E. Create the DBCSR Distribution and Initialize the Matrix
264 : CALL dbcsr_distribution_new(dist, template=dbcsr_dist_ks, &
265 0 : row_dist=row_dist, col_dist=col_dist)
266 :
267 : CALL dbcsr_create(mat_phi_mu_l, name="phi_val_sparse", dist=dist, &
268 : matrix_type=dbcsr_type_no_symmetry, &
269 0 : row_blk_size=r_blk_sizes, col_blk_size=sizes_AO)
270 :
271 : ! =========================================================================
272 : ! 2. STREAM DATA DIRECTLY INTO SPARSE MATRIX
273 : ! =========================================================================
274 : ! Iterate over the atoms assigned to this specific MPI rank
275 0 : DO iatom = para_env%mepos + 1, natom, para_env%num_pe
276 :
277 0 : c_size = sizes_AO(iatom)
278 :
279 : ! Allocate a temporary dense buffer just for this specific atom
280 0 : ALLOCATE (atom_col_buffer(bs_env%ri_rs%n_grid_points, c_size))
281 0 : atom_col_buffer = 0.0_dp
282 :
283 : ! Evaluate the basis functions on the grid. Skip grid points outside the spatial
284 : ! extent of the most diffuse AO Gaussian on iatom; beyond that radius the contribution
285 : ! is guaranteed below eps_filter.
286 : CALL fill_phi_for_atom(atom_col_buffer, ri_rs_grid_points, bs_env%ri_rs%n_grid_points, &
287 : iatom, bs_env, &
288 0 : r2_threshold=bs_env%ri_rs%radius_ao_per_atom(iatom)**2)
289 :
290 : ! Slice the dense column into chunks and insert into DBCSR
291 0 : DO i_blk = 1, num_grid_chunks
292 0 : r_start = (i_blk - 1)*chunk_size + 1
293 0 : r_end = MIN(i_blk*chunk_size, bs_env%ri_rs%n_grid_points)
294 :
295 : ! Apply dynamic sparsity filtering: Only store blocks with physical significance
296 0 : IF (MAXVAL(ABS(atom_col_buffer(r_start:r_end, 1:c_size))) > bs_env%eps_filter) THEN
297 : CALL dbcsr_put_block(mat_phi_mu_l, row=i_blk, col=iatom, &
298 0 : block=atom_col_buffer(r_start:r_end, 1:c_size))
299 : END IF
300 : END DO
301 :
302 0 : DEALLOCATE (atom_col_buffer)
303 :
304 : END DO
305 :
306 : ! Finalize triggers internal MPI communication to route blocks to their correct 2D process owners
307 0 : CALL dbcsr_finalize(mat_phi_mu_l)
308 :
309 0 : IF (bs_env%unit_nr > 0) THEN
310 0 : WRITE (bs_env%unit_nr, *) "Done with evaluation of phi"
311 : END IF
312 :
313 : ! -------------------------------------------------------------------------
314 : ! CLEANUP
315 : ! -------------------------------------------------------------------------
316 0 : DEALLOCATE (r_blk_sizes, row_dist, col_dist)
317 0 : CALL dbcsr_distribution_release(dist)
318 :
319 0 : CALL timestop(handle)
320 :
321 0 : END SUBROUTINE atomic_basis_at_grid_point
322 :
323 : ! **************************************************************************************************
324 : !> \brief Compute value of all basis functions for a single atom across all grid points.
325 : !> Sums contributions from periodic images of `iatom` (loop over (ix, iy, iz) cells gated
326 : !> by `cell%perd`). Each per-image squared distance is compared against `r2_threshold`
327 : !> (per-atom AO Gaussian extent²); images beyond that radius contribute below eps_filter
328 : !> and are skipped.
329 : !> \param phi_val ...
330 : !> \param ri_rs_grid ...
331 : !> \param npts ...
332 : !> \param iatom ...
333 : !> \param bs_env ...
334 : !> \param r2_threshold per-image squared-distance threshold; CYCLE if r² > r2_threshold. Pass
335 : !> HUGE(1.0_dp) to disable.
336 : ! **************************************************************************************************
337 :
338 0 : SUBROUTINE fill_phi_for_atom(phi_val, ri_rs_grid, npts, iatom, bs_env, r2_threshold)
339 :
340 : REAL(KIND=dp), INTENT(INOUT) :: phi_val(:, :)
341 : INTEGER, INTENT(IN) :: npts
342 : REAL(KIND=dp), INTENT(IN) :: ri_rs_grid(3, npts)
343 : INTEGER, INTENT(IN) :: iatom
344 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
345 : REAL(KIND=dp), INTENT(IN) :: r2_threshold
346 :
347 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_phi_for_atom'
348 :
349 : INTEGER :: first_sgf, handle, i_pt, ico, iend_co, ikind, ipgf, iset, isgf, ishell, &
350 : istart_co, ix, ix_max, ix_min, iy, iy_max, iy_min, iz, iz_max, iz_min, l, last_sgf, lx, &
351 : ly, lz, n_cart_total, row_idx
352 : REAL(KIND=dp) :: alpha, cell_vector(3), dist_vec(3), &
353 : dist_vec_raw(3), exp_val, poly, r2, &
354 : r_atom(3), weight
355 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
356 : TYPE(cell_type), POINTER :: cell
357 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
358 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
359 :
360 0 : CALL timeset(routineN, handle)
361 :
362 : ! Get Atom Info
363 0 : cell => bs_env%ri_rs%cell
364 0 : particle_set => bs_env%ri_rs%particle_set
365 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
366 0 : orb_basis_set => bs_env%basis_set_AO(ikind)%gto_basis_set
367 0 : CALL get_cell(cell=cell, h=hmat)
368 0 : CPASSERT(ASSOCIATED(orb_basis_set))
369 :
370 0 : IF (cell%perd(1) == 1) THEN; ix_min = -1; ix_max = 1; ELSE; ix_min = 0; ix_max = 0
371 : END IF
372 0 : IF (cell%perd(2) == 1) THEN; iy_min = -1; iy_max = 1; ELSE; iy_min = 0; iy_max = 0
373 : END IF
374 0 : IF (cell%perd(3) == 1) THEN; iz_min = -1; iz_max = 1; ELSE; iz_min = 0; iz_max = 0
375 : END IF
376 :
377 0 : r_atom = particle_set(iatom)%r
378 :
379 : !$OMP PARALLEL DO DEFAULT(NONE) &
380 : !$OMP SHARED(phi_val, ri_rs_grid, npts, orb_basis_set, r_atom, hmat, &
381 : !$OMP ncoset, indco, cell, ix_min, ix_max, &
382 : !$OMP iy_min, iy_max, iz_min, iz_max, r2_threshold) &
383 : !$OMP PRIVATE(i_pt, dist_vec_raw, ix, iy, iz, cell_vector, dist_vec, r2, iset, &
384 : !$OMP n_cart_total, ishell, l, istart_co, iend_co, first_sgf, last_sgf, &
385 : !$OMP ipgf, alpha, exp_val, isgf, ico, row_idx, weight, lx, ly, lz, poly) &
386 0 : !$OMP SCHEDULE(DYNAMIC)
387 :
388 : DO i_pt = 1, npts
389 :
390 : dist_vec_raw = ri_rs_grid(:, i_pt) - r_atom
391 :
392 : DO ix = ix_min, ix_max
393 : DO iy = iy_min, iy_max
394 : DO iz = iz_min, iz_max
395 :
396 : cell_vector(1:3) = MATMUL(hmat, REAL([ix, iy, iz], dp))
397 :
398 : dist_vec = dist_vec_raw - cell_vector
399 :
400 : r2 = DOT_PRODUCT(dist_vec, dist_vec)
401 :
402 : IF (r2 > r2_threshold) CYCLE
403 :
404 : DO iset = 1, orb_basis_set%nset
405 : n_cart_total = ncoset(orb_basis_set%lmax(iset))
406 :
407 : DO ishell = 1, orb_basis_set%nshell(iset)
408 : l = orb_basis_set%l(ishell, iset)
409 : istart_co = ncoset(l - 1) + 1
410 : iend_co = ncoset(l)
411 :
412 : first_sgf = orb_basis_set%first_sgf(ishell, iset)
413 : last_sgf = orb_basis_set%last_sgf(ishell, iset)
414 :
415 : DO ipgf = 1, orb_basis_set%npgf(iset)
416 : alpha = orb_basis_set%zet(ipgf, iset)
417 : exp_val = EXP(-alpha*r2)
418 :
419 : DO isgf = first_sgf, last_sgf
420 : DO ico = istart_co, iend_co
421 : row_idx = (ipgf - 1)*n_cart_total + ico
422 : weight = orb_basis_set%sphi(row_idx, isgf)
423 : lx = indco(1, ico)
424 : ly = indco(2, ico)
425 : lz = indco(3, ico)
426 : poly = (dist_vec(1)**lx)*(dist_vec(2)**ly)*(dist_vec(3)**lz)
427 :
428 : phi_val(i_pt, isgf) = phi_val(i_pt, isgf) + (weight*poly*exp_val)
429 :
430 : END DO
431 : END DO
432 : END DO
433 : END DO
434 : END DO
435 : END DO
436 : END DO
437 : END DO
438 : END DO
439 : !$OMP END PARALLEL DO
440 :
441 0 : CALL timestop(handle)
442 :
443 0 : END SUBROUTINE fill_phi_for_atom
444 :
445 : ! **************************************************************************************************
446 : !> \brief Compute RI-RS Coefficients (Z_lP)
447 : !> \param qs_env ...
448 : !> \param bs_env ...
449 : !> \param ri_rs_grid_points ...
450 : !> \param mat_phi_mu_l ...
451 : !> \param mat_Z_lP ...
452 : ! **************************************************************************************************
453 :
454 0 : SUBROUTINE compute_Z_lP(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
455 :
456 : ! Arguments
457 : TYPE(qs_environment_type), POINTER :: qs_env
458 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
459 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
460 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
461 : TYPE(dbcsr_type), INTENT(OUT) :: mat_Z_lP
462 :
463 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
464 : routineN = 'compute_Z_lP'
465 :
466 : INTEGER :: atom_P, col_end, col_start, group_handle, handle, handle_dpotrf, handle_dpotrs, &
467 : i_blk, ikind, info, j, l, loc_idx, max_ao_size, max_loc_ri, n_ao_total, n_grid_total, &
468 : n_loc_ri, n_local_grid, natom, nkind, npcol_phi, num_grid_chunks, P_loop_atom, ri_atom
469 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: local_grid_idx, row_offset
470 0 : INTEGER, DIMENSION(:), POINTER :: col_dist_ri, r_blk_sizes, ri_blk_sizes, &
471 0 : row_dist_grid
472 : REAL(KIND=dp) :: cutoff_ri, cutoff_ri_2, d_sP, dist2_min, &
473 : r2_threshold, r_c, t1, t2, t3
474 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cutoff_ri_per_atom, cutoff_ri_per_kind, &
475 0 : d_vec_local
476 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: D_local, d_lp_local, phi_local, &
477 0 : sphere_grid
478 : REAL(KIND=dp), DIMENSION(3) :: dist_vec_raw, pos_P
479 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
480 : TYPE(cell_type), POINTER :: cell
481 : TYPE(cp_logger_type), POINTER :: logger
482 : TYPE(dbcsr_distribution_type) :: dist_phi, dist_Z
483 0 : TYPE(gw_3c_ctx_type) :: ctx_3c
484 : TYPE(mp_para_env_type), POINTER :: para_env
485 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
486 : TYPE(section_vals_type), POINTER :: input
487 :
488 0 : CALL timeset(routineN, handle)
489 :
490 0 : t1 = m_walltime()
491 :
492 0 : CALL get_qs_env(qs_env, para_env=para_env, particle_set=particle_set, input=input, cell=cell)
493 0 : atomic_kind_set => bs_env%ri_rs%atomic_kind_set
494 :
495 0 : natom = bs_env%n_atom
496 0 : n_ao_total = bs_env%i_ao_end_from_atom(natom)
497 0 : n_grid_total = bs_env%ri_rs%n_grid_points
498 0 : CPASSERT(SIZE(ri_rs_grid_points, 2) == n_grid_total)
499 :
500 : ! =========================================================================
501 : ! 1. SETUP DBCSR TOPOLOGY & EXACT OFFSETS
502 : ! =========================================================================
503 0 : CALL dbcsr_get_info(mat_phi_mu_l, row_blk_size=r_blk_sizes, distribution=dist_phi)
504 : CALL dbcsr_distribution_get(dist_phi, row_dist=row_dist_grid, &
505 0 : group=group_handle, npcols=npcol_phi)
506 :
507 0 : num_grid_chunks = SIZE(r_blk_sizes)
508 :
509 0 : ALLOCATE (row_offset(num_grid_chunks))
510 0 : row_offset(1) = 0
511 0 : DO i_blk = 2, num_grid_chunks
512 0 : row_offset(i_blk) = row_offset(i_blk - 1) + r_blk_sizes(i_blk - 1)
513 : END DO
514 :
515 0 : ALLOCATE (ri_blk_sizes(natom), col_dist_ri(natom))
516 0 : DO atom_P = 1, natom
517 0 : ri_blk_sizes(atom_P) = bs_env%i_RI_end_from_atom(atom_P) - bs_env%i_RI_start_from_atom(atom_P) + 1
518 0 : col_dist_ri(atom_P) = MOD(atom_P - 1, npcol_phi)
519 : END DO
520 :
521 0 : CALL dbcsr_distribution_new(dist_Z, template=dist_phi, row_dist=row_dist_grid, col_dist=col_dist_ri)
522 :
523 0 : IF (bs_env%ri_rs%Z_lP_exists) THEN
524 : CALL dbcsr_binary_read(filepath=TRIM(bs_env%prefix)//"Z_lP.matrix", &
525 : distribution=dist_Z, &
526 0 : matrix_new=mat_Z_lP)
527 0 : IF (bs_env%unit_nr > 0) THEN
528 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
529 0 : 'Read Z_lP from file ', ' Execution time', m_walltime() - t1, ' s'
530 0 : WRITE (bs_env%unit_nr, '(A)') ' '
531 : END IF
532 : ELSE
533 :
534 : CALL dbcsr_create(mat_Z_lP, name="mat_Z_lP", dist=dist_Z, &
535 : matrix_type=dbcsr_type_no_symmetry, &
536 0 : row_blk_size=r_blk_sizes, col_blk_size=ri_blk_sizes)
537 :
538 0 : max_ao_size = 0
539 0 : DO j = 1, natom
540 0 : max_ao_size = MAX(max_ao_size, bs_env%i_ao_end_from_atom(j) - bs_env%i_ao_start_from_atom(j) + 1)
541 : END DO
542 0 : max_loc_ri = MAXVAL(ri_blk_sizes)
543 :
544 : ! Per-atom RI-RS integration sphere:
545 : ! cutoff_ri(P) = r_c + r_AO(P)
546 : ! where r_c is the truncated-Coulomb cutoff of the RI metric.
547 0 : nkind = SIZE(atomic_kind_set)
548 0 : ALLOCATE (cutoff_ri_per_atom(natom))
549 :
550 0 : r_c = bs_env%ri_metric%cutoff_radius
551 0 : DO P_loop_atom = 1, natom
552 : cutoff_ri_per_atom(P_loop_atom) = &
553 0 : r_c + bs_env%ri_rs%radius_ao_per_atom(P_loop_atom)
554 : END DO
555 :
556 0 : ALLOCATE (cutoff_ri_per_kind(nkind))
557 0 : cutoff_ri_per_kind(:) = 0.0_dp
558 0 : IF (bs_env%unit_nr > 0) THEN
559 0 : DO P_loop_atom = 1, natom
560 0 : ikind = particle_set(P_loop_atom)%atomic_kind%kind_number
561 : cutoff_ri_per_kind(ikind) = MAX(cutoff_ri_per_kind(ikind), &
562 0 : cutoff_ri_per_atom(P_loop_atom))
563 : END DO
564 0 : WRITE (bs_env%unit_nr, '(T2,A)') 'Per-kind maximum RI-RS sphere cutoff (Bohr):'
565 0 : WRITE (bs_env%unit_nr, '(T4,A4,A14)') 'Kind', 'max cutoff_ri'
566 0 : DO ikind = 1, nkind
567 : WRITE (bs_env%unit_nr, '(T4,A4,F14.4)') &
568 0 : atomic_kind_set(ikind)%element_symbol, &
569 0 : cutoff_ri_per_kind(ikind)
570 : END DO
571 0 : WRITE (bs_env%unit_nr, '(A)') ' '
572 0 : DEALLOCATE (cutoff_ri_per_kind)
573 : END IF
574 :
575 : ! Shared 3c-integral context: hoists libint / t_c_g0 / md_ftable / contracted
576 : ! sphi tables out of the per-triple call so compute_d_lp threads only allocate
577 : ! a lightweight per-thread workspace. MPI-collective; must be outside any
578 : ! OMP region.
579 : CALL gw_3c_ctx_create(ctx_3c, bs_env, bs_env%ri_metric, &
580 : basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, &
581 0 : basis_i=bs_env%basis_set_RI)
582 :
583 : ! =========================================================================
584 : ! 2. MPI LOOP OVER ATOMS (Fully independent, no MPI barriers inside)
585 : ! phi_local for each atom_P's cutoff sphere is built on the fly via
586 : ! fill_phi_for_atom — no dense replicated phi_global, no allreduce.
587 : ! =========================================================================
588 0 : DO atom_P = para_env%mepos + 1, natom, para_env%num_pe
589 :
590 0 : n_loc_ri = ri_blk_sizes(atom_P)
591 0 : pos_P(:) = particle_set(atom_P)%r(:)
592 :
593 0 : cutoff_ri = cutoff_ri_per_atom(atom_P)
594 0 : cutoff_ri_2 = cutoff_ri**2
595 :
596 : ! ---------------------------------------------------------------------
597 : ! A. Determine Local Grid Domain based on cutoff_ri (PBC distance)
598 : ! ---------------------------------------------------------------------
599 0 : n_local_grid = 0
600 0 : DO l = 1, n_grid_total
601 0 : dist_vec_raw = pbc(ri_rs_grid_points(1:3, l), pos_P(1:3), cell)
602 0 : dist2_min = DOT_PRODUCT(dist_vec_raw, dist_vec_raw)
603 0 : IF (dist2_min <= cutoff_ri_2) n_local_grid = n_local_grid + 1
604 : END DO
605 :
606 0 : ALLOCATE (local_grid_idx(n_local_grid))
607 :
608 0 : n_local_grid = 0
609 0 : DO l = 1, n_grid_total
610 0 : dist_vec_raw = pbc(ri_rs_grid_points(1:3, l), pos_P(1:3), cell)
611 0 : dist2_min = DOT_PRODUCT(dist_vec_raw, dist_vec_raw)
612 0 : IF (dist2_min <= cutoff_ri_2) THEN
613 0 : n_local_grid = n_local_grid + 1
614 0 : local_grid_idx(n_local_grid) = l
615 : END IF
616 : END DO
617 :
618 : ! ---------------------------------------------------------------------
619 : ! B. Build phi_local on the fly via fill_phi_for_atom.
620 : ! Only source atoms whose AO basis can reach the cutoff sphere of
621 : ! atom_P (MIC distance) contribute; the rest are pruned. The
622 : ! periodic fill_phi_for_atom sums over (ix, iy, iz) images of
623 : ! ri_atom internally.
624 : ! ---------------------------------------------------------------------
625 0 : ALLOCATE (sphere_grid(3, n_local_grid))
626 0 : DO loc_idx = 1, n_local_grid
627 0 : sphere_grid(:, loc_idx) = ri_rs_grid_points(:, local_grid_idx(loc_idx))
628 : END DO
629 :
630 0 : ALLOCATE (phi_local(n_local_grid, n_ao_total))
631 0 : phi_local = 0.0_dp
632 :
633 0 : DO ri_atom = 1, natom
634 0 : dist_vec_raw = pbc(particle_set(ri_atom)%r(:), pos_P(:), cell)
635 0 : d_sP = NORM2(dist_vec_raw)
636 0 : IF (d_sP > bs_env%ri_rs%radius_ao_per_atom(ri_atom) + cutoff_ri) CYCLE
637 :
638 0 : col_start = bs_env%i_ao_start_from_atom(ri_atom)
639 0 : col_end = bs_env%i_ao_end_from_atom(ri_atom)
640 0 : r2_threshold = bs_env%ri_rs%radius_ao_per_atom(ri_atom)**2
641 :
642 : CALL fill_phi_for_atom(phi_local(:, col_start:col_end), sphere_grid, &
643 0 : n_local_grid, ri_atom, bs_env, r2_threshold)
644 : END DO
645 :
646 0 : DEALLOCATE (sphere_grid)
647 :
648 : ! ---------------------------------------------------------------------
649 : ! C. Build Local RHS Matrix (d_lp_local). compute_d_lp does not depend
650 : ! on D_local or d_vec_local.
651 : ! ---------------------------------------------------------------------
652 0 : ALLOCATE (d_lp_local(n_local_grid, n_loc_ri))
653 0 : d_lp_local = 0.0_dp
654 :
655 0 : t2 = m_walltime()
656 :
657 : CALL compute_d_lp(bs_env, ctx_3c, phi_local, d_lp_local, n_local_grid, &
658 0 : n_loc_ri, atom_P, max_ao_size)
659 :
660 0 : t3 = m_walltime()
661 :
662 : ! ---------------------------------------------------------------------
663 : ! D. Build d_vec_local (Jacobi diagonal) + LHS
664 : ! ---------------------------------------------------------------------
665 0 : ALLOCATE (d_vec_local(n_local_grid))
666 :
667 : CALL build_gram_jacobi_blas(phi_local, n_local_grid, n_ao_total, &
668 0 : bs_env%ri_rs%tikhonov, D_local, d_vec_local)
669 :
670 : ! ---------------------------------------------------------------------
671 : ! E. Pre-scale d_lp by d_vec
672 : ! ---------------------------------------------------------------------
673 0 : CALL scale_rows_by_diag(d_lp_local, d_vec_local, n_local_grid, n_loc_ri)
674 :
675 : ! ---------------------------------------------------------------------
676 : ! F. Solve — LAPACK dpotrf/dpotrs
677 : ! ---------------------------------------------------------------------
678 0 : CALL timeset(routineN//"_dpotrf", handle_dpotrf)
679 0 : CALL dpotrf('L', n_local_grid, D_local, n_local_grid, info)
680 0 : CALL timestop(handle_dpotrf)
681 0 : IF (info /= 0) CPABORT("RI-RS Cholesky factorization failed")
682 0 : CALL timeset(routineN//"_dpotrs", handle_dpotrs)
683 : CALL dpotrs('L', n_local_grid, n_loc_ri, D_local, n_local_grid, &
684 0 : d_lp_local, n_local_grid, info)
685 0 : CALL timestop(handle_dpotrs)
686 0 : IF (info /= 0) CPABORT("RI-RS Cholesky solve failed")
687 0 : DEALLOCATE (D_local)
688 :
689 : ! ---------------------------------------------------------------------
690 : ! G. Post-scale solution by d_vec
691 : ! ---------------------------------------------------------------------
692 0 : CALL scale_rows_by_diag(d_lp_local, d_vec_local, n_local_grid, n_loc_ri)
693 :
694 : ! ---------------------------------------------------------------------
695 : ! H. Scatter Local Solution Back to Global DBCSR Matrix.
696 : ! DBCSR routes blocks to their global owner on finalize.
697 : ! local_grid_idx is ascending (built by the ordered scan above), so
698 : ! a single walking pointer over chunks works.
699 : ! ---------------------------------------------------------------------
700 : CALL store_Z_lP_columns(mat_Z_lP, d_lp_local, local_grid_idx, n_local_grid, &
701 : n_loc_ri, atom_P, r_blk_sizes, row_offset, &
702 0 : bs_env%eps_filter)
703 :
704 0 : DEALLOCATE (d_vec_local, d_lp_local)
705 0 : DEALLOCATE (local_grid_idx, phi_local)
706 :
707 : END DO
708 :
709 0 : DEALLOCATE (cutoff_ri_per_atom)
710 0 : CALL gw_3c_ctx_release(ctx_3c)
711 :
712 0 : CALL dbcsr_finalize(mat_Z_lP)
713 :
714 0 : IF (bs_env%unit_nr > 0) THEN
715 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
716 0 : 'Computed Z_lP ', ' Execution time', m_walltime() - t1, ' s'
717 0 : WRITE (bs_env%unit_nr, '(A)') ' '
718 : END IF
719 :
720 0 : logger => cp_get_default_logger()
721 :
722 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
723 0 : CALL dbcsr_binary_write(matrix=mat_Z_lP, filepath=TRIM(bs_env%prefix)//"Z_lP.matrix")
724 : END IF
725 :
726 : END IF
727 :
728 0 : DEALLOCATE (row_offset, ri_blk_sizes, col_dist_ri)
729 0 : CALL dbcsr_distribution_release(dist_Z)
730 :
731 0 : DEALLOCATE (ri_rs_grid_points)
732 :
733 0 : CALL timestop(handle)
734 :
735 0 : END SUBROUTINE compute_Z_lP
736 :
737 : ! **************************************************************************************************
738 : !> \brief Computes the dense localized RHS d_lp(l,P) = Σ_{μν,R,S} ϕ_μ(r_l)·ϕ_ν(r_l)·(μν|P) for one
739 : !> RI atom P. OMP-threaded over (atom_j, atom_k) AO-pair blocks: per thread, sweep all
740 : !> (cell_R, cell_S) periodic images of (atom_j, atom_k) about atom_P at cell (0,0,0); each
741 : !> 3c block is built by build_3c_integral_block_ctx (cached libint / sphi tables in ctx,
742 : !> kind-radius triangle screen → `screened` short-circuits negligible image triples), and
743 : !> grid-chunked pair densities are contracted into a private d_lp partial that is reduced
744 : !> into d_lp at the end of the parallel region.
745 : !> \param bs_env ...
746 : !> \param ctx shared 3c-integral context (gw_3c_ctx_create)
747 : !> \param phi_val ϕ_μ(r_l) on the local-sphere grid (n_grid_total × n_ao)
748 : !> \param d_lp output (n_grid_total × n_loc_ri), zeroed by the caller, accumulated here
749 : !> \param n_grid_total number of local-sphere grid rows
750 : !> \param n_loc_ri number of RI functions of atom_P
751 : !> \param atom_P RI atom (pinned to cell (0,0,0))
752 : !> \param max_ao_size ...
753 : ! **************************************************************************************************
754 :
755 0 : SUBROUTINE compute_d_lp(bs_env, ctx, phi_val, d_lp, n_grid_total, n_loc_ri, atom_P, &
756 : max_ao_size)
757 :
758 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
759 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
760 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val
761 : INTEGER, INTENT(IN) :: n_grid_total, n_loc_ri
762 : REAL(KIND=dp), INTENT(INOUT) :: d_lp(n_grid_total, n_loc_ri)
763 : INTEGER, INTENT(IN) :: atom_P, max_ao_size
764 :
765 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp'
766 : INTEGER, PARAMETER :: grid_chunk = 1024
767 :
768 : INTEGER :: atom_j, atom_k, c, handle, handle_dgemm, ix_max, ix_min, ix_R, ix_S, iy_max, &
769 : iy_min, iy_R, iy_S, iz_max, iz_min, iz_R, iz_S, j, jk_idx, jsize, jstart, k, ksize, &
770 : kstart, l, l0, natom, ri
771 : INTEGER, DIMENSION(3) :: cell_R_vec, cell_S_vec
772 : LOGICAL :: any_kept, screened
773 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_lp_prv, int_2d_prv, rho_chunk
774 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c_prv, int_3c_sum
775 0 : TYPE(gw_3c_ws_type) :: ws
776 :
777 0 : CALL timeset(routineN, handle)
778 :
779 0 : natom = bs_env%n_atom
780 :
781 0 : IF (bs_env%ri_rs%cell%perd(1) == 1) THEN; ix_min = -1; ix_max = 1; ELSE; ix_min = 0; ix_max = 0
782 : END IF
783 0 : IF (bs_env%ri_rs%cell%perd(2) == 1) THEN; iy_min = -1; iy_max = 1; ELSE; iy_min = 0; iy_max = 0
784 : END IF
785 0 : IF (bs_env%ri_rs%cell%perd(3) == 1) THEN; iz_min = -1; iz_max = 1; ELSE; iz_min = 0; iz_max = 0
786 : END IF
787 :
788 : !$OMP PARALLEL DEFAULT(NONE) &
789 : !$OMP SHARED(bs_env, ctx, phi_val, d_lp, n_grid_total, n_loc_ri, atom_P, max_ao_size, &
790 : !$OMP natom, ix_min, ix_max, iy_min, iy_max, iz_min, iz_max) &
791 : !$OMP PRIVATE(any_kept, atom_j, atom_k, c, handle_dgemm, j, jk_idx, jsize, jstart, k, &
792 : !$OMP ksize, kstart, l, l0, ri, ix_R, iy_R, iz_R, ix_S, iy_S, iz_S, cell_R_vec, &
793 0 : !$OMP cell_S_vec, screened, d_lp_prv, int_2d_prv, rho_chunk, int_3c_prv, int_3c_sum, ws)
794 :
795 : CALL gw_3c_ws_create(ws, ctx)
796 : ALLOCATE (int_3c_prv(max_ao_size, max_ao_size, n_loc_ri))
797 : ALLOCATE (int_3c_sum(max_ao_size, max_ao_size, n_loc_ri))
798 : ALLOCATE (int_2d_prv(max_ao_size*max_ao_size, n_loc_ri))
799 : ALLOCATE (rho_chunk(grid_chunk, max_ao_size*max_ao_size))
800 : ALLOCATE (d_lp_prv(n_grid_total, n_loc_ri))
801 : d_lp_prv(:, :) = 0.0_dp
802 :
803 : ! atom_P pinned at cell (0,0,0); enumerate (atom_j, cell_R) × (atom_k, cell_S). The ctx
804 : ! integral builder's kind_radius triangle screen sets screened=.TRUE. for the bulk of
805 : ! image triples (one or both AO atoms beyond the truncated-Coulomb reach of atom_P),
806 : ! so the 27 × 27 = 729 candidate cells collapse to "adjacent cells" in practice.
807 : ! The inner atom_k loop carries enough work for DYNAMIC.
808 : !$OMP DO SCHEDULE(DYNAMIC)
809 : DO atom_j = 1, natom
810 : DO atom_k = 1, natom
811 : jstart = bs_env%i_ao_start_from_atom(atom_j)
812 : jsize = bs_env%i_ao_end_from_atom(atom_j) - jstart + 1
813 : kstart = bs_env%i_ao_start_from_atom(atom_k)
814 : ksize = bs_env%i_ao_end_from_atom(atom_k) - kstart + 1
815 :
816 : int_3c_sum(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp
817 : any_kept = .FALSE.
818 :
819 : DO ix_R = ix_min, ix_max
820 : DO iy_R = iy_min, iy_max
821 : DO iz_R = iz_min, iz_max
822 : cell_R_vec = [ix_R, iy_R, iz_R]
823 : DO ix_S = ix_min, ix_max
824 : DO iy_S = iy_min, iy_max
825 : DO iz_S = iz_min, iz_max
826 : cell_S_vec = [ix_S, iy_S, iz_S]
827 :
828 : int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp
829 :
830 : CALL build_3c_integral_block_ctx(int_3c_prv( &
831 : 1:jsize, 1:ksize, 1:n_loc_ri), ctx, ws, &
832 : atom_j=atom_j, atom_k=atom_k, atom_i=atom_P, &
833 : cell_j=cell_R_vec, cell_k=cell_S_vec, cell_i=[0, 0, 0], &
834 : screened=screened)
835 : IF (screened) CYCLE
836 :
837 : any_kept = .TRUE.
838 : int_3c_sum(1:jsize, 1:ksize, 1:n_loc_ri) = &
839 : int_3c_sum(1:jsize, 1:ksize, 1:n_loc_ri) + &
840 : int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri)
841 : END DO
842 : END DO
843 : END DO
844 : END DO
845 : END DO
846 : END DO
847 :
848 : IF (.NOT. any_kept) CYCLE
849 :
850 : ! Flatten 3D B_{μν,P} → 2D B_{(μν),P}
851 : DO ri = 1, n_loc_ri
852 : DO k = 1, ksize
853 : DO j = 1, jsize
854 : jk_idx = (k - 1)*jsize + j
855 : int_2d_prv(jk_idx, ri) = int_3c_sum(j, k, ri)
856 : END DO
857 : END DO
858 : END DO
859 :
860 : ! Pair density ρ(l,μν) = ϕ_μ(r_l)ϕ_ν(r_l) in grid chunks, contracted on the fly:
861 : ! d_{l,P} += ρ(l,μν) B_{(μν),P} (dgemm runs serially inside the parallel region)
862 : DO l0 = 1, n_grid_total, grid_chunk
863 : c = MIN(grid_chunk, n_grid_total - l0 + 1)
864 : DO k = 1, ksize
865 : DO j = 1, jsize
866 : jk_idx = (k - 1)*jsize + j
867 : DO l = 1, c
868 : rho_chunk(l, jk_idx) = phi_val(l0 + l - 1, jstart + j - 1)* &
869 : phi_val(l0 + l - 1, kstart + k - 1)
870 : END DO
871 : END DO
872 : END DO
873 : CALL timeset(routineN//"_dgemm", handle_dgemm)
874 : CALL dgemm("N", "N", c, n_loc_ri, jsize*ksize, &
875 : 1.0_dp, rho_chunk, grid_chunk, &
876 : int_2d_prv, max_ao_size*max_ao_size, &
877 : 1.0_dp, d_lp_prv(l0, 1), n_grid_total)
878 : CALL timestop(handle_dgemm)
879 : END DO
880 : END DO
881 : END DO
882 : !$OMP END DO
883 :
884 : !$OMP CRITICAL (compute_d_lp_reduce)
885 : d_lp(1:n_grid_total, 1:n_loc_ri) = d_lp(1:n_grid_total, 1:n_loc_ri) + &
886 : d_lp_prv(1:n_grid_total, 1:n_loc_ri)
887 : !$OMP END CRITICAL (compute_d_lp_reduce)
888 :
889 : DEALLOCATE (int_3c_prv, int_3c_sum, int_2d_prv, rho_chunk, d_lp_prv)
890 : CALL gw_3c_ws_release(ws)
891 :
892 : !$OMP END PARALLEL
893 :
894 0 : CALL timestop(handle)
895 :
896 0 : END SUBROUTINE compute_d_lp
897 :
898 : ! **************************************************************************************************
899 : !> \brief Computes the χ(iτ, k=0) matrix
900 : !> \param bs_env ...
901 : !> \param mat_chi_Gamma_tau ...
902 : !> \param mat_phi_mu_l ...
903 : !> \param mat_Z_lP ...
904 : ! **************************************************************************************************
905 :
906 0 : SUBROUTINE get_mat_chi_Gamma_tau(bs_env, mat_chi_Gamma_tau, mat_phi_mu_l, mat_Z_lP)
907 :
908 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
909 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
910 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
911 :
912 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
913 :
914 : INTEGER :: handle, i, i_t, ispin, npcol
915 0 : INTEGER, DIMENSION(:), POINTER :: blk_ao, blk_grid, dist_col_grid, &
916 0 : dist_row_grid
917 : REAL(KIND=dp) :: t1, tau
918 : TYPE(dbcsr_distribution_type) :: dist_grid_grid, dist_phi
919 : TYPE(dbcsr_type) :: matrix_chi_grid, matrix_chi_grid_spin, &
920 : matrix_G_occ_grid, matrix_G_vir_grid
921 :
922 0 : CALL timeset(routineN, handle)
923 :
924 : ! =========================================================================
925 : ! 1. SETUP CORE TOPOLOGIES
926 : ! =========================================================================
927 0 : CALL dbcsr_get_info(mat_phi_mu_l, distribution=dist_phi, row_blk_size=blk_grid, col_blk_size=blk_ao)
928 0 : CALL dbcsr_distribution_get(dist_phi, row_dist=dist_row_grid, npcols=npcol)
929 :
930 : ! Build a perfectly safe column distribution for the Grid dimension
931 0 : ALLOCATE (dist_col_grid(SIZE(blk_grid)))
932 0 : DO i = 1, SIZE(blk_grid)
933 0 : dist_col_grid(i) = MOD(i - 1, npcol)
934 : END DO
935 :
936 : CALL dbcsr_distribution_new(dist_grid_grid, template=dist_phi, &
937 0 : row_dist=dist_row_grid, col_dist=dist_col_grid)
938 :
939 0 : CALL dbcsr_create(matrix_G_occ_grid, "G_occ_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
940 0 : CALL dbcsr_create(matrix_G_vir_grid, "G_vir_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
941 0 : CALL dbcsr_create(matrix_chi_grid, "chi_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
942 0 : CALL dbcsr_create(matrix_chi_grid_spin, "chi_grid_spin", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
943 :
944 : ! =========================================================================
945 : ! 2. MAIN IMAGINARY TIME LOOP
946 : ! =========================================================================
947 0 : DO i_t = 1, bs_env%num_time_freq_points
948 0 : t1 = m_walltime()
949 :
950 0 : tau = bs_env%time_frequency_grid%imaginary_time(i_t)
951 0 : CALL dbcsr_set(matrix_chi_grid, 0.0_dp)
952 :
953 : ! ----------------------------------------------------------------------
954 : ! A. SPIN LOOP (Allocations safely encapsulated in wrappers)
955 : ! ----------------------------------------------------------------------
956 0 : DO ispin = 1, bs_env%n_spin
957 :
958 : ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
959 : ! G^occ_ll'(i|τ|,k=0) = sum_µν ϕ_µ(r_l) G^occ_µν ϕ_ν(r_l')
960 : CALL build_G_grid(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, &
961 0 : matrix_G_occ_grid, bs_env%eps_filter)
962 :
963 : ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
964 : ! G^vir_ll'(i|τ|,k=0) = sum_µν ϕ_µ(r_l) G^vir_µν ϕ_ν(r_l')
965 : CALL build_G_grid(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, &
966 0 : matrix_G_vir_grid, bs_env%eps_filter)
967 :
968 : ! -------------------------------------------------------------------
969 : ! B. ELEMENT-WISE HADAMARD PRODUCT
970 : ! -------------------------------------------------------------------
971 : ! χ_ll'(iτ,k=0) = G^occ_ll'(i|τ|,k=0) * G^vir_ll'(i|τ|,k=0)
972 0 : CALL hadamard_product(matrix_G_occ_grid, matrix_G_vir_grid, matrix_chi_grid_spin, bs_env%spin_degeneracy)
973 :
974 : ! Accumulate spin contributions
975 0 : CALL dbcsr_add(matrix_chi_grid, matrix_chi_grid_spin, 1.0_dp, 1.0_dp)
976 :
977 : END DO ! ispin
978 :
979 : ! ----------------------------------------------------------------------
980 : ! C. TRANSFORM TO AUXILIARY BASIS & EXPORT DIRECTLY
981 : ! χ_aux = Z^T * χ_grid * Z
982 : ! χ_PQ(iτ,k=0) = sum_ll' Z_lP χ_ll'(iτ,k=0) Z_l'Q
983 : ! Result is dumped directly into the final array mat_chi_Gamma_tau!
984 : ! ----------------------------------------------------------------------
985 : CALL dbcsr_contract_ABA("T", "N", mat_Z_lP, matrix_chi_grid, &
986 0 : mat_chi_Gamma_tau(i_t)%matrix, bs_env%eps_filter)
987 :
988 0 : IF (bs_env%unit_nr > 0) THEN
989 : WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
990 0 : 'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
991 0 : ', Execution time', m_walltime() - t1, ' s'
992 : END IF
993 :
994 : END DO ! i_t
995 :
996 : ! =========================================================================
997 : ! 3. FINAL CLEANUP
998 : ! =========================================================================
999 0 : CALL dbcsr_release(matrix_G_occ_grid)
1000 0 : CALL dbcsr_release(matrix_G_vir_grid)
1001 0 : CALL dbcsr_release(matrix_chi_grid)
1002 0 : CALL dbcsr_release(matrix_chi_grid_spin)
1003 0 : CALL dbcsr_distribution_release(dist_grid_grid)
1004 0 : DEALLOCATE (dist_col_grid)
1005 :
1006 0 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1007 :
1008 0 : CALL timestop(handle)
1009 :
1010 0 : END SUBROUTINE get_mat_chi_Gamma_tau
1011 :
1012 : ! **************************************************************************************************
1013 : !> \brief Computes Green's Function in grid basis
1014 : !> \param bs_env ...
1015 : !> \param tau ...
1016 : !> \param ispin ...
1017 : !> \param occ ...
1018 : !> \param vir ...
1019 : !> \param mat_phi_mu_l ...
1020 : !> \param matrix_G_grid ...
1021 : !> \param eps_filter ...
1022 : ! **************************************************************************************************
1023 :
1024 0 : SUBROUTINE build_G_grid(bs_env, tau, ispin, occ, vir, mat_phi_mu_l, matrix_G_grid, eps_filter)
1025 :
1026 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1027 : REAL(KIND=dp), INTENT(IN) :: tau
1028 : INTEGER, INTENT(IN) :: ispin
1029 : LOGICAL, INTENT(IN) :: occ, vir
1030 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, matrix_G_grid
1031 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1032 :
1033 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_G_grid'
1034 :
1035 : INTEGER :: handle
1036 0 : INTEGER, DIMENSION(:), POINTER :: blk_ao, dist_row_ao
1037 : TYPE(cp_fm_type), POINTER :: fm_G
1038 : TYPE(dbcsr_distribution_type) :: dist_ao_ao
1039 : TYPE(dbcsr_type) :: matrix_G_ao
1040 :
1041 0 : CALL timeset(routineN, handle)
1042 :
1043 : ! 1. Select the correct FM matrix based on occ/vir flags
1044 0 : IF (occ) THEN
1045 0 : fm_G => bs_env%fm_Gocc
1046 : ELSE
1047 0 : fm_G => bs_env%fm_Gvir
1048 : END IF
1049 :
1050 : ! 2. Compute Dense FM Green's Function
1051 : ! G^occ/vir_µλ(i|τ|,k=0) = sum_G^occ/vir_µλn^occ/vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1052 0 : CALL G_occ_vir(bs_env, tau, fm_G, ispin, occ=occ, vir=vir)
1053 :
1054 : ! 3. Setup AO DBCSR Topology and Create Matrix dynamically
1055 0 : CALL setup_square_topology(mat_phi_mu_l, 'COL', dist_ao_ao, blk_ao, dist_row_ao)
1056 :
1057 : CALL dbcsr_create(matrix_G_ao, name="G_ao", dist=dist_ao_ao, &
1058 : matrix_type=dbcsr_type_no_symmetry, &
1059 0 : row_blk_size=blk_ao, col_blk_size=blk_ao)
1060 :
1061 : ! 4. Convert FM to Sparse DBCSR
1062 0 : CALL copy_fm_to_dbcsr(fm_G, matrix_G_ao, keep_sparsity=.FALSE.)
1063 :
1064 : ! 5. Transform to Grid Basis: G_grid = phi * G_ao * phi^T
1065 : ! G^occ/vir_ll'(i|τ|,k=0) = sum_µν ϕ_µ(r_l) G^occ/vir_µν ϕ_ν(r_l')
1066 : CALL dbcsr_contract_ABA("N", "T", mat_phi_mu_l, matrix_G_ao, matrix_G_grid, &
1067 0 : eps_filter)
1068 :
1069 : ! 6. Release AO matrix and topology
1070 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_ao_ao, mapped_dist=dist_row_ao, m1=matrix_G_ao)
1071 :
1072 0 : CALL timestop(handle)
1073 :
1074 0 : END SUBROUTINE build_G_grid
1075 :
1076 : ! **************************************************************************************************
1077 : !> \brief Computes the exact exchange part of the GW self-energy
1078 : !> \param bs_env ...
1079 : !> \param qs_env ...
1080 : !> \param mat_phi_mu_l ...
1081 : !> \param mat_Z_lP ...
1082 : !> \param fm_Sigma_x_Gamma ...
1083 : ! **************************************************************************************************
1084 :
1085 0 : SUBROUTINE compute_Sigma_x(bs_env, qs_env, mat_phi_mu_l, mat_Z_lP, fm_Sigma_x_Gamma)
1086 :
1087 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1088 : TYPE(qs_environment_type), POINTER :: qs_env
1089 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
1090 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
1091 :
1092 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_x'
1093 :
1094 : INTEGER :: handle, ispin
1095 0 : INTEGER, DIMENSION(:), POINTER :: blk_aux, blk_grid, dist_col_grid, &
1096 0 : dist_row_aux
1097 : REAL(KIND=dp) :: t1
1098 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
1099 : TYPE(dbcsr_distribution_type) :: dist_aux_aux, dist_grid_grid
1100 : TYPE(dbcsr_type) :: mat_Sigma_x_Gamma, matrix_D_grid, &
1101 : matrix_Sigma_x_grid, matrix_V_aux, &
1102 : matrix_V_grid
1103 :
1104 0 : CALL timeset(routineN, handle)
1105 :
1106 0 : t1 = m_walltime()
1107 :
1108 0 : ALLOCATE (fm_Sigma_x_Gamma(bs_env%n_spin))
1109 0 : DO ispin = 1, bs_env%n_spin
1110 0 : CALL cp_fm_create(fm_Sigma_x_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1111 : END DO
1112 :
1113 0 : CALL dbcsr_create(mat_Sigma_x_Gamma, template=bs_env%mat_ao_ao%matrix)
1114 :
1115 : ! =========================================================================
1116 : ! 1. SETUP CORE TOPOLOGIES
1117 : ! =========================================================================
1118 0 : CALL setup_square_topology(mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, dist_col_grid)
1119 0 : CALL setup_square_topology(mat_Z_lP, 'COL', dist_aux_aux, blk_aux, dist_row_aux)
1120 :
1121 : ! =========================================================================
1122 : ! 2. COMPUTE V^tr_ll'
1123 : ! =========================================================================
1124 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI%matrix_struct, bs_env%n_RI, &
1125 0 : bs_env%trunc_coulomb)
1126 :
1127 : ! M^-1(k=0) V^tr(τ) M^-1(k=0) -> fm_Vtr_Gamma
1128 0 : CALL fm_contract_ABA(bs_env%fm_Minv_Gamma, fm_Vtr_Gamma(:, 1))
1129 :
1130 0 : CALL dbcsr_create(matrix_V_aux, "V_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
1131 0 : CALL dbcsr_create(matrix_V_grid, "V_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1132 :
1133 0 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.FALSE.)
1134 :
1135 : ! V^tr_ll' = sum_PQ Z_lP V^trunc_PQ Z_l'Q
1136 : CALL dbcsr_contract_ABA("N", "T", mat_Z_lP, matrix_V_aux, matrix_V_grid, &
1137 0 : bs_env%eps_filter)
1138 0 : CALL dbcsr_release(matrix_V_aux)
1139 :
1140 : ! =========================================================================
1141 : ! 3. SPIN LOOP FOR EXACT EXCHANGE
1142 : ! =========================================================================
1143 0 : DO ispin = 1, bs_env%n_spin
1144 :
1145 : ! Density matrix on grid is essentially G_occ at tau = 0.0
1146 : ! D_µν = sum_n^occ C_µn(k=0) C_νn(k=0)
1147 : ! D_ll' = sum_µν ϕ_µ(r_l) D_µν ϕ_ν(r_l')
1148 0 : CALL dbcsr_create(matrix_D_grid, "D_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1149 0 : CALL build_G_grid(bs_env, 0.0_dp, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_D_grid, bs_env%eps_filter)
1150 :
1151 : ! Element-wise Hadamard product: Σ^x_grid = D_grid ◦ V_grid
1152 : ! Σ^x_ll' = D_ll' * V^tr_ll'
1153 0 : CALL dbcsr_create(matrix_Sigma_x_grid, template=matrix_V_grid)
1154 0 : CALL hadamard_product(matrix_D_grid, matrix_V_grid, matrix_Sigma_x_grid, 1.0_dp)
1155 :
1156 0 : CALL dbcsr_release(matrix_D_grid)
1157 :
1158 : ! Transform back to AO basis: Σ^x_ao = -1.0 * phi^T * Σ^x_grid * phi
1159 : ! Σ^x_λσ(k=0) = -sum_ll' ϕ_λ(r_l) Σ^x_ll' ϕ_σ(r_l')
1160 : CALL dbcsr_contract_ABA("T", "N", mat_phi_mu_l, matrix_Sigma_x_grid, &
1161 0 : mat_Sigma_x_Gamma, bs_env%eps_filter)
1162 0 : CALL dbcsr_scale(mat_Sigma_x_Gamma, -1.0_dp)
1163 :
1164 0 : CALL dbcsr_release(matrix_Sigma_x_grid)
1165 :
1166 : ! Data I/O and Export to CP2K Full Matrices
1167 0 : CALL copy_dbcsr_to_fm(mat_Sigma_x_Gamma, fm_Sigma_x_Gamma(ispin))
1168 :
1169 : END DO ! ispin
1170 :
1171 0 : IF (bs_env%unit_nr > 0) THEN
1172 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
1173 0 : 'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
1174 0 : WRITE (bs_env%unit_nr, '(A)') ' '
1175 : END IF
1176 :
1177 : ! =========================================================================
1178 : ! 4. CLEANUP
1179 : ! =========================================================================
1180 : CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid, &
1181 0 : m1=mat_Sigma_x_Gamma, m2=matrix_V_grid)
1182 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
1183 :
1184 0 : CALL cp_fm_release(fm_Vtr_Gamma)
1185 :
1186 0 : CALL timestop(handle)
1187 :
1188 0 : END SUBROUTINE compute_Sigma_x
1189 :
1190 : ! **************************************************************************************************
1191 : !> \brief Computes the correlation part of the GW self-energy
1192 : !> \param bs_env ...
1193 : !> \param fm_W_time ...
1194 : !> \param mat_phi_mu_l ...
1195 : !> \param mat_Z_lP ...
1196 : !> \param fm_Sigma_c_Gamma_time ...
1197 : ! **************************************************************************************************
1198 :
1199 0 : SUBROUTINE compute_Sigma_c(bs_env, fm_W_time, mat_phi_mu_l, mat_Z_lP, fm_Sigma_c_Gamma_time)
1200 :
1201 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1202 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
1203 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
1204 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
1205 :
1206 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c'
1207 :
1208 : INTEGER :: handle, i_t, ispin
1209 0 : INTEGER, DIMENSION(:), POINTER :: blk_aux, blk_grid, dist_col_grid, &
1210 0 : dist_row_aux
1211 : REAL(KIND=dp) :: t1, tau
1212 : TYPE(dbcsr_distribution_type) :: dist_aux_aux, dist_grid_grid
1213 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
1214 : TYPE(dbcsr_type) :: matrix_G_occ_grid, matrix_G_vir_grid, matrix_Sigma_neg_grid, &
1215 : matrix_Sigma_pos_grid, matrix_W_aux, matrix_W_grid
1216 :
1217 0 : CALL timeset(routineN, handle)
1218 :
1219 : ! =========================================================================
1220 : ! 1. SETUP CORE TOPOLOGIES AND PRE-ALLOCATE OUTPUT ARRAYS
1221 : ! =========================================================================
1222 0 : CALL setup_square_topology(mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, dist_col_grid)
1223 0 : CALL setup_square_topology(mat_Z_lP, 'COL', dist_aux_aux, blk_aux, dist_row_aux)
1224 :
1225 : ! Pre-allocate local DBCSR matrices to act as targets for final output
1226 0 : NULLIFY (mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1227 0 : ALLOCATE (mat_Sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1228 0 : ALLOCATE (mat_Sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1229 :
1230 0 : DO i_t = 1, bs_env%num_time_freq_points
1231 0 : DO ispin = 1, bs_env%n_spin
1232 0 : ALLOCATE (mat_Sigma_neg_tau(i_t, ispin)%matrix)
1233 0 : ALLOCATE (mat_Sigma_pos_tau(i_t, ispin)%matrix)
1234 0 : CALL dbcsr_create(mat_Sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1235 0 : CALL dbcsr_create(mat_Sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1236 : END DO
1237 : END DO
1238 :
1239 : ! =========================================================================
1240 : ! 2. MAIN IMAGINARY TIME LOOP
1241 : ! =========================================================================
1242 0 : DO i_t = 1, bs_env%num_time_freq_points
1243 0 : tau = bs_env%time_frequency_grid%imaginary_time(i_t)
1244 :
1245 : ! -------------------------------------------------------------------
1246 : ! Compute W_grid = Z * W_aux * Z^T
1247 : ! -------------------------------------------------------------------
1248 0 : CALL dbcsr_create(matrix_W_aux, "W_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
1249 0 : CALL dbcsr_create(matrix_W_grid, "W_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1250 :
1251 0 : CALL copy_fm_to_dbcsr(fm_W_time(i_t), matrix_W_aux, keep_sparsity=.FALSE.)
1252 :
1253 : ! W^MIC_ll'(iτ,k=0) = sum_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
1254 : CALL dbcsr_contract_ABA("N", "T", mat_Z_lP, matrix_W_aux, matrix_W_grid, &
1255 0 : bs_env%eps_filter)
1256 :
1257 0 : CALL dbcsr_release(matrix_W_aux) ! Clean up aux basis immediately
1258 :
1259 0 : DO ispin = 1, bs_env%n_spin
1260 0 : t1 = m_walltime()
1261 :
1262 : ! -------------------------------------------------------------------
1263 : ! A. Transform Green's Functions to the Grid
1264 : ! -------------------------------------------------------------------
1265 0 : CALL dbcsr_create(matrix_G_occ_grid, "G_occ_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1266 0 : CALL dbcsr_create(matrix_G_vir_grid, "G_vir_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1267 :
1268 : ! G^occ_µλ(i|τ|,k=0) = sum_G^occ_µλn^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1269 : ! G^occ_ll'(i|τ|,k=0) = sum_µν ϕ_µ(r_l) G^occ_µν ϕ_ν(r_l')
1270 0 : CALL build_G_grid(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_G_occ_grid, bs_env%eps_filter)
1271 :
1272 : ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1273 : ! G^vir_ll'(i|τ|,k=0) = sum_µν ϕ_µ(r_l) G^vir_µν ϕ_ν(r_l')
1274 0 : CALL build_G_grid(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, matrix_G_vir_grid, bs_env%eps_filter)
1275 :
1276 : ! -------------------------------------------------------------------
1277 : ! B. Element-wise Hadamard Products for Sigma_c on Grid
1278 : ! Σ_neg_grid = G_occ_grid ◦ W_grid
1279 : ! Σ_pos_grid = G_vir_grid ◦ W_grid
1280 : ! -------------------------------------------------------------------
1281 0 : CALL dbcsr_create(matrix_Sigma_neg_grid, template=matrix_W_grid)
1282 0 : CALL dbcsr_create(matrix_Sigma_pos_grid, template=matrix_W_grid)
1283 :
1284 : ! Σ^c_ll'(iτ,k=0) = -G^occ_ll'(i|τ|,k=0) * W^MIC_ll'(iτ,k=0), for τ < 0
1285 0 : CALL hadamard_product(matrix_G_occ_grid, matrix_W_grid, matrix_Sigma_neg_grid, 1.0_dp)
1286 :
1287 : ! Σ^c_ll'(iτ,k=0) = G^vir_ll'(i|τ|,k=0) * W^MIC_ll'(iτ,k=0), for τ > 0
1288 0 : CALL hadamard_product(matrix_G_vir_grid, matrix_W_grid, matrix_Sigma_pos_grid, 1.0_dp)
1289 :
1290 : ! Instantly purge massive G_grid arrays to save memory
1291 0 : CALL dbcsr_release(matrix_G_occ_grid)
1292 0 : CALL dbcsr_release(matrix_G_vir_grid)
1293 :
1294 : ! -------------------------------------------------------------------
1295 : ! C. Transform Sigma back to AO Basis
1296 : ! Σ_AO = phi^T * Σ_grid * phi
1297 : ! -------------------------------------------------------------------
1298 :
1299 : ! Σ^c_λσ(iτ,k=0) = sum_ll' ϕ_λ(r_l) Σ^c_ll'(iτ,k=0) ϕ_σ(r_l'), for τ < 0
1300 : CALL dbcsr_contract_ABA("T", "N", mat_phi_mu_l, matrix_Sigma_neg_grid, &
1301 0 : mat_Sigma_neg_tau(i_t, ispin)%matrix, bs_env%eps_filter)
1302 0 : CALL dbcsr_scale(mat_Sigma_neg_tau(i_t, ispin)%matrix, -1.0_dp)
1303 :
1304 : ! Σ^c_λσ(iτ,k=0) = sum_ll' ϕ_λ(r_l) Σ^c_ll'(iτ,k=0) ϕ_σ(r_l'), for τ > 0
1305 : CALL dbcsr_contract_ABA("T", "N", mat_phi_mu_l, matrix_Sigma_pos_grid, &
1306 0 : mat_Sigma_pos_tau(i_t, ispin)%matrix, bs_env%eps_filter)
1307 :
1308 : ! Purge Grid Sigma arrays
1309 0 : CALL dbcsr_release(matrix_Sigma_neg_grid)
1310 0 : CALL dbcsr_release(matrix_Sigma_pos_grid)
1311 :
1312 0 : IF (bs_env%unit_nr > 0) THEN
1313 : WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
1314 0 : 'Computed Σ^c(iτ,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1315 0 : ', Execution time', m_walltime() - t1, ' s'
1316 : END IF
1317 :
1318 : END DO ! ispin
1319 :
1320 : ! Release the W_grid for this time point
1321 0 : CALL dbcsr_release(matrix_W_grid)
1322 :
1323 : END DO ! i_t
1324 :
1325 0 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1326 :
1327 : ! -------------------------------------------------------------------------
1328 : ! 3. FINALIZE AND CLEANUP
1329 : ! -------------------------------------------------------------------------
1330 : CALL fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
1331 0 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
1332 :
1333 0 : CALL cp_fm_release(fm_W_time)
1334 :
1335 0 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_neg_tau)
1336 0 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_pos_tau)
1337 :
1338 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid)
1339 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
1340 :
1341 0 : CALL delete_unnecessary_files(bs_env)
1342 0 : CALL timestop(handle)
1343 :
1344 0 : END SUBROUTINE compute_Sigma_c
1345 :
1346 : ! **************************************************************************************************
1347 : !> \brief DBCSR Topology Generation
1348 : !> \param matrix_template ...
1349 : !> \param dim_type ...
1350 : !> \param square_dist ...
1351 : !> \param blk_sizes ...
1352 : !> \param mapped_dist ...
1353 : ! **************************************************************************************************
1354 :
1355 26880 : SUBROUTINE setup_square_topology(matrix_template, dim_type, square_dist, blk_sizes, mapped_dist)
1356 :
1357 : TYPE(dbcsr_type), INTENT(IN) :: matrix_template
1358 : CHARACTER(LEN=*), INTENT(IN) :: dim_type
1359 : TYPE(dbcsr_distribution_type), INTENT(OUT) :: square_dist
1360 : INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: blk_sizes, mapped_dist
1361 :
1362 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_square_topology'
1363 :
1364 : INTEGER :: handle, i, np, npcols, nprows
1365 8960 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist, row_blk, row_dist
1366 : TYPE(dbcsr_distribution_type) :: dist_template
1367 :
1368 8960 : CALL timeset(routineN, handle)
1369 :
1370 : CALL dbcsr_get_info(matrix_template, distribution=dist_template, &
1371 8960 : row_blk_size=row_blk, col_blk_size=col_blk)
1372 : CALL dbcsr_distribution_get(dist_template, row_dist=row_dist, col_dist=col_dist, &
1373 8960 : nprows=nprows, npcols=npcols)
1374 :
1375 8960 : IF (TRIM(dim_type) == 'ROW') THEN
1376 : ! Creates ROW x ROW (e.g., Grid x Grid from mat_phi_mu_l)
1377 600 : blk_sizes => row_blk
1378 600 : np = npcols
1379 1800 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
1380 6960 : DO i = 1, SIZE(blk_sizes)
1381 6960 : mapped_dist(i) = MOD(i - 1, np)
1382 : END DO
1383 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
1384 600 : row_dist=row_dist, col_dist=mapped_dist)
1385 :
1386 8360 : ELSE IF (TRIM(dim_type) == 'COL') THEN
1387 : ! Creates COL x COL (e.g., Aux x Aux from mat_Z_lP)
1388 8360 : blk_sizes => col_blk
1389 8360 : np = nprows
1390 25080 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
1391 25080 : DO i = 1, SIZE(blk_sizes)
1392 25080 : mapped_dist(i) = MOD(i - 1, np)
1393 : END DO
1394 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
1395 8360 : row_dist=mapped_dist, col_dist=col_dist)
1396 : END IF
1397 :
1398 8960 : CALL timestop(handle)
1399 :
1400 8960 : END SUBROUTINE setup_square_topology
1401 :
1402 : ! **************************************************************************************************
1403 : !> \brief DBCSR matrices deallocation
1404 : !> \param dist ...
1405 : !> \param mapped_dist ...
1406 : !> \param m1 ...
1407 : !> \param m2 ...
1408 : !> \param m3 ...
1409 : !> \param m4 ...
1410 : ! **************************************************************************************************
1411 :
1412 13560 : SUBROUTINE release_dbcsr_topology_and_matrices(dist, mapped_dist, m1, m2, m3, m4)
1413 :
1414 : TYPE(dbcsr_distribution_type), INTENT(INOUT), &
1415 : OPTIONAL :: dist
1416 : INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL, &
1417 : POINTER :: mapped_dist
1418 : TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: m1, m2, m3, m4
1419 :
1420 : CHARACTER(LEN=*), PARAMETER :: routineN = 'release_dbcsr_topology_and_matrices'
1421 :
1422 : INTEGER :: handle
1423 :
1424 13560 : CALL timeset(routineN, handle)
1425 :
1426 13560 : IF (PRESENT(dist)) CALL dbcsr_distribution_release(dist)
1427 13560 : IF (PRESENT(mapped_dist)) THEN
1428 8960 : IF (ASSOCIATED(mapped_dist)) THEN
1429 8960 : DEALLOCATE (mapped_dist)
1430 : NULLIFY (mapped_dist)
1431 : END IF
1432 : END IF
1433 13560 : IF (PRESENT(m1)) CALL dbcsr_release(m1)
1434 13560 : IF (PRESENT(m2)) CALL dbcsr_release(m2)
1435 13560 : IF (PRESENT(m3)) CALL dbcsr_release(m3)
1436 13560 : IF (PRESENT(m4)) CALL dbcsr_release(m4)
1437 :
1438 13560 : CALL timestop(handle)
1439 :
1440 13560 : END SUBROUTINE release_dbcsr_topology_and_matrices
1441 :
1442 : END MODULE gw_ri_rs_large_cell_Gamma
|
