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 :
14 : MODULE gw_non_periodic_ri_rs
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind_set
17 : USE cell_types, ONLY: cell_type
18 : USE constants_operator, ONLY: operator_coulomb
19 : USE cp_dbcsr_api, ONLY: &
20 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_get, &
21 : dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, &
22 : dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_data_size, dbcsr_get_info, &
23 : dbcsr_get_occupation, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
24 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
25 : dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
26 : dbcsr_type_no_symmetry
27 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
28 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
29 : copy_fm_to_dbcsr,&
30 : dbcsr_deallocate_matrix_set,&
31 : max_elements_per_block
32 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
33 : cp_fm_uplo_to_full
34 : USE cp_fm_diag, ONLY: cp_fm_geeig
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_diag,&
37 : cp_fm_get_info,&
38 : cp_fm_release,&
39 : cp_fm_set_all,&
40 : cp_fm_to_fm,&
41 : cp_fm_type
42 : USE gw_compute_Z_lP, ONLY: compute_Z_lP
43 : USE gw_fm_utils, ONLY: fm_contract_ABA,&
44 : fm_invert,&
45 : fm_sqrt
46 : USE gw_large_cell_gamma, ONLY: Fourier_transform_w_to_t,&
47 : G_occ_vir,&
48 : compute_fm_chi_Gamma_freq,&
49 : create_fm_W_MIC_time,&
50 : delete_unnecessary_files,&
51 : fill_fm_Sigma_c_Gamma_time,&
52 : fm_write
53 : USE gw_optimize_ri_rs_grid, ONLY: optimize_ri_rs_grid
54 : USE gw_utils, ONLY: analyt_conti_and_print,&
55 : de_init_bs_env,&
56 : time_to_freq
57 : USE input_constants, ONLY: G0W0,&
58 : evGW0,&
59 : rtp_method_bse
60 : USE kinds, ONLY: default_path_length,&
61 : default_string_length,&
62 : dp,&
63 : int_8,&
64 : max_line_length
65 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
66 : USE machine, ONLY: m_flush,&
67 : m_walltime
68 : USE message_passing, ONLY: mp_mem_avail_per_rank_GB,&
69 : mp_para_env_type,&
70 : mp_print_mem_per_rank
71 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
72 : USE parallel_gemm_api, ONLY: parallel_gemm
73 : USE particle_types, ONLY: particle_type
74 : USE physcon, ONLY: angstrom,&
75 : evolt
76 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type,&
77 : rirs_grid_type
78 : USE post_scf_bandstructure_utils, ONLY: get_all_VBM_CBM_bandgaps
79 : USE qs_environment_types, ONLY: get_qs_env,&
80 : qs_environment_type
81 : USE qs_kind_types, ONLY: get_qs_kind,&
82 : qs_kind_type
83 : USE rirs_grid_utils, ONLY: evaluate_ao_on_points,&
84 : get_rirs_grid_filepath,&
85 : initialize_rirs_grid,&
86 : read_rirs_grid_file
87 : USE util, ONLY: sort
88 : #include "./base/base_uses.f90"
89 :
90 : IMPLICIT NONE
91 :
92 : PRIVATE
93 :
94 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_non_periodic_ri_rs'
95 :
96 : PUBLIC :: gw_calc_non_periodic_ri_rs, ri_rs_grid_assembler, &
97 : get_basis_offsets, precompute_ri_rs_radii, &
98 : atomic_basis_at_grid_point, &
99 : reserve_blocks_within_radius
100 :
101 : CONTAINS
102 :
103 : ! **************************************************************************************************
104 : !> \brief GW calculation using RI-RS formalism for molecules
105 : !> \param qs_env ...
106 : !> \param bs_env ...
107 : ! **************************************************************************************************
108 38 : SUBROUTINE gw_calc_non_periodic_ri_rs(qs_env, bs_env)
109 :
110 : TYPE(qs_environment_type), POINTER :: qs_env
111 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
112 :
113 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_non_periodic_ri_rs'
114 :
115 : INTEGER :: handle
116 38 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_W_time
117 :
118 38 : CALL timeset(routineN, handle)
119 :
120 : ! ==============================================================================
121 : ! 0. Precompute AO and RI radii
122 : ! Per-atom cutoff radii from the most diffuse Gaussian primitives of
123 : ! the AO ("ORB") and RI auxiliary ("RI_AUX") basis sets:
124 : ! α_min,ao = min { ζ_ao | ζ_ao > 10⁻³ }, α_min,ri analogous
125 : ! r_ao = sqrt( -ln(ε) / α_min,ao ) (radius_ao_per_atom)
126 : ! r_ri = sqrt( -ln(ε) / α_min,ri ) (radius_ri_per_atom)
127 : ! ==============================================================================
128 38 : CALL precompute_ri_rs_radii(qs_env, bs_env)
129 :
130 : ! ==============================================================================
131 : ! 1. Grid generation for RI-RS
132 : ! Modified Lebedev atomic grids (Duchemin & Blase), one per atom,
133 : ! concatenated into a flat global list: r_l = R_A + r_l^(A)
134 : ! ==============================================================================
135 38 : CALL ri_rs_grid_assembler(qs_env, bs_env, bs_env%ri_rs%grid_points)
136 :
137 : ! ==============================================================================
138 : ! 2a. Atomic basis evaluation on the grid (grid x AO matrix)
139 : ! Φ_μl = Φ_μ(r_l) (mat_phi_mu_l)
140 : ! ==============================================================================
141 : CALL atomic_basis_at_grid_point(qs_env, bs_env, bs_env%ri_rs%grid_points, &
142 38 : bs_env%ri_rs%mat_phi_mu_l)
143 :
144 : ! ==============================================================================
145 : ! 2b. Print the memory estimate for the RI-RS calculation
146 : ! ==============================================================================
147 38 : CALL print_ri_rs_memory_estimate(qs_env, bs_env)
148 :
149 : ! ==============================================================================
150 : ! 3. RI-RS fitting coefficients Z_lP (grid x RI matrix)
151 : ! Per-atom regularized solve, restricted to grid points r_l within a
152 : ! cutoff distance of atom P:
153 : ! a. D_ll' = [ Σ_μ Φ_μ(r_l) Φ_μ(r_l') ]²
154 : ! b. D_lP = Σ_μν Φ_μ(r_l) Φ_ν(r_l) (μν|P)
155 : ! c. Jacobi conditioning with d_l = 1/sqrt(D_ll):
156 : ! D'_ll' = d_l D_ll' d_l' + λδ_ll' , D'_lP = d_l D_lP
157 : ! d. Solve Σ_l' D'_ll' Z'_l'P = D'_lP
158 : ! e. Rescale Z_lP = d_l Z'_lP (mat_Z_lP)
159 : ! ==============================================================================
160 : CALL compute_Z_lP(qs_env, bs_env, bs_env%ri_rs%grid_points, &
161 38 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
162 :
163 : ! flag the RI-RS grid as built so a subsequent RT-BSE run reuses Z_lP
164 : ! instead of rebuilding it
165 38 : bs_env%ri_rs%grid_built = .TRUE.
166 :
167 : CALL mp_print_mem_per_rank(bs_env%para_env, bs_env%unit_nr, &
168 38 : label='Memory per MPI process after computing Z_lP:')
169 :
170 : ! ==============================================================================
171 : ! 4. Polarizability matrix χ on the imaginary-time grid
172 : ! G^occ_µλ(i|τ|) = Σ_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
173 : ! G^vir_µλ(i|τ|) = Σ_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
174 : ! G^occ_ll'(i|τ|) = Σ_µν Φ_µ(r_l) G^occ_µν Φ_ν(r_l') (G^vir analogous)
175 : ! χ_ll'(iτ) = G^occ_ll'(i|τ|) ∘ G^vir_ll'(i|τ|) (element-wise)
176 : ! χ_PQ(iτ) = Σ_ll' Z_lP χ_ll'(iτ) Z_l'Q
177 : ! ==============================================================================
178 : CALL get_mat_chi_Gamma_tau(bs_env, bs_env%mat_chi_Gamma_tau, &
179 38 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
180 :
181 : ! ==============================================================================
182 : ! 5. Screened Coulomb interaction W (RI basis)
183 : ! χ_PQ(iτ) -> χ_PQ(iω) -> ε_PQ(iω) -> W_PQ(iω) -> W_PQ(iτ)
184 : ! ==============================================================================
185 38 : CALL compute_W(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_time)
186 :
187 : ! ==============================================================================
188 : ! 6. Exact-exchange self-energy Σ^x
189 : ! D_µν = Σ_n^occ C_µn C_νn (density matrix)
190 : ! D_ll' = Σ_µν Φ_µ(r_l) D_µν Φ_ν(r_l')
191 : ! V^tr_ll' = Σ_PQ Z_lP V^tr_PQ Z_l'Q (truncated Coulomb)
192 : ! Σ^x_ll' = D_ll' ∘ V^tr_ll'
193 : ! Σ^x_λσ(k=0) = -Σ_ll' Φ_λ(r_l) Σ^x_ll' Φ_σ(r_l')
194 : ! ==============================================================================
195 : CALL compute_Sigma_x(bs_env, qs_env, bs_env%ri_rs%mat_phi_mu_l, &
196 38 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_x_Gamma)
197 :
198 : ! ==============================================================================
199 : ! 7. Correlation self-energy Σ^c and quasiparticle energies, iterated
200 : ! until eigenvalue self-consistency if the &EVGW0 section is given and
201 : ! done in a single pass for G0W0.
202 : !
203 : ! (a) W_ll'(iτ) = Σ_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
204 : ! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) ∘ W_ll'(iτ), τ < 0
205 : ! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) ∘ W_ll'(iτ), τ > 0
206 : ! Σ^c_λσ(iτ) = Σ_ll' Φ_λ(r_l) Σ^c_ll'(iτ) Φ_σ(r_l')
207 : ! (b) Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
208 : ! ϵ_n^GW,(i) = ϵ_n^DFT + Σ^c_nn[G^(i-1),W](ϵ_n^GW,(i)) + Σ^x_nn - v^xc_nn
209 : ! ==============================================================================
210 38 : CALL compute_Sigma_c_and_QP_energies(bs_env, fm_W_time, fm_Sigma_x_Gamma)
211 :
212 38 : CALL de_init_bs_env(qs_env, bs_env)
213 :
214 38 : CALL timestop(handle)
215 :
216 38 : END SUBROUTINE gw_calc_non_periodic_ri_rs
217 :
218 : ! **************************************************************************************************
219 : !> \brief Correlation self-energy and quasiparticle energies, iterated to eigenvalue
220 : !> self-consistency in G (evGW0).
221 : !>
222 : !> For G0W0 this runs once with G^(0) built from the DFT eigenvalues. For evGW0 the
223 : !> Green's function is rebuilt from the quasiparticle energies of all states until the
224 : !> quasiparticle HOMO, LUMO and HOMO-LUMO gap change by less than EPS_ITER between two
225 : !> cycles, or MAX_ITER cycles are spent. W, Σ^x and everything computed before this
226 : !> routine stay frozen: they either do not depend on the eigenvalues at all (Σ^x is G
227 : !> at τ = 0, a pure density matrix) or are held fixed by construction in GW0. The
228 : !> imaginary-time grid is fixed as well, since W(iτ) lives on it.
229 : !>
230 : !> \param bs_env ...
231 : !> \param fm_W_time ...
232 : !> \param fm_Sigma_x_Gamma ...
233 : ! **************************************************************************************************
234 38 : SUBROUTINE compute_Sigma_c_and_QP_energies(bs_env, fm_W_time, fm_Sigma_x_Gamma)
235 :
236 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
237 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time, fm_Sigma_x_Gamma
238 :
239 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c_and_QP_energies'
240 :
241 : INTEGER :: handle, i_iter, n_iter
242 : LOGICAL :: converged
243 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_scf_Gamma_dft
244 : REAL(KIND=dp), DIMENSION(2) :: e_fermi_dft
245 : REAL(KIND=dp), DIMENSION(3, 2) :: band_prev
246 38 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
247 :
248 38 : CALL timeset(routineN, handle)
249 :
250 : converged = .FALSE.
251 :
252 38 : CALL init_evGW0(bs_env, n_iter, band_prev, eigenval_scf_Gamma_dft, e_fermi_dft)
253 :
254 : ! evGW0 self-consistency loop; for G0W0, loop is terminated after one iteration
255 52 : DO i_iter = 1, n_iter
256 :
257 52 : bs_env%ri_rs%evgw0_i_iter = i_iter
258 :
259 : ! W_ll'(iτ) = Σ_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
260 : ! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) ∘ W_ll'(iτ), τ < 0
261 : ! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) ∘ W_ll'(iτ), τ > 0
262 : ! Σ^c_λσ(iτ) = Σ_ll' Φ_λ(r_l) Σ^c_ll'(iτ) Φ_σ(r_l')
263 : CALL compute_Sigma_c(bs_env, fm_W_time, bs_env%ri_rs%mat_phi_mu_l, &
264 52 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_c_Gamma_time)
265 :
266 : ! Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
267 : ! ϵ_n^GW,(i) = ϵ_n^DFT + Σ^c_nn[G^(i-1),W](ϵ_n^GW,(i)) + Σ^x_nn - v^xc_nn
268 52 : CALL compute_QP_energies(bs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
269 :
270 52 : IF (bs_env%gw_flavour == G0W0) EXIT
271 :
272 18 : CALL print_evGW0_band_edges(bs_env, band_prev, i_iter, n_iter, converged)
273 :
274 18 : IF (i_iter == n_iter .OR. converged) EXIT
275 :
276 : ! eigenvalues ϵ_n to be updated in G:
277 : ! G^occ_µλ(i|τ|) = Σ_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
278 : ! G^vir_µλ(i|τ|) = Σ_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
279 70 : CALL update_eigenvalues_G(bs_env)
280 :
281 : END DO
282 :
283 38 : CALL cp_fm_release(fm_W_time)
284 :
285 38 : CALL reset_and_clean_bs_env(bs_env, eigenval_scf_Gamma_dft, e_fermi_dft, fm_Sigma_x_Gamma)
286 :
287 38 : CALL delete_unnecessary_files(bs_env)
288 :
289 38 : CALL timestop(handle)
290 :
291 38 : END SUBROUTINE compute_Sigma_c_and_QP_energies
292 :
293 : ! **************************************************************************************************
294 : !> \brief Sets up the evGW0 eigenvalue self-consistency loop: the cycle count, the eigenvalues
295 : !> the first Green's function is built from, and the DFT reference that the loop
296 : !> overwrites. A G0W0 run reduces to a single cycle and needs none of it.
297 : !> \param bs_env ...
298 : !> \param n_iter ...
299 : !> \param band_prev ...
300 : !> \param eigenval_scf_Gamma_dft ...
301 : !> \param e_fermi_dft ...
302 : ! **************************************************************************************************
303 38 : SUBROUTINE init_evGW0(bs_env, n_iter, band_prev, eigenval_scf_Gamma_dft, e_fermi_dft)
304 :
305 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
306 : INTEGER, INTENT(OUT) :: n_iter
307 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(OUT) :: band_prev
308 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
309 : INTENT(OUT) :: eigenval_scf_Gamma_dft
310 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: e_fermi_dft
311 :
312 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_evGW0'
313 :
314 : INTEGER :: handle
315 :
316 38 : CALL timeset(routineN, handle)
317 :
318 38 : n_iter = 1
319 38 : band_prev(:, :) = 0.0_dp
320 38 : e_fermi_dft(:) = 0.0_dp
321 :
322 38 : IF (bs_env%gw_flavour == evGW0) THEN
323 4 : n_iter = bs_env%ri_rs%evgw0_iter
324 :
325 : ! eigenvalues currently in G; the first cycle starts is a plain G0W0 step
326 70 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_scf(:, :, :)
327 : ! the loop overwrites these; post-GW printing expects the DFT values back
328 76 : ALLOCATE (eigenval_scf_Gamma_dft, SOURCE=bs_env%eigenval_scf_Gamma)
329 12 : e_fermi_dft(:) = bs_env%e_fermi(:)
330 : END IF
331 :
332 38 : CALL timestop(handle)
333 :
334 38 : END SUBROUTINE init_evGW0
335 :
336 : ! **************************************************************************************************
337 : !> \brief Feeds the quasiparticle energies of the current evGW0 cycle back into the Green's
338 : !> function used by the next one, and re-centres the Fermi level between the new
339 : !> quasiparticle HOMO and LUMO.
340 : !> \param bs_env ...
341 : ! **************************************************************************************************
342 14 : SUBROUTINE update_eigenvalues_G(bs_env)
343 :
344 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
345 :
346 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_eigenvalues_G'
347 :
348 : INTEGER :: handle, i_mo, ispin
349 :
350 14 : CALL timeset(routineN, handle)
351 :
352 : ! record this cycle's evGW0 result; it is also what the next Green's function is built from
353 230 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_GW(:, :, :)
354 :
355 34 : DO ispin = 1, bs_env%n_spin
356 : ! Update all physical states; exclude linear-dependency placeholders.
357 196 : DO i_mo = 1, bs_env%n_mo_retained
358 196 : bs_env%eigenval_scf_Gamma(i_mo, ispin) = bs_env%eigenval_GW(i_mo, 1, ispin)
359 : END DO
360 : bs_env%e_fermi(ispin) = &
361 : 0.5_dp*(bs_env%eigenval_GW(bs_env%n_occ(ispin), 1, ispin) + &
362 34 : bs_env%eigenval_GW(bs_env%n_occ(ispin) + 1, 1, ispin))
363 : END DO
364 :
365 14 : CALL timestop(handle)
366 :
367 14 : END SUBROUTINE update_eigenvalues_G
368 :
369 : ! **************************************************************************************************
370 : !> \brief Restores the DFT reference that the evGW0 loop overwrote and cleanup
371 : !> \param bs_env ...
372 : !> \param eigenval_scf_Gamma_dft ...
373 : !> \param e_fermi_dft ...
374 : !> \param fm_Sigma_x_Gamma ..
375 : ! **************************************************************************************************
376 38 : SUBROUTINE reset_and_clean_bs_env(bs_env, eigenval_scf_Gamma_dft, e_fermi_dft, fm_Sigma_x_Gamma)
377 :
378 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
379 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
380 : INTENT(INOUT) :: eigenval_scf_Gamma_dft
381 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: e_fermi_dft
382 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
383 :
384 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_and_clean_bs_env'
385 :
386 : INTEGER :: handle
387 :
388 38 : CALL timeset(routineN, handle)
389 :
390 38 : IF (bs_env%gw_flavour == evGW0) THEN
391 : ! Update array with final evGW0 result
392 70 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_GW(:, :, :)
393 : ! put the DFT reference back for the post-GW DOS/band-edge printing
394 64 : bs_env%eigenval_scf_Gamma(:, :) = eigenval_scf_Gamma_dft(:, :)
395 12 : bs_env%e_fermi(:) = e_fermi_dft(:)
396 4 : DEALLOCATE (eigenval_scf_Gamma_dft)
397 4 : CALL cp_fm_release(fm_Sigma_x_Gamma)
398 : END IF
399 :
400 38 : CALL timestop(handle)
401 :
402 38 : END SUBROUTINE reset_and_clean_bs_env
403 :
404 : ! **************************************************************************************************
405 : !> \brief Print the quasiparticle HOMO, LUMO and HOMO-LUMO gap of the current evGW0 cycle
406 : !> \param bs_env ...
407 : !> \param band_prev ...
408 : !> \param i_iter ...
409 : !> \param n_iter ...
410 : !> \param converged ...
411 : ! **************************************************************************************************
412 18 : SUBROUTINE print_evGW0_band_edges(bs_env, band_prev, i_iter, n_iter, converged)
413 :
414 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
415 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(INOUT) :: band_prev
416 : INTEGER, INTENT(IN) :: i_iter, n_iter
417 : LOGICAL, INTENT(OUT) :: converged
418 :
419 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_evGW0_band_edges'
420 :
421 : INTEGER :: handle, homo, ispin, u
422 : REAL(KIND=dp) :: max_delta
423 : REAL(KIND=dp), DIMENSION(3) :: band
424 :
425 18 : CALL timeset(routineN, handle)
426 :
427 18 : u = bs_env%unit_nr
428 :
429 18 : converged = (i_iter > 1)
430 18 : max_delta = 0.0_dp
431 :
432 18 : IF (u > 0) THEN
433 9 : WRITE (u, '(A)') ' '
434 9 : WRITE (u, '(T2,A)') REPEAT('-', 79)
435 9 : WRITE (u, '(T2,A,I4,A,I4)') 'evGW0 cycle', i_iter, ' /', n_iter
436 9 : WRITE (u, '(T2,A)') REPEAT('-', 79)
437 : END IF
438 :
439 44 : DO ispin = 1, bs_env%n_spin
440 :
441 26 : homo = bs_env%n_occ(ispin)
442 26 : band(1) = bs_env%eigenval_GW(homo, 1, ispin)
443 26 : band(2) = bs_env%eigenval_GW(homo + 1, 1, ispin)
444 26 : band(3) = band(2) - band(1)
445 :
446 26 : IF (i_iter > 1) THEN
447 80 : max_delta = MAX(max_delta, MAXVAL(ABS(band(:) - band_prev(:, ispin))))
448 : END IF
449 :
450 26 : IF (u > 0) THEN
451 13 : IF (bs_env%n_spin == 2) WRITE (u, '(T2,A,I0)') 'Spin ', ispin
452 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 HOMO (eV)', band(1)*evolt
453 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 LUMO (eV)', band(2)*evolt
454 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 HOMO-LUMO gap (eV)', band(3)*evolt
455 : END IF
456 :
457 122 : band_prev(:, ispin) = band(:)
458 :
459 : END DO
460 :
461 18 : IF (i_iter > 1) THEN
462 14 : converged = (max_delta < bs_env%ri_rs%evgw0_eps_iter)
463 21 : IF (u > 0) WRITE (u, '(T2,A,T61,F20.6)') 'Max. change to previous cycle (eV)', &
464 14 : max_delta*evolt
465 : END IF
466 :
467 11 : IF (u > 0) WRITE (u, '(T2,A)') REPEAT('-', 79)
468 :
469 18 : IF (converged) THEN
470 4 : IF (u > 0) THEN
471 2 : WRITE (u, '(A)') ' '
472 : WRITE (u, '(T2,A,I4,A)') &
473 2 : 'evGW0 eigenvalue self-consistency reached in', i_iter, ' cycles.'
474 2 : WRITE (u, '(A)') ' '
475 : END IF
476 14 : ELSE IF (i_iter == n_iter) THEN
477 : CALL cp_warn(__LOCATION__, &
478 : "The evGW0 eigenvalue self-consistency cycle did not converge "// &
479 : "within MAX_ITER cycles. The reported quasiparticle energies are "// &
480 0 : "those of the last cycle.")
481 : END IF
482 :
483 18 : CALL timestop(handle)
484 :
485 18 : END SUBROUTINE print_evGW0_band_edges
486 :
487 : ! **************************************************************************************************
488 : !> \brief Compute per-atom AO and RI basis radii from the most diffuse Gaussian
489 : !> primitive in the AO ("ORB") and RI auxiliary ("RI_AUX") basis sets.
490 : !> Stores results in bs_env%ri_rs%radius_ao_per_atom(:) and
491 : !> bs_env%ri_rs%radius_ri_per_atom(:) and prints a per-atom table.
492 : !> Radius: r_kind = sqrt(-log(eps) / alpha_min_kind)
493 : !> with eps = eps_filter.
494 : !> \param qs_env ...
495 : !> \param bs_env ...
496 : ! **************************************************************************************************
497 46 : SUBROUTINE precompute_ri_rs_radii(qs_env, bs_env)
498 :
499 : TYPE(qs_environment_type), POINTER :: qs_env
500 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
501 :
502 : CHARACTER(LEN=*), PARAMETER :: routineN = 'precompute_ri_rs_radii'
503 :
504 : INTEGER :: handle, i, iatom, ikind, j, natom, nkind
505 46 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
506 : REAL(KIND=dp) :: eps
507 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha_min_ao_kind, alpha_min_ri_kind
508 46 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet_ao, zet_ri
509 46 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
510 46 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
511 :
512 46 : CALL timeset(routineN, handle)
513 :
514 : CALL get_qs_env(qs_env, nkind=nkind, atomic_kind_set=atomic_kind_set, &
515 46 : particle_set=particle_set)
516 46 : natom = bs_env%n_atom
517 :
518 46 : eps = bs_env%eps_filter
519 :
520 184 : ALLOCATE (alpha_min_ao_kind(nkind), alpha_min_ri_kind(nkind))
521 112 : alpha_min_ao_kind = HUGE(1.0_dp)
522 112 : alpha_min_ri_kind = HUGE(1.0_dp)
523 :
524 112 : DO ikind = 1, nkind
525 66 : zet_ao => bs_env%basis_set_AO(ikind)%gto_basis_set%zet
526 66 : zet_ri => bs_env%basis_set_RI(ikind)%gto_basis_set%zet
527 :
528 220 : DO i = 1, SIZE(zet_ao, 1)
529 512 : DO j = 1, SIZE(zet_ao, 2)
530 446 : IF (zet_ao(i, j) > 1.0E-3_dp) THEN
531 292 : alpha_min_ao_kind(ikind) = MIN(alpha_min_ao_kind(ikind), zet_ao(i, j))
532 : END IF
533 : END DO
534 : END DO
535 178 : DO i = 1, SIZE(zet_ri, 1)
536 580 : DO j = 1, SIZE(zet_ri, 2)
537 514 : IF (zet_ri(i, j) > 1.0E-3_dp) THEN
538 448 : alpha_min_ri_kind(ikind) = MIN(alpha_min_ri_kind(ikind), zet_ri(i, j))
539 : END IF
540 : END DO
541 : END DO
542 : END DO
543 :
544 46 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
545 :
546 138 : ALLOCATE (bs_env%ri_rs%radius_ao_per_atom(natom))
547 92 : ALLOCATE (bs_env%ri_rs%radius_ri_per_atom(natom))
548 158 : DO iatom = 1, natom
549 112 : ikind = kind_of(iatom)
550 112 : bs_env%ri_rs%radius_ao_per_atom(iatom) = SQRT(-LOG(eps)/alpha_min_ao_kind(ikind))
551 158 : bs_env%ri_rs%radius_ri_per_atom(iatom) = SQRT(-LOG(eps)/alpha_min_ri_kind(ikind))
552 : END DO
553 :
554 46 : IF (bs_env%unit_nr > 0) THEN
555 23 : WRITE (bs_env%unit_nr, '(T2,A)') 'RI-RS basis radii (Å):'
556 23 : WRITE (bs_env%unit_nr, '(T4,A6,2X,A4,2A14)') 'Kind', 'Elem', 'r_AO (Å)', 'r_RI (Å)'
557 56 : DO ikind = 1, nkind
558 : WRITE (bs_env%unit_nr, '(T4,I6,2X,A4,2F14.4)') &
559 33 : ikind, &
560 33 : atomic_kind_set(ikind)%element_symbol, &
561 33 : SQRT(-LOG(eps)/alpha_min_ao_kind(ikind))*angstrom, &
562 89 : SQRT(-LOG(eps)/alpha_min_ri_kind(ikind))*angstrom
563 : END DO
564 23 : WRITE (bs_env%unit_nr, '(A)') ' '
565 : END IF
566 :
567 46 : DEALLOCATE (alpha_min_ao_kind, alpha_min_ri_kind, kind_of)
568 :
569 46 : CALL timestop(handle)
570 :
571 46 : END SUBROUTINE precompute_ri_rs_radii
572 :
573 : ! **************************************************************************************************
574 : !> \brief Spreads the low 21 bits of a into every third bit (bits 0,3,6,...,60): the 1-D helper
575 : !> for a 3-D Morton (Z-order) code. Standard 64-bit magic-mask implementation.
576 : !> \param a value in [0, 2^21)
577 : !> \param x a with two zero bits inserted between consecutive input bits
578 : ! **************************************************************************************************
579 336 : SUBROUTINE morton_split3(a, x)
580 : INTEGER(KIND=int_8), INTENT(IN) :: a
581 : INTEGER(KIND=int_8), INTENT(OUT) :: x
582 :
583 336 : x = IAND(a, INT(z'1FFFFF', int_8))
584 336 : x = IAND(IOR(x, ISHFT(x, 32)), INT(z'1F00000000FFFF', int_8))
585 336 : x = IAND(IOR(x, ISHFT(x, 16)), INT(z'1F0000FF0000FF', int_8))
586 336 : x = IAND(IOR(x, ISHFT(x, 8)), INT(z'100F00F00F00F00F', int_8))
587 336 : x = IAND(IOR(x, ISHFT(x, 4)), INT(z'10C30C30C30C30C3', int_8))
588 336 : x = IAND(IOR(x, ISHFT(x, 2)), INT(z'1249249249249249', int_8))
589 336 : END SUBROUTINE morton_split3
590 :
591 : ! **************************************************************************************************
592 : !> \brief Returns a permutation of atom indices in Morton (Z-order) space-filling order of their
593 : !> Cartesian centers, so consecutive atoms are spatial neighbors. The RI-RS grid rows are
594 : !> laid down in this order, so a contiguous grid panel maps to a compact spatial region and
595 : !> the CUTOFF_RADIUS_RL_W neighborhood of every panel shrinks. The grid row index
596 : !> is a summed contraction index, so ANY permutation is result-preserving; this one is
597 : !> chosen purely to improve locality. Coordinates are normalized to the atom bounding box
598 : !> and quantized to 21 bits per axis (sub-picometre for any real cell).
599 : !> \param particle_set ...
600 : !> \param order order(i) = atom index placed at layout position i
601 : ! **************************************************************************************************
602 46 : SUBROUTINE spatial_atom_order(particle_set, order)
603 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
604 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: order
605 :
606 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spatial_atom_order'
607 : INTEGER, PARAMETER :: nbits = 21
608 :
609 : INTEGER :: handle, ia, k, natom
610 : INTEGER(KIND=int_8) :: cmax, ic(3), m1, m2, m3
611 46 : INTEGER(KIND=int_8), ALLOCATABLE :: mcode(:)
612 : REAL(KIND=dp) :: hi(3), lo(3), span(3)
613 :
614 46 : CALL timeset(routineN, handle)
615 :
616 46 : natom = SIZE(particle_set)
617 230 : ALLOCATE (order(natom), mcode(natom))
618 : cmax = ISHFT(1_int_8, nbits) - 1_int_8
619 :
620 184 : lo(:) = HUGE(1.0_dp)
621 184 : hi(:) = -HUGE(1.0_dp)
622 158 : DO ia = 1, natom
623 494 : DO k = 1, 3
624 336 : lo(k) = MIN(lo(k), particle_set(ia)%r(k))
625 448 : hi(k) = MAX(hi(k), particle_set(ia)%r(k))
626 : END DO
627 : END DO
628 184 : span(:) = hi(:) - lo(:)
629 184 : DO k = 1, 3
630 184 : IF (span(k) <= 0.0_dp) span(k) = 1.0_dp
631 : END DO
632 :
633 158 : DO ia = 1, natom
634 448 : DO k = 1, 3
635 336 : ic(k) = INT(((particle_set(ia)%r(k) - lo(k))/span(k))*REAL(cmax, dp), int_8)
636 448 : ic(k) = MIN(cmax, MAX(0_int_8, ic(k)))
637 : END DO
638 112 : CALL morton_split3(ic(1), m1)
639 112 : CALL morton_split3(ic(2), m2)
640 112 : CALL morton_split3(ic(3), m3)
641 158 : mcode(ia) = IOR(IOR(m1, ISHFT(m2, 1)), ISHFT(m3, 2))
642 : END DO
643 :
644 : ! sort(mcode, natom, order): order(i) = original atom index with the i-th smallest code
645 46 : CALL sort(mcode, natom, order)
646 :
647 46 : DEALLOCATE (mcode)
648 :
649 46 : CALL timestop(handle)
650 :
651 46 : END SUBROUTINE spatial_atom_order
652 :
653 : ! **************************************************************************************************
654 : !> \brief Compute grid points for RI-RS, two options:
655 : !> 1) Read tabulated atomic grids by Duchemin & Blase:
656 : !> JCP 150, 174120 (2019), JCTC 17, 2383 (2021)
657 : !> 2) Optimize grid for specific molecule / unit cell, start from tabulated grid
658 : !> \param qs_env ...
659 : !> \param bs_env ...
660 : !> \param ri_rs_grid_points ...
661 : ! **************************************************************************************************
662 46 : SUBROUTINE ri_rs_grid_assembler(qs_env, bs_env, ri_rs_grid_points)
663 :
664 : TYPE(qs_environment_type), POINTER :: qs_env
665 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
666 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: ri_rs_grid_points(:, :)
667 :
668 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ri_rs_grid_assembler'
669 :
670 : INTEGER :: handle
671 46 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
672 46 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
673 :
674 46 : CALL timeset(routineN, handle)
675 :
676 : !! Get the information about the atoms in the system
677 46 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, particle_set=particle_set)
678 :
679 : !! 1. Read the tabulated atomic grids into cache
680 46 : CALL build_grid_cache(bs_env, atomic_kind_set)
681 :
682 : !! 2. Optimize RIRS grid points
683 46 : IF (bs_env%ri_rs%grid_opt%enabled) THEN
684 8 : CALL optimize_ri_rs_grid(qs_env, bs_env)
685 : END IF
686 :
687 : !! 3. Assemble and store the global molecular grid
688 46 : CALL assemble_RI_RS_grid_points(bs_env, particle_set, ri_rs_grid_points)
689 :
690 : !! 4. Release the atom-relative grid cache
691 46 : CALL release_RI_RS_grid_cache(bs_env)
692 :
693 46 : CALL timestop(handle)
694 :
695 46 : END SUBROUTINE ri_rs_grid_assembler
696 :
697 : ! **************************************************************************************************
698 : !> \brief Assemble the global RI-RS grid in spatial atom order from atom-relative cached grids.
699 : !> For point l on atom A, r_l = R_A + Δr_l^(A).
700 : !> \param bs_env GW environment containing the cached grids and receiving the atom boundaries.
701 : !> \param particle_set Atomic positions R_A.
702 : !> \param ri_rs_grid_points Global Cartesian grid coordinates in spatial atom order.
703 : ! **************************************************************************************************
704 46 : SUBROUTINE assemble_RI_RS_grid_points(bs_env, particle_set, ri_rs_grid_points)
705 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
706 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
707 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: ri_rs_grid_points(:, :)
708 :
709 : INTEGER :: atom, atom_grid_end, atom_grid_start, &
710 : layout_position, natom, total_grid_npts
711 46 : INTEGER, ALLOCATABLE :: atom_grid_offsets(:), atom_order(:)
712 : REAL(KIND=dp) :: atom_center(3)
713 :
714 46 : natom = bs_env%n_atom
715 46 : CALL spatial_atom_order(particle_set, atom_order)
716 :
717 : ! grid_atom_boundaries(k) starts the k-th atom's run in spatial layout order.
718 : ! atom_grid_offsets(A) locates atom A in the resulting global coordinate array.
719 230 : ALLOCATE (bs_env%ri_rs%grid_atom_boundaries(natom + 1), atom_grid_offsets(natom))
720 46 : total_grid_npts = 0
721 158 : DO layout_position = 1, natom
722 112 : atom = atom_order(layout_position)
723 112 : atom_grid_offsets(atom) = total_grid_npts + 1
724 112 : bs_env%ri_rs%grid_atom_boundaries(layout_position) = total_grid_npts + 1
725 158 : total_grid_npts = total_grid_npts + bs_env%ri_rs%grid_cache(atom)%npts
726 : END DO
727 46 : bs_env%ri_rs%grid_atom_boundaries(natom + 1) = total_grid_npts + 1
728 :
729 46 : IF (bs_env%unit_nr > 0) THEN
730 : WRITE (bs_env%unit_nr, FMT="(T2,A,T69,I12)") &
731 23 : 'Total grid points used for RI-RS:', total_grid_npts
732 23 : WRITE (bs_env%unit_nr, "(A)") ' '
733 : END IF
734 :
735 138 : ALLOCATE (ri_rs_grid_points(3, total_grid_npts))
736 : !$OMP PARALLEL DO DEFAULT(NONE) &
737 : !$OMP SHARED(ri_rs_grid_points, atom_grid_offsets, particle_set, bs_env, natom) &
738 : !$OMP PRIVATE(atom, atom_center, atom_grid_start, atom_grid_end) &
739 46 : !$OMP SCHEDULE(DYNAMIC, 1)
740 : DO atom = 1, natom
741 : atom_center(:) = particle_set(atom)%r(:)
742 : atom_grid_start = atom_grid_offsets(atom)
743 : atom_grid_end = atom_grid_start + bs_env%ri_rs%grid_cache(atom)%npts - 1
744 :
745 : ! Translate the cached atom-relative points Δr_l^(A) to r_l = R_A + Δr_l^(A).
746 : ri_rs_grid_points(1, atom_grid_start:atom_grid_end) = &
747 : bs_env%ri_rs%grid_cache(atom)%raw_points(1, :) + atom_center(1)
748 : ri_rs_grid_points(2, atom_grid_start:atom_grid_end) = &
749 : bs_env%ri_rs%grid_cache(atom)%raw_points(2, :) + atom_center(2)
750 : ri_rs_grid_points(3, atom_grid_start:atom_grid_end) = &
751 : bs_env%ri_rs%grid_cache(atom)%raw_points(3, :) + atom_center(3)
752 : END DO
753 : !$OMP END PARALLEL DO
754 46 : END SUBROUTINE assemble_RI_RS_grid_points
755 :
756 : ! **************************************************************************************************
757 : !> \brief Release the atom-relative RI-RS grid cache after assembling the global molecular grid.
758 : !> \param bs_env GW environment containing the grid cache.
759 : ! **************************************************************************************************
760 46 : SUBROUTINE release_RI_RS_grid_cache(bs_env)
761 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
762 :
763 : INTEGER :: atom
764 :
765 158 : DO atom = 1, SIZE(bs_env%ri_rs%grid_cache)
766 158 : DEALLOCATE (bs_env%ri_rs%grid_cache(atom)%raw_points)
767 : END DO
768 158 : DEALLOCATE (bs_env%ri_rs%grid_cache)
769 46 : END SUBROUTINE release_RI_RS_grid_cache
770 :
771 : ! **************************************************************************************************
772 : !> \brief Reads grids from .ion files and stores them in memory based on grid_select
773 : !> \param bs_env ...
774 : !> \param atomic_kind_set ...
775 : ! **************************************************************************************************
776 46 : SUBROUTINE build_grid_cache(bs_env, atomic_kind_set)
777 :
778 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
779 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
780 :
781 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_grid_cache'
782 :
783 : CHARACTER(LEN=default_path_length) :: full_path, warning_msg
784 : CHARACTER(LEN=default_string_length) :: atom_sym
785 : INTEGER :: handle, iatom, ikind, &
786 : n_voronoi_candidates, natom, nkind, &
787 : requested_npts
788 46 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
789 46 : TYPE(rirs_grid_type), ALLOCATABLE :: kind_grid(:)
790 :
791 46 : CALL timeset(routineN, handle)
792 :
793 46 : nkind = SIZE(atomic_kind_set)
794 46 : particle_set => bs_env%ri_rs%particle_set
795 46 : CPASSERT(ASSOCIATED(particle_set))
796 46 : natom = bs_env%n_atom
797 204 : ALLOCATE (kind_grid(nkind))
798 :
799 112 : DO ikind = 1, nkind
800 66 : atom_sym = TRIM(atomic_kind_set(ikind)%element_symbol)
801 : CALL get_rirs_grid_filepath(atom_sym, bs_env%ri_rs%grid_select, &
802 66 : bs_env%ri_rs%grid_file_suffix, full_path)
803 66 : CALL read_rirs_grid_file(TRIM(full_path), kind_grid(ikind)%raw_points)
804 112 : kind_grid(ikind)%npts = SIZE(kind_grid(ikind)%raw_points, 2)
805 : END DO
806 :
807 : ! Make every atom own an independent copy. The optimizer can therefore adapt symmetry-
808 : ! equivalent atoms to their individual local environments, and GW consumes those copies
809 : ! directly without writing intermediate .ion files.
810 46 : IF (ALLOCATED(bs_env%ri_rs%grid_cache)) DEALLOCATE (bs_env%ri_rs%grid_cache)
811 250 : ALLOCATE (bs_env%ri_rs%grid_cache(natom))
812 158 : DO iatom = 1, natom
813 112 : ikind = particle_set(iatom)%atomic_kind%kind_number
814 112 : bs_env%ri_rs%grid_cache(iatom) = kind_grid(ikind)
815 158 : IF (bs_env%ri_rs%grid_opt%enabled) THEN
816 24 : atom_sym = TRIM(atomic_kind_set(ikind)%element_symbol)
817 24 : requested_npts = bs_env%ri_rs%grid_opt%grid_npoints(ikind)
818 : CALL initialize_rirs_grid(bs_env%ri_rs%grid_cache(iatom)%raw_points, &
819 : requested_npts, iatom, particle_set, &
820 24 : n_voronoi_candidates)
821 24 : IF (n_voronoi_candidates < requested_npts) THEN
822 : WRITE (warning_msg, '(A,I0,A,A,A,I0,A,I0,A)') &
823 0 : "Atom ", iatom, " (", TRIM(atom_sym), "): Grid selection requests ", requested_npts, &
824 0 : " points, but only ", n_voronoi_candidates, &
825 0 : " are inside its Voronoi volume. Selecting from the full RI-RS source grid."
826 0 : CALL cp_warn(__LOCATION__, TRIM(warning_msg))
827 : END IF
828 24 : bs_env%ri_rs%grid_cache(iatom)%npts = requested_npts
829 : END IF
830 : END DO
831 112 : DEALLOCATE (kind_grid)
832 :
833 46 : CALL timestop(handle)
834 :
835 46 : END SUBROUTINE build_grid_cache
836 :
837 : ! **************************************************************************************************
838 : !> \brief Evaluates the AO basis on the RI-RS grid and stores it as the sparse DBCSR matrix
839 : !> Φ_μl = Φ_μ(r_l) (rows = grid points in atom-aligned blocks of at most
840 : !> max_elements_per_block points, columns = one block per atom's full AO set).
841 : !> Grid points outside the reach of an atom's most
842 : !> diffuse Gaussian (or the CUTOFF_RADIUS_RL_AO) are skipped, and only blocks
843 : !> with at least one element > eps_filter are stored. This locality is the source of
844 : !> ALL grid-dimension sparsity used downstream. Also caches the atom centers and the
845 : !> per-chunk centroids needed by the optional CUTOFF_RADIUS_G_W / CUTOFF_RADIUS_RL_W
846 : !> operator truncations.
847 : !> \param qs_env ...
848 : !> \param bs_env ...
849 : !> \param ri_rs_grid_points ...
850 : !> \param mat_phi_mu_l ...
851 : ! **************************************************************************************************
852 46 : SUBROUTINE atomic_basis_at_grid_point(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l)
853 :
854 : TYPE(qs_environment_type), POINTER :: qs_env
855 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
856 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
857 : TYPE(dbcsr_type), INTENT(OUT) :: mat_phi_mu_l
858 :
859 : CHARACTER(LEN=*), PARAMETER :: routineN = 'atomic_basis_at_grid_point'
860 :
861 : INTEGER :: bs_eff, c_size, dimen_ORB, handle, i, i_blk, ia, iatom, natom, npcol, nprow, &
862 : num_grid_chunks, r_end, r_start, remaining, run, safe_max, total_grid_npts
863 46 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_row_start, first_sgf
864 46 : INTEGER, DIMENSION(:), POINTER :: c_blk_sizes, col_dist, r_blk_sizes, &
865 46 : row_dist
866 : REAL(KIND=dp) :: r2_threshold
867 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atom_col_buffer
868 46 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
869 : TYPE(cell_type), POINTER :: cell
870 : TYPE(dbcsr_distribution_type) :: dist
871 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist_ks
872 : TYPE(mp_para_env_type), POINTER :: para_env
873 46 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
874 46 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
875 :
876 46 : CALL timeset(routineN, handle)
877 :
878 : ! Extract environment variables
879 : CALL get_qs_env(qs_env, cell=cell, atomic_kind_set=atomic_kind_set, &
880 : qs_kind_set=qs_kind_set, particle_set=particle_set, &
881 46 : para_env=para_env)
882 :
883 46 : natom = bs_env%n_atom
884 46 : total_grid_npts = SIZE(ri_rs_grid_points, 2)
885 :
886 : ! Map the starting indices of spherical gaussian functions (SGF) for each atom
887 138 : ALLOCATE (first_sgf(natom + 1))
888 46 : CALL get_basis_offsets(particle_set, qs_kind_set, first_sgf, dimen_ORB)
889 :
890 : ! =========================================================================
891 : ! 1. SETUP DBCSR MATRIX TOPOLOGY
892 : ! =========================================================================
893 :
894 : ! A. Define Column Block Sizes (1 Block = 1 Atom's full basis set)
895 138 : ALLOCATE (c_blk_sizes(natom))
896 158 : DO iatom = 1, natom
897 158 : c_blk_sizes(iatom) = first_sgf(iatom + 1) - first_sgf(iatom)
898 : END DO
899 :
900 : ! B. Define Row Block Sizes: atom-aligned blocks (a block never spans two atoms' grid
901 : ! runs), each atom's run subdivided into blocks of at most bs_eff points.
902 :
903 : ! Fetch CP2K's default process grid configuration
904 46 : CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist_ks)
905 46 : CALL dbcsr_distribution_get(dbcsr_dist_ks, nprows=nprow, npcols=npcol)
906 :
907 : ! Overflow-safe upper bound on the block size (see bs_env%dbcsr_msg_elem_limit).
908 : safe_max = INT(0.5_dp*REAL(bs_env%dbcsr_msg_elem_limit, dp)* &
909 : REAL(MAX(MIN(nprow, npcol), 1), dp)/ &
910 46 : REAL(total_grid_npts, dp))
911 46 : safe_max = MAX(1, safe_max)
912 : ! Block size = CP2K's global max_elements_per_block (GLOBAL/DBCSR input; default 32),
913 : ! overflow-capped.
914 46 : bs_eff = MAX(1, MIN(max_elements_per_block, safe_max))
915 :
916 : ! Count the atom-aligned blocks, then fill r_blk_sizes and each block's starting grid row.
917 46 : num_grid_chunks = 0
918 158 : DO ia = 1, natom
919 112 : run = bs_env%ri_rs%grid_atom_boundaries(ia + 1) - bs_env%ri_rs%grid_atom_boundaries(ia)
920 158 : IF (run > 0) num_grid_chunks = num_grid_chunks + (run + bs_eff - 1)/bs_eff
921 : END DO
922 184 : ALLOCATE (r_blk_sizes(num_grid_chunks), blk_row_start(num_grid_chunks))
923 46 : i_blk = 0
924 46 : r_start = 1
925 158 : DO ia = 1, natom
926 112 : remaining = bs_env%ri_rs%grid_atom_boundaries(ia + 1) - bs_env%ri_rs%grid_atom_boundaries(ia)
927 760 : DO WHILE (remaining > 0)
928 602 : i_blk = i_blk + 1
929 602 : r_blk_sizes(i_blk) = MIN(bs_eff, remaining)
930 602 : blk_row_start(i_blk) = r_start
931 602 : r_start = r_start + r_blk_sizes(i_blk)
932 714 : remaining = remaining - r_blk_sizes(i_blk)
933 : END DO
934 : END DO
935 :
936 46 : IF (bs_env%unit_nr > 0) THEN
937 : ! T71 compensates for the two two-byte Greek characters in the label.
938 : WRITE (bs_env%unit_nr, '(T2,A,T71,I12)') &
939 23 : 'RI-RS grid row-blocks of Φ_μ(r_l)', num_grid_chunks
940 23 : WRITE (bs_env%unit_nr, '(T2,A,T69,I12)') 'RI-RS grid points per block (max)', bs_eff
941 : END IF
942 :
943 : ! Cache atomic positions: AO and RI blocks are one-block-per-atom, so these are the
944 : ! block centers used by the optional CUTOFF_RADIUS_G_W atom-pair truncation.
945 46 : IF (ALLOCATED(bs_env%ri_rs%atom_centers)) DEALLOCATE (bs_env%ri_rs%atom_centers)
946 138 : ALLOCATE (bs_env%ri_rs%atom_centers(3, natom))
947 158 : DO iatom = 1, natom
948 942 : bs_env%ri_rs%atom_centers(1:3, iatom) = particle_set(iatom)%r(1:3)
949 : END DO
950 :
951 : ! Cache per-chunk centroids for the optional CUTOFF_RADIUS_RL_W / CUTOFF_RADIUS_W0 block
952 : ! truncations. Left unallocated otherwise, so PRESENT(centroids) stays .FALSE. at the
953 : ! contract_grid_panels* calls.
954 46 : IF (bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp .OR. &
955 : bs_env%ri_rs%cutoff_radius_w0 > 0.0_dp) THEN
956 8 : IF (ALLOCATED(bs_env%ri_rs%chunk_centroids)) DEALLOCATE (bs_env%ri_rs%chunk_centroids)
957 24 : ALLOCATE (bs_env%ri_rs%chunk_centroids(3, num_grid_chunks))
958 110 : DO i_blk = 1, num_grid_chunks
959 102 : r_start = blk_row_start(i_blk)
960 102 : r_end = r_start + r_blk_sizes(i_blk) - 1
961 : bs_env%ri_rs%chunk_centroids(1, i_blk) = &
962 2950 : SUM(ri_rs_grid_points(1, r_start:r_end))/REAL(r_blk_sizes(i_blk), dp)
963 : bs_env%ri_rs%chunk_centroids(2, i_blk) = &
964 2950 : SUM(ri_rs_grid_points(2, r_start:r_end))/REAL(r_blk_sizes(i_blk), dp)
965 : bs_env%ri_rs%chunk_centroids(3, i_blk) = &
966 2996 : SUM(ri_rs_grid_points(3, r_start:r_end))/REAL(r_blk_sizes(i_blk), dp)
967 : END DO
968 : END IF
969 :
970 : ! C. Build Custom Mappings using Round-Robin across the 2D process grid
971 :
972 92 : ALLOCATE (row_dist(num_grid_chunks))
973 648 : DO i = 1, num_grid_chunks
974 648 : row_dist(i) = MOD(i - 1, nprow)
975 : END DO
976 :
977 92 : ALLOCATE (col_dist(natom))
978 158 : DO i = 1, natom
979 158 : col_dist(i) = MOD(i - 1, npcol)
980 : END DO
981 :
982 : ! E. Create the DBCSR Distribution and Initialize the Matrix
983 : CALL dbcsr_distribution_new(dist, template=dbcsr_dist_ks, &
984 46 : row_dist=row_dist, col_dist=col_dist)
985 :
986 : CALL dbcsr_create(mat_phi_mu_l, name="phi_val_sparse", dist=dist, &
987 : matrix_type=dbcsr_type_no_symmetry, &
988 46 : row_blk_size=r_blk_sizes, col_blk_size=c_blk_sizes)
989 :
990 : ! =========================================================================
991 : ! 2. STREAM DATA DIRECTLY INTO SPARSE MATRIX
992 : ! =========================================================================
993 : ! Iterate over the atoms assigned to this specific MPI rank
994 46 : DO iatom = para_env%mepos + 1, natom, para_env%num_pe
995 :
996 56 : c_size = c_blk_sizes(iatom)
997 :
998 : ! Allocate a temporary dense buffer just for this specific atom
999 224 : ALLOCATE (atom_col_buffer(total_grid_npts, c_size))
1000 56 : atom_col_buffer = 0.0_dp
1001 :
1002 : ! Evaluate the basis functions on the grid. Skip grid points outside
1003 : ! the spatial extent of the most diffuse AO Gaussian on iatom; beyond
1004 : ! that radius the contribution is guaranteed below eps_filter. A positive
1005 : ! CUTOFF_RADIUS_RL_AO overrides this with a user-defined hard cutoff.
1006 56 : IF (bs_env%ri_rs%cutoff_radius_ri_ao > 0.0_dp) THEN
1007 3 : r2_threshold = bs_env%ri_rs%cutoff_radius_ri_ao**2
1008 : ELSE
1009 53 : r2_threshold = bs_env%ri_rs%radius_ao_per_atom(iatom)**2
1010 : END IF
1011 : CALL fill_phi_for_atom(atom_col_buffer, ri_rs_grid_points, total_grid_npts, &
1012 56 : iatom, particle_set, qs_kind_set, cell, r2_threshold)
1013 :
1014 : ! Slice the dense column into the atom-aligned grid row-blocks and insert into DBCSR
1015 801 : DO i_blk = 1, num_grid_chunks
1016 745 : r_start = blk_row_start(i_blk)
1017 745 : r_end = r_start + r_blk_sizes(i_blk) - 1
1018 :
1019 : ! Apply dynamic sparsity filtering: Only store blocks with physical significance
1020 82176 : IF (MAXVAL(ABS(atom_col_buffer(r_start:r_end, 1:c_size))) > bs_env%eps_filter) THEN
1021 : CALL dbcsr_put_block(mat_phi_mu_l, row=i_blk, col=iatom, &
1022 733 : block=atom_col_buffer(r_start:r_end, 1:c_size))
1023 : END IF
1024 : END DO
1025 :
1026 56 : DEALLOCATE (atom_col_buffer)
1027 :
1028 : END DO
1029 :
1030 46 : CALL dbcsr_finalize(mat_phi_mu_l)
1031 :
1032 46 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1033 46 : CALL print_matrix_occupation(mat_phi_mu_l, 'Φ_μ(r_l)', bs_env)
1034 :
1035 : ! -------------------------------------------------------------------------
1036 : ! CLEANUP
1037 : ! -------------------------------------------------------------------------
1038 46 : DEALLOCATE (first_sgf, r_blk_sizes, c_blk_sizes, row_dist, col_dist, blk_row_start)
1039 46 : CALL dbcsr_distribution_release(dist)
1040 :
1041 46 : CALL timestop(handle)
1042 :
1043 230 : END SUBROUTINE atomic_basis_at_grid_point
1044 :
1045 : ! **************************************************************************************************
1046 : !> \brief Computes the AO basis offsets: first_sgf(iatom) is the global index of the first
1047 : !> spherical Gaussian function (SGF) of iatom, first_sgf(natom+1) = total_sgf + 1,
1048 : !> and total_sgf is the total number of AO basis functions.
1049 : !> \param particle_set ...
1050 : !> \param qs_kind_set ...
1051 : !> \param first_sgf ...
1052 : !> \param total_sgf ...
1053 : ! **************************************************************************************************
1054 46 : SUBROUTINE get_basis_offsets(particle_set, qs_kind_set, first_sgf, total_sgf)
1055 :
1056 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1057 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1058 : INTEGER, INTENT(OUT) :: first_sgf(:), total_sgf
1059 :
1060 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_basis_offsets'
1061 :
1062 : INTEGER :: handle, iatom, ikind, nsgf
1063 :
1064 46 : CALL timeset(routineN, handle)
1065 :
1066 46 : total_sgf = 0
1067 158 : DO iatom = 1, SIZE(particle_set)
1068 112 : first_sgf(iatom) = total_sgf + 1
1069 112 : ikind = particle_set(iatom)%atomic_kind%kind_number
1070 112 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf, basis_type="ORB")
1071 158 : total_sgf = total_sgf + nsgf
1072 : END DO
1073 46 : first_sgf(SIZE(particle_set) + 1) = total_sgf + 1
1074 :
1075 46 : CALL timestop(handle)
1076 :
1077 46 : END SUBROUTINE get_basis_offsets
1078 :
1079 : ! **************************************************************************************************
1080 : !> \brief Computes χ_PQ(iτ) from the occupied and virtual Green's functions.
1081 : !> \param bs_env GW calculation environment and imaginary-time grid
1082 : !> \param mat_chi_Gamma_tau Response matrices χ_PQ(iτ) to be computed
1083 : !> \param mat_phi_mu_l AO values Φ_μ(r_l) on the RI-RS grid
1084 : !> \param mat_Z_lP RI-RS fitting matrix Z_lP
1085 : ! **************************************************************************************************
1086 38 : SUBROUTINE get_mat_chi_Gamma_tau(bs_env, mat_chi_Gamma_tau, mat_phi_mu_l, mat_Z_lP)
1087 :
1088 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1089 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1090 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
1091 :
1092 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
1093 :
1094 : INTEGER :: handle, i_t, ispin, n_panels
1095 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
1096 : REAL(KIND=dp) :: grid_occ, t1, tau
1097 : TYPE(dbcsr_type) :: matrix_G_occ_ao, matrix_G_vir_ao
1098 :
1099 38 : CALL timeset(routineN, handle)
1100 :
1101 : ! Panel boundaries for the grid-streaming contraction.
1102 : ! The panels are identical for χ, Σ^x and Σ^c, so the count is reported once here for all three stages.
1103 :
1104 38 : CALL resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
1105 :
1106 38 : n_panels = SIZE(pan_first)
1107 38 : IF (bs_env%unit_nr > 0) THEN
1108 : WRITE (bs_env%unit_nr, '(T2,A,T74,I9)') &
1109 19 : 'Number of batches for χ, Σ matrices', n_panels
1110 19 : WRITE (bs_env%unit_nr, '(A)') ' '
1111 19 : CALL m_flush(bs_env%unit_nr)
1112 : END IF
1113 :
1114 : ! =========================================================================
1115 : ! IMAGINARY TIME LOOP
1116 : ! χ_PQ(iτ) = Σ_s g_s · Z^T ( (φ G^occ_s φ^T) ∘ (φ G^vir_s φ^T) ) Z
1117 : ! (g_s = spin degeneracy)
1118 : ! =========================================================================
1119 568 : DO i_t = 1, bs_env%num_time_freq_points
1120 530 : t1 = m_walltime()
1121 530 : tau = bs_env%time_frequency_grid%imaginary_time(i_t)
1122 :
1123 1140 : DO ispin = 1, bs_env%n_spin
1124 :
1125 : ! AO-space Green's functions G^occ_µν, G^vir_µν (dense AO x AO, small)
1126 610 : CALL build_G_ao(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_G_occ_ao)
1127 610 : CALL build_G_ao(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, matrix_G_vir_ao)
1128 :
1129 : ! χ_PQ += g_s · Z^T ( (φ G^occ φ^T) ∘ (φ G^vir φ^T) ) Z
1130 : CALL contract_grid_panels(L_A=mat_phi_mu_l, M_A=matrix_G_occ_ao, &
1131 : L_B=mat_phi_mu_l, M_B=matrix_G_vir_ao, &
1132 : L_out=mat_Z_lP, mat_out=mat_chi_Gamma_tau(i_t)%matrix, &
1133 : scale=bs_env%spin_degeneracy, eps=bs_env%eps_filter, &
1134 : para_env=bs_env%para_env, &
1135 : pan_first=pan_first, pan_last=pan_last, &
1136 : lb_eq_la=.TRUE., lout_eq_la=.FALSE., &
1137 : zero_out=(ispin == 1), &
1138 : keep_sparsity=bs_env%ri_rs%keep_sparsity_rirs, &
1139 : centroids=bs_env%ri_rs%chunk_centroids, &
1140 : cutoff=bs_env%ri_rs%cutoff_radius_v_w, &
1141 610 : grid_occupation=grid_occ)
1142 :
1143 610 : CALL dbcsr_release(matrix_G_occ_ao)
1144 1140 : CALL dbcsr_release(matrix_G_vir_ao)
1145 :
1146 : END DO ! ispin
1147 :
1148 : ! Sparsity reports
1149 530 : IF (i_t == 1) THEN
1150 38 : CALL print_matrix_occupation(mat_Z_lP, 'Z_lP', bs_env)
1151 38 : IF (bs_env%unit_nr > 0) THEN
1152 : WRITE (bs_env%unit_nr, '(T2,A,T73,F7.2,A)') &
1153 19 : 'Percentage of non-zero matrix elements in G_ll'', χ_ll'', W_ll''', &
1154 38 : grid_occ*100.0_dp, ' %'
1155 19 : CALL m_flush(bs_env%unit_nr)
1156 : END IF
1157 : CALL print_matrix_occupation(mat_chi_Gamma_tau(i_t)%matrix, 'χ_PQ', bs_env, &
1158 38 : suffix=' for time point 1')
1159 38 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1160 : END IF
1161 :
1162 568 : IF (bs_env%unit_nr > 0) THEN
1163 : WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
1164 265 : 'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
1165 530 : ', Execution time', m_walltime() - t1, ' s'
1166 : END IF
1167 :
1168 : END DO ! i_t
1169 :
1170 38 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1171 :
1172 38 : CALL timestop(handle)
1173 :
1174 76 : END SUBROUTINE get_mat_chi_Gamma_tau
1175 :
1176 : ! **************************************************************************************************
1177 : !> \brief Marks the grid blocks whose centroid lies within cutoff of the bounding box of the
1178 : !> panel [blk0, blk1]'s chunk centroids.
1179 : !> \param centroids ...
1180 : !> \param blk0 ...
1181 : !> \param blk1 ...
1182 : !> \param cutoff ...
1183 : !> \param used ...
1184 : ! **************************************************************************************************
1185 92 : SUBROUTINE mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, used)
1186 :
1187 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
1188 : INTEGER, INTENT(IN) :: blk0, blk1
1189 : REAL(KIND=dp), INTENT(IN) :: cutoff
1190 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: used
1191 :
1192 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mask_grid_blocks_near_panel'
1193 :
1194 : INTEGER :: c, handle, k
1195 : REAL(KIND=dp) :: cutoff2, d2, dx
1196 : REAL(KIND=dp), DIMENSION(3) :: hi, lo
1197 :
1198 92 : CALL timeset(routineN, handle)
1199 :
1200 92 : cutoff2 = cutoff**2
1201 92 : lo(:) = MINVAL(centroids(:, blk0:blk1), DIM=2)
1202 92 : hi(:) = MAXVAL(centroids(:, blk0:blk1), DIM=2)
1203 :
1204 276 : ALLOCATE (used(SIZE(centroids, 2)))
1205 2024 : DO c = 1, SIZE(centroids, 2)
1206 : d2 = 0.0_dp
1207 7728 : DO k = 1, 3
1208 5796 : dx = MAX(0.0_dp, lo(k) - centroids(k, c), centroids(k, c) - hi(k))
1209 7728 : d2 = d2 + dx*dx
1210 : END DO
1211 2024 : used(c) = (d2 <= cutoff2)
1212 : END DO
1213 :
1214 92 : CALL timestop(handle)
1215 :
1216 92 : END SUBROUTINE mask_grid_blocks_near_panel
1217 :
1218 : ! **************************************************************************************************
1219 : !> \brief Exact allocated-element count of the geo template of panel [blk0, blk1]: the very same
1220 : !> per-block-pair centroid test as build_geo_template_panel, so this is the true DBCSR
1221 : !> data size of A_pan/B_pan/C_pan (DBCSR stores whole blocks).
1222 : !> \param r_blk_sizes ...
1223 : !> \param centroids ...
1224 : !> \param used ...
1225 : !> \param blk0 ...
1226 : !> \param blk1 ...
1227 : !> \param cutoff ...
1228 : !> \param nze_tmpl ...
1229 : ! **************************************************************************************************
1230 8 : SUBROUTINE panel_template_elems(r_blk_sizes, centroids, used, blk0, blk1, cutoff, nze_tmpl)
1231 :
1232 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes
1233 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
1234 : LOGICAL, DIMENSION(:), INTENT(IN) :: used
1235 : INTEGER, INTENT(IN) :: blk0, blk1
1236 : REAL(KIND=dp), INTENT(IN) :: cutoff
1237 : INTEGER(KIND=int_8), INTENT(OUT) :: nze_tmpl
1238 :
1239 : CHARACTER(LEN=*), PARAMETER :: routineN = 'panel_template_elems'
1240 :
1241 : INTEGER :: c, handle, ib, n_used
1242 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: used_idx
1243 : REAL(KIND=dp) :: cutoff2
1244 :
1245 8 : CALL timeset(routineN, handle)
1246 :
1247 : ! Compress the near mask once so the pair loop only visits candidate columns.
1248 176 : n_used = COUNT(used)
1249 24 : ALLOCATE (used_idx(n_used))
1250 8 : n_used = 0
1251 176 : DO c = 1, SIZE(used)
1252 176 : IF (used(c)) THEN
1253 168 : n_used = n_used + 1
1254 168 : used_idx(n_used) = c
1255 : END IF
1256 : END DO
1257 :
1258 8 : cutoff2 = cutoff**2
1259 8 : nze_tmpl = 0_int_8
1260 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(blk0, blk1, n_used, used_idx, centroids, cutoff2, &
1261 8 : !$OMP r_blk_sizes) PRIVATE(ib, c) REDUCTION(+:nze_tmpl)
1262 : DO ib = blk0, blk1
1263 : DO c = 1, n_used
1264 : IF (SUM((centroids(:, ib) - centroids(:, used_idx(c)))**2) <= cutoff2) THEN
1265 : nze_tmpl = nze_tmpl + INT(r_blk_sizes(ib), int_8)*INT(r_blk_sizes(used_idx(c)), int_8)
1266 : END IF
1267 : END DO
1268 : END DO
1269 : !$OMP END PARALLEL DO
1270 :
1271 8 : CALL timestop(handle)
1272 :
1273 16 : END SUBROUTINE panel_template_elems
1274 :
1275 : ! **************************************************************************************************
1276 : !> \brief Per-rank peak memory (GB) of one panel step of the neighborhood-restricted
1277 : !> contractions: three grid x grid panels of the template size (A_pan, B_pan, C_pan)
1278 : !> plus the grid x RI intermediates (tmp2 and the accumulation operand) and the
1279 : !> grid x AO intermediate (tmpA), whose column support is the panel's geometric
1280 : !> neighborhood fraction f_near = width/n_grid. Shared by the panel planner and
1281 : !> \param nze_tmpl ...
1282 : !> \param pan_rows ...
1283 : !> \param width ...
1284 : !> \param n_grid_total ...
1285 : !> \param n_RI ...
1286 : !> \param n_ao ...
1287 : !> \param n_procs ...
1288 : !> \param mem_GB ...
1289 : ! **************************************************************************************************
1290 8 : SUBROUTINE panel_mem_estimate_GB(nze_tmpl, pan_rows, width, n_grid_total, n_RI, n_ao, &
1291 : n_procs, mem_GB)
1292 :
1293 : INTEGER(KIND=int_8), INTENT(IN) :: nze_tmpl
1294 : INTEGER, INTENT(IN) :: pan_rows, width, n_grid_total, n_RI, &
1295 : n_ao, n_procs
1296 : REAL(KIND=dp), INTENT(OUT) :: mem_GB
1297 :
1298 : REAL(KIND=dp) :: f_near
1299 :
1300 8 : f_near = REAL(width, dp)/REAL(MAX(n_grid_total, 1), dp)
1301 : mem_GB = (3.0_dp*REAL(nze_tmpl, dp) + &
1302 : REAL(pan_rows, dp)*f_near*(2.0_dp*REAL(n_RI, dp) + REAL(n_ao, dp)))* &
1303 8 : 8.0_dp/REAL(MAX(n_procs, 1), dp)*1.0E-9_dp
1304 :
1305 8 : END SUBROUTINE panel_mem_estimate_GB
1306 :
1307 : ! **************************************************************************************************
1308 : !> \brief Plans the panel boundaries for the streaming contractions. Panels grow by whole grid
1309 : !> row-blocks towards ~panel_size rows. When the neighborhood restriction is active
1310 : !> (centroids+cutoff), each candidate panel is additionally checked against
1311 : !> (a) the 32-bit message bound with the panel's TRUE occupancy
1312 : !> (b) the per-rank memory budget: panel_mem_estimate_GB <= mem_budget_GB.
1313 : !> \param bs_env ...
1314 : !> \param r_blk_sizes ...
1315 : !> \param panel_size ...
1316 : !> \param min_dim ...
1317 : !> \param pan_first ...
1318 : !> \param pan_last ...
1319 : !> \param centroids ...
1320 : !> \param cutoff ...
1321 : !> \param n_RI ...
1322 : !> \param n_ao ...
1323 : !> \param n_procs ...
1324 : !> \param mem_budget_GB ...
1325 : !> \param honor_exact ...
1326 : !> \param unsafe ...
1327 : ! **************************************************************************************************
1328 38 : SUBROUTINE plan_grid_panels(bs_env, r_blk_sizes, panel_size, min_dim, pan_first, pan_last, &
1329 38 : centroids, cutoff, n_RI, n_ao, n_procs, mem_budget_GB, &
1330 : honor_exact, unsafe)
1331 :
1332 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1333 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes
1334 : INTEGER, INTENT(IN) :: panel_size, min_dim
1335 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: pan_first, pan_last
1336 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
1337 : OPTIONAL :: centroids
1338 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
1339 : INTEGER, INTENT(IN), OPTIONAL :: n_RI, n_ao, n_procs
1340 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: mem_budget_GB
1341 : LOGICAL, INTENT(IN), OPTIONAL :: honor_exact
1342 : LOGICAL, INTENT(OUT), OPTIONAL :: unsafe
1343 :
1344 : CHARACTER(LEN=*), PARAMETER :: routineN = 'plan_grid_panels'
1345 :
1346 : INTEGER :: blk0, blk1, handle, ib, n_grid_blocks, &
1347 : n_grid_total, n_panels, rows_acc, &
1348 : TARGET, width
1349 : INTEGER(KIND=int_8) :: msg, nze_tmpl, side
1350 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_first, tmp_last
1351 : LOGICAL :: fits, my_honor_exact, my_unsafe, &
1352 : use_cutoff
1353 38 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1354 : REAL(KIND=dp) :: f_near, mem_GB
1355 :
1356 38 : CALL timeset(routineN, handle)
1357 :
1358 38 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
1359 38 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
1360 2 : IF (use_cutoff) THEN
1361 2 : CPASSERT(PRESENT(n_RI) .AND. PRESENT(n_ao) .AND. PRESENT(n_procs))
1362 : END IF
1363 :
1364 : ! honor_exact: use N_PANELS as requested -- do NOT split a panel further even if it trips
1365 : ! the message-overflow / memory-budget check; instead flag `unsafe` so the caller can warn.
1366 38 : my_honor_exact = .FALSE.
1367 38 : IF (PRESENT(honor_exact)) my_honor_exact = honor_exact
1368 38 : my_unsafe = .FALSE.
1369 :
1370 38 : n_grid_blocks = SIZE(r_blk_sizes)
1371 556 : n_grid_total = SUM(r_blk_sizes)
1372 152 : ALLOCATE (tmp_first(n_grid_blocks), tmp_last(n_grid_blocks))
1373 :
1374 38 : n_panels = 0
1375 38 : blk0 = 1
1376 78 : DO WHILE (blk0 <= n_grid_blocks)
1377 40 : TARGET = panel_size
1378 0 : DO
1379 40 : rows_acc = 0
1380 40 : blk1 = blk0
1381 520 : DO ib = blk0, n_grid_blocks
1382 518 : rows_acc = rows_acc + r_blk_sizes(ib)
1383 518 : blk1 = ib
1384 520 : IF (rows_acc >= TARGET) EXIT
1385 : END DO
1386 40 : IF (.NOT. use_cutoff .OR. blk1 == blk0) EXIT
1387 4 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, used)
1388 88 : width = SUM(r_blk_sizes, MASK=used)
1389 4 : CALL panel_template_elems(r_blk_sizes, centroids, used, blk0, blk1, cutoff, nze_tmpl)
1390 4 : f_near = REAL(width, dp)/REAL(MAX(n_grid_total, 1), dp)
1391 4 : side = INT(REAL(rows_acc, dp)*f_near*REAL(MAX(n_RI, n_ao), dp), int_8)
1392 4 : msg = MAX(nze_tmpl, side)/INT(MAX(min_dim, 1), int_8)
1393 4 : fits = (msg <= bs_env%dbcsr_msg_elem_limit/4)
1394 4 : IF (fits .AND. PRESENT(mem_budget_GB)) THEN
1395 4 : IF (mem_budget_GB > 0.0_dp) THEN
1396 : CALL panel_mem_estimate_GB(nze_tmpl, rows_acc, width, n_grid_total, &
1397 4 : n_RI, n_ao, n_procs, mem_GB)
1398 4 : fits = (mem_GB <= mem_budget_GB)
1399 : END IF
1400 : END IF
1401 : ! Panel size is bounded only by the message-overflow and memory checks above; there is
1402 : ! no neighborhood-width (f_near) cap. mp_waitall is dominated by the NUMBER of panel
1403 : ! multiplies, so fewer/larger panels are cheaper here -- panel count is driven DOWN by
1404 : ! the N_PANELS keyword (panel_size), not split up by a width heuristic.
1405 4 : IF (my_honor_exact) THEN
1406 : ! Keep exactly the requested grouping; just record if it exceeds a safety limit.
1407 4 : IF (.NOT. fits) my_unsafe = .TRUE.
1408 : EXIT
1409 : END IF
1410 0 : IF (fits) EXIT
1411 36 : TARGET = MAX(1, MIN(TARGET, rows_acc)/2)
1412 : END DO
1413 40 : n_panels = n_panels + 1
1414 40 : tmp_first(n_panels) = blk0
1415 40 : tmp_last(n_panels) = blk1
1416 40 : blk0 = blk1 + 1
1417 : END DO
1418 :
1419 152 : ALLOCATE (pan_first(n_panels), pan_last(n_panels))
1420 78 : pan_first(:) = tmp_first(1:n_panels)
1421 78 : pan_last(:) = tmp_last(1:n_panels)
1422 38 : DEALLOCATE (tmp_first, tmp_last)
1423 :
1424 38 : IF (PRESENT(unsafe)) unsafe = my_unsafe
1425 :
1426 38 : CALL timestop(handle)
1427 :
1428 76 : END SUBROUTINE plan_grid_panels
1429 :
1430 : ! **************************************************************************************************
1431 : !> \brief Resolves the panel boundaries for the streaming contractions from the bs_env settings:
1432 : !> \param bs_env ...
1433 : !> \param mat_phi_mu_l ...
1434 : !> \param pan_first ...
1435 : !> \param pan_last ...
1436 : ! **************************************************************************************************
1437 166 : SUBROUTINE resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
1438 :
1439 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1440 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
1441 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: pan_first, pan_last
1442 :
1443 : CHARACTER(LEN=*), PARAMETER :: routineN = 'resolve_grid_panels'
1444 :
1445 : CHARACTER(LEN=max_line_length) :: msg
1446 : INTEGER :: handle, min_dim, n_grid_total, &
1447 : n_panels_req, npcols, nprows, &
1448 : panel_size, safe_max
1449 166 : INTEGER, DIMENSION(:), POINTER :: r_blk_sizes
1450 : LOGICAL :: honor_exact, panels_unsafe, use_cutoff
1451 : REAL(KIND=dp) :: mem_avail_GB, mem_budget_GB
1452 : TYPE(dbcsr_distribution_type) :: dist
1453 :
1454 166 : CALL timeset(routineN, handle)
1455 :
1456 166 : IF (ALLOCATED(bs_env%ri_rs%pan_first)) THEN
1457 518 : ALLOCATE (pan_first, SOURCE=bs_env%ri_rs%pan_first)
1458 518 : ALLOCATE (pan_last, SOURCE=bs_env%ri_rs%pan_last)
1459 128 : CALL timestop(handle)
1460 : RETURN
1461 : END IF
1462 :
1463 : use_cutoff = bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp .AND. &
1464 38 : ALLOCATED(bs_env%ri_rs%chunk_centroids)
1465 :
1466 : ! MIN(nprows, npcols) is the divisor that bounds the worst-rank Cannon message: a
1467 : ! P x n_grid panel is replicated into block row strips (P/nprows x n_grid) or column
1468 : ! strips (P x n_grid/npcols) during multiply_cannon, so the largest single-rank
1469 : ! message is ~ P*n_grid / MIN(nprows,npcols) elements.
1470 : CALL dbcsr_get_info(mat_phi_mu_l, nfullrows_total=n_grid_total, row_blk_size=r_blk_sizes, &
1471 38 : distribution=dist)
1472 38 : CALL dbcsr_distribution_get(dist, nprows=nprows, npcols=npcols)
1473 38 : min_dim = MAX(MIN(nprows, npcols), 1)
1474 :
1475 : ! Panel height such that NO per-rank DBCSR message can overflow the 32-bit length field
1476 : ! (see bs_env%dbcsr_msg_elem_limit): requiring the worst-rank message to stay under
1477 : ! 0.5 * HUGE(int_4) gives the safe height P_safe = 0.5 * HUGE(int_4) * min_dim / n_grid.
1478 38 : IF (use_cutoff) THEN
1479 2 : safe_max = n_grid_total
1480 : ELSE
1481 : safe_max = INT(0.5_dp*REAL(bs_env%dbcsr_msg_elem_limit, dp)*REAL(min_dim, dp)/ &
1482 36 : REAL(n_grid_total, dp))
1483 36 : safe_max = MAX(1, MIN(safe_max, n_grid_total))
1484 : END IF
1485 :
1486 : ! A user-set N_PANELS ( > 1 ) is honored EXACTLY: the planner produces that many panels
1487 : ! (up to grid-block granularity) and never force-splits them for the message/memory safety
1488 : ! limits -- if a limit is tripped it warns instead of silently changing the count.
1489 38 : n_panels_req = bs_env%ri_rs%n_panels
1490 38 : honor_exact = (n_panels_req > 1)
1491 38 : panels_unsafe = .FALSE.
1492 38 : IF (n_panels_req > 1) THEN
1493 : ! ceil(n_grid/n_panels_req) rows per panel => exactly n_panels_req panels. With the
1494 : ! cutoff active safe_max = n_grid_total (no clamp, honored exactly); without it, safe_max
1495 : ! is the int32-overflow ceiling and MUST still bound the panel (the non-cutoff planner
1496 : ! loop has no in-loop message-size check).
1497 2 : panel_size = MIN((n_grid_total + n_panels_req - 1)/n_panels_req, safe_max)
1498 : ELSE
1499 : ! Default (<= 1): a single whole-grid panel, clamped to the overflow-safe ceiling.
1500 36 : panel_size = safe_max
1501 : END IF
1502 38 : panel_size = MAX(1, panel_size)
1503 :
1504 38 : IF (use_cutoff) THEN
1505 : ! Half of the measured free memory as panel budget.
1506 2 : CALL mp_mem_avail_per_rank_GB(bs_env%para_env, mem_avail_GB)
1507 2 : mem_budget_GB = 0.5_dp*mem_avail_GB
1508 : CALL plan_grid_panels(bs_env, r_blk_sizes, panel_size, min_dim, pan_first, pan_last, &
1509 : centroids=bs_env%ri_rs%chunk_centroids, &
1510 : cutoff=bs_env%ri_rs%cutoff_radius_v_w, &
1511 : n_RI=bs_env%n_RI, n_ao=bs_env%n_ao, &
1512 : n_procs=bs_env%para_env%num_pe, mem_budget_GB=mem_budget_GB, &
1513 2 : honor_exact=honor_exact, unsafe=panels_unsafe)
1514 : ELSE
1515 36 : CALL plan_grid_panels(bs_env, r_blk_sizes, panel_size, min_dim, pan_first, pan_last)
1516 : END IF
1517 :
1518 38 : IF (honor_exact .AND. panels_unsafe) THEN
1519 : WRITE (msg, '(A,I0,A)') &
1520 0 : "N_PANELS = ", n_panels_req, " is used as requested, but one or more panels "// &
1521 : "exceed the DBCSR 32-bit message length or the memory budget. The run may abort "// &
1522 0 : "or swap; increase N_PANELS if it does."
1523 0 : CPWARN(TRIM(msg))
1524 : END IF
1525 :
1526 154 : ALLOCATE (bs_env%ri_rs%pan_first, SOURCE=pan_first)
1527 154 : ALLOCATE (bs_env%ri_rs%pan_last, SOURCE=pan_last)
1528 :
1529 38 : CALL timestop(handle)
1530 :
1531 204 : END SUBROUTINE resolve_grid_panels
1532 :
1533 : ! **************************************************************************************************
1534 : !> \brief Estimates and prints per-process memory requirements for the RI-RS GW calculation.
1535 : !> \param qs_env ...
1536 : !> \param bs_env ...
1537 : ! **************************************************************************************************
1538 38 : SUBROUTINE print_ri_rs_memory_estimate(qs_env, bs_env)
1539 :
1540 : !$ USE OMP_LIB, ONLY: omp_get_max_threads
1541 :
1542 : TYPE(qs_environment_type), POINTER :: qs_env
1543 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1544 :
1545 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_ri_rs_memory_estimate'
1546 :
1547 : CHARACTER(LEN=max_line_length) :: msg
1548 : INTEGER :: handle, iatom, ipan, l, &
1549 : max_n_ao_used, max_n_local_grid, &
1550 : n_ao_used_atom, n_grid_total, &
1551 : n_local_grid, n_loc_ri_max, n_procs, &
1552 : n_procs_per_atom, n_RI, n_threads, &
1553 : natom, pan_rows, pan_width
1554 : INTEGER(KIND=int_8) :: nze_tmpl
1555 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
1556 38 : INTEGER, DIMENSION(:), POINTER :: r_blk_sizes
1557 : LOGICAL :: use_cutoff
1558 38 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used
1559 : REAL(KIND=dp) :: cutoff_ri, mem_avail_GB, mem_D_local_GB, &
1560 : mem_dlp_GB, mem_pan_GB, mem_panels_GB, &
1561 : mem_phi_local_GB, mem_Z_lP_GB, &
1562 : mem_Zlp_peak_GB, pos_P(3)
1563 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1564 :
1565 38 : CALL timeset(routineN, handle)
1566 :
1567 : CALL dbcsr_get_info(bs_env%ri_rs%mat_phi_mu_l, nfullrows_total=n_grid_total, &
1568 38 : row_blk_size=r_blk_sizes)
1569 38 : n_RI = bs_env%n_RI
1570 38 : n_procs = bs_env%para_env%num_pe
1571 :
1572 : ! Z_lP upper bound: dense n_grid × n_RI, distributed evenly across all ranks.
1573 : ! The actual sparse Z_lP is smaller due to the per-atom locality cutoff.
1574 : mem_Z_lP_GB = REAL(n_grid_total, dp)*REAL(n_RI, dp)*8.0_dp/ &
1575 38 : REAL(n_procs, dp)*1.0E-9_dp
1576 :
1577 : ! Peak panel memory during Σ^c: two G panels (A_occ, A_vir) + one W panel plus the
1578 : ! grid × RI / grid × AO intermediates. With the CUTOFF_RADIUS_RL_W restriction the panel
1579 : ! matrices only allocate the geo-template blocks, so use the same nze-aware model as the
1580 : ! panel planner (panel_mem_estimate_GB); without the cutoff, dense panel_rows × n_grid.
1581 : ! Plus the n_RI × n_RI W_aux matrix. Distributed over n_procs ranks.
1582 : use_cutoff = bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp .AND. &
1583 38 : ALLOCATED(bs_env%ri_rs%chunk_centroids)
1584 38 : CALL resolve_grid_panels(bs_env, bs_env%ri_rs%mat_phi_mu_l, pan_first, pan_last)
1585 38 : mem_panels_GB = 0.0_dp
1586 78 : DO ipan = 1, SIZE(pan_first)
1587 558 : pan_rows = SUM(r_blk_sizes(pan_first(ipan):pan_last(ipan)))
1588 40 : IF (use_cutoff) THEN
1589 : CALL mask_grid_blocks_near_panel(bs_env%ri_rs%chunk_centroids, pan_first(ipan), &
1590 : pan_last(ipan), bs_env%ri_rs%cutoff_radius_v_w, &
1591 4 : grid_used)
1592 88 : pan_width = SUM(r_blk_sizes, MASK=grid_used)
1593 : CALL panel_template_elems(r_blk_sizes, bs_env%ri_rs%chunk_centroids, &
1594 : grid_used, pan_first(ipan), pan_last(ipan), &
1595 4 : bs_env%ri_rs%cutoff_radius_v_w, nze_tmpl)
1596 : CALL panel_mem_estimate_GB(nze_tmpl, pan_rows, pan_width, n_grid_total, &
1597 4 : n_RI, bs_env%n_ao, n_procs, mem_pan_GB)
1598 : ELSE
1599 36 : pan_width = n_grid_total
1600 : mem_pan_GB = (3.0_dp*REAL(pan_rows, dp)*REAL(pan_width, dp) + &
1601 : 2.0_dp*REAL(pan_rows, dp)*REAL(n_RI, dp))* &
1602 36 : 8.0_dp/REAL(n_procs, dp)*1.0E-9_dp
1603 : END IF
1604 78 : mem_panels_GB = MAX(mem_panels_GB, mem_pan_GB)
1605 : END DO
1606 : mem_panels_GB = mem_panels_GB + &
1607 38 : REAL(n_RI, dp)*REAL(n_RI, dp)*8.0_dp/REAL(n_procs, dp)*1.0E-9_dp
1608 :
1609 : ! Z_lP SOLVE peak (compute_Z_lP). For the atom P with the largest integration
1610 : ! sphere, one rank holds simultaneously:
1611 : ! D_local : n_local_grid x n_local_grid (dense Gram, BLAS path only; O(n_local_grid^2))
1612 : ! phi_local: n_local_grid x n_ao_used (AOs reaching into the sphere only)
1613 : ! d_lp : n_local_grid x n_loc_ri, replicated once + one private copy per OMP thread
1614 : ! n_local_grid = # grid points within cutoff_ri(P) = CUTOFF_RADIUS_RL_RI (if > 0) else
1615 : ! r_c(RI metric) + r_RI(P). This is NOT evenly distributed: n_local_grid depends on the
1616 : ! local density of atoms/grid, so the rank owning the densest atom peaks well above the
1617 : ! average. We report the worst-case (max over atoms) as a per-rank upper bound.
1618 38 : CALL get_qs_env(qs_env, particle_set=particle_set)
1619 38 : natom = bs_env%n_atom
1620 :
1621 38 : max_n_local_grid = 0
1622 38 : n_loc_ri_max = 0
1623 38 : max_n_ao_used = 0
1624 134 : DO iatom = 1, natom
1625 96 : IF (bs_env%ri_rs%cutoff_radius_ri_rs > 0.0_dp) THEN
1626 6 : cutoff_ri = bs_env%ri_rs%cutoff_radius_ri_rs
1627 : ELSE
1628 90 : cutoff_ri = bs_env%ri_metric%cutoff_radius + bs_env%ri_rs%radius_ri_per_atom(iatom)
1629 : END IF
1630 384 : pos_P(:) = particle_set(iatom)%r(:)
1631 : n_local_grid = 0
1632 37108 : DO l = 1, n_grid_total
1633 148144 : IF (SUM((bs_env%ri_rs%grid_points(1:3, l) - pos_P(1:3))**2) <= cutoff_ri**2) THEN
1634 35626 : n_local_grid = n_local_grid + 1
1635 : END IF
1636 : END DO
1637 96 : max_n_local_grid = MAX(max_n_local_grid, n_local_grid)
1638 96 : CALL get_n_ao_in_sphere(bs_env, particle_set, iatom, cutoff_ri, n_ao_used_atom)
1639 96 : max_n_ao_used = MAX(max_n_ao_used, n_ao_used_atom)
1640 : n_loc_ri_max = MAX(n_loc_ri_max, &
1641 134 : bs_env%i_RI_end_from_atom(iatom) - bs_env%i_RI_start_from_atom(iatom) + 1)
1642 : END DO
1643 :
1644 38 : n_procs_per_atom = MIN(MAX(bs_env%ri_rs%n_procs_per_atom_z_lp, 1), n_procs)
1645 38 : n_threads = 1
1646 38 : !$ n_threads = omp_get_max_threads()
1647 :
1648 : ! D_local: dense on one rank for the BLAS path; block-cyclic over the subgroup (=> /G) for
1649 : ! the ScaLAPACK path (N_PROCS_PER_ATOM_Z_LP = G > 1). phi_local/d_lp stay per-rank either way.
1650 38 : IF (n_procs_per_atom > 1) THEN
1651 0 : mem_D_local_GB = REAL(max_n_local_grid, dp)**2*8.0_dp/REAL(n_procs_per_atom, dp)*1.0E-9_dp
1652 : ELSE
1653 38 : mem_D_local_GB = REAL(max_n_local_grid, dp)**2*8.0_dp*1.0E-9_dp
1654 : END IF
1655 38 : mem_phi_local_GB = REAL(max_n_local_grid, dp)*REAL(max_n_ao_used, dp)*8.0_dp*1.0E-9_dp
1656 : mem_dlp_GB = REAL(max_n_local_grid, dp)*REAL(n_loc_ri_max, dp)*8.0_dp* &
1657 38 : REAL(1 + n_threads, dp)*1.0E-9_dp
1658 38 : mem_Zlp_peak_GB = mem_D_local_GB + mem_phi_local_GB + mem_dlp_GB
1659 :
1660 : ! Available memory per process = node MemLikelyFree / ranks-per-node, min across ranks
1661 : ! (0 on non-Linux => warnings suppressed below).
1662 38 : CALL mp_mem_avail_per_rank_GB(bs_env%para_env, mem_avail_GB)
1663 :
1664 38 : IF (bs_env%unit_nr > 0) THEN
1665 19 : WRITE (bs_env%unit_nr, '(A)') ' '
1666 19 : WRITE (bs_env%unit_nr, '(T2,A)') 'RI-RS memory estimate per MPI process:'
1667 : WRITE (bs_env%unit_nr, '(T4,A,F37.2,A)') &
1668 19 : 'Available memory per process (system)', mem_avail_GB, ' GB'
1669 : WRITE (bs_env%unit_nr, '(T4,A,F18.2,A)') &
1670 19 : 'Required for Z_lP (dense upper bound; actual is sparser)', mem_Z_lP_GB, ' GB'
1671 : WRITE (bs_env%unit_nr, '(T4,A,F25.2,A)') &
1672 19 : 'Required for χ, W, Σ panels (peak per panel step)', mem_panels_GB, ' GB'
1673 : WRITE (bs_env%unit_nr, '(T4,A,F17.2,A)') &
1674 19 : 'Required for Z_lP solve peak (D_local+φ, worst-case atom)', mem_Zlp_peak_GB, ' GB'
1675 : WRITE (bs_env%unit_nr, '(T4,A,T69,I12)') &
1676 19 : 'Worst-case local-grid number of grid points:', max_n_local_grid
1677 : WRITE (bs_env%unit_nr, '(T4,A,T69,F9.2,A)') &
1678 19 : 'Worst-case memory D_local:', mem_D_local_GB, ' GB'
1679 :
1680 : END IF
1681 :
1682 38 : IF (mem_avail_GB > 0.0_dp .AND. mem_Z_lP_GB > mem_avail_GB) THEN
1683 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
1684 0 : "The estimated memory for Z_lP, ", mem_Z_lP_GB, " GB per process, exceeds the "// &
1685 0 : "available ", mem_avail_GB, " GB. Z_lP (n_grid x n_RI) is distributed across all "// &
1686 : "MPI ranks, so add nodes, use fewer MPI ranks per node, or raise "// &
1687 : "N_PROCS_PER_ATOM_Z_LP to distribute each atom block via ScaLAPACK, which reduces "// &
1688 0 : "the per-rank memory roughly by the number of ranks per atom."
1689 0 : CPWARN(TRIM(msg))
1690 : END IF
1691 :
1692 38 : IF (mem_avail_GB > 0.0_dp .AND. mem_panels_GB > mem_avail_GB) THEN
1693 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
1694 0 : "The estimated peak memory of the chi/W/Sigma panels, ", mem_panels_GB, &
1695 0 : " GB per process, exceeds the available ", mem_avail_GB, " GB. Panel memory "// &
1696 : "scales roughly as 3*panel_size*n_grid/n_procs, so add nodes, use fewer MPI ranks "// &
1697 0 : "per node, or raise N_PANELS for more but smaller panels."
1698 0 : CPWARN(TRIM(msg))
1699 : END IF
1700 :
1701 38 : IF (mem_avail_GB > 0.0_dp .AND. mem_Zlp_peak_GB > mem_avail_GB) THEN
1702 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
1703 0 : "The estimated peak memory of the Z_lP solve, ", mem_Zlp_peak_GB, &
1704 0 : " GB per process, exceeds the available ", mem_avail_GB, &
1705 : " GB. The per-atom Gram matrix D_local(n_local_grid, n_local_grid) dominates and "// &
1706 : "scales as n_local_grid^2, and it is not balanced across ranks: the rank owning "// &
1707 : "the atom with the largest integration sphere peaks well above the average. "// &
1708 : "Either raise N_PROCS_PER_ATOM_Z_LP to distribute D_local block-cyclic via "// &
1709 : "ScaLAPACK, which reduces that term roughly by the number of ranks per atom at no "// &
1710 : "loss of accuracy, or lower CUTOFF_RADIUS_RL_RI, which shrinks D_local as "// &
1711 : "n_local_grid^2 but trades accuracy, or use fewer MPI ranks per node so that each "// &
1712 0 : "rank has more memory for the peak atom."
1713 0 : CPWARN(TRIM(msg))
1714 : END IF
1715 :
1716 38 : CALL timestop(handle)
1717 :
1718 76 : END SUBROUTINE print_ri_rs_memory_estimate
1719 :
1720 : ! **************************************************************************************************
1721 : !> \brief Creates an empty (panel_chunks x neighborhood_chunks) DBCSR matrix with zero blocks
1722 : !> pre-allocated only where |centroid(panel_row r) - centroid(column c)| <= cutoff.
1723 : !> Used with retain_sparsity=.TRUE. in the subsequent dbcsr_multiply so distant blocks
1724 : !> of the grid-basis panels (φ G φ^T, Z W Z^T, ...) are never computed at all.
1725 : !> \param L_pan ...
1726 : !> \param L_full ...
1727 : !> \param centroids ...
1728 : !> \param cutoff ...
1729 : !> \param blk0 ...
1730 : !> \param A_template ...
1731 : !> \param col_map ...
1732 : ! **************************************************************************************************
1733 124 : SUBROUTINE build_geo_template_panel(L_pan, L_full, centroids, cutoff, blk0, A_template, col_map)
1734 : TYPE(dbcsr_type), INTENT(IN) :: L_pan, L_full
1735 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
1736 : REAL(KIND=dp), INTENT(IN) :: cutoff
1737 : INTEGER, INTENT(IN) :: blk0
1738 : TYPE(dbcsr_type), INTENT(OUT) :: A_template
1739 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: col_map
1740 :
1741 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_geo_template_panel'
1742 :
1743 : INTEGER :: c, cg, cs, handle, my_pcol, my_prow, &
1744 : n_grid_blks, n_pan_blks, npcols, &
1745 : nprows, r, rs
1746 124 : INTEGER, DIMENSION(:), POINTER :: grid_blk_sizes, pan_blk_sizes
1747 : REAL(KIND=dp) :: cutoff2
1748 124 : REAL(KIND=dp), ALLOCATABLE :: zero_blk(:, :)
1749 : TYPE(dbcsr_distribution_type) :: dist
1750 :
1751 124 : CALL timeset(routineN, handle)
1752 :
1753 124 : cutoff2 = cutoff**2
1754 124 : CALL dbcsr_get_info(L_pan, nblkrows_total=n_pan_blks, row_blk_size=pan_blk_sizes)
1755 124 : CALL dbcsr_get_info(L_full, nblkrows_total=n_grid_blks, row_blk_size=grid_blk_sizes)
1756 :
1757 : ! create_product_matrix assigns row r to process MOD(r-1,nprows) and
1758 : ! col c to MOD(c-1,npcols), so we can determine local ownership analytically.
1759 124 : CALL create_product_matrix(L_pan, L_full, 'N', 'T', A_template)
1760 124 : CALL dbcsr_get_info(A_template, distribution=dist)
1761 : CALL dbcsr_distribution_get(dist, nprows=nprows, npcols=npcols, &
1762 124 : myprow=my_prow, mypcol=my_pcol)
1763 :
1764 : ALLOCATE (zero_blk(MAXVAL(pan_blk_sizes(1:n_pan_blks)), &
1765 4402 : MAXVAL(grid_blk_sizes(1:n_grid_blks))))
1766 124 : zero_blk(:, :) = 0.0_dp
1767 :
1768 1426 : DO r = 1, n_pan_blks
1769 1302 : IF (MOD(r - 1, nprows) /= my_prow) CYCLE
1770 651 : rs = pan_blk_sizes(r)
1771 14446 : DO c = 1, n_grid_blks
1772 13671 : IF (MOD(c - 1, npcols) /= my_pcol) CYCLE
1773 13671 : cg = c
1774 13671 : IF (PRESENT(col_map)) cg = col_map(c)
1775 13671 : IF ((centroids(1, blk0 + r - 1) - centroids(1, cg))**2 + &
1776 : (centroids(2, blk0 + r - 1) - centroids(2, cg))**2 + &
1777 1302 : (centroids(3, blk0 + r - 1) - centroids(3, cg))**2 <= cutoff2) THEN
1778 12741 : cs = grid_blk_sizes(c)
1779 12741 : CALL dbcsr_put_block(A_template, r, c, zero_blk(1:rs, 1:cs))
1780 : END IF
1781 : END DO
1782 : END DO
1783 124 : CALL dbcsr_finalize(A_template)
1784 :
1785 124 : DEALLOCATE (zero_blk)
1786 124 : CALL timestop(handle)
1787 :
1788 124 : END SUBROUTINE build_geo_template_panel
1789 :
1790 : ! **************************************************************************************************
1791 : !> \brief Slices a contiguous range of grid row-blocks [blk0, blk1] out of a (grid x n) DBCSR
1792 : !> matrix into a new (P x n) panel matrix: iterate the source's local blocks, put the
1793 : !> in-range ones into the panel with a remapped row-block index, then finalize. Row-block
1794 : !> index i of the panel corresponds to source row-block blk0+i-1.
1795 : !> \param mat_full ...
1796 : !> \param blk0 ...
1797 : !> \param blk1 ...
1798 : !> \param mat_panel ...
1799 : ! **************************************************************************************************
1800 3012 : SUBROUTINE extract_grid_panel(mat_full, blk0, blk1, mat_panel)
1801 :
1802 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_full
1803 : INTEGER, INTENT(IN) :: blk0, blk1
1804 : TYPE(dbcsr_type), INTENT(OUT) :: mat_panel
1805 :
1806 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extract_grid_panel'
1807 :
1808 : INTEGER :: handle, ib, jb, npb
1809 3012 : INTEGER, DIMENSION(:), POINTER :: col_blk_full, col_dist_full, &
1810 3012 : row_blk_full, row_blk_pan, &
1811 3012 : row_dist_full, row_dist_pan
1812 3012 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
1813 : TYPE(dbcsr_distribution_type) :: dist_full, dist_pan
1814 : TYPE(dbcsr_iterator_type) :: iter
1815 :
1816 3012 : CALL timeset(routineN, handle)
1817 :
1818 : CALL dbcsr_get_info(mat_full, distribution=dist_full, &
1819 3012 : row_blk_size=row_blk_full, col_blk_size=col_blk_full)
1820 3012 : CALL dbcsr_distribution_get(dist_full, row_dist=row_dist_full, col_dist=col_dist_full)
1821 :
1822 3012 : npb = blk1 - blk0 + 1
1823 12048 : ALLOCATE (row_dist_pan(npb), row_blk_pan(npb))
1824 91344 : row_dist_pan(:) = row_dist_full(blk0:blk1)
1825 91344 : row_blk_pan(:) = row_blk_full(blk0:blk1)
1826 :
1827 : CALL dbcsr_distribution_new(dist_pan, template=dist_full, &
1828 3012 : row_dist=row_dist_pan, col_dist=col_dist_full)
1829 : CALL dbcsr_create(mat_panel, name="grid_panel", dist=dist_pan, &
1830 : matrix_type=dbcsr_type_no_symmetry, &
1831 3012 : row_blk_size=row_blk_pan, col_blk_size=col_blk_full)
1832 :
1833 3012 : CALL dbcsr_iterator_start(iter, mat_full)
1834 55416 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1835 52404 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
1836 52404 : IF (ib < blk0 .OR. ib > blk1) CYCLE
1837 55416 : CALL dbcsr_put_block(mat_panel, ib - blk0 + 1, jb, blk)
1838 : END DO
1839 3012 : CALL dbcsr_iterator_stop(iter)
1840 3012 : CALL dbcsr_finalize(mat_panel)
1841 :
1842 3012 : CALL dbcsr_distribution_release(dist_pan)
1843 3012 : DEALLOCATE (row_dist_pan, row_blk_pan)
1844 :
1845 3012 : CALL timestop(handle)
1846 :
1847 9036 : END SUBROUTINE extract_grid_panel
1848 :
1849 : ! **************************************************************************************************
1850 : !> \brief Marks which column blocks of a DBCSR matrix carry at least one non-zero block anywhere
1851 : !> (global union). Used to restrict the inner index of the panel multiplies to the
1852 : !> AO/RI atoms that actually touch the panel (exact: dropped rows only meet zeros).
1853 : !> \param matrix ...
1854 : !> \param para_env ...
1855 : !> \param used ...
1856 : ! **************************************************************************************************
1857 2382 : SUBROUTINE collect_used_col_blocks(matrix, para_env, used)
1858 :
1859 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1860 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
1861 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: used
1862 :
1863 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collect_used_col_blocks'
1864 :
1865 : INTEGER :: handle, ib, jb, nblkcols
1866 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iused
1867 2382 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
1868 : TYPE(dbcsr_iterator_type) :: iter
1869 :
1870 2382 : CALL timeset(routineN, handle)
1871 :
1872 2382 : CALL dbcsr_get_info(matrix, nblkcols_total=nblkcols)
1873 7146 : ALLOCATE (iused(nblkcols))
1874 2382 : iused(:) = 0
1875 :
1876 2382 : CALL dbcsr_iterator_start(iter, matrix)
1877 43039 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1878 40657 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
1879 40657 : iused(jb) = 1
1880 : END DO
1881 2382 : CALL dbcsr_iterator_stop(iter)
1882 :
1883 2382 : CALL para_env%sum(iused)
1884 :
1885 4764 : ALLOCATE (used(nblkcols))
1886 8040 : used(:) = (iused(:) > 0)
1887 2382 : DEALLOCATE (iused)
1888 :
1889 2382 : CALL timestop(handle)
1890 :
1891 4764 : END SUBROUTINE collect_used_col_blocks
1892 :
1893 : ! **************************************************************************************************
1894 : !> \brief Copies the flagged block rows (compress_rows=.TRUE.) or block columns (.FALSE.) of a
1895 : !> DBCSR matrix into a compressed matrix. The subset keeps the parent's process assignment
1896 : !> along the compressed dimension, so every block stays on its owning rank: the extraction
1897 : !> is purely local (zero communication), like extract_grid_panel.
1898 : !> \param mat_full ...
1899 : !> \param used ...
1900 : !> \param mat_out ...
1901 : !> \param compress_rows ...
1902 : !> \param blk_map ...
1903 : ! **************************************************************************************************
1904 6392 : SUBROUTINE extract_masked_blocks(mat_full, used, mat_out, compress_rows, blk_map)
1905 :
1906 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_full
1907 : LOGICAL, DIMENSION(:), INTENT(IN) :: used
1908 : TYPE(dbcsr_type), INTENT(OUT) :: mat_out
1909 : LOGICAL, INTENT(IN) :: compress_rows
1910 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT), &
1911 : OPTIONAL :: blk_map
1912 :
1913 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extract_masked_blocks'
1914 :
1915 : INTEGER :: handle, ib, jb, n_blk, n_sub, r
1916 6392 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inv_map
1917 6392 : INTEGER, DIMENSION(:), POINTER :: blk_full, blk_sub, col_blk_full, &
1918 6392 : col_dist_full, dist_full_1d, &
1919 6392 : dist_sub_1d, row_blk_full, &
1920 6392 : row_dist_full
1921 6392 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
1922 : TYPE(dbcsr_distribution_type) :: dist_full, dist_sub
1923 : TYPE(dbcsr_iterator_type) :: iter
1924 :
1925 6392 : CALL timeset(routineN, handle)
1926 :
1927 : CALL dbcsr_get_info(mat_full, distribution=dist_full, &
1928 6392 : row_blk_size=row_blk_full, col_blk_size=col_blk_full)
1929 6392 : CALL dbcsr_distribution_get(dist_full, row_dist=row_dist_full, col_dist=col_dist_full)
1930 :
1931 6392 : IF (compress_rows) THEN
1932 4010 : blk_full => row_blk_full
1933 4010 : dist_full_1d => row_dist_full
1934 : ELSE
1935 2382 : blk_full => col_blk_full
1936 2382 : dist_full_1d => col_dist_full
1937 : END IF
1938 6392 : n_blk = SIZE(blk_full)
1939 :
1940 24696 : n_sub = COUNT(used)
1941 6392 : CPASSERT(n_sub > 0)
1942 38352 : ALLOCATE (inv_map(n_blk), blk_sub(n_sub), dist_sub_1d(n_sub))
1943 6476 : IF (PRESENT(blk_map)) ALLOCATE (blk_map(n_sub))
1944 6392 : inv_map(:) = 0
1945 6392 : r = 0
1946 24696 : DO ib = 1, n_blk
1947 24696 : IF (used(ib)) THEN
1948 18304 : r = r + 1
1949 18304 : inv_map(ib) = r
1950 18304 : blk_sub(r) = blk_full(ib)
1951 18304 : dist_sub_1d(r) = dist_full_1d(ib)
1952 18304 : IF (PRESENT(blk_map)) blk_map(r) = ib
1953 : END IF
1954 : END DO
1955 :
1956 6392 : IF (compress_rows) THEN
1957 : CALL dbcsr_distribution_new(dist_sub, template=dist_full, &
1958 4010 : row_dist=dist_sub_1d, col_dist=col_dist_full)
1959 : CALL dbcsr_create(mat_out, name="row_subset", dist=dist_sub, &
1960 : matrix_type=dbcsr_type_no_symmetry, &
1961 4010 : row_blk_size=blk_sub, col_blk_size=col_blk_full)
1962 : ELSE
1963 : CALL dbcsr_distribution_new(dist_sub, template=dist_full, &
1964 2382 : row_dist=row_dist_full, col_dist=dist_sub_1d)
1965 : CALL dbcsr_create(mat_out, name="col_subset", dist=dist_sub, &
1966 : matrix_type=dbcsr_type_no_symmetry, &
1967 2382 : row_blk_size=row_blk_full, col_blk_size=blk_sub)
1968 : END IF
1969 :
1970 6392 : CALL dbcsr_iterator_start(iter, mat_full)
1971 62602 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1972 56210 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
1973 62602 : IF (compress_rows) THEN
1974 15553 : IF (inv_map(ib) > 0) CALL dbcsr_put_block(mat_out, inv_map(ib), jb, blk)
1975 : ELSE
1976 40657 : IF (inv_map(jb) > 0) CALL dbcsr_put_block(mat_out, ib, inv_map(jb), blk)
1977 : END IF
1978 : END DO
1979 6392 : CALL dbcsr_iterator_stop(iter)
1980 6392 : CALL dbcsr_finalize(mat_out)
1981 :
1982 6392 : CALL dbcsr_distribution_release(dist_sub)
1983 6392 : DEALLOCATE (inv_map, blk_sub, dist_sub_1d)
1984 :
1985 6392 : CALL timestop(handle)
1986 :
1987 12784 : END SUBROUTINE extract_masked_blocks
1988 :
1989 : ! **************************************************************************************************
1990 : !> \brief Pre-seeds a square blocked DBCSR matrix with zero blocks only for block pairs whose
1991 : !> centers lie within radius, for use with copy_fm_to_dbcsr(keep_sparsity=T) or
1992 : !> dbcsr_multiply(retain_sparsity=T). Consumers: the CUTOFF_RADIUS_G_W operator truncation
1993 : !> (atom-blocked, centers = atom_centers) and the RT-BSE CUTOFF_RADIUS_W0 truncation of the
1994 : !> grid-basis W^0 (grid-blocked, centers = chunk_centroids).
1995 : !> \param matrix ...
1996 : !> \param centers block positions, one column per block row/column of matrix
1997 : !> \param radius truncation radius, same units as centers (bohr)
1998 : ! **************************************************************************************************
1999 6 : SUBROUTINE reserve_blocks_within_radius(matrix, centers, radius)
2000 :
2001 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
2002 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centers
2003 : REAL(KIND=dp), INTENT(IN) :: radius
2004 :
2005 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reserve_blocks_within_radius'
2006 :
2007 : INTEGER :: handle, i, j, my_pcol, my_prow, &
2008 : nblkcols, nblkrows
2009 6 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist, row_blk, row_dist
2010 : REAL(KIND=dp) :: radius2
2011 6 : REAL(KIND=dp), ALLOCATABLE :: zero_blk(:, :)
2012 : TYPE(dbcsr_distribution_type) :: dist
2013 :
2014 6 : CALL timeset(routineN, handle)
2015 :
2016 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows, nblkcols_total=nblkcols, &
2017 6 : row_blk_size=row_blk, col_blk_size=col_blk, distribution=dist)
2018 : CALL dbcsr_distribution_get(dist, row_dist=row_dist, col_dist=col_dist, &
2019 6 : myprow=my_prow, mypcol=my_pcol)
2020 6 : CPASSERT(nblkrows == SIZE(centers, 2))
2021 6 : CPASSERT(nblkcols == SIZE(centers, 2))
2022 :
2023 6 : radius2 = radius**2
2024 144 : ALLOCATE (zero_blk(MAXVAL(row_blk(1:nblkrows)), MAXVAL(col_blk(1:nblkcols))))
2025 6 : zero_blk(:, :) = 0.0_dp
2026 :
2027 66 : DO i = 1, nblkrows
2028 60 : IF (row_dist(i) /= my_prow) CYCLE
2029 336 : DO j = 1, nblkcols
2030 300 : IF (col_dist(j) /= my_pcol) CYCLE
2031 300 : IF ((centers(1, i) - centers(1, j))**2 + (centers(2, i) - centers(2, j))**2 + &
2032 60 : (centers(3, i) - centers(3, j))**2 <= radius2) THEN
2033 168 : CALL dbcsr_put_block(matrix, i, j, zero_blk(1:row_blk(i), 1:col_blk(j)))
2034 : END IF
2035 : END DO
2036 : END DO
2037 6 : CALL dbcsr_finalize(matrix)
2038 :
2039 6 : DEALLOCATE (zero_blk)
2040 6 : CALL timestop(handle)
2041 :
2042 6 : END SUBROUTINE reserve_blocks_within_radius
2043 :
2044 : ! **************************************************************************************************
2045 : !> \brief Creates the (empty) result matrix of op(mat_left) * op(mat_right) with the correct block
2046 : !> structure and a distribution on the shared process grid, ready to be filled by
2047 : !> dbcsr_multiply. Row structure comes from op(left), column structure from op(right).
2048 : !> \param mat_left ...
2049 : !> \param mat_right ...
2050 : !> \param transa 'N' or 'T' applied to mat_left
2051 : !> \param transb 'N' or 'T' applied to mat_right
2052 : !> \param mat_out ...
2053 : ! **************************************************************************************************
2054 8514 : SUBROUTINE create_product_matrix(mat_left, mat_right, transa, transb, mat_out)
2055 :
2056 : TYPE(dbcsr_type), INTENT(IN) :: mat_left, mat_right
2057 : CHARACTER(LEN=1), INTENT(IN) :: transa, transb
2058 : TYPE(dbcsr_type), INTENT(OUT) :: mat_out
2059 :
2060 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_product_matrix'
2061 :
2062 : INTEGER :: handle, i, npcols, nprows
2063 8514 : INTEGER, DIMENSION(:), POINTER :: col_blk_l, col_blk_r, out_col_blk, &
2064 8514 : out_col_dist, out_row_blk, &
2065 8514 : out_row_dist, row_blk_l, row_blk_r
2066 : TYPE(dbcsr_distribution_type) :: dist_l, dist_out
2067 :
2068 8514 : CALL timeset(routineN, handle)
2069 :
2070 8514 : CALL dbcsr_get_info(mat_left, distribution=dist_l, row_blk_size=row_blk_l, col_blk_size=col_blk_l)
2071 8514 : CALL dbcsr_get_info(mat_right, row_blk_size=row_blk_r, col_blk_size=col_blk_r)
2072 8514 : CALL dbcsr_distribution_get(dist_l, nprows=nprows, npcols=npcols)
2073 :
2074 : ! block SIZES follow op(left)/op(right); DISTRIBUTIONS are freshly round-robined onto the
2075 : ! shared process grid (a transposed operand's row-dist is NOT a valid col-dist on a
2076 : ! non-square grid). dbcsr_multiply redistributes internally, so any valid mapping works.
2077 8514 : IF (transa == 'N') THEN
2078 8514 : out_row_blk => row_blk_l
2079 : ELSE
2080 0 : out_row_blk => col_blk_l
2081 : END IF
2082 8514 : IF (transb == 'N') THEN
2083 6178 : out_col_blk => col_blk_r
2084 : ELSE
2085 2336 : out_col_blk => row_blk_r
2086 : END IF
2087 :
2088 42570 : ALLOCATE (out_row_dist(SIZE(out_row_blk)), out_col_dist(SIZE(out_col_blk)))
2089 131364 : DO i = 1, SIZE(out_row_blk)
2090 131364 : out_row_dist(i) = MOD(i - 1, nprows)
2091 : END DO
2092 58318 : DO i = 1, SIZE(out_col_blk)
2093 58318 : out_col_dist(i) = MOD(i - 1, npcols)
2094 : END DO
2095 :
2096 : CALL dbcsr_distribution_new(dist_out, template=dist_l, &
2097 8514 : row_dist=out_row_dist, col_dist=out_col_dist)
2098 : CALL dbcsr_create(mat_out, name="panel_product", dist=dist_out, &
2099 : matrix_type=dbcsr_type_no_symmetry, &
2100 8514 : row_blk_size=out_row_blk, col_blk_size=out_col_blk)
2101 8514 : CALL dbcsr_distribution_release(dist_out)
2102 8514 : DEALLOCATE (out_row_dist, out_col_dist)
2103 :
2104 8514 : CALL timestop(handle)
2105 :
2106 25542 : END SUBROUTINE create_product_matrix
2107 :
2108 : ! **************************************************************************************************
2109 : !> \brief Builds the AO-space Green's function operator G^occ/vir_µν (AO x AO DBCSR)
2110 : !> \param bs_env ...
2111 : !> \param tau ...
2112 : !> \param ispin ...
2113 : !> \param occ ...
2114 : !> \param vir ...
2115 : !> \param template ...
2116 : !> \param matrix_G_ao ...
2117 : ! **************************************************************************************************
2118 5768 : SUBROUTINE build_G_ao(bs_env, tau, ispin, occ, vir, template, matrix_G_ao)
2119 :
2120 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2121 : REAL(KIND=dp), INTENT(IN) :: tau
2122 : INTEGER, INTENT(IN) :: ispin
2123 : LOGICAL, INTENT(IN) :: occ, vir
2124 : TYPE(dbcsr_type), INTENT(INOUT) :: template
2125 : TYPE(dbcsr_type), INTENT(OUT) :: matrix_G_ao
2126 :
2127 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_G_ao'
2128 :
2129 : INTEGER :: handle
2130 2884 : INTEGER, DIMENSION(:), POINTER :: blk_ao, dist_row_ao
2131 : TYPE(cp_fm_type), POINTER :: fm_G
2132 : TYPE(dbcsr_distribution_type) :: dist_ao_ao
2133 :
2134 2884 : CALL timeset(routineN, handle)
2135 :
2136 2884 : IF (occ) THEN
2137 1464 : fm_G => bs_env%fm_Gocc
2138 : ELSE
2139 1420 : fm_G => bs_env%fm_Gvir
2140 : END IF
2141 :
2142 2884 : CALL G_occ_vir(bs_env, tau, fm_G, ispin, occ=occ, vir=vir)
2143 :
2144 2884 : CALL setup_square_topology(template, dist_ao_ao, blk_ao, dist_row_ao)
2145 : CALL dbcsr_create(matrix_G_ao, name="G_ao", dist=dist_ao_ao, &
2146 : matrix_type=dbcsr_type_no_symmetry, &
2147 2884 : row_blk_size=blk_ao, col_blk_size=blk_ao)
2148 :
2149 : ! Optional CUTOFF_RADIUS_G_W operator truncation: only atom-pair blocks within the radius
2150 : ! are reserved and filled.
2151 2884 : IF (bs_env%ri_rs%cutoff_radius_g_w > 0.0_dp .AND. &
2152 : ALLOCATED(bs_env%ri_rs%atom_centers)) THEN
2153 : CALL reserve_blocks_within_radius(matrix_G_ao, bs_env%ri_rs%atom_centers, &
2154 0 : bs_env%ri_rs%cutoff_radius_g_w)
2155 0 : CALL copy_fm_to_dbcsr(fm_G, matrix_G_ao, keep_sparsity=.TRUE.)
2156 : ELSE
2157 2884 : CALL copy_fm_to_dbcsr(fm_G, matrix_G_ao, keep_sparsity=.FALSE.)
2158 : END IF
2159 2884 : CALL dbcsr_filter(matrix_G_ao, bs_env%eps_filter)
2160 :
2161 : ! release only the topology; keep matrix_G_ao for the caller
2162 2884 : CALL release_square_topology(dist=dist_ao_ao, mapped_dist=dist_row_ao)
2163 :
2164 2884 : CALL timestop(handle)
2165 :
2166 2884 : END SUBROUTINE build_G_ao
2167 :
2168 : ! **************************************************************************************************
2169 : !> \brief Panel-streaming evaluation of out += scale * L_out^T (A_grid ∘ B_grid) L_out,
2170 : !> where A_grid = L_A M_A L_A^T and B_grid = L_B M_B L_B^T, WITHOUT ever forming the full
2171 : !> grid x grid objects. The grid (row) index is processed in panels of ~panel_size rows; for
2172 : !> each panel only P x grid slabs are built, Hadamard-multiplied, and contracted into the
2173 : !> (small) output. Algebraically identical to L_out^T (A_grid ∘ B_grid) L_out summed over
2174 : !> grid rows, so the result matches the non-streamed path to eps_filter.
2175 : !>
2176 : !> Mapping (L in {phi (grid x AO), Z (grid x RI)}, M the AO/RI-space operator):
2177 : !> chi : L_A=L_B=phi, M_A=G_occ_ao, M_B=G_vir_ao, L_out=Z -> RI x RI
2178 : !> Sig : L_A=phi (M_A=D/G), L_B=Z (M_B=V/W), L_out=phi -> AO x AO
2179 : !> \param L_A ...
2180 : !> \param M_A ...
2181 : !> \param L_B ...
2182 : !> \param M_B ...
2183 : !> \param L_out ...
2184 : !> \param mat_out ...
2185 : !> \param scale ...
2186 : !> \param eps ...
2187 : !> \param para_env ...
2188 : !> \param pan_first ...
2189 : !> \param pan_last ...
2190 : !> \param lb_eq_la ...
2191 : !> \param lout_eq_la ...
2192 : !> \param zero_out ...
2193 : !> \param keep_sparsity ...
2194 : !> \param centroids ...
2195 : !> \param cutoff ...
2196 : !> \param grid_occupation ...
2197 : ! **************************************************************************************************
2198 654 : SUBROUTINE contract_grid_panels(L_A, M_A, L_B, M_B, L_out, mat_out, scale, eps, para_env, &
2199 654 : pan_first, pan_last, lb_eq_la, lout_eq_la, zero_out, &
2200 654 : keep_sparsity, centroids, cutoff, grid_occupation)
2201 :
2202 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_A
2203 : TYPE(dbcsr_type), INTENT(INOUT) :: M_A
2204 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_B
2205 : TYPE(dbcsr_type), INTENT(INOUT) :: M_B
2206 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_out
2207 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_out
2208 : REAL(KIND=dp), INTENT(IN) :: scale, eps
2209 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
2210 : INTEGER, DIMENSION(:), INTENT(IN) :: pan_first, pan_last
2211 : LOGICAL, INTENT(IN) :: lb_eq_la, lout_eq_la, zero_out
2212 : LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
2213 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
2214 : OPTIONAL :: centroids
2215 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
2216 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: grid_occupation
2217 :
2218 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_grid_panels'
2219 :
2220 : INTEGER :: blk0, blk1, handle, ipan, n_grid_total, &
2221 : ncols_pan, nrows_pan
2222 654 : INTEGER, ALLOCATABLE, DIMENSION(:) :: gmap
2223 : LOGICAL :: my_keep_sparsity, use_cutoff
2224 654 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used, usedA, usedB
2225 : TYPE(dbcsr_type) :: A_pan, B_pan, C_pan, LA_pan, LA_panC, &
2226 : LB_pan, LB_panC, Lout_pan, MA_sub, &
2227 : MB_sub, tmp2, tmpA, tmpB
2228 : TYPE(dbcsr_type), POINTER :: RB_A, RB_B, RB_out
2229 : TYPE(dbcsr_type), TARGET :: LA_near, LB_near, Lout_near
2230 :
2231 654 : CALL timeset(routineN, handle)
2232 :
2233 654 : my_keep_sparsity = .FALSE.
2234 654 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
2235 654 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
2236 654 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
2237 654 : IF (PRESENT(grid_occupation)) grid_occupation = 0.0_dp
2238 :
2239 654 : CALL dbcsr_get_info(L_A, nfullrows_total=n_grid_total)
2240 :
2241 654 : IF (zero_out) CALL dbcsr_set(mat_out, 0.0_dp)
2242 :
2243 1330 : DO ipan = 1, SIZE(pan_first)
2244 676 : blk0 = pan_first(ipan)
2245 676 : blk1 = pan_last(ipan)
2246 :
2247 : ! phi/Z panel slices (P x n)
2248 676 : CALL extract_grid_panel(L_A, blk0, blk1, LA_pan)
2249 676 : IF (.NOT. lb_eq_la) CALL extract_grid_panel(L_B, blk0, blk1, LB_pan)
2250 676 : IF (.NOT. lout_eq_la) CALL extract_grid_panel(L_out, blk0, blk1, Lout_pan)
2251 :
2252 : ! Which AO/RI atoms (column blocks) actually touch this panel: the inner index of
2253 : ! every multiply below is restricted to them, so only the matching rows of the
2254 : ! system-wide operators M_A/M_B ever enter Cannon (exact: dropped rows meet zeros).
2255 676 : CALL collect_used_col_blocks(LA_pan, para_env, usedA)
2256 676 : IF (.NOT. lb_eq_la) THEN
2257 46 : CALL collect_used_col_blocks(LB_pan, para_env, usedB)
2258 : ELSE
2259 630 : IF (ALLOCATED(usedB)) DEALLOCATE (usedB)
2260 3380 : ALLOCATE (usedB, SOURCE=usedA)
2261 : END IF
2262 676 : IF (.NOT. (ANY(usedA) .AND. ANY(usedB))) THEN
2263 : ! empty panel slice: its Hadamard contribution is exactly zero on all ranks
2264 0 : CALL dbcsr_release(LA_pan)
2265 0 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_pan)
2266 0 : IF (.NOT. lout_eq_la) CALL dbcsr_release(Lout_pan)
2267 : CYCLE
2268 : END IF
2269 :
2270 : ! Grid rows within reach of the panel: with the CUTOFF_RADIUS_RL_W truncation only
2271 : ! they can appear as columns of the panel products / inner rows of the L_out multiply,
2272 : ! so the system-wide phi/Z right operands are cut down to this neighborhood slice
2273 : ! (local extraction, zero communication; exact w.r.t. the geo template).
2274 676 : IF (use_cutoff) THEN
2275 44 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, grid_used)
2276 44 : CALL extract_masked_blocks(L_A, grid_used, LA_near, compress_rows=.TRUE., blk_map=gmap)
2277 44 : IF (.NOT. lb_eq_la) CALL extract_masked_blocks(L_B, grid_used, LB_near, compress_rows=.TRUE.)
2278 44 : IF (.NOT. lout_eq_la) CALL extract_masked_blocks(L_out, grid_used, Lout_near, compress_rows=.TRUE.)
2279 : RB_A => LA_near
2280 : ELSE
2281 : RB_A => L_A
2282 : END IF
2283 676 : IF (lb_eq_la) THEN
2284 : RB_B => RB_A
2285 46 : ELSE IF (use_cutoff) THEN
2286 : RB_B => LB_near
2287 : ELSE
2288 42 : RB_B => L_B
2289 : END IF
2290 676 : IF (lout_eq_la) THEN
2291 : RB_out => RB_A
2292 630 : ELSE IF (use_cutoff) THEN
2293 : RB_out => Lout_near
2294 : ELSE
2295 590 : RB_out => L_out
2296 : END IF
2297 :
2298 : ! A_pan = LA_pan * M_A * L_A^T (P x grid_near).
2299 : ! When cutoff is active, A_pan is pre-seeded with only nearby blocks via
2300 : ! build_geo_template_panel, and the multiply uses retain_sparsity to skip
2301 : ! computing distant blocks entirely (exact: they are zero by locality).
2302 676 : CALL extract_masked_blocks(LA_pan, usedA, LA_panC, compress_rows=.FALSE.)
2303 676 : CALL extract_masked_blocks(M_A, usedA, MA_sub, compress_rows=.TRUE.)
2304 676 : CALL create_product_matrix(LA_panC, MA_sub, 'N', 'N', tmpA)
2305 676 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LA_panC, MA_sub, 0.0_dp, tmpA, filter_eps=eps)
2306 676 : CALL dbcsr_release(MA_sub)
2307 676 : IF (use_cutoff) THEN
2308 : CALL build_geo_template_panel(LA_pan, LA_near, centroids, cutoff, blk0, A_pan, &
2309 44 : col_map=gmap)
2310 : ELSE
2311 632 : CALL create_product_matrix(tmpA, RB_A, 'N', 'T', A_pan)
2312 : END IF
2313 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpA, RB_A, 0.0_dp, A_pan, &
2314 676 : filter_eps=eps, retain_sparsity=use_cutoff)
2315 676 : CALL dbcsr_release(tmpA)
2316 :
2317 : ! Grid-basis occupation of A_pan = φ G φ^T, accumulated over ALL panels into the
2318 : ! occupation of the (never formed) full grid x grid object:
2319 : ! sum_panels nnz(A_pan) / n_grid^2, with nnz = occ * pan_rows * pan_cols.
2320 : ! Panel-independent by construction -- a single-panel sample would instead report the
2321 : ! local neighbor count of whichever region happens to land in that panel.
2322 676 : IF (PRESENT(grid_occupation)) THEN
2323 630 : CALL dbcsr_get_info(A_pan, nfullrows_total=nrows_pan, nfullcols_total=ncols_pan)
2324 : grid_occupation = grid_occupation + dbcsr_get_occupation(A_pan)* &
2325 : REAL(ncols_pan, dp)*REAL(nrows_pan, dp)/ &
2326 630 : (REAL(n_grid_total, dp)*REAL(n_grid_total, dp))
2327 : END IF
2328 :
2329 : ! B_pan = LB_pan * M_B * L_B^T (P x grid_near); reuse the L_A slices when L_B == L_A.
2330 : ! With keep_sparsity, B_pan is pre-populated with A_pan's block structure so that
2331 : ! retain_sparsity forces the final multiply to fill only those blocks (exact for ∘).
2332 676 : IF (lb_eq_la) THEN
2333 630 : CALL extract_masked_blocks(M_B, usedA, MB_sub, compress_rows=.TRUE.)
2334 630 : CALL create_product_matrix(LA_panC, MB_sub, 'N', 'N', tmpB)
2335 630 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LA_panC, MB_sub, 0.0_dp, tmpB, filter_eps=eps)
2336 : ELSE
2337 46 : CALL extract_masked_blocks(LB_pan, usedB, LB_panC, compress_rows=.FALSE.)
2338 46 : CALL extract_masked_blocks(M_B, usedB, MB_sub, compress_rows=.TRUE.)
2339 46 : CALL create_product_matrix(LB_panC, MB_sub, 'N', 'N', tmpB)
2340 46 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LB_panC, MB_sub, 0.0_dp, tmpB, filter_eps=eps)
2341 46 : CALL dbcsr_release(LB_panC)
2342 : END IF
2343 676 : CALL dbcsr_release(MB_sub)
2344 676 : IF (my_keep_sparsity) THEN
2345 676 : CALL dbcsr_create(B_pan, template=A_pan)
2346 676 : CALL dbcsr_copy(B_pan, A_pan)
2347 676 : CALL dbcsr_set(B_pan, 0.0_dp)
2348 : ! The output pattern is already fixed. Omitting redundant on-the-fly filtering also
2349 : ! avoids overflowing DBCSR's single-precision screening norms for conditioned Z fits.
2350 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_B, 0.0_dp, B_pan, &
2351 676 : retain_sparsity=.TRUE.)
2352 : ELSE
2353 0 : CALL create_product_matrix(tmpB, RB_B, 'N', 'T', B_pan)
2354 0 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_B, 0.0_dp, B_pan, filter_eps=eps)
2355 : END IF
2356 676 : CALL dbcsr_release(tmpB)
2357 676 : CALL dbcsr_release(LA_panC)
2358 :
2359 : ! C_pan = scale * (A_pan ∘ B_pan) (P x grid_near)
2360 676 : CALL dbcsr_create(C_pan, template=A_pan)
2361 676 : CALL hadamard_product(A_pan, B_pan, C_pan, scale)
2362 676 : CALL dbcsr_release(A_pan)
2363 676 : CALL dbcsr_release(B_pan)
2364 :
2365 : ! tmp2 = C_pan * L_out (P x n_out; inner index restricted to the neighborhood)
2366 676 : CALL create_product_matrix(C_pan, RB_out, 'N', 'N', tmp2)
2367 676 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_pan, RB_out, 0.0_dp, tmp2, filter_eps=eps)
2368 676 : CALL dbcsr_release(C_pan)
2369 :
2370 : ! mat_out += L_out_pan^T * tmp2 (accumulate: beta = 1)
2371 676 : IF (lout_eq_la) THEN
2372 46 : CALL dbcsr_multiply('T', 'N', 1.0_dp, LA_pan, tmp2, 1.0_dp, mat_out, filter_eps=eps)
2373 : ELSE
2374 630 : CALL dbcsr_multiply('T', 'N', 1.0_dp, Lout_pan, tmp2, 1.0_dp, mat_out, filter_eps=eps)
2375 630 : CALL dbcsr_release(Lout_pan)
2376 : END IF
2377 676 : CALL dbcsr_release(tmp2)
2378 676 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_pan)
2379 676 : CALL dbcsr_release(LA_pan)
2380 3358 : IF (use_cutoff) THEN
2381 44 : CALL dbcsr_release(LA_near)
2382 44 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_near)
2383 44 : IF (.NOT. lout_eq_la) CALL dbcsr_release(Lout_near)
2384 : END IF
2385 :
2386 : END DO
2387 :
2388 654 : CALL timestop(handle)
2389 :
2390 1308 : END SUBROUTINE contract_grid_panels
2391 :
2392 : ! **************************************************************************************************
2393 : !> \brief Σ^c-specific panel loop: computes both the occupied (neg) and virtual (pos) contributions
2394 : !> in a single pass over grid panels, forming W_pan = Z_panel × W_aux × Z^T only ONCE per
2395 : !> panel and reusing it for both the G^occ and G^vir Hadamard contractions.
2396 : !>
2397 : !> Computes:
2398 : !> mat_Sigma_neg = φ^T ( (φ G^occ φ^T) ∘ (Z W^MIC Z^T) ) φ
2399 : !> mat_Sigma_pos = φ^T ( (φ G^vir φ^T) ∘ (Z W^MIC Z^T) ) φ
2400 : !>
2401 : !> \param mat_phi ...
2402 : !> \param mat_Z ...
2403 : !> \param mat_G_occ_ao ...
2404 : !> \param mat_G_vir_ao ...
2405 : !> \param mat_W_aux ...
2406 : !> \param mat_Sigma_neg ...
2407 : !> \param mat_Sigma_pos ...
2408 : !> \param eps ...
2409 : !> \param para_env ...
2410 : !> \param pan_first ...
2411 : !> \param pan_last ...
2412 : !> \param keep_sparsity ...
2413 : !> \param centroids ...
2414 : !> \param cutoff ...
2415 : ! **************************************************************************************************
2416 810 : SUBROUTINE contract_grid_panels_sigma_c(mat_phi, mat_Z, mat_G_occ_ao, mat_G_vir_ao, &
2417 : mat_W_aux, mat_Sigma_neg, mat_Sigma_pos, eps, para_env, &
2418 810 : pan_first, pan_last, keep_sparsity, centroids, cutoff)
2419 :
2420 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: mat_phi, mat_Z
2421 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_G_occ_ao, mat_G_vir_ao, mat_W_aux, &
2422 : mat_Sigma_neg, mat_Sigma_pos
2423 : REAL(KIND=dp), INTENT(IN) :: eps
2424 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
2425 : INTEGER, DIMENSION(:), INTENT(IN) :: pan_first, pan_last
2426 : LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
2427 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
2428 : OPTIONAL :: centroids
2429 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
2430 :
2431 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_grid_panels_sigma_c'
2432 :
2433 : INTEGER :: blk0, blk1, handle, ipan
2434 810 : INTEGER, ALLOCATABLE, DIMENSION(:) :: gmap
2435 : LOGICAL :: my_keep_sparsity, use_cutoff
2436 810 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used, used_ao, used_ri
2437 : TYPE(dbcsr_type) :: A_occ, A_vir, C_pan, G_occ_sub, &
2438 : G_vir_sub, phi_pan, phi_panC, tmp2, &
2439 : tmpA, tmpB, W_pan, W_sub, Z_pan, Z_panC
2440 : TYPE(dbcsr_type), POINTER :: RB_phi, RB_Z
2441 : TYPE(dbcsr_type), TARGET :: phi_near, Z_near
2442 :
2443 810 : CALL timeset(routineN, handle)
2444 :
2445 810 : my_keep_sparsity = .FALSE.
2446 810 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
2447 810 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
2448 810 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
2449 :
2450 810 : CALL dbcsr_set(mat_Sigma_neg, 0.0_dp)
2451 810 : CALL dbcsr_set(mat_Sigma_pos, 0.0_dp)
2452 :
2453 1640 : DO ipan = 1, SIZE(pan_first)
2454 830 : blk0 = pan_first(ipan)
2455 830 : blk1 = pan_last(ipan)
2456 :
2457 830 : CALL extract_grid_panel(mat_phi, blk0, blk1, phi_pan)
2458 830 : CALL extract_grid_panel(mat_Z, blk0, blk1, Z_pan)
2459 :
2460 : ! AO/RI atoms touching this panel: only the matching rows of G_occ/G_vir/W ever
2461 : ! enter the multiplies below (exact: dropped rows meet zero columns of the panel).
2462 830 : CALL collect_used_col_blocks(phi_pan, para_env, used_ao)
2463 830 : CALL collect_used_col_blocks(Z_pan, para_env, used_ri)
2464 830 : IF (.NOT. (ANY(used_ao) .AND. ANY(used_ri))) THEN
2465 0 : CALL dbcsr_release(phi_pan)
2466 0 : CALL dbcsr_release(Z_pan)
2467 0 : CYCLE
2468 : END IF
2469 :
2470 : ! Neighborhood slices of phi/Z (grid rows within cutoff of the panel): they replace
2471 : ! the system-wide right operands in every multiply (local extraction, zero comm).
2472 830 : IF (use_cutoff) THEN
2473 40 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, grid_used)
2474 40 : CALL extract_masked_blocks(mat_phi, grid_used, phi_near, compress_rows=.TRUE., blk_map=gmap)
2475 40 : CALL extract_masked_blocks(mat_Z, grid_used, Z_near, compress_rows=.TRUE.)
2476 40 : RB_phi => phi_near
2477 40 : RB_Z => Z_near
2478 : ELSE
2479 : RB_phi => mat_phi
2480 : RB_Z => mat_Z
2481 : END IF
2482 :
2483 830 : CALL extract_masked_blocks(phi_pan, used_ao, phi_panC, compress_rows=.FALSE.)
2484 830 : CALL extract_masked_blocks(Z_pan, used_ri, Z_panC, compress_rows=.FALSE.)
2485 830 : CALL extract_masked_blocks(mat_G_occ_ao, used_ao, G_occ_sub, compress_rows=.TRUE.)
2486 830 : CALL extract_masked_blocks(mat_G_vir_ao, used_ao, G_vir_sub, compress_rows=.TRUE.)
2487 830 : CALL extract_masked_blocks(mat_W_aux, used_ri, W_sub, compress_rows=.TRUE.)
2488 :
2489 : ! A_occ = phi_pan × G_occ × phi^T (built first so W_pan can inherit its pattern).
2490 : ! With cutoff active, A_occ is pre-seeded with geo-local blocks so that the
2491 : ! phi^T multiply uses retain_sparsity and never computes distant blocks.
2492 830 : CALL create_product_matrix(phi_panC, G_occ_sub, 'N', 'N', tmpA)
2493 830 : CALL dbcsr_multiply('N', 'N', 1.0_dp, phi_panC, G_occ_sub, 0.0_dp, tmpA, filter_eps=eps)
2494 830 : IF (use_cutoff) THEN
2495 : CALL build_geo_template_panel(phi_pan, phi_near, centroids, cutoff, blk0, A_occ, &
2496 40 : col_map=gmap)
2497 : ELSE
2498 790 : CALL create_product_matrix(tmpA, RB_phi, 'N', 'T', A_occ)
2499 : END IF
2500 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpA, RB_phi, 0.0_dp, A_occ, &
2501 830 : filter_eps=eps, retain_sparsity=use_cutoff)
2502 830 : CALL dbcsr_release(tmpA)
2503 830 : CALL dbcsr_release(G_occ_sub)
2504 :
2505 : ! A_vir = phi_pan × G_vir × phi^T (same pre-screen as A_occ)
2506 830 : CALL create_product_matrix(phi_panC, G_vir_sub, 'N', 'N', tmpA)
2507 830 : CALL dbcsr_multiply('N', 'N', 1.0_dp, phi_panC, G_vir_sub, 0.0_dp, tmpA, filter_eps=eps)
2508 830 : IF (use_cutoff) THEN
2509 : CALL build_geo_template_panel(phi_pan, phi_near, centroids, cutoff, blk0, A_vir, &
2510 40 : col_map=gmap)
2511 : ELSE
2512 790 : CALL create_product_matrix(tmpA, RB_phi, 'N', 'T', A_vir)
2513 : END IF
2514 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpA, RB_phi, 0.0_dp, A_vir, &
2515 830 : filter_eps=eps, retain_sparsity=use_cutoff)
2516 830 : CALL dbcsr_release(tmpA)
2517 830 : CALL dbcsr_release(G_vir_sub)
2518 :
2519 : ! W_pan = Z_pan × W_aux × Z^T (computed once, reused for both Σ^c terms).
2520 : ! With keep_sparsity, W_pan is pre-seeded with the union of A_occ and A_vir block
2521 : ! patterns so that retain_sparsity forces the multiply to fill only those blocks:
2522 : ! exact since W outside G_occ∪G_vir is multiplied by zero in the Hadamard.
2523 830 : CALL create_product_matrix(Z_panC, W_sub, 'N', 'N', tmpB)
2524 830 : CALL dbcsr_multiply('N', 'N', 1.0_dp, Z_panC, W_sub, 0.0_dp, tmpB, filter_eps=eps)
2525 830 : IF (my_keep_sparsity) THEN
2526 830 : CALL dbcsr_create(W_pan, template=A_occ)
2527 830 : CALL dbcsr_copy(W_pan, A_occ)
2528 830 : CALL dbcsr_add(W_pan, A_vir, 1.0_dp, 1.0_dp)
2529 830 : CALL dbcsr_set(W_pan, 0.0_dp)
2530 : ! The retained pattern is exact for the following Hadamard products; screening it
2531 : ! again is redundant and may overflow DBCSR's single-precision block-norm product.
2532 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_Z, 0.0_dp, W_pan, &
2533 830 : retain_sparsity=.TRUE.)
2534 : ELSE
2535 0 : CALL create_product_matrix(tmpB, RB_Z, 'N', 'T', W_pan)
2536 0 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_Z, 0.0_dp, W_pan, filter_eps=eps)
2537 : END IF
2538 830 : CALL dbcsr_release(tmpB)
2539 830 : CALL dbcsr_release(W_sub)
2540 830 : CALL dbcsr_release(phi_panC)
2541 830 : CALL dbcsr_release(Z_panC)
2542 :
2543 : ! Σ^c_neg: φ^T ( A_occ ∘ W_pan ) φ
2544 830 : CALL dbcsr_create(C_pan, template=A_occ)
2545 830 : CALL hadamard_product(A_occ, W_pan, C_pan, 1.0_dp)
2546 830 : CALL dbcsr_release(A_occ)
2547 830 : CALL create_product_matrix(C_pan, RB_phi, 'N', 'N', tmp2)
2548 830 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_pan, RB_phi, 0.0_dp, tmp2, filter_eps=eps)
2549 830 : CALL dbcsr_release(C_pan)
2550 830 : CALL dbcsr_multiply('T', 'N', 1.0_dp, phi_pan, tmp2, 1.0_dp, mat_Sigma_neg, filter_eps=eps)
2551 830 : CALL dbcsr_release(tmp2)
2552 :
2553 : ! Σ^c_pos: φ^T ( A_vir ∘ W_pan ) φ — W_pan reused
2554 830 : CALL dbcsr_create(C_pan, template=A_vir)
2555 830 : CALL hadamard_product(A_vir, W_pan, C_pan, 1.0_dp)
2556 830 : CALL dbcsr_release(A_vir)
2557 830 : CALL create_product_matrix(C_pan, RB_phi, 'N', 'N', tmp2)
2558 830 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_pan, RB_phi, 0.0_dp, tmp2, filter_eps=eps)
2559 830 : CALL dbcsr_release(C_pan)
2560 830 : CALL dbcsr_multiply('T', 'N', 1.0_dp, phi_pan, tmp2, 1.0_dp, mat_Sigma_pos, filter_eps=eps)
2561 830 : CALL dbcsr_release(tmp2)
2562 :
2563 830 : CALL dbcsr_release(W_pan)
2564 830 : CALL dbcsr_release(Z_pan)
2565 830 : CALL dbcsr_release(phi_pan)
2566 6620 : IF (use_cutoff) THEN
2567 40 : CALL dbcsr_release(phi_near)
2568 40 : CALL dbcsr_release(Z_near)
2569 : END IF
2570 :
2571 : END DO
2572 :
2573 810 : CALL timestop(handle)
2574 :
2575 1620 : END SUBROUTINE contract_grid_panels_sigma_c
2576 :
2577 : ! **************************************************************************************************
2578 : !> \brief Computes C = A ◦ B (Element-wise Hadamard product) for sparse DBCSR matrices.
2579 : !> \param matrix_A ...
2580 : !> \param matrix_B ...
2581 : !> \param matrix_C ...
2582 : !> \param fac (Scaling factor applied to the product)
2583 : ! **************************************************************************************************
2584 4672 : SUBROUTINE hadamard_product(matrix_A, matrix_B, matrix_C, fac)
2585 :
2586 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B, matrix_C
2587 : REAL(KIND=dp), INTENT(IN) :: fac
2588 :
2589 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hadamard_product'
2590 :
2591 : INTEGER :: col, handle, row
2592 : LOGICAL :: found
2593 2336 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk_B, blk_C
2594 : TYPE(dbcsr_iterator_type) :: iter
2595 :
2596 2336 : CALL timeset(routineN, handle)
2597 :
2598 2336 : CALL dbcsr_copy(matrix_C, matrix_A)
2599 :
2600 2336 : CALL dbcsr_iterator_start(iter, matrix_C)
2601 298394 : DO WHILE (dbcsr_iterator_blocks_left(iter))
2602 296058 : CALL dbcsr_iterator_next_block(iter, row, col, blk_C)
2603 :
2604 296058 : CALL dbcsr_get_block_p(matrix_B, row, col, blk_B, found)
2605 :
2606 298394 : IF (found) THEN
2607 527574254 : blk_C(:, :) = fac*blk_C(:, :)*blk_B(:, :)
2608 : ELSE
2609 : ! If B is sparse here, the product is zero
2610 0 : blk_C(:, :) = 0.0_dp
2611 : END IF
2612 : END DO
2613 2336 : CALL dbcsr_iterator_stop(iter)
2614 :
2615 2336 : CALL timestop(handle)
2616 :
2617 2336 : END SUBROUTINE hadamard_product
2618 :
2619 : ! **************************************************************************************************
2620 : !> \brief Computes the screened Coulomb interaction on the imaginary-time grid, entirely in the
2621 : !> RI auxiliary (PQ) basis:
2622 : !> χ_PQ(iω) = Σ_τ w(ω,τ) cos(ωτ) χ_PQ(iτ) (cosine transform)
2623 : !> ε(iω) = Id - V^0.5 M^-1 χ(iω) M^-1 V^0.5 (dielectric function)
2624 : !> W(iω) = V^0.5 ( ε^-1(iω) - Id ) V^0.5 (correlation part only)
2625 : !> W(iτ) = Σ_ω w̃(τ,ω) cos(ωτ) W(iω) (back transform)
2626 : !> W(iτ) <- M^-1 W(iτ) M^-1 (fold in the RI metric)
2627 : !> where V is the bare Coulomb matrix and M the RI metric.
2628 : !> \param bs_env ...
2629 : !> \param qs_env ...
2630 : !> \param mat_chi_Gamma_tau ...
2631 : !> \param fm_W_time ...
2632 : ! **************************************************************************************************
2633 152 : SUBROUTINE compute_W(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_time)
2634 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2635 : TYPE(qs_environment_type), POINTER :: qs_env
2636 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
2637 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
2638 :
2639 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_W'
2640 :
2641 : INTEGER :: handle, i_t, j_w
2642 : REAL(KIND=dp) :: t1
2643 : TYPE(cp_fm_type) :: fm_M_inv_V_sqrt, fm_V, fm_V_sqrt
2644 :
2645 38 : CALL timeset(routineN, handle)
2646 :
2647 38 : t1 = m_walltime()
2648 :
2649 38 : CALL create_fm_W_MIC_time(bs_env, fm_W_time)
2650 :
2651 : ! 1. Allocate V and M matrices
2652 38 : CALL cp_fm_create(fm_V, bs_env%fm_RI_RI%matrix_struct)
2653 38 : CALL cp_fm_create(fm_V_sqrt, bs_env%fm_RI_RI%matrix_struct)
2654 38 : CALL cp_fm_create(fm_M_inv_V_sqrt, bs_env%fm_RI_RI%matrix_struct)
2655 :
2656 : ! Compute V and M^-1 * V^0.5
2657 38 : CALL compute_V_MinvVsqrt(bs_env, qs_env, fm_V, fm_V_sqrt, fm_M_inv_V_sqrt)
2658 :
2659 : ! 2. Loop over frequencies
2660 568 : DO j_w = 1, bs_env%num_time_freq_points
2661 : ! Fourier transformation of χ_PQ(iτ) to χ_PQ(iω_j)
2662 530 : CALL compute_fm_chi_Gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
2663 :
2664 : ! ε(iω_j) = Id - V^0.5*M^-1*χ(iω_j)*M^-1*V^0.5
2665 : ! W(iω_j) = V^0.5*(ε^-1(iω_j)-Id)*V^0.5
2666 : CALL compute_fm_W_freq(bs_env, bs_env%fm_chi_Gamma_freq, fm_V_sqrt, &
2667 530 : fm_M_inv_V_sqrt, bs_env%fm_W_MIC_freq)
2668 :
2669 : ! Fourier transform from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
2670 568 : CALL Fourier_transform_w_to_t(bs_env, fm_W_time, bs_env%fm_W_MIC_freq, j_w)
2671 : END DO
2672 :
2673 : ! M^-1(k=0) W^MIC(iτ) M^-1(k=0) -> fm_W_time
2674 38 : CALL fm_contract_ABA(bs_env%fm_Minv_Gamma, fm_W_time)
2675 :
2676 38 : IF (bs_env%unit_nr > 0) THEN
2677 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
2678 19 : 'Computed W(iτ),', ' Execution time', m_walltime() - t1, ' s'
2679 : END IF
2680 :
2681 38 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
2682 :
2683 : ! Cleanup
2684 38 : CALL cp_fm_release(fm_V)
2685 38 : CALL cp_fm_release(fm_V_sqrt)
2686 38 : CALL cp_fm_release(fm_M_inv_V_sqrt)
2687 :
2688 : ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
2689 38 : IF (bs_env%rtp_method == rtp_method_bse) THEN
2690 14 : t1 = m_walltime()
2691 14 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
2692 : ! Set to zero
2693 14 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
2694 : ! Sum over all times
2695 294 : DO i_t = 1, bs_env%num_time_freq_points
2696 : ! Add the relevant structure with correct weight
2697 : CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
2698 294 : bs_env%time_frequency_grid%time_weights_at_zero_frequency(i_t), fm_W_time(i_t))
2699 : END DO
2700 : ! Done, save to file
2701 14 : CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
2702 : ! Report calculation
2703 14 : IF (bs_env%unit_nr > 0) THEN
2704 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
2705 7 : 'Computed W(0),', ' Execution time', m_walltime() - t1, ' s'
2706 : END IF
2707 : END IF
2708 :
2709 38 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
2710 :
2711 38 : CALL timestop(handle)
2712 :
2713 38 : END SUBROUTINE compute_W
2714 :
2715 : ! **************************************************************************************************
2716 : !> \brief Computes V, V^0.5, and M^-1 V^0.5 for the RI-RS dielectric function.
2717 : !> The Coulomb matrix V is constructed by the RI-RS k-point path. The inverse metric
2718 : !> M^-1(k=0) is precomputed once in gw_utils and stored in bs_env.
2719 : !> \param bs_env GW environment containing M^-1(k=0) and the RI basis
2720 : !> \param qs_env electronic-structure environment used to evaluate V(k=0)
2721 : !> \param fm_V Coulomb matrix V(k=0)
2722 : !> \param fm_V_sqrt symmetric factor V^0.5
2723 : !> \param fm_Minv_Vsqrt product M^-1(k=0) V^0.5
2724 : ! **************************************************************************************************
2725 38 : SUBROUTINE compute_V_MinvVsqrt(bs_env, qs_env, fm_V, fm_V_sqrt, fm_Minv_Vsqrt)
2726 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2727 : TYPE(qs_environment_type), POINTER :: qs_env
2728 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_V, fm_V_sqrt, fm_Minv_Vsqrt
2729 :
2730 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_MinvVsqrt'
2731 :
2732 : INTEGER :: handle
2733 38 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2734 : TYPE(cell_type), POINTER :: cell
2735 38 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_kp
2736 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2737 38 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2738 :
2739 38 : CALL timeset(routineN, handle)
2740 :
2741 38 : IF (bs_env%auto_ri%enabled) THEN
2742 : ! -------------------------------------------------------------------
2743 : ! 1a. The optimized AB functions span two atoms, so their previously transformed
2744 : ! Coulomb matrix cannot be rebuilt by the atom-local k-point integral routine.
2745 : ! -------------------------------------------------------------------
2746 2 : CALL cp_fm_to_fm(bs_env%auto_ri%V_pq, fm_V)
2747 : ELSE
2748 : ! -------------------------------------------------------------------
2749 : ! 1b. Build Coulomb Matrix V(k=0) using the kp-routine but only for ikp=1.
2750 : ! -------------------------------------------------------------------
2751 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, &
2752 36 : qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
2753 :
2754 180 : ALLOCATE (mat_V_kp(1:1, 1:2))
2755 36 : NULLIFY (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
2756 36 : ALLOCATE (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
2757 36 : CALL dbcsr_create(mat_V_kp(1, 1)%matrix, template=bs_env%mat_RI_RI%matrix)
2758 36 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 1)%matrix)
2759 36 : CALL dbcsr_set(mat_V_kp(1, 1)%matrix, 0.0_dp)
2760 36 : CALL dbcsr_create(mat_V_kp(1, 2)%matrix, template=bs_env%mat_RI_RI%matrix)
2761 36 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 2)%matrix)
2762 : ! The dummy imaginary part is required only by the k-point routine interface.
2763 36 : CALL dbcsr_set(mat_V_kp(1, 2)%matrix, 0.0_dp)
2764 :
2765 144 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
2766 : CALL build_2c_coulomb_matrix_kp(mat_V_kp, bs_env%kpoints_chi_eps_W, "RI_AUX", cell, &
2767 : particle_set, qs_kind_set, atomic_kind_set, &
2768 36 : bs_env%size_lattice_sum_V, operator_coulomb, 1, 1)
2769 36 : CALL copy_dbcsr_to_fm(mat_V_kp(1, 1)%matrix, fm_V)
2770 :
2771 36 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 1)%matrix)
2772 36 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 2)%matrix)
2773 36 : DEALLOCATE (mat_V_kp)
2774 : END IF
2775 :
2776 : ! -----------------------------------------------------------------------
2777 : ! 2. V -> V^0.5.
2778 : ! -----------------------------------------------------------------------
2779 38 : CALL fm_sqrt(fm_V, fm_V_sqrt, bs_env%eps_eigval_mat_RI, bs_env%unit_nr)
2780 :
2781 : ! -----------------------------------------------------------------------
2782 : ! 3. M^-1(k=0) V^0.5.
2783 : ! -----------------------------------------------------------------------
2784 : CALL parallel_gemm("N", "T", bs_env%n_RI, bs_env%n_RI, bs_env%n_RI, 1.0_dp, &
2785 38 : bs_env%fm_Minv_Gamma, fm_V_sqrt, 0.0_dp, fm_Minv_Vsqrt)
2786 :
2787 38 : CALL timestop(handle)
2788 :
2789 38 : END SUBROUTINE compute_V_MinvVsqrt
2790 :
2791 : ! **************************************************************************************************
2792 : !> \brief Computes the screened interaction at one imaginary frequency:
2793 : !> ε(iω_j) = Id - (M^-1 V^0.5)^T χ(iω_j) (M^-1 V^0.5)
2794 : !> W(iω_j) = V^0.5^T ( ε^-1(iω_j) - Id ) V^0.5
2795 : !> ε is inverted via Cholesky; if that fails due to conditioning, via
2796 : !> eigendecomposition (cp_fm_power) with eigenvalue filtering.
2797 : !> \param bs_env ...
2798 : !> \param fm_chi_freq_j ...
2799 : !> \param fm_V_sqrt ...
2800 : !> \param fm_Minv_Vsqrt ...
2801 : !> \param fm_W_freq_j ...
2802 : ! **************************************************************************************************
2803 1590 : SUBROUTINE compute_fm_W_freq(bs_env, fm_chi_freq_j, fm_V_sqrt, fm_Minv_Vsqrt, fm_W_freq_j)
2804 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2805 : TYPE(cp_fm_type), INTENT(IN) :: fm_chi_freq_j, fm_V_sqrt, fm_Minv_Vsqrt
2806 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_W_freq_j
2807 :
2808 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_freq'
2809 :
2810 : INTEGER :: handle, n_RI
2811 : TYPE(cp_fm_type) :: fm_eps_freq_j, fm_work
2812 :
2813 530 : CALL timeset(routineN, handle)
2814 :
2815 530 : n_RI = bs_env%n_RI
2816 :
2817 530 : CALL cp_fm_create(fm_eps_freq_j, fm_chi_freq_j%matrix_struct)
2818 530 : CALL cp_fm_create(fm_work, fm_chi_freq_j%matrix_struct)
2819 :
2820 : ! -----------------------------------------------------------------------
2821 : ! 1. ε(iω_j) = Id - (M^-1 * V^0.5)^T * χ(iω_j) * (M^-1 * V^0.5)
2822 : ! -----------------------------------------------------------------------
2823 : ! work = χ(iω_j) * (M^-1 * V^0.5)
2824 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
2825 530 : fm_chi_freq_j, fm_Minv_Vsqrt, 0.0_dp, fm_work)
2826 :
2827 : ! eps_work = (M^-1 * V^0.5)^T * work
2828 : CALL parallel_gemm('T', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
2829 530 : fm_Minv_Vsqrt, fm_work, 0.0_dp, fm_eps_freq_j)
2830 :
2831 : ! ε(iω_j) = Id - eps_work --> -eps_work + Id
2832 530 : CALL fm_add_on_diag(fm_eps_freq_j, 1.0_dp)
2833 :
2834 : ! Force perfect symmetry before Cholesky to avoid info != 0 due to GEMM noise
2835 530 : CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
2836 :
2837 : ! -----------------------------------------------------------------------
2838 : ! 2. W(iω_j) = V^0.5^T * (ε^-1(iω_j) - Id) * V^0.5
2839 : ! -----------------------------------------------------------------------
2840 :
2841 : ! a) Invert ε by Cholesky decomposition or, if that fails, by diagonalization.
2842 530 : CALL fm_invert(fm_eps_freq_j, bs_env%eps_eigval_mat_RI, bs_env%unit_nr)
2843 :
2844 : ! b) ε^-1(iω_j) - Id
2845 530 : CALL fm_add_on_diag(fm_eps_freq_j, -1.0_dp)
2846 :
2847 : ! c) W(iω_j) = V^0.5^T * (ε^-1(iω_j) - Id) * V^0.5
2848 530 : CALL fm_contract_ABA(fm_V_sqrt, fm_eps_freq_j, fm_W_freq_j)
2849 :
2850 : ! Cleanup
2851 530 : CALL cp_fm_release(fm_work)
2852 530 : CALL cp_fm_release(fm_eps_freq_j)
2853 :
2854 530 : CALL timestop(handle)
2855 :
2856 530 : END SUBROUTINE compute_fm_W_freq
2857 :
2858 : ! **************************************************************************************************
2859 : !> \brief Adds a real scalar value to the diagonal of a real full matrix
2860 : !> \param fm ...
2861 : !> \param alpha ...
2862 : ! **************************************************************************************************
2863 2120 : SUBROUTINE fm_add_on_diag(fm, alpha)
2864 : TYPE(cp_fm_type), INTENT(INOUT) :: fm
2865 : REAL(KIND=dp), INTENT(IN) :: alpha
2866 :
2867 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_add_on_diag'
2868 :
2869 : INTEGER :: handle, i_global, i_row, j_col, &
2870 : j_global, ncol_local, nrow_local
2871 1060 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2872 :
2873 1060 : CALL timeset(routineN, handle)
2874 :
2875 : CALL cp_fm_get_info(matrix=fm, &
2876 : nrow_local=nrow_local, &
2877 : ncol_local=ncol_local, &
2878 : row_indices=row_indices, &
2879 1060 : col_indices=col_indices)
2880 :
2881 37280 : DO j_col = 1, ncol_local
2882 36220 : j_global = col_indices(j_col)
2883 739390 : DO i_row = 1, nrow_local
2884 702110 : i_global = row_indices(i_row)
2885 738330 : IF (j_global == i_global) THEN
2886 18110 : fm%local_data(i_row, j_col) = fm%local_data(i_row, j_col) + alpha
2887 : END IF
2888 : END DO
2889 : END DO
2890 :
2891 1060 : CALL timestop(handle)
2892 :
2893 1060 : END SUBROUTINE fm_add_on_diag
2894 :
2895 : ! **************************************************************************************************
2896 : !> \brief Computes the exact-exchange part of the GW self-energy:
2897 : !> D_μν = Σ_n^occ C_μn C_νn (density matrix = G^occ at τ=0)
2898 : !> V^tr_PQ = M^-1 (P|Q)_trunc M^-1 (truncated Coulomb, RI basis)
2899 : !> Σ^x_λσ(k=0) = -Σ_ll' Φ_λ(r_l) [ (φ D φ^T)_ll' ∘ (Z V^tr Z^T)_ll' ] Φ_σ(r_l')
2900 : !> \param bs_env ...
2901 : !> \param qs_env ...
2902 : !> \param mat_phi_mu_l ...
2903 : !> \param mat_Z_lP ...
2904 : !> \param fm_Sigma_x_Gamma ...
2905 : ! **************************************************************************************************
2906 38 : SUBROUTINE compute_Sigma_x(bs_env, qs_env, mat_phi_mu_l, mat_Z_lP, fm_Sigma_x_Gamma)
2907 :
2908 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2909 : TYPE(qs_environment_type), POINTER :: qs_env
2910 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
2911 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
2912 :
2913 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_x'
2914 :
2915 : INTEGER :: handle, ispin
2916 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
2917 38 : INTEGER, DIMENSION(:), POINTER :: blk_aux, dist_row_aux
2918 : REAL(KIND=dp) :: t1
2919 38 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
2920 : TYPE(dbcsr_distribution_type) :: dist_aux_aux
2921 : TYPE(dbcsr_type) :: mat_Sigma_x_Gamma, matrix_D_ao, &
2922 : matrix_V_aux
2923 :
2924 38 : CALL timeset(routineN, handle)
2925 :
2926 38 : t1 = m_walltime()
2927 :
2928 158 : ALLOCATE (fm_Sigma_x_Gamma(bs_env%n_spin))
2929 82 : DO ispin = 1, bs_env%n_spin
2930 82 : CALL cp_fm_create(fm_Sigma_x_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
2931 : END DO
2932 :
2933 38 : CALL dbcsr_create(mat_Sigma_x_Gamma, template=bs_env%mat_ao_ao%matrix)
2934 :
2935 38 : CALL resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
2936 :
2937 : ! =========================================================================
2938 : ! 1. COMPUTE V^tr_PQ (RI x RI)
2939 : ! =========================================================================
2940 38 : CALL setup_square_topology(mat_Z_lP, dist_aux_aux, blk_aux, dist_row_aux)
2941 :
2942 38 : IF (bs_env%auto_ri%enabled) THEN
2943 6 : ALLOCATE (fm_Vtr_Gamma(1, 1))
2944 2 : CALL cp_fm_create(fm_Vtr_Gamma(1, 1), bs_env%fm_RI_RI%matrix_struct)
2945 2 : CALL cp_fm_to_fm(bs_env%auto_ri%V_pq, fm_Vtr_Gamma(1, 1))
2946 : ELSE
2947 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
2948 36 : bs_env%trunc_coulomb, do_kpoints=.FALSE.)
2949 : END IF
2950 :
2951 : ! M^-1(k=0) V^tr(τ) M^-1(k=0) -> fm_Vtr_Gamma
2952 38 : CALL fm_contract_ABA(bs_env%fm_Minv_Gamma, fm_Vtr_Gamma(:, 1))
2953 :
2954 38 : CALL dbcsr_create(matrix_V_aux, "V_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
2955 : ! Optional CUTOFF_RADIUS_G_W operator truncation + filter (see build_G_ao)
2956 38 : IF (bs_env%ri_rs%cutoff_radius_g_w > 0.0_dp .AND. ALLOCATED(bs_env%ri_rs%atom_centers)) THEN
2957 : CALL reserve_blocks_within_radius(matrix_V_aux, bs_env%ri_rs%atom_centers, &
2958 0 : bs_env%ri_rs%cutoff_radius_g_w)
2959 0 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.TRUE.)
2960 : ELSE
2961 38 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.FALSE.)
2962 : END IF
2963 38 : CALL dbcsr_filter(matrix_V_aux, bs_env%eps_filter)
2964 :
2965 : ! =========================================================================
2966 : ! 2. SPIN LOOP FOR EXACT EXCHANGE
2967 : ! Σ^x_λσ = -Σ_ll' Φ_λ(r_l) ( D_ll' V^tr_ll' ) Φ_σ(r_l')
2968 : ! = -φ^T ( (φ D φ^T) ∘ (Z V^tr Z^T) ) φ
2969 : ! =========================================================================
2970 82 : DO ispin = 1, bs_env%n_spin
2971 :
2972 : ! AO-space density matrix D_µν = G^occ at τ = 0
2973 44 : CALL build_G_ao(bs_env, 0.0_dp, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_D_ao)
2974 :
2975 : CALL contract_grid_panels(L_A=mat_phi_mu_l, M_A=matrix_D_ao, &
2976 : L_B=mat_Z_lP, M_B=matrix_V_aux, &
2977 : L_out=mat_phi_mu_l, mat_out=mat_Sigma_x_Gamma, &
2978 : scale=1.0_dp, eps=bs_env%eps_filter, &
2979 : para_env=bs_env%para_env, &
2980 : pan_first=pan_first, pan_last=pan_last, &
2981 : lb_eq_la=.FALSE., lout_eq_la=.TRUE., zero_out=.TRUE., &
2982 : keep_sparsity=bs_env%ri_rs%keep_sparsity_rirs, &
2983 : centroids=bs_env%ri_rs%chunk_centroids, &
2984 44 : cutoff=bs_env%ri_rs%cutoff_radius_v_w)
2985 44 : CALL dbcsr_scale(mat_Sigma_x_Gamma, -1.0_dp)
2986 :
2987 44 : CALL dbcsr_release(matrix_D_ao)
2988 :
2989 : ! Data I/O and Export to CP2K Full Matrices
2990 82 : CALL copy_dbcsr_to_fm(mat_Sigma_x_Gamma, fm_Sigma_x_Gamma(ispin))
2991 :
2992 : END DO ! ispin
2993 :
2994 38 : IF (bs_env%unit_nr > 0) THEN
2995 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
2996 19 : 'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
2997 19 : WRITE (bs_env%unit_nr, '(A)') ' '
2998 : END IF
2999 :
3000 : ! =========================================================================
3001 : ! 3. CLEANUP
3002 : ! =========================================================================
3003 38 : CALL dbcsr_release(matrix_V_aux)
3004 38 : CALL dbcsr_release(mat_Sigma_x_Gamma)
3005 38 : CALL release_square_topology(dist=dist_aux_aux, mapped_dist=dist_row_aux)
3006 :
3007 38 : CALL cp_fm_release(fm_Vtr_Gamma)
3008 :
3009 38 : CALL timestop(handle)
3010 :
3011 114 : END SUBROUTINE compute_Sigma_x
3012 :
3013 : ! **************************************************************************************************
3014 : !> \brief Computes the correlation part of the GW self-energy on the imaginary-time grid:
3015 : !> Σ^c_λσ(iτ<0) = -Σ_ll' Φ_λ(r_l) [ (φ G^occ φ^T)_ll' ∘ (Z W^MIC Z^T)_ll' ] Φ_σ(r_l')
3016 : !> Σ^c_λσ(iτ>0) = +Σ_ll' Φ_λ(r_l) [ (φ G^vir φ^T)_ll' ∘ (Z W^MIC Z^T)_ll' ] Φ_σ(r_l')
3017 : !> \param bs_env ...
3018 : !> \param fm_W_time ...
3019 : !> \param mat_phi_mu_l ...
3020 : !> \param mat_Z_lP ...
3021 : !> \param fm_Sigma_c_Gamma_time ...
3022 : ! **************************************************************************************************
3023 52 : SUBROUTINE compute_Sigma_c(bs_env, fm_W_time, mat_phi_mu_l, mat_Z_lP, fm_Sigma_c_Gamma_time)
3024 :
3025 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3026 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
3027 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
3028 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
3029 :
3030 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c'
3031 :
3032 : INTEGER :: handle, i_t, ispin
3033 52 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
3034 52 : INTEGER, DIMENSION(:), POINTER :: blk_aux, dist_row_aux
3035 : REAL(KIND=dp) :: t1, tau
3036 : TYPE(dbcsr_distribution_type) :: dist_aux_aux
3037 52 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
3038 : TYPE(dbcsr_type) :: matrix_G_occ_ao, matrix_G_vir_ao, &
3039 : matrix_W_aux
3040 :
3041 52 : CALL timeset(routineN, handle)
3042 :
3043 : ! =========================================================================
3044 : ! 1. SETUP AUXILIARY TOPOLOGY AND PRE-ALLOCATE OUTPUT ARRAYS
3045 : ! =========================================================================
3046 52 : CALL setup_square_topology(mat_Z_lP, dist_aux_aux, blk_aux, dist_row_aux)
3047 :
3048 52 : CALL resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
3049 :
3050 : ! Pre-allocate local DBCSR matrices to act as targets for final output
3051 52 : NULLIFY (mat_Sigma_neg_tau, mat_Sigma_pos_tau)
3052 1082 : ALLOCATE (mat_Sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
3053 1082 : ALLOCATE (mat_Sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
3054 :
3055 722 : DO i_t = 1, bs_env%num_time_freq_points
3056 1532 : DO ispin = 1, bs_env%n_spin
3057 810 : ALLOCATE (mat_Sigma_neg_tau(i_t, ispin)%matrix)
3058 810 : ALLOCATE (mat_Sigma_pos_tau(i_t, ispin)%matrix)
3059 810 : CALL dbcsr_create(mat_Sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
3060 1480 : CALL dbcsr_create(mat_Sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
3061 : END DO
3062 : END DO
3063 :
3064 : ! =========================================================================
3065 : ! 2. IMAGINARY TIME LOOP
3066 : ! Σ^c_neg_λσ(iτ) = -φ^T ( (φ G^occ φ^T) ∘ (Z W^MIC Z^T) ) φ
3067 : ! Σ^c_pos_λσ(iτ) = φ^T ( (φ G^vir φ^T) ∘ (Z W^MIC Z^T) ) φ
3068 : ! =========================================================================
3069 722 : DO i_t = 1, bs_env%num_time_freq_points
3070 670 : tau = bs_env%time_frequency_grid%imaginary_time(i_t)
3071 :
3072 : CALL dbcsr_create(matrix_W_aux, "W_aux", dist_aux_aux, dbcsr_type_no_symmetry, &
3073 670 : blk_aux, blk_aux)
3074 670 : IF (bs_env%ri_rs%cutoff_radius_g_w > 0.0_dp .AND. ALLOCATED(bs_env%ri_rs%atom_centers)) THEN
3075 : CALL reserve_blocks_within_radius(matrix_W_aux, bs_env%ri_rs%atom_centers, &
3076 0 : bs_env%ri_rs%cutoff_radius_g_w)
3077 0 : CALL copy_fm_to_dbcsr(fm_W_time(i_t), matrix_W_aux, keep_sparsity=.TRUE.)
3078 : ELSE
3079 670 : CALL copy_fm_to_dbcsr(fm_W_time(i_t), matrix_W_aux, keep_sparsity=.FALSE.)
3080 : END IF
3081 670 : CALL dbcsr_filter(matrix_W_aux, bs_env%eps_filter)
3082 :
3083 1480 : DO ispin = 1, bs_env%n_spin
3084 810 : t1 = m_walltime()
3085 :
3086 : ! AO-space Green's functions G^occ_µν, G^vir_µν (dense AO x AO, small)
3087 810 : CALL build_G_ao(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_G_occ_ao)
3088 810 : CALL build_G_ao(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, matrix_G_vir_ao)
3089 :
3090 : ! Σ^c_neg and Σ^c_pos in a single panel loop: W_pan = Z_panel × W × Z^T built once
3091 : CALL contract_grid_panels_sigma_c(mat_phi=mat_phi_mu_l, mat_Z=mat_Z_lP, &
3092 : mat_G_occ_ao=matrix_G_occ_ao, &
3093 : mat_G_vir_ao=matrix_G_vir_ao, &
3094 : mat_W_aux=matrix_W_aux, &
3095 : mat_Sigma_neg=mat_Sigma_neg_tau(i_t, ispin)%matrix, &
3096 : mat_Sigma_pos=mat_Sigma_pos_tau(i_t, ispin)%matrix, &
3097 : eps=bs_env%eps_filter, &
3098 : para_env=bs_env%para_env, &
3099 : pan_first=pan_first, pan_last=pan_last, &
3100 : keep_sparsity=bs_env%ri_rs%keep_sparsity_rirs, &
3101 : centroids=bs_env%ri_rs%chunk_centroids, &
3102 810 : cutoff=bs_env%ri_rs%cutoff_radius_v_w)
3103 810 : CALL dbcsr_scale(mat_Sigma_neg_tau(i_t, ispin)%matrix, -1.0_dp)
3104 :
3105 810 : CALL dbcsr_release(matrix_G_occ_ao)
3106 810 : CALL dbcsr_release(matrix_G_vir_ao)
3107 :
3108 1480 : IF (bs_env%unit_nr > 0) THEN
3109 : WRITE (bs_env%unit_nr, '(T2,A,I15,A,I3,A,F7.1,A)') &
3110 405 : 'Computed Σ^c(iτ) for time point', i_t, ' /', bs_env%num_time_freq_points, &
3111 810 : ', Execution time', m_walltime() - t1, ' s'
3112 : END IF
3113 :
3114 : END DO ! ispin
3115 :
3116 722 : CALL dbcsr_release(matrix_W_aux)
3117 :
3118 : END DO ! i_t
3119 :
3120 : ! -------------------------------------------------------------------------
3121 : ! 3. FINALIZE AND CLEANUP
3122 : ! -------------------------------------------------------------------------
3123 : CALL fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
3124 52 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
3125 :
3126 : ! fm_W_time and the scratch files are released by the caller: in an evGW0 cycle this
3127 : ! routine is entered once per iteration and both have to survive until it is done.
3128 52 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_neg_tau)
3129 52 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_pos_tau)
3130 :
3131 52 : CALL release_square_topology(dist=dist_aux_aux, mapped_dist=dist_row_aux)
3132 :
3133 52 : CALL timestop(handle)
3134 :
3135 104 : END SUBROUTINE compute_Sigma_c
3136 :
3137 : ! **************************************************************************************************
3138 : !> \brief Builds the DBCSR distribution.
3139 : !> \param matrix_template ...
3140 : !> \param square_dist ...
3141 : !> \param blk_sizes ...
3142 : !> \param mapped_dist ...
3143 : ! **************************************************************************************************
3144 2974 : SUBROUTINE setup_square_topology(matrix_template, square_dist, blk_sizes, mapped_dist)
3145 :
3146 : TYPE(dbcsr_type), INTENT(IN) :: matrix_template
3147 : TYPE(dbcsr_distribution_type), INTENT(OUT) :: square_dist
3148 : INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: blk_sizes, mapped_dist
3149 :
3150 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_square_topology'
3151 :
3152 : INTEGER :: handle, i, nprows
3153 2974 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist
3154 : TYPE(dbcsr_distribution_type) :: dist_template
3155 :
3156 2974 : CALL timeset(routineN, handle)
3157 :
3158 2974 : CALL dbcsr_get_info(matrix_template, distribution=dist_template, col_blk_size=col_blk)
3159 2974 : CALL dbcsr_distribution_get(dist_template, col_dist=col_dist, nprows=nprows)
3160 :
3161 2974 : blk_sizes => col_blk
3162 8922 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
3163 9990 : DO i = 1, SIZE(blk_sizes)
3164 9990 : mapped_dist(i) = MOD(i - 1, nprows)
3165 : END DO
3166 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
3167 2974 : row_dist=mapped_dist, col_dist=col_dist)
3168 :
3169 2974 : CALL timestop(handle)
3170 :
3171 5948 : END SUBROUTINE setup_square_topology
3172 :
3173 : ! **************************************************************************************************
3174 : !> \brief Releases a distribution created by setup_square_topology.
3175 : !> \param dist ...
3176 : !> \param mapped_dist ...
3177 : ! **************************************************************************************************
3178 2974 : SUBROUTINE release_square_topology(dist, mapped_dist)
3179 :
3180 : TYPE(dbcsr_distribution_type), INTENT(INOUT) :: dist
3181 : INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: mapped_dist
3182 :
3183 2974 : CALL dbcsr_distribution_release(dist)
3184 2974 : IF (ASSOCIATED(mapped_dist)) THEN
3185 2974 : DEALLOCATE (mapped_dist)
3186 : NULLIFY (mapped_dist)
3187 : END IF
3188 :
3189 2974 : END SUBROUTINE release_square_topology
3190 :
3191 : ! **************************************************************************************************
3192 : !> \brief Σ^c_λσ(iτ) -> Σ^c_nn(ϵ) and the quasi-particle levels of the non-periodic RI-RS path,
3193 : !> ϵ_n^GW = ϵ_n^DFT + Σ^c_nn(ϵ_n^GW) + Σ^x_nn - v^xc_nn.
3194 : !> \param bs_env ...
3195 : !> \param fm_Sigma_x_Gamma ...
3196 : !> \param fm_Sigma_c_Gamma_time ...
3197 : ! **************************************************************************************************
3198 52 : SUBROUTINE compute_QP_energies(bs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
3199 :
3200 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3201 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
3202 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
3203 :
3204 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
3205 :
3206 : INTEGER :: handle, ispin, j_t
3207 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Sigma_x_n, V_xc_n
3208 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Sigma_c_n_freq, Sigma_c_n_time
3209 : TYPE(cp_fm_type) :: fm_ks, fm_mos, fm_s, fm_work
3210 :
3211 52 : CALL timeset(routineN, handle)
3212 :
3213 52 : CALL cp_fm_create(fm_ks, bs_env%fm_s_Gamma%matrix_struct)
3214 52 : CALL cp_fm_create(fm_s, bs_env%fm_s_Gamma%matrix_struct)
3215 52 : CALL cp_fm_create(fm_mos, bs_env%fm_s_Gamma%matrix_struct)
3216 52 : CALL cp_fm_create(fm_work, bs_env%fm_s_Gamma%matrix_struct)
3217 :
3218 208 : ALLOCATE (V_xc_n(bs_env%n_ao), Sigma_x_n(bs_env%n_ao))
3219 260 : ALLOCATE (Sigma_c_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
3220 156 : ALLOCATE (Sigma_c_n_freq(bs_env%n_ao, bs_env%num_time_freq_points, 2))
3221 :
3222 116 : DO ispin = 1, bs_env%n_spin
3223 :
3224 : ! 1. Roothaan-Hall H^KS_µν C_νn = S_µν C_νn ϵ_n
3225 64 : CALL cp_fm_to_fm(bs_env%fm_ks_Gamma(ispin), fm_ks)
3226 64 : CALL cp_fm_to_fm(bs_env%fm_s_Gamma, fm_s)
3227 64 : CALL cp_fm_geeig(fm_ks, fm_s, fm_mos, bs_env%eigenval_scf(:, 1, ispin), fm_work)
3228 :
3229 : ! 2. v^xc_µν -> v^xc_nn and Σ^x_µν -> Σ^x_nn
3230 64 : CALL to_Gamma_and_mo_real(V_xc_n, bs_env%fm_V_xc_Gamma(ispin), fm_mos)
3231 64 : CALL to_Gamma_and_mo_real(Sigma_x_n, fm_Sigma_x_Gamma(ispin), fm_mos)
3232 :
3233 : ! 3. Σ^c_µν(+/-i|τ_j|) -> Σ^c_nn(+/-i|τ_j|)
3234 874 : DO j_t = 1, bs_env%num_time_freq_points
3235 : CALL to_Gamma_and_mo_real(Sigma_c_n_time(:, j_t, 1), &
3236 810 : fm_Sigma_c_Gamma_time(j_t, 1, ispin), fm_mos)
3237 : CALL to_Gamma_and_mo_real(Sigma_c_n_time(:, j_t, 2), &
3238 874 : fm_Sigma_c_Gamma_time(j_t, 2, ispin), fm_mos)
3239 : END DO
3240 :
3241 : ! 4. Σ^c_nn(iτ) -> Σ^c_nn(iω)
3242 64 : CALL time_to_freq(bs_env, Sigma_c_n_time, Sigma_c_n_freq, ispin)
3243 :
3244 : ! 5. Analytic continuation Σ^c_nn(iω) -> Σ^c_nn(ϵ) and the QP levels
3245 : CALL analyt_conti_and_print(bs_env, Sigma_c_n_freq, Sigma_x_n, V_xc_n, &
3246 116 : bs_env%eigenval_scf(:, 1, ispin), 1, ispin)
3247 :
3248 : END DO ! ispin
3249 :
3250 52 : CALL get_all_VBM_CBM_bandgaps(bs_env)
3251 :
3252 52 : IF (bs_env%gw_flavour == G0W0) CALL cp_fm_release(fm_Sigma_x_Gamma)
3253 52 : CALL cp_fm_release(fm_Sigma_c_Gamma_time)
3254 :
3255 52 : CALL cp_fm_release(fm_ks)
3256 52 : CALL cp_fm_release(fm_s)
3257 52 : CALL cp_fm_release(fm_mos)
3258 52 : CALL cp_fm_release(fm_work)
3259 :
3260 52 : CALL timestop(handle)
3261 :
3262 104 : END SUBROUTINE compute_QP_energies
3263 :
3264 : ! **************************************************************************************************
3265 : !> \brief AO -> MO transform of a Γ-point matrix
3266 : !> \param array_n ...
3267 : !> \param fm_Gamma ...
3268 : !> \param fm_mos ...
3269 : ! **************************************************************************************************
3270 1748 : SUBROUTINE to_Gamma_and_mo_real(array_n, fm_Gamma, fm_mos)
3271 :
3272 : REAL(KIND=dp), DIMENSION(:) :: array_n
3273 : TYPE(cp_fm_type) :: fm_Gamma, fm_mos
3274 :
3275 : CHARACTER(LEN=*), PARAMETER :: routineN = 'to_Gamma_and_mo_real'
3276 :
3277 : INTEGER :: handle
3278 : TYPE(cp_fm_type) :: fm_mo
3279 :
3280 1748 : CALL timeset(routineN, handle)
3281 :
3282 1748 : CALL cp_fm_create(fm_mo, fm_Gamma%matrix_struct)
3283 :
3284 : ! A_nn' = Σ_μν C_μn A_μν C_νn'
3285 1748 : CALL fm_contract_ABA(fm_mos, fm_Gamma, fm_mo)
3286 :
3287 1748 : CALL cp_fm_get_diag(fm_mo, array_n)
3288 :
3289 1748 : CALL cp_fm_release(fm_mo)
3290 :
3291 1748 : CALL timestop(handle)
3292 :
3293 1748 : END SUBROUTINE to_Gamma_and_mo_real
3294 :
3295 : ! **************************************************************************************************
3296 : !> \brief Evaluates the atom-local AO functions Φ_μ(r_l) on the RI-RS grid points.
3297 : !> \param phi_val Values Φ_μ(r_l) for the requested atom
3298 : !> \param ri_rs_grid Cartesian RI-RS grid points r_l
3299 : !> \param npts Number of RI-RS grid points
3300 : !> \param iatom Atom whose AO functions are evaluated
3301 : !> \param particle_set Atomic positions and kind assignments
3302 : !> \param qs_kind_set Atomic-kind basis information
3303 : !> \param cell Simulation cell and boundary conditions
3304 : !> \param r2_threshold Squared radial cutoff for neglecting AO values
3305 : ! **************************************************************************************************
3306 56 : SUBROUTINE fill_phi_for_atom(phi_val, ri_rs_grid, npts, iatom, &
3307 : particle_set, qs_kind_set, cell, r2_threshold)
3308 :
3309 : REAL(KIND=dp), INTENT(INOUT) :: phi_val(:, :)
3310 : INTEGER, INTENT(IN) :: npts
3311 : REAL(KIND=dp), INTENT(IN) :: ri_rs_grid(3, npts)
3312 : INTEGER, INTENT(IN) :: iatom
3313 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3314 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3315 : TYPE(cell_type), POINTER :: cell
3316 : REAL(KIND=dp), INTENT(IN) :: r2_threshold
3317 :
3318 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_phi_for_atom'
3319 :
3320 : INTEGER :: handle
3321 :
3322 56 : CALL timeset(routineN, handle)
3323 :
3324 : CALL evaluate_ao_on_points(phi_val, ri_rs_grid, iatom, particle_set, qs_kind_set, cell, &
3325 56 : cutoff_squared=r2_threshold)
3326 :
3327 56 : CALL timestop(handle)
3328 :
3329 56 : END SUBROUTINE fill_phi_for_atom
3330 :
3331 : ! **************************************************************************************************
3332 : !> \brief Counts the AO functions whose radial support intersects an RI fitting sphere.
3333 : !> \param bs_env GW calculation environment and atom-to-AO offsets
3334 : !> \param particle_set Atomic positions
3335 : !> \param atom_P Atom at the center of the RI fitting sphere
3336 : !> \param cutoff_ri Radius of the RI fitting sphere
3337 : !> \param n_ao_used Number of intersecting AO functions
3338 : ! **************************************************************************************************
3339 96 : SUBROUTINE get_n_ao_in_sphere(bs_env, particle_set, atom_P, cutoff_ri, n_ao_used)
3340 :
3341 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3342 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3343 : INTEGER, INTENT(IN) :: atom_P
3344 : REAL(KIND=dp), INTENT(IN) :: cutoff_ri
3345 : INTEGER, INTENT(OUT) :: n_ao_used
3346 :
3347 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_n_ao_in_sphere'
3348 :
3349 : INTEGER :: handle, ri_atom
3350 :
3351 96 : CALL timeset(routineN, handle)
3352 :
3353 96 : n_ao_used = 0
3354 348 : DO ri_atom = 1, bs_env%n_atom
3355 1008 : IF (NORM2(particle_set(ri_atom)%r(:) - particle_set(atom_P)%r(:)) > &
3356 : bs_env%ri_rs%radius_ao_per_atom(ri_atom) + cutoff_ri) CYCLE
3357 : n_ao_used = n_ao_used + bs_env%i_ao_end_from_atom(ri_atom) - &
3358 348 : bs_env%i_ao_start_from_atom(ri_atom) + 1
3359 : END DO
3360 :
3361 96 : CALL timestop(handle)
3362 :
3363 96 : END SUBROUTINE get_n_ao_in_sphere
3364 :
3365 : ! **************************************************************************************************
3366 : !> \brief Prints the percentage of non-zero elements in a distributed RI-RS matrix.
3367 : !> \param matrix Distributed matrix whose occupation is reported
3368 : !> \param label Mathematical matrix label used in the output
3369 : !> \param bs_env GW calculation environment and output unit
3370 : !> \param suffix Optional text appended to the matrix label
3371 : ! **************************************************************************************************
3372 122 : SUBROUTINE print_matrix_occupation(matrix, label, bs_env, suffix)
3373 :
3374 : TYPE(dbcsr_type), INTENT(IN) :: matrix
3375 : CHARACTER(LEN=*), INTENT(IN) :: label
3376 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3377 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: suffix
3378 :
3379 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_matrix_occupation'
3380 :
3381 : CHARACTER(LEN=32) :: output_format
3382 : CHARACTER(LEN=max_line_length) :: msg, output_label
3383 : INTEGER :: handle, i, unicode_shift
3384 : REAL(KIND=dp) :: frac_2p31, max_loc, occ
3385 :
3386 122 : CALL timeset(routineN, handle)
3387 :
3388 122 : occ = dbcsr_get_occupation(matrix)
3389 122 : max_loc = REAL(dbcsr_get_data_size(matrix), dp)
3390 122 : CALL bs_env%para_env%max(max_loc)
3391 :
3392 122 : IF (bs_env%unit_nr > 0) THEN
3393 61 : frac_2p31 = max_loc/REAL(bs_env%dbcsr_msg_elem_limit, dp)
3394 61 : output_label = 'Percentage of non-zero matrix elements in '//TRIM(label)
3395 61 : IF (PRESENT(suffix)) output_label = TRIM(output_label)//TRIM(suffix)
3396 : ! Fortran counts UTF-8 bytes, whereas the terminal displays each Greek letter in one
3397 : ! column. Shift the absolute output tab once for every continuation byte.
3398 61 : unicode_shift = 0
3399 3347 : DO i = 1, LEN_TRIM(output_label)
3400 3347 : IF (IAND(IACHAR(output_label(i:i)), 192) == 128) unicode_shift = unicode_shift + 1
3401 : END DO
3402 61 : WRITE (output_format, '(A,I0,A)') '(T2,A,T', 72 + unicode_shift, ',F7.2,A)'
3403 61 : WRITE (bs_env%unit_nr, output_format) TRIM(output_label), occ*100.0_dp, ' %'
3404 61 : IF (frac_2p31 > 0.5_dp) THEN
3405 : WRITE (msg, '(3A,F0.2,A)') &
3406 0 : "The largest per-rank message of ", TRIM(label), " reaches ", frac_2p31, &
3407 : " of the 32-bit limit that DBCSR uses for its message length. Beyond it the "// &
3408 : "length overflows and multiply_cannon fails. Reduce the per-rank block size, "// &
3409 0 : "for instance with more MPI ranks or a larger N_PANELS."
3410 0 : CPWARN(TRIM(msg))
3411 : END IF
3412 61 : CALL m_flush(bs_env%unit_nr)
3413 : END IF
3414 :
3415 122 : CALL timestop(handle)
3416 :
3417 122 : END SUBROUTINE print_matrix_occupation
3418 :
3419 : END MODULE gw_non_periodic_ri_rs
|
