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