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