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