Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief GW using RI-RS Approximation for molecules
10 : !> \par History
11 : !> 04.2026 created [Ritaj Tyagi]
12 : ! **************************************************************************************************
13 :
14 : MODULE gw_non_periodic_ri_rs
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind_set
17 : USE cell_types, ONLY: cell_type
18 : USE constants_operator, ONLY: operator_coulomb
19 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
20 : cp_blacs_env_release,&
21 : cp_blacs_env_type
22 : USE cp_dbcsr_api, ONLY: &
23 : dbcsr_add, dbcsr_binary_read, dbcsr_binary_write, dbcsr_copy, dbcsr_create, &
24 : dbcsr_deallocate_matrix, dbcsr_distribution_get, dbcsr_distribution_new, &
25 : dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
26 : dbcsr_get_block_p, dbcsr_get_data_size, dbcsr_get_info, dbcsr_get_occupation, &
27 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
28 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_put_block, &
29 : dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
30 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 : copy_fm_to_dbcsr,&
33 : dbcsr_deallocate_matrix_set,&
34 : max_elements_per_block
35 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
36 : cp_fm_uplo_to_full
37 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
38 : cp_fm_cholesky_invert,&
39 : cp_fm_cholesky_solve
40 : USE cp_fm_diag, ONLY: cp_fm_geeig,&
41 : cp_fm_power
42 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
43 : cp_fm_struct_release,&
44 : cp_fm_struct_type
45 : USE cp_fm_types, ONLY: &
46 : cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, &
47 : cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
48 : USE cp_log_handling, ONLY: cp_get_default_logger,&
49 : cp_logger_type
50 : USE cp_output_handling, ONLY: cp_p_file,&
51 : cp_print_key_should_output
52 : USE gw_integrals, ONLY: build_3c_integral_block_ctx,&
53 : gw_3c_ctx_create,&
54 : gw_3c_ctx_release,&
55 : gw_3c_ctx_type,&
56 : gw_3c_ws_create,&
57 : gw_3c_ws_release,&
58 : gw_3c_ws_type
59 : USE gw_large_cell_gamma, ONLY: Fourier_transform_w_to_t,&
60 : G_occ_vir,&
61 : compute_fm_chi_Gamma_freq,&
62 : create_fm_W_MIC_time,&
63 : delete_unnecessary_files,&
64 : fill_fm_Sigma_c_Gamma_time,&
65 : fm_write,&
66 : multiply_fm_W_MIC_time_with_Minv_Gamma
67 : USE gw_optimize_ri_rs_grid, ONLY: optimize_ri_rs_grid
68 : USE gw_utils, ONLY: analyt_conti_and_print,&
69 : de_init_bs_env,&
70 : time_to_freq
71 : USE input_constants, ONLY: G0W0,&
72 : evGW0,&
73 : rtp_method_bse
74 : USE input_section_types, ONLY: section_vals_type
75 : USE kinds, ONLY: default_path_length,&
76 : default_string_length,&
77 : dp,&
78 : int_8,&
79 : max_line_length
80 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
81 : USE machine, ONLY: m_flush,&
82 : m_hostnm,&
83 : m_memory_details,&
84 : m_walltime
85 : USE message_passing, ONLY: mp_para_env_type
86 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
87 : USE parallel_gemm_api, ONLY: parallel_gemm
88 : USE particle_types, ONLY: particle_type
89 : USE physcon, ONLY: angstrom,&
90 : evolt
91 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type,&
92 : rirs_grid_type
93 : USE post_scf_bandstructure_utils, ONLY: get_all_VBM_CBM_bandgaps
94 : USE qs_environment_types, ONLY: get_qs_env,&
95 : qs_environment_type
96 : USE qs_kind_types, ONLY: get_qs_kind,&
97 : qs_kind_type
98 : USE rirs_grid_utils, ONLY: evaluate_ao_on_points,&
99 : get_rirs_grid_filepath,&
100 : initialize_rirs_grid,&
101 : read_rirs_grid_file
102 : USE util, ONLY: sort
103 : #include "./base/base_uses.f90"
104 :
105 : IMPLICIT NONE
106 :
107 : PRIVATE
108 :
109 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_non_periodic_ri_rs'
110 :
111 : PUBLIC :: gw_calc_non_periodic_ri_rs, ri_rs_grid_assembler, &
112 : get_basis_offsets, precompute_ri_rs_radii, solve_D_lp_distributed, &
113 : atomic_basis_at_grid_point, compute_coeff_Z_lP, &
114 : reserve_blocks_within_radius
115 :
116 : CONTAINS
117 :
118 : ! **************************************************************************************************
119 : !> \brief GW calculation using RI-RS formalism for molecules
120 : !> \param qs_env ...
121 : !> \param bs_env ...
122 : ! **************************************************************************************************
123 36 : SUBROUTINE gw_calc_non_periodic_ri_rs(qs_env, bs_env)
124 :
125 : TYPE(qs_environment_type), POINTER :: qs_env
126 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
127 :
128 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_non_periodic_ri_rs'
129 :
130 : INTEGER :: handle
131 36 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_W_time
132 :
133 36 : CALL timeset(routineN, handle)
134 :
135 : ! ==============================================================================
136 : ! 0. Precompute AO and RI radii
137 : ! Per-atom cutoff radii from the most diffuse Gaussian primitives of
138 : ! the AO ("ORB") and RI auxiliary ("RI_AUX") basis sets:
139 : ! α_min,ao = min { ζ_ao | ζ_ao > 10⁻³ }, α_min,ri analogous
140 : ! r_ao = sqrt( -ln(ε) / α_min,ao ) (radius_ao_per_atom)
141 : ! r_ri = sqrt( -ln(ε) / α_min,ri ) (radius_ri_per_atom)
142 : ! ==============================================================================
143 36 : CALL precompute_ri_rs_radii(qs_env, bs_env)
144 :
145 : ! ==============================================================================
146 : ! 1. Grid generation for RI-RS
147 : ! Modified Lebedev atomic grids (Duchemin & Blase), one per atom,
148 : ! concatenated into a flat global list: r_l = R_A + r_l^(A)
149 : ! ==============================================================================
150 36 : CALL ri_rs_grid_assembler(qs_env, bs_env, bs_env%ri_rs%grid_points)
151 :
152 : ! ==============================================================================
153 : ! 2a. Atomic basis evaluation on the grid (grid x AO matrix)
154 : ! Φ_μl = Φ_μ(r_l) (mat_phi_mu_l)
155 : ! ==============================================================================
156 : CALL atomic_basis_at_grid_point(qs_env, bs_env, bs_env%ri_rs%grid_points, &
157 36 : bs_env%ri_rs%mat_phi_mu_l)
158 :
159 : ! ==============================================================================
160 : ! 2b. Print the memory estimate for the RI-RS calculation
161 : ! ==============================================================================
162 36 : CALL print_ri_rs_memory_estimate(qs_env, bs_env)
163 :
164 : ! ==============================================================================
165 : ! 3. RI-RS fitting coefficients Z_lP (grid x RI matrix)
166 : ! Per-atom regularized solve, restricted to grid points r_l within a
167 : ! cutoff distance of atom P:
168 : ! a. D_ll' = [ Σ_μ Φ_μ(r_l) Φ_μ(r_l') ]²
169 : ! b. D_lP = Σ_μν Φ_μ(r_l) Φ_ν(r_l) (μν|P)
170 : ! c. Jacobi conditioning with d_l = 1/sqrt(D_ll):
171 : ! D'_ll' = d_l D_ll' d_l' + λδ_ll' , D'_lP = d_l D_lP
172 : ! d. Solve Σ_l' D'_ll' Z'_l'P = D'_lP
173 : ! e. Rescale Z_lP = d_l Z'_lP (mat_Z_lP)
174 : ! ==============================================================================
175 : CALL compute_coeff_Z_lP(qs_env, bs_env, bs_env%ri_rs%grid_points, &
176 36 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
177 :
178 : ! flag the RI-RS grid as built so a subsequent RT-BSE run reuses Z_lP
179 : ! instead of rebuilding it
180 36 : bs_env%ri_rs%grid_built = .TRUE.
181 :
182 : ! ==============================================================================
183 : ! 4. Polarizability matrix χ on the imaginary-time grid
184 : ! G^occ_µλ(i|τ|) = Σ_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
185 : ! G^vir_µλ(i|τ|) = Σ_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
186 : ! G^occ_ll'(i|τ|) = Σ_µν Φ_µ(r_l) G^occ_µν Φ_ν(r_l') (G^vir analogous)
187 : ! χ_ll'(iτ) = G^occ_ll'(i|τ|) ∘ G^vir_ll'(i|τ|) (element-wise)
188 : ! χ_PQ(iτ) = Σ_ll' Z_lP χ_ll'(iτ) Z_l'Q
189 : ! ==============================================================================
190 : CALL get_mat_chi_Gamma_tau(bs_env, bs_env%mat_chi_Gamma_tau, &
191 36 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
192 :
193 : ! ==============================================================================
194 : ! 5. Screened Coulomb interaction W (RI basis)
195 : ! χ_PQ(iτ) -> χ_PQ(iω) -> ε_PQ(iω) -> W_PQ(iω) -> W_PQ(iτ)
196 : ! ==============================================================================
197 36 : CALL compute_W(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_time)
198 :
199 : ! ==============================================================================
200 : ! 6. Exact-exchange self-energy Σ^x
201 : ! D_µν = Σ_n^occ C_µn C_νn (density matrix)
202 : ! D_ll' = Σ_µν Φ_µ(r_l) D_µν Φ_ν(r_l')
203 : ! V^tr_ll' = Σ_PQ Z_lP V^tr_PQ Z_l'Q (truncated Coulomb)
204 : ! Σ^x_ll' = D_ll' ∘ V^tr_ll'
205 : ! Σ^x_λσ(k=0) = -Σ_ll' Φ_λ(r_l) Σ^x_ll' Φ_σ(r_l')
206 : ! ==============================================================================
207 : CALL compute_Sigma_x(bs_env, qs_env, bs_env%ri_rs%mat_phi_mu_l, &
208 36 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_x_Gamma)
209 :
210 : ! ==============================================================================
211 : ! 7. Correlation self-energy Σ^c and quasiparticle energies, iterated
212 : ! until eigenvalue self-consistency if the &EVGW0 section is given and
213 : ! done in a single pass for G0W0.
214 : !
215 : ! (a) W_ll'(iτ) = Σ_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
216 : ! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) ∘ W_ll'(iτ), τ < 0
217 : ! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) ∘ W_ll'(iτ), τ > 0
218 : ! Σ^c_λσ(iτ) = Σ_ll' Φ_λ(r_l) Σ^c_ll'(iτ) Φ_σ(r_l')
219 : ! (b) Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
220 : ! ϵ_n^GW,(i) = ϵ_n^DFT + Σ^c_nn[G^(i-1),W](ϵ_n^GW,(i)) + Σ^x_nn - v^xc_nn
221 : ! ==============================================================================
222 36 : CALL compute_Sigma_c_and_QP_energies(bs_env, fm_W_time, fm_Sigma_x_Gamma)
223 :
224 36 : CALL de_init_bs_env(qs_env, bs_env)
225 :
226 36 : CALL timestop(handle)
227 :
228 36 : END SUBROUTINE gw_calc_non_periodic_ri_rs
229 :
230 : ! **************************************************************************************************
231 : !> \brief Correlation self-energy and quasiparticle energies, iterated to eigenvalue
232 : !> self-consistency in G (evGW0).
233 : !>
234 : !> For G0W0 this runs once with G^(0) built from the DFT eigenvalues. For evGW0 the
235 : !> Green's function is rebuilt from the quasiparticle energies of all states until the
236 : !> quasiparticle HOMO, LUMO and HOMO-LUMO gap change by less than EPS_ITER between two
237 : !> cycles, or MAX_ITER cycles are spent. W, Σ^x and everything computed before this
238 : !> routine stay frozen: they either do not depend on the eigenvalues at all (Σ^x is G
239 : !> at τ = 0, a pure density matrix) or are held fixed by construction in GW0. The
240 : !> imaginary-time grid is fixed as well, since W(iτ) lives on it.
241 : !>
242 : !> \param bs_env ...
243 : !> \param fm_W_time ...
244 : !> \param fm_Sigma_x_Gamma ...
245 : ! **************************************************************************************************
246 36 : SUBROUTINE compute_Sigma_c_and_QP_energies(bs_env, fm_W_time, fm_Sigma_x_Gamma)
247 :
248 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
249 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time, fm_Sigma_x_Gamma
250 :
251 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c_and_QP_energies'
252 :
253 : INTEGER :: handle, i_iter, n_iter
254 : LOGICAL :: converged
255 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_scf_Gamma_dft
256 : REAL(KIND=dp), DIMENSION(2) :: e_fermi_dft
257 : REAL(KIND=dp), DIMENSION(3, 2) :: band_prev
258 36 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
259 :
260 36 : CALL timeset(routineN, handle)
261 :
262 : converged = .FALSE.
263 :
264 36 : CALL init_evGW0(bs_env, n_iter, band_prev, eigenval_scf_Gamma_dft, e_fermi_dft)
265 :
266 : ! evGW0 self-consistency loop; for G0W0, loop is terminated after one iteration
267 50 : DO i_iter = 1, n_iter
268 :
269 50 : bs_env%ri_rs%evgw0_i_iter = i_iter
270 :
271 : ! W_ll'(iτ) = Σ_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
272 : ! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) ∘ W_ll'(iτ), τ < 0
273 : ! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) ∘ W_ll'(iτ), τ > 0
274 : ! Σ^c_λσ(iτ) = Σ_ll' Φ_λ(r_l) Σ^c_ll'(iτ) Φ_σ(r_l')
275 : CALL compute_Sigma_c(bs_env, fm_W_time, bs_env%ri_rs%mat_phi_mu_l, &
276 50 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_c_Gamma_time)
277 :
278 : ! Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
279 : ! ϵ_n^GW,(i) = ϵ_n^DFT + Σ^c_nn[G^(i-1),W](ϵ_n^GW,(i)) + Σ^x_nn - v^xc_nn
280 50 : CALL compute_QP_energies(bs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
281 :
282 50 : IF (bs_env%gw_flavour == G0W0) EXIT
283 :
284 18 : CALL print_evGW0_band_edges(bs_env, band_prev, i_iter, n_iter, converged)
285 :
286 18 : IF (i_iter == n_iter .OR. converged) EXIT
287 :
288 : ! eigenvalues ϵ_n to be updated in G:
289 : ! G^occ_µλ(i|τ|) = Σ_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
290 : ! G^vir_µλ(i|τ|) = Σ_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
291 68 : CALL update_eigenvalues_G(bs_env)
292 :
293 : END DO
294 :
295 36 : CALL cp_fm_release(fm_W_time)
296 :
297 36 : CALL reset_and_clean_bs_env(bs_env, eigenval_scf_Gamma_dft, e_fermi_dft, fm_Sigma_x_Gamma)
298 :
299 36 : CALL delete_unnecessary_files(bs_env)
300 :
301 36 : CALL timestop(handle)
302 :
303 36 : END SUBROUTINE compute_Sigma_c_and_QP_energies
304 :
305 : ! **************************************************************************************************
306 : !> \brief Sets up the evGW0 eigenvalue self-consistency loop: the cycle count, the eigenvalues
307 : !> the first Green's function is built from, and the DFT reference that the loop
308 : !> overwrites. A G0W0 run reduces to a single cycle and needs none of it.
309 : !> \param bs_env ...
310 : !> \param n_iter ...
311 : !> \param band_prev ...
312 : !> \param eigenval_scf_Gamma_dft ...
313 : !> \param e_fermi_dft ...
314 : ! **************************************************************************************************
315 36 : SUBROUTINE init_evGW0(bs_env, n_iter, band_prev, eigenval_scf_Gamma_dft, e_fermi_dft)
316 :
317 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
318 : INTEGER, INTENT(OUT) :: n_iter
319 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(OUT) :: band_prev
320 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
321 : INTENT(OUT) :: eigenval_scf_Gamma_dft
322 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: e_fermi_dft
323 :
324 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_evGW0'
325 :
326 : INTEGER :: handle
327 :
328 36 : CALL timeset(routineN, handle)
329 :
330 36 : n_iter = 1
331 36 : band_prev(:, :) = 0.0_dp
332 36 : e_fermi_dft(:) = 0.0_dp
333 :
334 36 : IF (bs_env%gw_flavour == evGW0) THEN
335 4 : n_iter = bs_env%ri_rs%evgw0_iter
336 :
337 : ! eigenvalues currently in G; the first cycle starts is a plain G0W0 step
338 70 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_scf(:, :, :)
339 : ! the loop overwrites these; post-GW printing expects the DFT values back
340 76 : ALLOCATE (eigenval_scf_Gamma_dft, SOURCE=bs_env%eigenval_scf_Gamma)
341 12 : e_fermi_dft(:) = bs_env%e_fermi(:)
342 : END IF
343 :
344 36 : CALL timestop(handle)
345 :
346 36 : END SUBROUTINE init_evGW0
347 :
348 : ! **************************************************************************************************
349 : !> \brief Feeds the quasiparticle energies of the current evGW0 cycle back into the Green's
350 : !> function used by the next one, and re-centres the Fermi level between the new
351 : !> quasiparticle HOMO and LUMO.
352 : !> \param bs_env ...
353 : ! **************************************************************************************************
354 14 : SUBROUTINE update_eigenvalues_G(bs_env)
355 :
356 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
357 :
358 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_eigenvalues_G'
359 :
360 : INTEGER :: handle, i_mo, ispin
361 :
362 14 : CALL timeset(routineN, handle)
363 :
364 : ! record this cycle's evGW0 result; it is also what the next Green's function is built from
365 230 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_GW(:, :, :)
366 :
367 34 : DO ispin = 1, bs_env%n_spin
368 : ! Update all physical states; exclude linear-dependency placeholders.
369 196 : DO i_mo = 1, bs_env%n_mo_retained
370 196 : bs_env%eigenval_scf_Gamma(i_mo, ispin) = bs_env%eigenval_GW(i_mo, 1, ispin)
371 : END DO
372 : bs_env%e_fermi(ispin) = &
373 : 0.5_dp*(bs_env%eigenval_GW(bs_env%n_occ(ispin), 1, ispin) + &
374 34 : bs_env%eigenval_GW(bs_env%n_occ(ispin) + 1, 1, ispin))
375 : END DO
376 :
377 14 : CALL timestop(handle)
378 :
379 14 : END SUBROUTINE update_eigenvalues_G
380 :
381 : ! **************************************************************************************************
382 : !> \brief Restores the DFT reference that the evGW0 loop overwrote and cleanup
383 : !> \param bs_env ...
384 : !> \param eigenval_scf_Gamma_dft ...
385 : !> \param e_fermi_dft ...
386 : !> \param fm_Sigma_x_Gamma ..
387 : ! **************************************************************************************************
388 36 : SUBROUTINE reset_and_clean_bs_env(bs_env, eigenval_scf_Gamma_dft, e_fermi_dft, fm_Sigma_x_Gamma)
389 :
390 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
391 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
392 : INTENT(INOUT) :: eigenval_scf_Gamma_dft
393 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: e_fermi_dft
394 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
395 :
396 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_and_clean_bs_env'
397 :
398 : INTEGER :: handle
399 :
400 36 : CALL timeset(routineN, handle)
401 :
402 36 : IF (bs_env%gw_flavour == evGW0) THEN
403 : ! Update array with final evGW0 result
404 70 : bs_env%eigenval_evGW0(:, :, :) = bs_env%eigenval_GW(:, :, :)
405 : ! put the DFT reference back for the post-GW DOS/band-edge printing
406 64 : bs_env%eigenval_scf_Gamma(:, :) = eigenval_scf_Gamma_dft(:, :)
407 12 : bs_env%e_fermi(:) = e_fermi_dft(:)
408 4 : DEALLOCATE (eigenval_scf_Gamma_dft)
409 4 : CALL cp_fm_release(fm_Sigma_x_Gamma)
410 : END IF
411 :
412 36 : CALL timestop(handle)
413 :
414 36 : END SUBROUTINE reset_and_clean_bs_env
415 :
416 : ! **************************************************************************************************
417 : !> \brief Print the quasiparticle HOMO, LUMO and HOMO-LUMO gap of the current evGW0 cycle
418 : !> \param bs_env ...
419 : !> \param band_prev ...
420 : !> \param i_iter ...
421 : !> \param n_iter ...
422 : !> \param converged ...
423 : ! **************************************************************************************************
424 18 : SUBROUTINE print_evGW0_band_edges(bs_env, band_prev, i_iter, n_iter, converged)
425 :
426 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
427 : REAL(KIND=dp), DIMENSION(3, 2), INTENT(INOUT) :: band_prev
428 : INTEGER, INTENT(IN) :: i_iter, n_iter
429 : LOGICAL, INTENT(OUT) :: converged
430 :
431 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_evGW0_band_edges'
432 :
433 : INTEGER :: handle, homo, ispin, u
434 : REAL(KIND=dp) :: max_delta
435 : REAL(KIND=dp), DIMENSION(3) :: band
436 :
437 18 : CALL timeset(routineN, handle)
438 :
439 18 : u = bs_env%unit_nr
440 :
441 18 : converged = (i_iter > 1)
442 18 : max_delta = 0.0_dp
443 :
444 18 : IF (u > 0) THEN
445 9 : WRITE (u, '(A)') ' '
446 9 : WRITE (u, '(T2,A)') REPEAT('-', 79)
447 9 : WRITE (u, '(T2,A,I4,A,I4)') 'evGW0 cycle', i_iter, ' /', n_iter
448 9 : WRITE (u, '(T2,A)') REPEAT('-', 79)
449 : END IF
450 :
451 44 : DO ispin = 1, bs_env%n_spin
452 :
453 26 : homo = bs_env%n_occ(ispin)
454 26 : band(1) = bs_env%eigenval_GW(homo, 1, ispin)
455 26 : band(2) = bs_env%eigenval_GW(homo + 1, 1, ispin)
456 26 : band(3) = band(2) - band(1)
457 :
458 26 : IF (i_iter > 1) THEN
459 80 : max_delta = MAX(max_delta, MAXVAL(ABS(band(:) - band_prev(:, ispin))))
460 : END IF
461 :
462 26 : IF (u > 0) THEN
463 13 : IF (bs_env%n_spin == 2) WRITE (u, '(T2,A,I0)') 'Spin ', ispin
464 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 HOMO (eV)', band(1)*evolt
465 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 LUMO (eV)', band(2)*evolt
466 13 : WRITE (u, '(T2,A,T61,F20.3)') 'evGW0 HOMO-LUMO gap (eV)', band(3)*evolt
467 : END IF
468 :
469 122 : band_prev(:, ispin) = band(:)
470 :
471 : END DO
472 :
473 18 : IF (i_iter > 1) THEN
474 14 : converged = (max_delta < bs_env%ri_rs%evgw0_eps_iter)
475 21 : IF (u > 0) WRITE (u, '(T2,A,T61,F20.6)') 'Max. change to previous cycle (eV)', &
476 14 : max_delta*evolt
477 : END IF
478 :
479 11 : IF (u > 0) WRITE (u, '(T2,A)') REPEAT('-', 79)
480 :
481 18 : IF (converged) THEN
482 4 : IF (u > 0) THEN
483 2 : WRITE (u, '(A)') ' '
484 : WRITE (u, '(T2,A,I4,A)') &
485 2 : 'evGW0 eigenvalue self-consistency reached in', i_iter, ' cycles.'
486 2 : WRITE (u, '(A)') ' '
487 : END IF
488 14 : ELSE IF (i_iter == n_iter) THEN
489 : CALL cp_warn(__LOCATION__, &
490 : "The evGW0 eigenvalue self-consistency cycle did not converge "// &
491 : "within MAX_ITER cycles. The reported quasiparticle energies are "// &
492 0 : "those of the last cycle.")
493 : END IF
494 :
495 18 : CALL timestop(handle)
496 :
497 18 : END SUBROUTINE print_evGW0_band_edges
498 :
499 : ! **************************************************************************************************
500 : !> \brief Compute per-atom AO and RI basis radii from the most diffuse Gaussian
501 : !> primitive in the AO ("ORB") and RI auxiliary ("RI_AUX") basis sets.
502 : !> Stores results in bs_env%ri_rs%radius_ao_per_atom(:) and
503 : !> bs_env%ri_rs%radius_ri_per_atom(:) and prints a per-atom table.
504 : !> Radius: r_kind = sqrt(-log(eps) / alpha_min_kind)
505 : !> with eps = eps_filter.
506 : !> \param qs_env ...
507 : !> \param bs_env ...
508 : ! **************************************************************************************************
509 44 : SUBROUTINE precompute_ri_rs_radii(qs_env, bs_env)
510 :
511 : TYPE(qs_environment_type), POINTER :: qs_env
512 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
513 :
514 : CHARACTER(LEN=*), PARAMETER :: routineN = 'precompute_ri_rs_radii'
515 :
516 : INTEGER :: handle, i, iatom, ikind, j, natom, nkind
517 44 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
518 : REAL(KIND=dp) :: eps
519 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha_min_ao_kind, alpha_min_ri_kind
520 44 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet_ao, zet_ri
521 44 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
522 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
523 :
524 44 : CALL timeset(routineN, handle)
525 :
526 : CALL get_qs_env(qs_env, nkind=nkind, atomic_kind_set=atomic_kind_set, &
527 44 : particle_set=particle_set)
528 44 : natom = SIZE(particle_set)
529 :
530 44 : eps = bs_env%eps_filter
531 :
532 176 : ALLOCATE (alpha_min_ao_kind(nkind), alpha_min_ri_kind(nkind))
533 106 : alpha_min_ao_kind = HUGE(1.0_dp)
534 106 : alpha_min_ri_kind = HUGE(1.0_dp)
535 :
536 106 : DO ikind = 1, nkind
537 62 : zet_ao => bs_env%basis_set_AO(ikind)%gto_basis_set%zet
538 62 : zet_ri => bs_env%basis_set_RI(ikind)%gto_basis_set%zet
539 :
540 204 : DO i = 1, SIZE(zet_ao, 1)
541 478 : DO j = 1, SIZE(zet_ao, 2)
542 416 : IF (zet_ao(i, j) > 1.0E-3_dp) THEN
543 274 : alpha_min_ao_kind(ikind) = MIN(alpha_min_ao_kind(ikind), zet_ao(i, j))
544 : END IF
545 : END DO
546 : END DO
547 168 : DO i = 1, SIZE(zet_ri, 1)
548 486 : DO j = 1, SIZE(zet_ri, 2)
549 424 : IF (zet_ri(i, j) > 1.0E-3_dp) THEN
550 362 : alpha_min_ri_kind(ikind) = MIN(alpha_min_ri_kind(ikind), zet_ri(i, j))
551 : END IF
552 : END DO
553 : END DO
554 : END DO
555 :
556 44 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
557 :
558 132 : ALLOCATE (bs_env%ri_rs%radius_ao_per_atom(natom))
559 88 : ALLOCATE (bs_env%ri_rs%radius_ri_per_atom(natom))
560 150 : DO iatom = 1, natom
561 106 : ikind = kind_of(iatom)
562 106 : bs_env%ri_rs%radius_ao_per_atom(iatom) = SQRT(-LOG(eps)/alpha_min_ao_kind(ikind))
563 150 : bs_env%ri_rs%radius_ri_per_atom(iatom) = SQRT(-LOG(eps)/alpha_min_ri_kind(ikind))
564 : END DO
565 :
566 44 : IF (bs_env%unit_nr > 0) THEN
567 22 : WRITE (bs_env%unit_nr, '(T2,A)') 'Per-kind RI-RS basis radii (Å):'
568 22 : WRITE (bs_env%unit_nr, '(T4,A6,2X,A4,2A14)') 'Kind', 'Elem', 'r_AO (Å)', 'r_RI (Å)'
569 53 : DO ikind = 1, nkind
570 : WRITE (bs_env%unit_nr, '(T4,I6,2X,A4,2F14.4)') &
571 31 : ikind, &
572 31 : atomic_kind_set(ikind)%element_symbol, &
573 31 : SQRT(-LOG(eps)/alpha_min_ao_kind(ikind))*angstrom, &
574 84 : SQRT(-LOG(eps)/alpha_min_ri_kind(ikind))*angstrom
575 : END DO
576 22 : WRITE (bs_env%unit_nr, '(A)') ' '
577 : END IF
578 :
579 44 : DEALLOCATE (alpha_min_ao_kind, alpha_min_ri_kind, kind_of)
580 :
581 44 : CALL timestop(handle)
582 :
583 44 : END SUBROUTINE precompute_ri_rs_radii
584 :
585 : ! **************************************************************************************************
586 : !> \brief Spreads the low 21 bits of a into every third bit (bits 0,3,6,...,60): the 1-D helper
587 : !> for a 3-D Morton (Z-order) code. Standard 64-bit magic-mask implementation.
588 : !> \param a value in [0, 2^21)
589 : !> \param x a with two zero bits inserted between consecutive input bits
590 : ! **************************************************************************************************
591 318 : SUBROUTINE morton_split3(a, x)
592 : INTEGER(KIND=int_8), INTENT(IN) :: a
593 : INTEGER(KIND=int_8), INTENT(OUT) :: x
594 :
595 318 : x = IAND(a, INT(z'1FFFFF', int_8))
596 318 : x = IAND(IOR(x, ISHFT(x, 32)), INT(z'1F00000000FFFF', int_8))
597 318 : x = IAND(IOR(x, ISHFT(x, 16)), INT(z'1F0000FF0000FF', int_8))
598 318 : x = IAND(IOR(x, ISHFT(x, 8)), INT(z'100F00F00F00F00F', int_8))
599 318 : x = IAND(IOR(x, ISHFT(x, 4)), INT(z'10C30C30C30C30C3', int_8))
600 318 : x = IAND(IOR(x, ISHFT(x, 2)), INT(z'1249249249249249', int_8))
601 318 : END SUBROUTINE morton_split3
602 :
603 : ! **************************************************************************************************
604 : !> \brief Returns a permutation of atom indices in Morton (Z-order) space-filling order of their
605 : !> Cartesian centers, so consecutive atoms are spatial neighbors. The RI-RS grid rows are
606 : !> laid down in this order, so a contiguous grid panel maps to a compact spatial region and
607 : !> the CUTOFF_RADIUS_RL_W neighborhood of every panel shrinks. The grid row index
608 : !> is a summed contraction index, so ANY permutation is result-preserving; this one is
609 : !> chosen purely to improve locality. Coordinates are normalized to the atom bounding box
610 : !> and quantized to 21 bits per axis (sub-picometre for any real cell).
611 : !> \param particle_set ...
612 : !> \param order order(i) = atom index placed at layout position i
613 : ! **************************************************************************************************
614 44 : SUBROUTINE spatial_atom_order(particle_set, order)
615 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
616 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: order
617 :
618 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spatial_atom_order'
619 : INTEGER, PARAMETER :: nbits = 21
620 :
621 : INTEGER :: handle, ia, k, natom
622 : INTEGER(KIND=int_8) :: cmax, ic(3), m1, m2, m3
623 44 : INTEGER(KIND=int_8), ALLOCATABLE :: mcode(:)
624 : REAL(KIND=dp) :: hi(3), lo(3), span(3)
625 :
626 44 : CALL timeset(routineN, handle)
627 :
628 44 : natom = SIZE(particle_set)
629 220 : ALLOCATE (order(natom), mcode(natom))
630 : cmax = ISHFT(1_int_8, nbits) - 1_int_8
631 :
632 176 : lo(:) = HUGE(1.0_dp)
633 176 : hi(:) = -HUGE(1.0_dp)
634 150 : DO ia = 1, natom
635 468 : DO k = 1, 3
636 318 : lo(k) = MIN(lo(k), particle_set(ia)%r(k))
637 424 : hi(k) = MAX(hi(k), particle_set(ia)%r(k))
638 : END DO
639 : END DO
640 176 : span(:) = hi(:) - lo(:)
641 176 : DO k = 1, 3
642 176 : IF (span(k) <= 0.0_dp) span(k) = 1.0_dp
643 : END DO
644 :
645 150 : DO ia = 1, natom
646 424 : DO k = 1, 3
647 318 : ic(k) = INT(((particle_set(ia)%r(k) - lo(k))/span(k))*REAL(cmax, dp), int_8)
648 424 : ic(k) = MIN(cmax, MAX(0_int_8, ic(k)))
649 : END DO
650 106 : CALL morton_split3(ic(1), m1)
651 106 : CALL morton_split3(ic(2), m2)
652 106 : CALL morton_split3(ic(3), m3)
653 150 : mcode(ia) = IOR(IOR(m1, ISHFT(m2, 1)), ISHFT(m3, 2))
654 : END DO
655 :
656 : ! sort(mcode, natom, order): order(i) = original atom index with the i-th smallest code
657 44 : CALL sort(mcode, natom, order)
658 :
659 44 : DEALLOCATE (mcode)
660 :
661 44 : CALL timestop(handle)
662 :
663 44 : END SUBROUTINE spatial_atom_order
664 :
665 : ! **************************************************************************************************
666 : !> \brief Compute grid points for RI-RS
667 : !> Right now based on Ivan and Xavier implementation
668 : !> JCP 150, 174120 (2019), JCTC 17, 2383 (2021)
669 : !> \param qs_env ...
670 : !> \param bs_env ...
671 : !> \param ri_rs_grid_points ...
672 : ! **************************************************************************************************
673 44 : SUBROUTINE ri_rs_grid_assembler(qs_env, bs_env, ri_rs_grid_points)
674 :
675 : TYPE(qs_environment_type), POINTER :: qs_env
676 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
677 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: ri_rs_grid_points(:, :)
678 :
679 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ri_rs_grid_assembler'
680 :
681 : INTEGER :: atom_idx, end_idx, handle, i_layout, &
682 : natom, start_idx, total_grid_npts
683 44 : INTEGER, ALLOCATABLE :: atom_order(:), ri_rs_grid_offsets(:)
684 : REAL(KIND=dp) :: atomic_center(3)
685 44 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
686 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
687 :
688 44 : CALL timeset(routineN, handle)
689 :
690 : !! Get the information about the atoms in the system
691 44 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, particle_set=particle_set)
692 :
693 44 : natom = SIZE(particle_set)
694 :
695 : !! 1. Build the grid cache
696 44 : CALL build_grid_cache(bs_env, atomic_kind_set)
697 :
698 : !! 2. Optimize RIRS grid points
699 44 : IF (bs_env%ri_rs%grid_opt%enabled) THEN
700 6 : CALL optimize_ri_rs_grid(qs_env, bs_env)
701 : END IF
702 :
703 : !! 3. Calculate grid counts and offsets
704 132 : ALLOCATE (ri_rs_grid_offsets(natom + 1))
705 :
706 : ! grid_atom_boundaries(k) = 1-based start of the k-th atom's grid run in LAYOUT order
707 : ! (the order in which the points are laid down below). Used later to build atom-aligned
708 : ! DBCSR grid row-blocks in atomic_basis_at_grid_point.
709 44 : IF (ALLOCATED(bs_env%ri_rs%grid_atom_boundaries)) DEALLOCATE (bs_env%ri_rs%grid_atom_boundaries)
710 88 : ALLOCATE (bs_env%ri_rs%grid_atom_boundaries(natom + 1))
711 :
712 : ! Lay the grid rows down in Morton (spatial) atom order, so a
713 : ! contiguous grid maps to a compact spatial region.
714 44 : CALL spatial_atom_order(particle_set, atom_order)
715 :
716 44 : total_grid_npts = 0
717 150 : DO i_layout = 1, natom
718 106 : atom_idx = atom_order(i_layout)
719 106 : ri_rs_grid_offsets(atom_idx) = total_grid_npts + 1
720 106 : bs_env%ri_rs%grid_atom_boundaries(i_layout) = total_grid_npts + 1
721 150 : total_grid_npts = total_grid_npts + bs_env%ri_rs%grid_cache(atom_idx)%npts
722 : END DO
723 :
724 44 : ri_rs_grid_offsets(natom + 1) = total_grid_npts + 1
725 44 : bs_env%ri_rs%grid_atom_boundaries(natom + 1) = total_grid_npts + 1
726 :
727 44 : IF (bs_env%unit_nr > 0) THEN
728 : WRITE (bs_env%unit_nr, FMT="(T2,A,T69,I12)") &
729 22 : 'Total grid points used for RI-RS:', total_grid_npts
730 22 : WRITE (bs_env%unit_nr, "(A)") ' '
731 : END IF
732 :
733 : !! 4. Allocate the global ri_rs_grid arrays
734 132 : ALLOCATE (ri_rs_grid_points(3, total_grid_npts))
735 :
736 : !! 5. Parallelize the grid generation loop
737 : !$OMP PARALLEL DO DEFAULT(NONE) &
738 : !$OMP SHARED(ri_rs_grid_points, ri_rs_grid_offsets, &
739 : !$OMP particle_set, bs_env, natom) &
740 : !$OMP PRIVATE(atom_idx, atomic_center, start_idx, end_idx) &
741 44 : !$OMP SCHEDULE(DYNAMIC, 1)
742 : DO atom_idx = 1, natom
743 : atomic_center(:) = particle_set(atom_idx)%r(:)
744 :
745 : start_idx = ri_rs_grid_offsets(atom_idx)
746 : end_idx = start_idx + bs_env%ri_rs%grid_cache(atom_idx)%npts - 1
747 :
748 : !! Shift the cached origin grid by the atom's center
749 : ri_rs_grid_points(1, start_idx:end_idx) = bs_env%ri_rs%grid_cache(atom_idx)%raw_points(1, :) + atomic_center(1)
750 : ri_rs_grid_points(2, start_idx:end_idx) = bs_env%ri_rs%grid_cache(atom_idx)%raw_points(2, :) + atomic_center(2)
751 : ri_rs_grid_points(3, start_idx:end_idx) = bs_env%ri_rs%grid_cache(atom_idx)%raw_points(3, :) + atomic_center(3)
752 :
753 : END DO
754 : !$OMP END PARALLEL DO
755 :
756 : !! 6. Cleanup memory
757 44 : IF (ALLOCATED(bs_env%ri_rs%grid_cache)) THEN
758 150 : DO atom_idx = 1, natom
759 150 : IF (ALLOCATED(bs_env%ri_rs%grid_cache(atom_idx)%raw_points)) THEN
760 106 : DEALLOCATE (bs_env%ri_rs%grid_cache(atom_idx)%raw_points)
761 : END IF
762 : END DO
763 150 : DEALLOCATE (bs_env%ri_rs%grid_cache)
764 : END IF
765 :
766 44 : DEALLOCATE (atom_order, ri_rs_grid_offsets)
767 :
768 44 : CALL timestop(handle)
769 :
770 88 : END SUBROUTINE ri_rs_grid_assembler
771 :
772 : ! **************************************************************************************************
773 : !> \brief Reads grids from .ion files and stores them in memory based on grid_select
774 : !> \param bs_env ...
775 : !> \param atomic_kind_set ...
776 : ! **************************************************************************************************
777 44 : SUBROUTINE build_grid_cache(bs_env, atomic_kind_set)
778 :
779 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
780 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
781 :
782 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_grid_cache'
783 :
784 : CHARACTER(LEN=default_path_length) :: full_path, warning_msg
785 : CHARACTER(LEN=default_string_length) :: atom_sym
786 : INTEGER :: handle, iatom, ikind, &
787 : n_voronoi_candidates, natom, nkind, &
788 : requested_npts
789 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
790 44 : TYPE(rirs_grid_type), ALLOCATABLE :: kind_grid(:)
791 :
792 44 : CALL timeset(routineN, handle)
793 :
794 44 : nkind = SIZE(atomic_kind_set)
795 44 : particle_set => bs_env%ri_rs%particle_set
796 44 : CPASSERT(ASSOCIATED(particle_set))
797 44 : natom = SIZE(particle_set)
798 194 : ALLOCATE (kind_grid(nkind))
799 :
800 106 : DO ikind = 1, nkind
801 62 : atom_sym = TRIM(atomic_kind_set(ikind)%element_symbol)
802 : CALL get_rirs_grid_filepath( &
803 62 : atom_sym, bs_env%ri_rs%grid_select, bs_env%ri_rs%grid_file_suffix, full_path)
804 62 : CALL read_rirs_grid_file(TRIM(full_path), kind_grid(ikind)%raw_points)
805 106 : kind_grid(ikind)%npts = SIZE(kind_grid(ikind)%raw_points, 2)
806 : END DO
807 :
808 : ! Make every atom own an independent copy. The optimizer can therefore adapt symmetry-
809 : ! equivalent atoms to their individual local environments, and GW consumes those copies
810 : ! directly without writing intermediate .ion files.
811 44 : IF (ALLOCATED(bs_env%ri_rs%grid_cache)) DEALLOCATE (bs_env%ri_rs%grid_cache)
812 238 : ALLOCATE (bs_env%ri_rs%grid_cache(natom))
813 150 : DO iatom = 1, natom
814 106 : ikind = particle_set(iatom)%atomic_kind%kind_number
815 106 : bs_env%ri_rs%grid_cache(iatom) = kind_grid(ikind)
816 150 : IF (bs_env%ri_rs%grid_opt%enabled) THEN
817 18 : atom_sym = TRIM(atomic_kind_set(ikind)%element_symbol)
818 18 : requested_npts = bs_env%ri_rs%grid_opt%grid_npoints(ikind)
819 : CALL initialize_rirs_grid( &
820 : bs_env%ri_rs%grid_cache(iatom)%raw_points, requested_npts, iatom, particle_set, &
821 18 : n_voronoi_candidates)
822 18 : IF (n_voronoi_candidates < requested_npts) THEN
823 : WRITE (warning_msg, '(A,I0,A,A,A,I0,A,I0,A)') &
824 0 : "Atom ", iatom, " (", TRIM(atom_sym), "): Grid selection requests ", requested_npts, &
825 0 : " points, but only ", n_voronoi_candidates, &
826 0 : " are inside its Voronoi volume. Selecting from the full RI-RS source grid."
827 0 : CALL cp_warn(__LOCATION__, TRIM(warning_msg))
828 : END IF
829 18 : bs_env%ri_rs%grid_cache(iatom)%npts = requested_npts
830 : END IF
831 : END DO
832 106 : DEALLOCATE (kind_grid)
833 :
834 44 : CALL timestop(handle)
835 :
836 44 : END SUBROUTINE build_grid_cache
837 :
838 : ! **************************************************************************************************
839 : !> \brief Evaluates the AO basis on the RI-RS grid and stores it as the sparse DBCSR matrix
840 : !> Φ_μl = Φ_μ(r_l) (rows = grid points in atom-aligned blocks of at most
841 : !> max_elements_per_block points, columns = one block per atom's full AO set).
842 : !> Grid points outside the reach of an atom's most
843 : !> diffuse Gaussian (or the CUTOFF_RADIUS_RL_AO) are skipped, and only blocks
844 : !> with at least one element > eps_filter are stored. This locality is the source of
845 : !> ALL grid-dimension sparsity used downstream. Also caches the atom centers and the
846 : !> per-chunk centroids needed by the optional CUTOFF_RADIUS_G_W / CUTOFF_RADIUS_RL_W
847 : !> operator truncations.
848 : !> \param qs_env ...
849 : !> \param bs_env ...
850 : !> \param ri_rs_grid_points ...
851 : !> \param mat_phi_mu_l ...
852 : ! **************************************************************************************************
853 44 : SUBROUTINE atomic_basis_at_grid_point(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l)
854 :
855 : TYPE(qs_environment_type), POINTER :: qs_env
856 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
857 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
858 : TYPE(dbcsr_type), INTENT(OUT) :: mat_phi_mu_l
859 :
860 : CHARACTER(LEN=*), PARAMETER :: routineN = 'atomic_basis_at_grid_point'
861 :
862 : INTEGER :: bs_eff, c_size, dimen_ORB, handle, i, i_blk, ia, iatom, natom, npcol, nprow, &
863 : num_grid_chunks, r_end, r_start, remaining, run, safe_max, total_grid_npts
864 44 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_row_start, first_sgf
865 44 : INTEGER, DIMENSION(:), POINTER :: c_blk_sizes, col_dist, r_blk_sizes, &
866 44 : row_dist
867 : REAL(KIND=dp) :: r2_threshold
868 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atom_col_buffer
869 44 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
870 : TYPE(cell_type), POINTER :: cell
871 : TYPE(dbcsr_distribution_type) :: dist
872 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist_ks
873 : TYPE(mp_para_env_type), POINTER :: para_env
874 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
875 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
876 :
877 44 : CALL timeset(routineN, handle)
878 :
879 : ! Extract environment variables
880 : CALL get_qs_env(qs_env, cell=cell, atomic_kind_set=atomic_kind_set, &
881 : qs_kind_set=qs_kind_set, particle_set=particle_set, &
882 44 : para_env=para_env)
883 :
884 44 : natom = SIZE(particle_set)
885 44 : total_grid_npts = SIZE(ri_rs_grid_points, 2)
886 :
887 : ! Map the starting indices of spherical gaussian functions (SGF) for each atom
888 132 : ALLOCATE (first_sgf(natom + 1))
889 44 : CALL get_basis_offsets(particle_set, qs_kind_set, first_sgf, dimen_ORB)
890 :
891 : ! =========================================================================
892 : ! 1. SETUP DBCSR MATRIX TOPOLOGY
893 : ! =========================================================================
894 :
895 : ! A. Define Column Block Sizes (1 Block = 1 Atom's full basis set)
896 132 : ALLOCATE (c_blk_sizes(natom))
897 150 : DO iatom = 1, natom
898 150 : c_blk_sizes(iatom) = first_sgf(iatom + 1) - first_sgf(iatom)
899 : END DO
900 :
901 : ! B. Define Row Block Sizes: atom-aligned blocks (a block never spans two atoms' grid
902 : ! runs), each atom's run subdivided into blocks of at most bs_eff points.
903 :
904 : ! Fetch CP2K's default process grid configuration
905 44 : CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist_ks)
906 44 : CALL dbcsr_distribution_get(dbcsr_dist_ks, nprows=nprow, npcols=npcol)
907 :
908 : ! Overflow-safe upper bound on the block size (see bs_env%dbcsr_msg_elem_limit).
909 : safe_max = INT(0.5_dp*REAL(bs_env%dbcsr_msg_elem_limit, dp)* &
910 : REAL(MAX(MIN(nprow, npcol), 1), dp)/ &
911 44 : REAL(total_grid_npts, dp))
912 44 : safe_max = MAX(1, safe_max)
913 : ! Block size = CP2K's global max_elements_per_block (GLOBAL/DBCSR input; default 32),
914 : ! overflow-capped.
915 44 : bs_eff = MAX(1, MIN(max_elements_per_block, safe_max))
916 :
917 : ! Count the atom-aligned blocks, then fill r_blk_sizes and each block's starting grid row.
918 44 : num_grid_chunks = 0
919 150 : DO ia = 1, natom
920 106 : run = bs_env%ri_rs%grid_atom_boundaries(ia + 1) - bs_env%ri_rs%grid_atom_boundaries(ia)
921 150 : IF (run > 0) num_grid_chunks = num_grid_chunks + (run + bs_eff - 1)/bs_eff
922 : END DO
923 176 : ALLOCATE (r_blk_sizes(num_grid_chunks), blk_row_start(num_grid_chunks))
924 44 : i_blk = 0
925 44 : r_start = 1
926 150 : DO ia = 1, natom
927 106 : remaining = bs_env%ri_rs%grid_atom_boundaries(ia + 1) - bs_env%ri_rs%grid_atom_boundaries(ia)
928 744 : DO WHILE (remaining > 0)
929 594 : i_blk = i_blk + 1
930 594 : r_blk_sizes(i_blk) = MIN(bs_eff, remaining)
931 594 : blk_row_start(i_blk) = r_start
932 594 : r_start = r_start + r_blk_sizes(i_blk)
933 700 : remaining = remaining - r_blk_sizes(i_blk)
934 : END DO
935 : END DO
936 :
937 44 : IF (bs_env%unit_nr > 0) THEN
938 22 : WRITE (bs_env%unit_nr, '(T2,A,T69,I12)') 'RI-RS grid row-blocks of phi(mu,l)', num_grid_chunks
939 22 : WRITE (bs_env%unit_nr, '(T2,A,T69,I12)') 'RI-RS grid points per block (max)', bs_eff
940 : END IF
941 :
942 : ! Cache atomic positions: AO and RI blocks are one-block-per-atom, so these are the
943 : ! block centers used by the optional CUTOFF_RADIUS_G_W atom-pair truncation.
944 44 : IF (ALLOCATED(bs_env%ri_rs%atom_centers)) DEALLOCATE (bs_env%ri_rs%atom_centers)
945 132 : ALLOCATE (bs_env%ri_rs%atom_centers(3, natom))
946 150 : DO iatom = 1, natom
947 892 : bs_env%ri_rs%atom_centers(1:3, iatom) = particle_set(iatom)%r(1:3)
948 : END DO
949 :
950 : ! Cache per-chunk centroids for the optional CUTOFF_RADIUS_RL_W / CUTOFF_RADIUS_W0 block
951 : ! truncations. Left unallocated otherwise, so PRESENT(centroids) stays .FALSE. at the
952 : ! contract_grid_panels* calls.
953 44 : IF (bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp .OR. &
954 : bs_env%ri_rs%cutoff_radius_w0 > 0.0_dp) THEN
955 8 : IF (ALLOCATED(bs_env%ri_rs%chunk_centroids)) DEALLOCATE (bs_env%ri_rs%chunk_centroids)
956 24 : ALLOCATE (bs_env%ri_rs%chunk_centroids(3, num_grid_chunks))
957 110 : DO i_blk = 1, num_grid_chunks
958 102 : r_start = blk_row_start(i_blk)
959 102 : r_end = r_start + r_blk_sizes(i_blk) - 1
960 : bs_env%ri_rs%chunk_centroids(1, i_blk) = &
961 2950 : SUM(ri_rs_grid_points(1, r_start:r_end))/REAL(r_blk_sizes(i_blk), dp)
962 : bs_env%ri_rs%chunk_centroids(2, i_blk) = &
963 2950 : SUM(ri_rs_grid_points(2, r_start:r_end))/REAL(r_blk_sizes(i_blk), dp)
964 : bs_env%ri_rs%chunk_centroids(3, i_blk) = &
965 2994 : SUM(ri_rs_grid_points(3, r_start:r_end))/REAL(r_blk_sizes(i_blk), dp)
966 : END DO
967 : END IF
968 :
969 : ! C. Build Custom Mappings using Round-Robin across the 2D process grid
970 :
971 88 : ALLOCATE (row_dist(num_grid_chunks))
972 638 : DO i = 1, num_grid_chunks
973 638 : row_dist(i) = MOD(i - 1, nprow)
974 : END DO
975 :
976 88 : ALLOCATE (col_dist(natom))
977 150 : DO i = 1, natom
978 150 : col_dist(i) = MOD(i - 1, npcol)
979 : END DO
980 :
981 : ! E. Create the DBCSR Distribution and Initialize the Matrix
982 : CALL dbcsr_distribution_new(dist, template=dbcsr_dist_ks, &
983 44 : row_dist=row_dist, col_dist=col_dist)
984 :
985 : CALL dbcsr_create(mat_phi_mu_l, name="phi_val_sparse", dist=dist, &
986 : matrix_type=dbcsr_type_no_symmetry, &
987 44 : row_blk_size=r_blk_sizes, col_blk_size=c_blk_sizes)
988 :
989 : ! =========================================================================
990 : ! 2. STREAM DATA DIRECTLY INTO SPARSE MATRIX
991 : ! =========================================================================
992 : ! Iterate over the atoms assigned to this specific MPI rank
993 44 : DO iatom = para_env%mepos + 1, natom, para_env%num_pe
994 :
995 53 : c_size = c_blk_sizes(iatom)
996 :
997 : ! Allocate a temporary dense buffer just for this specific atom
998 212 : ALLOCATE (atom_col_buffer(total_grid_npts, c_size))
999 53 : atom_col_buffer = 0.0_dp
1000 :
1001 : ! Evaluate the basis functions on the grid. Skip grid points outside
1002 : ! the spatial extent of the most diffuse AO Gaussian on iatom; beyond
1003 : ! that radius the contribution is guaranteed below eps_filter. A positive
1004 : ! CUTOFF_RADIUS_RL_AO overrides this with a user-defined hard cutoff.
1005 53 : IF (bs_env%ri_rs%cutoff_radius_ri_ao > 0.0_dp) THEN
1006 3 : r2_threshold = bs_env%ri_rs%cutoff_radius_ri_ao**2
1007 : ELSE
1008 50 : r2_threshold = bs_env%ri_rs%radius_ao_per_atom(iatom)**2
1009 : END IF
1010 : CALL fill_phi_for_atom(atom_col_buffer, ri_rs_grid_points, total_grid_npts, &
1011 53 : iatom, particle_set, qs_kind_set, cell, r2_threshold)
1012 :
1013 : ! Slice the dense column into the atom-aligned grid row-blocks and insert into DBCSR
1014 786 : DO i_blk = 1, num_grid_chunks
1015 733 : r_start = blk_row_start(i_blk)
1016 733 : r_end = r_start + r_blk_sizes(i_blk) - 1
1017 :
1018 : ! Apply dynamic sparsity filtering: Only store blocks with physical significance
1019 81503 : IF (MAXVAL(ABS(atom_col_buffer(r_start:r_end, 1:c_size))) > bs_env%eps_filter) THEN
1020 : CALL dbcsr_put_block(mat_phi_mu_l, row=i_blk, col=iatom, &
1021 721 : block=atom_col_buffer(r_start:r_end, 1:c_size))
1022 : END IF
1023 : END DO
1024 :
1025 53 : DEALLOCATE (atom_col_buffer)
1026 :
1027 : END DO
1028 :
1029 44 : CALL dbcsr_finalize(mat_phi_mu_l)
1030 :
1031 44 : CALL print_matrix_occupation(mat_phi_mu_l, 'φ(μ,l)', bs_env)
1032 :
1033 : ! -------------------------------------------------------------------------
1034 : ! CLEANUP
1035 : ! -------------------------------------------------------------------------
1036 44 : DEALLOCATE (first_sgf, r_blk_sizes, c_blk_sizes, row_dist, col_dist, blk_row_start)
1037 44 : CALL dbcsr_distribution_release(dist)
1038 :
1039 44 : CALL timestop(handle)
1040 :
1041 220 : END SUBROUTINE atomic_basis_at_grid_point
1042 :
1043 : ! **************************************************************************************************
1044 : !> \brief Evaluates all spherical AO basis functions of one atom on a set of grid points and
1045 : !> ACCUMULATES them into phi_val (+=). For each point within the cutoff radius,
1046 : !> Φ_μ(r) = Σ_pgf Σ_cart sphi(cart,μ) · (x-X_A)^lx (y-Y_A)^ly (z-Z_A)^lz · e^(-ζ_pgf |r-R_A|²)
1047 : !> i.e. contracted Cartesian Gaussians transformed to the spherical basis via the sphi
1048 : !> coefficients. Distances are minimum-image wrapped (pbc); points with
1049 : !> |r - R_A|² > r2_threshold are skipped since every primitive is below eps there.
1050 : !> OMP-parallel over grid points.
1051 : !> \param phi_val ...
1052 : !> \param ri_rs_grid ...
1053 : !> \param npts ...
1054 : !> \param iatom ...
1055 : !> \param particle_set ...
1056 : !> \param qs_kind_set ...
1057 : !> \param cell ...
1058 : !> \param r2_threshold ...
1059 : ! **************************************************************************************************
1060 177 : SUBROUTINE fill_phi_for_atom(phi_val, ri_rs_grid, npts, iatom, &
1061 : particle_set, qs_kind_set, cell, r2_threshold)
1062 :
1063 : REAL(KIND=dp), INTENT(INOUT) :: phi_val(:, :)
1064 : INTEGER, INTENT(IN) :: npts
1065 : REAL(KIND=dp), INTENT(IN) :: ri_rs_grid(3, npts)
1066 : INTEGER, INTENT(IN) :: iatom
1067 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1068 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1069 : TYPE(cell_type), POINTER :: cell
1070 : REAL(KIND=dp), INTENT(IN) :: r2_threshold
1071 :
1072 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_phi_for_atom'
1073 :
1074 : INTEGER :: handle
1075 :
1076 177 : CALL timeset(routineN, handle)
1077 :
1078 : CALL evaluate_ao_on_points(phi_val, ri_rs_grid, iatom, particle_set, qs_kind_set, cell, &
1079 177 : cutoff_squared=r2_threshold)
1080 :
1081 177 : CALL timestop(handle)
1082 :
1083 177 : END SUBROUTINE fill_phi_for_atom
1084 :
1085 : ! **************************************************************************************************
1086 : !> \brief Computes the AO basis offsets: first_sgf(iatom) is the global index of the first
1087 : !> spherical Gaussian function (SGF) of iatom, first_sgf(natom+1) = total_sgf + 1,
1088 : !> and total_sgf is the total number of AO basis functions.
1089 : !> \param particle_set ...
1090 : !> \param qs_kind_set ...
1091 : !> \param first_sgf ...
1092 : !> \param total_sgf ...
1093 : ! **************************************************************************************************
1094 44 : SUBROUTINE get_basis_offsets(particle_set, qs_kind_set, first_sgf, total_sgf)
1095 :
1096 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1097 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1098 : INTEGER, INTENT(OUT) :: first_sgf(:), total_sgf
1099 :
1100 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_basis_offsets'
1101 :
1102 : INTEGER :: handle, iatom, ikind, nsgf
1103 :
1104 44 : CALL timeset(routineN, handle)
1105 :
1106 44 : total_sgf = 0
1107 150 : DO iatom = 1, SIZE(particle_set)
1108 106 : first_sgf(iatom) = total_sgf + 1
1109 106 : ikind = particle_set(iatom)%atomic_kind%kind_number
1110 106 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf, basis_type="ORB")
1111 150 : total_sgf = total_sgf + nsgf
1112 : END DO
1113 44 : first_sgf(SIZE(particle_set) + 1) = total_sgf + 1
1114 :
1115 44 : CALL timestop(handle)
1116 :
1117 44 : END SUBROUTINE get_basis_offsets
1118 :
1119 : ! **************************************************************************************************
1120 : !> \brief Computes the RI-RS fitting coefficients Z_lP by solving, independently for every RI
1121 : !> atom P, a Jacobi-conditioned, Tikhonov-regularized linear system restricted to the
1122 : !> grid points r_l inside P's integration sphere |r_l - R_P| <= cutoff_ri(P):
1123 : !> D_ll' = [ Σ_μ Φ_μ(r_l) Φ_μ(r_l') ]² (squared grid Gram matrix, Eq. 13)
1124 : !> D_lP = Σ_μν Φ_μ(r_l) Φ_ν(r_l) (μν|P) (grid-RI right-hand side, Eq. 15)
1125 : !> d_l = 1 / sqrt(D_ll) (Jacobi conditioning vector)
1126 : !> D'_ll' = d_l D_ll' d_l' + λ δ_ll' (λ = TIKHONOV_SIGMA regularization)
1127 : !> Σ_l' D'_ll' Z'_l'P = d_l D_lP (Cholesky solve, Eq. 14)
1128 : !> Z_lP = d_l Z'_l'P (undo the conditioning)
1129 : !> Work is distributed over atoms in two phases (planned by classify_z_lp_atoms and
1130 : !> lpt_assign_atoms): Phase A solves "small" atoms with single-rank LAPACK
1131 : !> (dpotrf/dpotrs); Phase B solves "big" atoms, whose dense Gram matrix would exceed one
1132 : !> rank's memory, with ScaLAPACK (pdpotrf/pdpotrs) over rank subgroups of size G.
1133 : !> The solved Z columns are scattered into the sparse global mat_Z_lP.
1134 : !> If a Z_lP restart file exists, it is read instead and the solve is skipped entirely.
1135 : !> \param qs_env ...
1136 : !> \param bs_env ...
1137 : !> \param ri_rs_grid_points ...
1138 : !> \param mat_phi_mu_l ...
1139 : !> \param mat_Z_lP ...
1140 : ! **************************************************************************************************
1141 88 : SUBROUTINE compute_coeff_Z_lP(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
1142 :
1143 : TYPE(qs_environment_type), POINTER :: qs_env
1144 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1145 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
1146 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
1147 : TYPE(dbcsr_type), INTENT(OUT) :: mat_Z_lP
1148 :
1149 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
1150 : routineN = 'compute_coeff_Z_lP'
1151 :
1152 : INTEGER :: atom_j_mepos, atom_j_stride, atom_P, G, handle, handle_dpotrf, handle_dpotrs, &
1153 : i_blk, idx, info, iphase, j, max_ao_size, my_group, n_ao_total, n_ao_used, n_big, n_done, &
1154 : n_groups, n_loc_ri, n_local_grid, n_my_atoms, n_small, natom, next_pct, npcol_phi, &
1155 : num_grid_chunks, P_loop_atom, phase_hi
1156 88 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_col_map, big_list, local_grid_idx, &
1157 44 : my_atoms_A, my_atoms_B, &
1158 44 : n_local_grid_atom, row_offset, &
1159 44 : small_list
1160 88 : INTEGER, DIMENSION(:), POINTER :: col_dist_ri, r_blk_sizes, ri_blk_sizes, &
1161 44 : row_dist_grid
1162 : LOGICAL :: do_scatter, use_dist
1163 : REAL(KIND=dp) :: balance_A, balance_B, cutoff_ri, r_c, t1
1164 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cutoff_ri_per_atom, d_vec_local
1165 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: D_local, d_lp_local, phi_local
1166 44 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1167 : TYPE(cell_type), POINTER :: cell
1168 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
1169 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_b, fm_struct_D
1170 : TYPE(cp_fm_type) :: fm_b, fm_D
1171 : TYPE(cp_logger_type), POINTER :: logger
1172 : TYPE(dbcsr_distribution_type) :: dist_phi, dist_Z
1173 572 : TYPE(gw_3c_ctx_type) :: ctx_3c
1174 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
1175 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1176 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1177 : TYPE(section_vals_type), POINTER :: input
1178 :
1179 44 : CALL timeset(routineN, handle)
1180 :
1181 44 : t1 = m_walltime()
1182 :
1183 : CALL get_qs_env(qs_env, para_env=para_env, particle_set=particle_set, input=input, &
1184 44 : qs_kind_set=qs_kind_set, cell=cell, atomic_kind_set=atomic_kind_set)
1185 :
1186 44 : NULLIFY (para_env_sub, blacs_env_sub)
1187 :
1188 44 : natom = SIZE(bs_env%i_RI_start_from_atom)
1189 44 : n_ao_total = bs_env%i_ao_end_from_atom(natom)
1190 :
1191 : ! =========================================================================
1192 : ! 1. SETUP DBCSR TOPOLOGY & EXACT OFFSETS
1193 : ! mat_Z_lP inherits the grid row blocking (and row distribution) of
1194 : ! mat_phi_mu_l; its columns are one block per RI atom.
1195 : ! =========================================================================
1196 44 : CALL dbcsr_get_info(mat_phi_mu_l, row_blk_size=r_blk_sizes, distribution=dist_phi)
1197 44 : CALL dbcsr_distribution_get(dist_phi, row_dist=row_dist_grid, npcols=npcol_phi)
1198 :
1199 44 : num_grid_chunks = SIZE(r_blk_sizes)
1200 :
1201 132 : ALLOCATE (row_offset(num_grid_chunks))
1202 44 : row_offset(1) = 0
1203 594 : DO i_blk = 2, num_grid_chunks
1204 594 : row_offset(i_blk) = row_offset(i_blk - 1) + r_blk_sizes(i_blk - 1)
1205 : END DO
1206 :
1207 176 : ALLOCATE (ri_blk_sizes(natom), col_dist_ri(natom))
1208 150 : DO atom_P = 1, natom
1209 106 : ri_blk_sizes(atom_P) = bs_env%i_RI_end_from_atom(atom_P) - bs_env%i_RI_start_from_atom(atom_P) + 1
1210 150 : col_dist_ri(atom_P) = MOD(atom_P - 1, npcol_phi)
1211 : END DO
1212 :
1213 44 : CALL dbcsr_distribution_new(dist_Z, template=dist_phi, row_dist=row_dist_grid, col_dist=col_dist_ri)
1214 :
1215 44 : IF (bs_env%ri_rs%Z_lP_exists) THEN
1216 : CALL dbcsr_binary_read(filepath=TRIM(bs_env%prefix)//"Z_lP.matrix", &
1217 : distribution=dist_Z, &
1218 2 : matrix_new=mat_Z_lP)
1219 2 : IF (bs_env%unit_nr > 0) THEN
1220 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
1221 1 : 'Read Z_lP from file ', ' Execution time', m_walltime() - t1, ' s'
1222 : ! The grid rows are laid out in Morton (spatial) order (spatial_atom_order); a Z_lP.matrix
1223 : ! written by an older build with a different grid ordering would be silently misread into
1224 : ! the current row order. Delete stale Z_lP.matrix files and recompute if in doubt.
1225 : WRITE (bs_env%unit_nr, '(T2,A)') &
1226 1 : '*** NOTE: Z_lP restart must match the current (spatial) grid row ordering ***'
1227 1 : WRITE (bs_env%unit_nr, '(A)') ' '
1228 : END IF
1229 : ELSE
1230 :
1231 : CALL dbcsr_create(mat_Z_lP, name="mat_Z_lP", dist=dist_Z, &
1232 : matrix_type=dbcsr_type_no_symmetry, &
1233 42 : row_blk_size=r_blk_sizes, col_blk_size=ri_blk_sizes)
1234 :
1235 : ! Largest per-atom AO block, needed to size the 3c-integral work buffers.
1236 42 : max_ao_size = 0
1237 142 : DO j = 1, SIZE(bs_env%i_ao_start_from_atom)
1238 142 : max_ao_size = MAX(max_ao_size, bs_env%i_ao_end_from_atom(j) - bs_env%i_ao_start_from_atom(j) + 1)
1239 : END DO
1240 :
1241 : ! Per-atom RI-RS integration sphere:
1242 : ! cutoff_ri(P) = r_c + r_RI(P)
1243 : ! where r_c is the truncated-Coulomb cutoff of the RI metric and r_RI the radius of
1244 : ! the most diffuse RI auxiliary Gaussian on P. The CUTOFF_RADIUS_RL_RI keyword
1245 : ! (when > 0) overrides the entire cutoff calculation.
1246 126 : ALLOCATE (cutoff_ri_per_atom(natom))
1247 :
1248 42 : IF (bs_env%ri_rs%cutoff_radius_ri_rs > 0.0_dp) THEN
1249 8 : cutoff_ri_per_atom(:) = bs_env%ri_rs%cutoff_radius_ri_rs
1250 : ELSE
1251 40 : r_c = bs_env%ri_metric%cutoff_radius
1252 134 : DO P_loop_atom = 1, natom
1253 : cutoff_ri_per_atom(P_loop_atom) = &
1254 134 : r_c + bs_env%ri_rs%radius_ri_per_atom(P_loop_atom)
1255 : END DO
1256 : END IF
1257 :
1258 42 : CALL print_sphere_cutoff_table(bs_env, atomic_kind_set, particle_set, cutoff_ri_per_atom)
1259 :
1260 : ! =========================================================================
1261 : ! 2. PER-ATOM SOLVER CLASSIFICATION
1262 : ! Split the atoms into "small" (single-rank LAPACK, Phase A) and "big"
1263 : ! (distributed ScaLAPACK over subgroups of G ranks, Phase B) by comparing
1264 : ! each atom's estimated solve peak memory against the measured budget.
1265 : ! =========================================================================
1266 : CALL classify_z_lp_atoms(bs_env, para_env, ri_rs_grid_points, particle_set, &
1267 : cutoff_ri_per_atom, ri_blk_sizes, &
1268 42 : n_local_grid_atom, small_list, n_small, big_list, n_big, G)
1269 :
1270 : ! LPT scheduling: sort the atoms of each phase by estimated solve cost
1271 : ! (n_local_grid^3, Cholesky-dominated) and greedily assign to the least-loaded
1272 : ! rank (Phase A) / subgroup (Phase B).
1273 : CALL lpt_assign_atoms(small_list, n_small, n_local_grid_atom, para_env%num_pe, &
1274 42 : para_env%mepos, my_atoms_A, balance_A)
1275 42 : IF (n_big > 0) THEN
1276 0 : n_groups = para_env%num_pe/G
1277 0 : my_group = MIN(para_env%mepos/G, n_groups - 1)
1278 : CALL lpt_assign_atoms(big_list, n_big, n_local_grid_atom, n_groups, my_group, &
1279 0 : my_atoms_B, balance_B)
1280 : ELSE
1281 42 : ALLOCATE (my_atoms_B(0))
1282 42 : balance_B = 1.0_dp
1283 : END IF
1284 :
1285 : ! Atoms this rank will process across both phases for rank-0 progress
1286 42 : n_my_atoms = SIZE(my_atoms_A) + SIZE(my_atoms_B)
1287 42 : n_done = 0
1288 42 : next_pct = 25
1289 :
1290 42 : IF (bs_env%unit_nr > 0) THEN
1291 : WRITE (bs_env%unit_nr, '(T2,A,I7,A,I8,A)') &
1292 21 : 'RI-RS Z_lP solver: ', n_small, ' atoms single-rank (BLAS), ', n_big, &
1293 42 : ' atoms distributed'
1294 21 : IF (n_small > 0) WRITE (bs_env%unit_nr, '(T4,A,F18.2)') &
1295 21 : 'estimated single-rank load balance (max/mean cost per rank)', balance_A
1296 21 : IF (n_big > 0) THEN
1297 : WRITE (bs_env%unit_nr, '(T4,A,I44,A)') &
1298 0 : 'distributed subgroup size G', G, ' ranks'
1299 : WRITE (bs_env%unit_nr, '(T4,A,F17.2)') &
1300 0 : 'estimated distributed load balance (max/mean cost per group)', balance_B
1301 : END IF
1302 21 : WRITE (bs_env%unit_nr, '(A)') ' '
1303 : END IF
1304 :
1305 : ! Shared context for the three-center integrals (μν|P) of the RHS build
1306 : CALL gw_3c_ctx_create(ctx_3c, qs_env, bs_env%ri_metric, &
1307 : basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, &
1308 42 : basis_i=bs_env%basis_set_RI)
1309 :
1310 : ! =========================================================================
1311 : ! 3. TWO-PHASE LOOP OVER ATOMS
1312 : ! Phase A processes the "small" atoms with the single-rank BLAS path
1313 : ! Phase B processes the "big" atoms with the distributed ScaLAPACK path
1314 : ! over rank subgroups of size G. phi_local for each atom's cutoff sphere
1315 : ! is built on the fly to avoid replicating a global grid x AO matrix.
1316 : ! =========================================================================
1317 126 : DO iphase = 1, 2
1318 84 : IF (iphase == 1) THEN
1319 42 : use_dist = .FALSE.
1320 42 : atom_j_mepos = 0
1321 42 : atom_j_stride = 1
1322 42 : phase_hi = SIZE(my_atoms_A)
1323 : ELSE
1324 42 : IF (n_big == 0) CYCLE
1325 0 : use_dist = .TRUE.
1326 0 : n_groups = para_env%num_pe/G
1327 0 : my_group = MIN(para_env%mepos/G, n_groups - 1)
1328 0 : ALLOCATE (para_env_sub)
1329 0 : CALL para_env_sub%from_split(para_env, my_group)
1330 0 : CALL cp_blacs_env_create(blacs_env=blacs_env_sub, para_env=para_env_sub)
1331 0 : atom_j_mepos = para_env_sub%mepos
1332 0 : atom_j_stride = para_env_sub%num_pe
1333 : ! All ranks of a subgroup share my_group, hence the identical my_atoms_B list
1334 : ! (the per-atom ScaLAPACK solve is collective over the subgroup).
1335 0 : phase_hi = SIZE(my_atoms_B)
1336 : END IF
1337 :
1338 92 : DO idx = 1, phase_hi
1339 50 : IF (iphase == 1) THEN
1340 50 : atom_P = my_atoms_A(idx)
1341 : ELSE
1342 0 : atom_P = my_atoms_B(idx)
1343 : END IF
1344 :
1345 50 : n_loc_ri = ri_blk_sizes(atom_P)
1346 50 : cutoff_ri = cutoff_ri_per_atom(atom_P)
1347 :
1348 : ! ---------------------------------------------------------------------
1349 : ! A. Sphere-local AO matrix Φ_μ(r_l): select the grid points with
1350 : ! |r_l - R_P| <= cutoff_ri(P), evaluate every AO on them, and drop
1351 : ! points whose largest AO amplitude is below EPS_FILTER.
1352 : ! ---------------------------------------------------------------------
1353 : CALL build_phi_on_sphere(bs_env, particle_set, qs_kind_set, cell, &
1354 : ri_rs_grid_points, atom_P, cutoff_ri, n_ao_total, &
1355 : local_grid_idx, n_local_grid, phi_local, &
1356 50 : ao_col_map, n_ao_used)
1357 :
1358 : ! ---------------------------------------------------------------------
1359 : ! B. Right-hand side D_lP = Σ_μν Φ_μ(r_l) Φ_ν(r_l) (μν|P)
1360 : ! ---------------------------------------------------------------------
1361 200 : ALLOCATE (d_lp_local(n_local_grid, n_loc_ri))
1362 50 : d_lp_local = 0.0_dp
1363 :
1364 : CALL compute_d_lp(bs_env, ctx_3c, phi_local, ao_col_map, d_lp_local, n_local_grid, &
1365 50 : n_loc_ri, atom_P, max_ao_size, atom_j_mepos, atom_j_stride)
1366 :
1367 : ! Reduce per-subgroup-rank partials into the replicated d_lp_local.
1368 : ! Skipped for BLAS path: each rank has the full sum locally.
1369 50 : IF (use_dist) THEN
1370 0 : CALL para_env_sub%sum(d_lp_local)
1371 : END IF
1372 :
1373 : ! ---------------------------------------------------------------------
1374 : ! C. Jacobi conditioning vector d_l = 1/sqrt(D_ll) and, on the BLAS path,
1375 : ! the dense conditioned Gram matrix D'_ll' = d_l D_ll' d_l' + λδ_ll'.
1376 : ! ---------------------------------------------------------------------
1377 150 : ALLOCATE (d_vec_local(n_local_grid))
1378 :
1379 50 : IF (.NOT. use_dist) THEN
1380 : CALL build_gram_jacobi_blas(phi_local, n_local_grid, n_ao_used, &
1381 50 : bs_env%ri_rs%tikhonov, D_local, d_vec_local)
1382 : ELSE
1383 : ! ScaLAPACK path: only d_vec is needed here (= 1/||phi(r_l)||^2);
1384 : ! solve_D_lp_distributed builds its block-cyclic slice of D' internally
1385 : ! with the squared+scaled values, so no dense D_local on this rank.
1386 : CALL build_jacobi_diag_from_phi(phi_local, n_local_grid, n_ao_used, &
1387 0 : d_vec_local)
1388 : END IF
1389 :
1390 : ! ---------------------------------------------------------------------
1391 : ! D. Pre-scale the RHS: D'_lP = d_l * D_lP
1392 : ! ---------------------------------------------------------------------
1393 50 : CALL scale_rows_by_diag(d_lp_local, d_vec_local, n_local_grid, n_loc_ri)
1394 :
1395 : ! ---------------------------------------------------------------------
1396 : ! E. Cholesky solve Σ_l' D'_ll' Z'_l'P = D'_lP
1397 : ! (BLAS dpotrf/dpotrs or ScaLAPACK pdpotrf/pdpotrs)
1398 : ! ---------------------------------------------------------------------
1399 50 : IF (.NOT. use_dist) THEN
1400 50 : CALL timeset(routineN//"_dpotrf", handle_dpotrf)
1401 50 : CALL dpotrf('L', n_local_grid, D_local, n_local_grid, info)
1402 50 : CALL timestop(handle_dpotrf)
1403 50 : CALL timeset(routineN//"_dpotrs", handle_dpotrs)
1404 : CALL dpotrs('L', n_local_grid, n_loc_ri, D_local, n_local_grid, &
1405 50 : d_lp_local, n_local_grid, info)
1406 50 : CALL timestop(handle_dpotrs)
1407 50 : DEALLOCATE (D_local)
1408 : ELSE
1409 : CALL solve_D_lp_distributed(phi_local, d_vec_local, d_lp_local, &
1410 : n_local_grid, n_ao_used, n_loc_ri, &
1411 : bs_env%ri_rs%tikhonov, &
1412 : para_env_sub, blacs_env_sub, &
1413 0 : fm_struct_D, fm_struct_b, fm_D, fm_b, info)
1414 : END IF
1415 :
1416 : ! ---------------------------------------------------------------------
1417 : ! F. Undo the conditioning: Z_lP = d_l * Z'_lP
1418 : ! ---------------------------------------------------------------------
1419 50 : CALL scale_rows_by_diag(d_lp_local, d_vec_local, n_local_grid, n_loc_ri)
1420 :
1421 : ! ---------------------------------------------------------------------
1422 : ! G. Scatter the solved Z columns back into the global sparse mat_Z_lP.
1423 : ! ---------------------------------------------------------------------
1424 50 : do_scatter = .TRUE.
1425 50 : IF (use_dist) do_scatter = (para_env_sub%mepos == 0)
1426 0 : IF (do_scatter) THEN
1427 : CALL scatter_z_columns(mat_Z_lP, d_lp_local, local_grid_idx, n_local_grid, &
1428 : n_loc_ri, atom_P, r_blk_sizes, row_offset, &
1429 50 : bs_env%eps_filter)
1430 : END IF
1431 :
1432 50 : DEALLOCATE (d_vec_local, d_lp_local)
1433 50 : DEALLOCATE (local_grid_idx, phi_local, ao_col_map)
1434 :
1435 : ! Progress based on rank 0 only
1436 50 : n_done = n_done + 1
1437 92 : IF (bs_env%unit_nr > 0 .AND. n_my_atoms > 0) THEN
1438 112 : DO WHILE (next_pct <= 100 .AND. n_done*100 >= next_pct*n_my_atoms)
1439 84 : WRITE (bs_env%unit_nr, '(T2,A,I57,A)') 'Computing Z_lP:', next_pct, ' % done'
1440 91 : next_pct = next_pct + 25
1441 : END DO
1442 : END IF
1443 : END DO ! idx: atoms of this phase owned by this rank / subgroup
1444 :
1445 : ! Tear down the Phase-B subgroup (all ranks created it collectively).
1446 84 : IF (iphase == 2) THEN
1447 0 : CALL cp_blacs_env_release(blacs_env_sub)
1448 0 : CALL para_env_sub%free()
1449 0 : DEALLOCATE (para_env_sub)
1450 : END IF
1451 : END DO ! iphase
1452 :
1453 42 : DEALLOCATE (cutoff_ri_per_atom)
1454 42 : DEALLOCATE (small_list, big_list)
1455 :
1456 42 : CALL gw_3c_ctx_release(ctx_3c)
1457 :
1458 42 : CALL dbcsr_finalize(mat_Z_lP)
1459 :
1460 42 : CALL print_matrix_occupation(mat_Z_lP, 'Z(l,P)', bs_env)
1461 :
1462 42 : IF (bs_env%unit_nr > 0) THEN
1463 21 : WRITE (bs_env%unit_nr, '(A)') ' '
1464 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
1465 21 : 'Computed Z_lP ', ' Execution time', m_walltime() - t1, ' s'
1466 21 : WRITE (bs_env%unit_nr, '(A)') ' '
1467 : END IF
1468 :
1469 42 : logger => cp_get_default_logger()
1470 :
1471 42 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
1472 6 : CALL dbcsr_binary_write(matrix=mat_Z_lP, filepath=TRIM(bs_env%prefix)//"Z_lP.matrix")
1473 : END IF
1474 :
1475 : END IF
1476 :
1477 44 : DEALLOCATE (row_offset, ri_blk_sizes, col_dist_ri)
1478 44 : CALL dbcsr_distribution_release(dist_Z)
1479 :
1480 44 : DEALLOCATE (ri_rs_grid_points)
1481 :
1482 44 : CALL timestop(handle)
1483 :
1484 220 : END SUBROUTINE compute_coeff_Z_lP
1485 :
1486 : ! **************************************************************************************************
1487 : !> \brief Prints the per-kind maximum RI-RS integration-sphere cutoff table.
1488 : !> \param bs_env ...
1489 : !> \param atomic_kind_set ...
1490 : !> \param particle_set ...
1491 : !> \param cutoff_ri_per_atom ...
1492 : ! **************************************************************************************************
1493 42 : SUBROUTINE print_sphere_cutoff_table(bs_env, atomic_kind_set, particle_set, cutoff_ri_per_atom)
1494 :
1495 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1496 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1497 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1498 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutoff_ri_per_atom
1499 :
1500 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_sphere_cutoff_table'
1501 :
1502 : INTEGER :: handle, iatom, ikind, nkind
1503 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cutoff_ri_per_kind
1504 :
1505 42 : CALL timeset(routineN, handle)
1506 :
1507 42 : IF (bs_env%unit_nr <= 0) THEN
1508 21 : CALL timestop(handle)
1509 : RETURN
1510 : END IF
1511 :
1512 21 : nkind = SIZE(atomic_kind_set)
1513 63 : ALLOCATE (cutoff_ri_per_kind(nkind))
1514 21 : cutoff_ri_per_kind(:) = 0.0_dp
1515 :
1516 71 : DO iatom = 1, SIZE(particle_set)
1517 50 : ikind = particle_set(iatom)%atomic_kind%kind_number
1518 71 : cutoff_ri_per_kind(ikind) = MAX(cutoff_ri_per_kind(ikind), cutoff_ri_per_atom(iatom))
1519 : END DO
1520 :
1521 21 : WRITE (bs_env%unit_nr, '(T2,A)') 'Per-kind maximum RI-RS sphere cutoff (Å):'
1522 21 : WRITE (bs_env%unit_nr, '(T4,A4,A14)') 'Kind', 'cutoff (Å)'
1523 50 : DO ikind = 1, nkind
1524 : WRITE (bs_env%unit_nr, '(T4,A4,F14.4)') &
1525 29 : atomic_kind_set(ikind)%element_symbol, &
1526 79 : cutoff_ri_per_kind(ikind)*angstrom
1527 : END DO
1528 21 : WRITE (bs_env%unit_nr, '(A)') ' '
1529 :
1530 21 : DEALLOCATE (cutoff_ri_per_kind)
1531 :
1532 21 : CALL timestop(handle)
1533 :
1534 42 : END SUBROUTINE print_sphere_cutoff_table
1535 :
1536 : ! **************************************************************************************************
1537 : !> \brief Splits the atoms of the Z_lP solve into a single-rank list ("small", Phase A: LAPACK
1538 : !> dpotrf/dpotrs on one rank) and a distributed list ("big", Phase B: ScaLAPACK
1539 : !> pdpotrf/pdpotrs over a rank subgroup of size G), and sizes G.
1540 : !> AUTO mode (N_PROCS_PER_ATOM_Z_LP <= 0, the default): estimate each atom's single-rank
1541 : !> peak memory
1542 : !> peak(P) = 8*n_local_grid(P)^2 (dense Gram matrix D_local)
1543 : !> + 8*n_local_grid(P)*n_ao_used(P) (phi_local)
1544 : !> + 8*n_local_grid(P)*n_RI(P)*(1+n_threads) (d_lp + OMP partials)
1545 : !> and send atoms whose peak exceeds mem_safety * available-memory-per-proc to the
1546 : !> distributed path; G is auto-sized so the biggest atom's distributed D_local (/G)
1547 : !> fits alongside the replicated phi_local + d_lp.
1548 : !> MANUAL mode (> 0): 1 forces the single-rank path for every atom; > 1 keeps the
1549 : !> memory-based classification but forces that fixed subgroup size G.
1550 : !> In every mode G is floored by the ScaLAPACK 32-bit index limit (a local block-cyclic
1551 : !> slice of ~n_local_grid^2/G elements must stay below 2^31 or pdpotrf segfaults).
1552 : !> \param bs_env ...
1553 : !> \param para_env ...
1554 : !> \param ri_rs_grid_points ...
1555 : !> \param particle_set ...
1556 : !> \param cutoff_ri_per_atom ...
1557 : !> \param ri_blk_sizes per-atom ...
1558 : !> \param n_local_grid_atom ...
1559 : !> \param small_list ...
1560 : !> \param n_small ...
1561 : !> \param big_list ...
1562 : !> \param n_big ...
1563 : !> \param G ...
1564 : ! **************************************************************************************************
1565 42 : SUBROUTINE classify_z_lp_atoms(bs_env, para_env, ri_rs_grid_points, particle_set, &
1566 42 : cutoff_ri_per_atom, ri_blk_sizes, &
1567 : n_local_grid_atom, small_list, n_small, big_list, n_big, G)
1568 :
1569 : !$ USE OMP_LIB, ONLY: omp_get_max_threads
1570 :
1571 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1572 : TYPE(mp_para_env_type), POINTER :: para_env
1573 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ri_rs_grid_points
1574 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1575 :
1576 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutoff_ri_per_atom
1577 : INTEGER, DIMENSION(:), INTENT(IN) :: ri_blk_sizes
1578 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: n_local_grid_atom, small_list, big_list
1579 : INTEGER, INTENT(OUT) :: n_small, n_big, G
1580 : CHARACTER(LEN=*), PARAMETER :: routineN = 'classify_z_lp_atoms'
1581 :
1582 : INTEGER :: handle
1583 :
1584 : ! Conservative fraction of measured available memory usable per rank for the Z_lP
1585 : REAL(KIND=dp), PARAMETER :: mem_safety = 0.8_dp
1586 :
1587 : ! ScaLAPACK/BLACS index the per-rank local block-cyclic slice (~n_local_grid^2/G
1588 : ! elements) with 32-bit integers; keep it safely below 2^31 or pdpotrf segfaults.
1589 : REAL(KIND=dp), PARAMETER :: scalapack_loc_limit = 2.0E9_dp
1590 :
1591 : INTEGER :: G_atom, G_int32, G_int32_max, l, &
1592 : n_ao_used_atom, n_grid_total, &
1593 : n_local_grid, natom, nthreads_cls, &
1594 : P_loop_atom
1595 : LOGICAL :: auto_mode
1596 : REAL(KIND=dp) :: budget_bytes, cutoff_ri, dlp_bytes, &
1597 : mem_avail_GB, ng, nri, peak_bytes, &
1598 : phi_bytes
1599 : REAL(KIND=dp), DIMENSION(3) :: pos_P
1600 :
1601 42 : CALL timeset(routineN, handle)
1602 :
1603 42 : natom = SIZE(particle_set)
1604 42 : n_grid_total = SIZE(ri_rs_grid_points, 2)
1605 :
1606 : ! Per-atom sphere size: n_local_grid(P) = number of l with |r_l - R_P| <= cutoff_ri(P).
1607 : ! It sets both the memory footprint (D_local is n_local_grid^2) and the solve
1608 : ! cost (~n_local_grid^3), so it drives classification and the LPT load balancing.
1609 210 : ALLOCATE (n_local_grid_atom(natom), small_list(natom), big_list(natom))
1610 142 : DO P_loop_atom = 1, natom
1611 400 : pos_P(:) = particle_set(P_loop_atom)%r(:)
1612 100 : cutoff_ri = cutoff_ri_per_atom(P_loop_atom)
1613 100 : n_local_grid = 0
1614 37524 : DO l = 1, n_grid_total
1615 149796 : IF (SUM((ri_rs_grid_points(1:3, l) - pos_P(1:3))**2) <= cutoff_ri**2) THEN
1616 36038 : n_local_grid = n_local_grid + 1
1617 : END IF
1618 : END DO
1619 142 : n_local_grid_atom(P_loop_atom) = n_local_grid
1620 : END DO
1621 :
1622 42 : nthreads_cls = 1
1623 42 : !$ nthreads_cls = omp_get_max_threads()
1624 : ! N_PROCS_PER_ATOM_Z_LP: -1 (default) = AUTO (classify by memory, auto-size G);
1625 : ! 1 = force single-rank BLAS for every atom; >1 = classify by memory but use this
1626 : ! fixed subgroup size G for the big atoms.
1627 42 : auto_mode = (bs_env%ri_rs%n_procs_per_atom_z_lp <= 0)
1628 42 : CALL ri_rs_mem_avail_per_proc_GB(bs_env, mem_avail_GB) ! collective over all ranks
1629 42 : budget_bytes = mem_safety*mem_avail_GB*1.0E9_dp
1630 :
1631 42 : n_small = 0
1632 42 : n_big = 0
1633 42 : G = 1
1634 42 : G_atom = 1 ! max G a big atom needs (memory + ScaLAPACK int32 floor)
1635 42 : G_int32_max = 1 ! max ScaLAPACK-int32 floor over the distributed atoms
1636 42 : IF (bs_env%ri_rs%n_procs_per_atom_z_lp == 1) THEN
1637 : ! Force single-rank BLAS for every atom.
1638 0 : DO P_loop_atom = 1, natom
1639 0 : n_small = n_small + 1
1640 0 : small_list(n_small) = P_loop_atom
1641 : END DO
1642 42 : ELSE IF (mem_avail_GB <= 0.0_dp) THEN
1643 : ! No /proc/meminfo => cannot size by memory.
1644 0 : IF (auto_mode) THEN
1645 0 : IF (bs_env%unit_nr > 0) THEN
1646 0 : CPWARN("RI-RS Z_lP: no meminfo; single-rank solve for all atoms")
1647 : END IF
1648 0 : DO P_loop_atom = 1, natom
1649 0 : n_small = n_small + 1
1650 0 : small_list(n_small) = P_loop_atom
1651 : END DO
1652 : ELSE
1653 : ! Fixed G, no meminfo: distribute all atoms; still floor G by the int32 limit.
1654 0 : DO P_loop_atom = 1, natom
1655 0 : ng = REAL(n_local_grid_atom(P_loop_atom), dp)
1656 0 : G_int32_max = MAX(G_int32_max, CEILING(ng*ng/scalapack_loc_limit))
1657 0 : n_big = n_big + 1
1658 0 : big_list(n_big) = P_loop_atom
1659 : END DO
1660 0 : G = MIN(bs_env%ri_rs%n_procs_per_atom_z_lp, para_env%num_pe)
1661 0 : IF (G < G_int32_max) THEN
1662 0 : G = MIN(G_int32_max, para_env%num_pe)
1663 0 : IF (bs_env%unit_nr > 0) THEN
1664 0 : CPWARN("RI-RS Z_lP: raised G to avoid ScaLAPACK overflow")
1665 : END IF
1666 : END IF
1667 : END IF
1668 : ELSE
1669 : ! Classify by memory: peak (D_local + phi_local + d_lp) vs budget. Small -> BLAS,
1670 : ! big -> distributed. Same classification for AUTO and fixed-G modes.
1671 142 : DO P_loop_atom = 1, natom
1672 100 : ng = REAL(n_local_grid_atom(P_loop_atom), dp)
1673 100 : nri = REAL(ri_blk_sizes(P_loop_atom), dp)
1674 : CALL get_n_ao_in_sphere(bs_env, particle_set, P_loop_atom, &
1675 100 : cutoff_ri_per_atom(P_loop_atom), n_ao_used_atom)
1676 100 : phi_bytes = 8.0_dp*ng*REAL(n_ao_used_atom, dp)
1677 100 : dlp_bytes = 8.0_dp*ng*nri*REAL(1 + nthreads_cls, dp)
1678 100 : peak_bytes = 8.0_dp*ng*ng + phi_bytes + dlp_bytes
1679 142 : IF (peak_bytes <= budget_bytes) THEN
1680 100 : n_small = n_small + 1
1681 100 : small_list(n_small) = P_loop_atom
1682 : ELSE
1683 0 : n_big = n_big + 1
1684 0 : big_list(n_big) = P_loop_atom
1685 : ! G must satisfy BOTH: (a) memory — distributed D_local (/G) fits next to the
1686 : ! replicated phi_local + d_lp; (b) ScaLAPACK — local ~ng^2/G below the int32 limit.
1687 0 : G_int32 = CEILING(ng*ng/scalapack_loc_limit)
1688 0 : G_int32_max = MAX(G_int32_max, G_int32)
1689 : G_atom = MAX(G_atom, G_int32, &
1690 0 : CEILING(8.0_dp*ng*ng/MAX(budget_bytes - phi_bytes - dlp_bytes, 1.0_dp)))
1691 : END IF
1692 : END DO
1693 42 : IF (n_big > 0) THEN
1694 0 : IF (auto_mode) THEN
1695 : ! Auto-size G from the most demanding big atom.
1696 0 : IF (G_atom > para_env%num_pe) THEN
1697 : CALL cp_abort(__LOCATION__, &
1698 : "RI-RS Z_lP: an atom is too large to fit even when "// &
1699 : "distributed over all ranks. Add nodes, use fewer MPI ranks "// &
1700 : "per node, lower CUTOFF_RADIUS_RL_RI, or raise EPS_FILTER "// &
1701 0 : "for more grid screening.")
1702 : END IF
1703 0 : G = MIN(MAX(G_atom, 2), para_env%num_pe)
1704 : ELSE
1705 : ! Fixed G from the keyword. Hard-floor by the ScaLAPACK int32 limit (below it
1706 : ! pdpotrf segfaults); warn if it is still below the memory recommendation.
1707 0 : G = MIN(bs_env%ri_rs%n_procs_per_atom_z_lp, para_env%num_pe)
1708 0 : IF (G < G_int32_max) THEN
1709 0 : G = MIN(G_int32_max, para_env%num_pe)
1710 0 : IF (bs_env%unit_nr > 0) THEN
1711 0 : CPWARN("RI-RS Z_lP: raised G to avoid ScaLAPACK overflow")
1712 : END IF
1713 0 : ELSE IF (G < G_atom .AND. bs_env%unit_nr > 0) THEN
1714 0 : CPWARN("RI-RS Z_lP: N_PROCS_PER_ATOM_Z_LP too small for the largest atom")
1715 : END IF
1716 : END IF
1717 : END IF
1718 : END IF
1719 :
1720 42 : CALL timestop(handle)
1721 :
1722 84 : END SUBROUTINE classify_z_lp_atoms
1723 :
1724 : ! **************************************************************************************************
1725 : !> \brief Number of AO basis functions that can be non-zero inside the RI-RS integration sphere
1726 : !> \param bs_env ...
1727 : !> \param particle_set ...
1728 : !> \param atom_P ...
1729 : !> \param cutoff_ri ...
1730 : !> \param n_ao_used ...
1731 : ! **************************************************************************************************
1732 190 : SUBROUTINE get_n_ao_in_sphere(bs_env, particle_set, atom_P, cutoff_ri, n_ao_used)
1733 :
1734 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1735 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1736 : INTEGER, INTENT(IN) :: atom_P
1737 : REAL(KIND=dp), INTENT(IN) :: cutoff_ri
1738 : INTEGER, INTENT(OUT) :: n_ao_used
1739 :
1740 : INTEGER :: source_atom
1741 :
1742 190 : n_ao_used = 0
1743 672 : DO source_atom = 1, SIZE(particle_set)
1744 1928 : IF (NORM2(particle_set(source_atom)%r(:) - particle_set(atom_P)%r(:)) > &
1745 : bs_env%ri_rs%radius_ao_per_atom(source_atom) + cutoff_ri) CYCLE
1746 : n_ao_used = n_ao_used + bs_env%i_ao_end_from_atom(source_atom) - &
1747 672 : bs_env%i_ao_start_from_atom(source_atom) + 1
1748 : END DO
1749 :
1750 190 : END SUBROUTINE get_n_ao_in_sphere
1751 :
1752 : ! **************************************************************************************************
1753 : !> \brief Builds the sphere-local AO matrix phi_local(l, μ) = Φ_μ(r_l) for one RI atom P
1754 : !> \param bs_env ...
1755 : !> \param particle_set ...
1756 : !> \param qs_kind_set ...
1757 : !> \param cell ...
1758 : !> \param ri_rs_grid_points ...
1759 : !> \param atom_P ...
1760 : !> \param cutoff_ri ...
1761 : !> \param n_ao_total ...
1762 : !> \param local_grid_idx ...
1763 : !> \param n_local_grid ...
1764 : !> \param phi_local ...
1765 : !> \param ao_col_map ...
1766 : !> \param n_ao_used ...
1767 : ! **************************************************************************************************
1768 50 : SUBROUTINE build_phi_on_sphere(bs_env, particle_set, qs_kind_set, cell, ri_rs_grid_points, &
1769 : atom_P, cutoff_ri, n_ao_total, local_grid_idx, n_local_grid, &
1770 : phi_local, ao_col_map, n_ao_used)
1771 :
1772 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1773 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1774 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1775 : TYPE(cell_type), POINTER :: cell
1776 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ri_rs_grid_points
1777 : INTEGER, INTENT(IN) :: atom_P
1778 : REAL(KIND=dp), INTENT(IN) :: cutoff_ri
1779 : INTEGER, INTENT(IN) :: n_ao_total
1780 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: local_grid_idx
1781 : INTEGER, INTENT(OUT) :: n_local_grid
1782 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1783 : INTENT(OUT) :: phi_local
1784 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_col_map
1785 : INTEGER, INTENT(OUT) :: n_ao_used
1786 :
1787 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_phi_on_sphere'
1788 :
1789 : INTEGER :: col_end, col_start, handle, j, k, l, &
1790 : loc_idx, n_grid_total, n_keep, &
1791 : source_atom
1792 : REAL(KIND=dp) :: d_sP, dist, r2_threshold
1793 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w_pt
1794 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: phi_keep, sphere_grid
1795 : REAL(KIND=dp), DIMENSION(3) :: pos_P
1796 :
1797 50 : CALL timeset(routineN, handle)
1798 :
1799 50 : n_grid_total = SIZE(ri_rs_grid_points, 2)
1800 200 : pos_P(:) = particle_set(atom_P)%r(:)
1801 :
1802 50 : n_local_grid = 0
1803 18762 : DO l = 1, n_grid_total
1804 74848 : dist = SQRT(SUM((ri_rs_grid_points(1:3, l) - pos_P(1:3))**2))
1805 18762 : IF (dist <= cutoff_ri) n_local_grid = n_local_grid + 1
1806 : END DO
1807 :
1808 150 : ALLOCATE (local_grid_idx(n_local_grid))
1809 :
1810 50 : n_local_grid = 0
1811 18762 : DO l = 1, n_grid_total
1812 74848 : dist = SQRT(SUM((ri_rs_grid_points(1:3, l) - pos_P(1:3))**2))
1813 18762 : IF (dist <= cutoff_ri) THEN
1814 18019 : n_local_grid = n_local_grid + 1
1815 18019 : local_grid_idx(n_local_grid) = l
1816 : END IF
1817 : END DO
1818 :
1819 150 : ALLOCATE (sphere_grid(3, n_local_grid))
1820 18069 : DO loc_idx = 1, n_local_grid
1821 72126 : sphere_grid(:, loc_idx) = ri_rs_grid_points(:, local_grid_idx(loc_idx))
1822 : END DO
1823 :
1824 : ! Only AOs on atoms that reach into the sphere can be non-zero here
1825 150 : ALLOCATE (ao_col_map(n_ao_total))
1826 50 : ao_col_map(:) = 0
1827 50 : n_ao_used = 0
1828 174 : DO source_atom = 1, SIZE(particle_set)
1829 496 : d_sP = NORM2(particle_set(source_atom)%r(:) - pos_P(:))
1830 124 : IF (d_sP > bs_env%ri_rs%radius_ao_per_atom(source_atom) + cutoff_ri) CYCLE
1831 :
1832 646 : DO j = bs_env%i_ao_start_from_atom(source_atom), bs_env%i_ao_end_from_atom(source_atom)
1833 472 : n_ao_used = n_ao_used + 1
1834 596 : ao_col_map(j) = n_ao_used
1835 : END DO
1836 : END DO
1837 :
1838 200 : ALLOCATE (phi_local(n_local_grid, n_ao_used))
1839 50 : phi_local = 0.0_dp
1840 :
1841 174 : DO source_atom = 1, SIZE(particle_set)
1842 496 : d_sP = NORM2(particle_set(source_atom)%r(:) - pos_P(:))
1843 124 : IF (d_sP > bs_env%ri_rs%radius_ao_per_atom(source_atom) + cutoff_ri) CYCLE
1844 :
1845 124 : col_start = ao_col_map(bs_env%i_ao_start_from_atom(source_atom))
1846 124 : col_end = ao_col_map(bs_env%i_ao_end_from_atom(source_atom))
1847 : ! A positive CUTOFF_RADIUS_RI_AO overrides the per-atom Gaussian radius
1848 : ! with a user-defined hard cutoff.
1849 124 : IF (bs_env%ri_rs%cutoff_radius_ri_ao > 0.0_dp) THEN
1850 9 : r2_threshold = bs_env%ri_rs%cutoff_radius_ri_ao**2
1851 : ELSE
1852 115 : r2_threshold = bs_env%ri_rs%radius_ao_per_atom(source_atom)**2
1853 : END IF
1854 :
1855 : CALL fill_phi_for_atom(phi_local(:, col_start:col_end), sphere_grid, &
1856 : n_local_grid, source_atom, particle_set, qs_kind_set, &
1857 174 : cell, r2_threshold)
1858 : END DO
1859 :
1860 50 : DEALLOCATE (sphere_grid)
1861 :
1862 50 : IF (n_local_grid > 0) THEN
1863 150 : ALLOCATE (w_pt(n_local_grid))
1864 : !$OMP PARALLEL DO DEFAULT(NONE) &
1865 : !$OMP SHARED(n_local_grid, n_ao_used, phi_local, w_pt) &
1866 50 : !$OMP PRIVATE(l, j) SCHEDULE(STATIC)
1867 : DO l = 1, n_local_grid
1868 : w_pt(l) = 0.0_dp
1869 : DO j = 1, n_ao_used
1870 : w_pt(l) = MAX(w_pt(l), ABS(phi_local(l, j)))
1871 : END DO
1872 : END DO
1873 : !$OMP END PARALLEL DO
1874 18069 : n_keep = COUNT(w_pt > bs_env%eps_filter)
1875 50 : IF (n_keep < n_local_grid) THEN
1876 0 : ALLOCATE (phi_keep(n_keep, n_ao_used))
1877 0 : k = 0
1878 0 : DO l = 1, n_local_grid
1879 0 : IF (w_pt(l) > bs_env%eps_filter) THEN
1880 0 : k = k + 1
1881 0 : phi_keep(k, :) = phi_local(l, :)
1882 0 : local_grid_idx(k) = local_grid_idx(l)
1883 : END IF
1884 : END DO
1885 0 : CALL MOVE_ALLOC(phi_keep, phi_local)
1886 0 : n_local_grid = n_keep
1887 : END IF
1888 50 : DEALLOCATE (w_pt)
1889 : END IF
1890 :
1891 50 : CALL timestop(handle)
1892 :
1893 100 : END SUBROUTINE build_phi_on_sphere
1894 :
1895 : ! **************************************************************************************************
1896 : !> \brief Builds the dense Jacobi-conditioned Gram matrix and the conditioning vector for the
1897 : !> single-rank (BLAS/LAPACK) Z_lP solve:
1898 : !> D_ll' = [ Σ_μ Φ_μ(r_l) Φ_μ(r_l') ]² (dsyrk of phi_local, then squared)
1899 : !> d_l = 1 / sqrt(D_ll)
1900 : !> D'_ll' = d_l D_ll' d_l' + λ δ_ll'
1901 : !> Only the lower triangle is referenced by the subsequent dpotrf('L')
1902 : !> \param phi_local ...
1903 : !> \param n_local_grid ...
1904 : !> \param n_ao_used ...
1905 : !> \param tikhonov ...
1906 : !> \param D_local ...
1907 : !> \param d_vec_local ...
1908 : ! **************************************************************************************************
1909 50 : SUBROUTINE build_gram_jacobi_blas(phi_local, n_local_grid, n_ao_used, tikhonov, D_local, &
1910 50 : d_vec_local)
1911 :
1912 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_local
1913 : INTEGER, INTENT(IN) :: n_local_grid, n_ao_used
1914 : REAL(KIND=dp), INTENT(IN) :: tikhonov
1915 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1916 : INTENT(OUT) :: D_local
1917 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d_vec_local
1918 :
1919 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gram_jacobi_blas'
1920 :
1921 : INTEGER :: handle, handle_dsyrk, i, j
1922 :
1923 50 : CALL timeset(routineN, handle)
1924 :
1925 200 : ALLOCATE (D_local(n_local_grid, n_local_grid))
1926 50 : D_local = 0.0_dp
1927 :
1928 : ! D_ll' = Σ_μ Φ_μ(r_l) Φ_μ(r_l') (lower triangle only)
1929 50 : CALL timeset(routineN//"_dsyrk", handle_dsyrk)
1930 : CALL dsyrk("L", "N", n_local_grid, n_ao_used, 1.0_dp, phi_local, &
1931 50 : n_local_grid, 0.0_dp, D_local, n_local_grid)
1932 50 : CALL timestop(handle_dsyrk)
1933 :
1934 : ! Diagonal: square, derive d_l = 1/sqrt(D_ll), scale, add Tikhonov λ
1935 : !$OMP PARALLEL DO DEFAULT(NONE) &
1936 : !$OMP SHARED(n_local_grid, D_local, d_vec_local, tikhonov) &
1937 : !$OMP PRIVATE(i) &
1938 50 : !$OMP SCHEDULE(STATIC)
1939 : DO i = 1, n_local_grid
1940 : D_local(i, i) = D_local(i, i)**2
1941 : d_vec_local(i) = 1.0_dp/SQRT(MAX(D_local(i, i), 1.0E-16_dp))
1942 : D_local(i, i) = (D_local(i, i)*d_vec_local(i)**2) + tikhonov
1943 : END DO
1944 : !$OMP END PARALLEL DO
1945 :
1946 : ! Off-diagonal: D'_ll' = d_l D_ll'^2 d_l' (mirror to the upper triangle)
1947 : !$OMP PARALLEL DO DEFAULT(NONE) &
1948 : !$OMP SHARED(n_local_grid, D_local, d_vec_local) &
1949 : !$OMP PRIVATE(j, i) &
1950 50 : !$OMP SCHEDULE(DYNAMIC)
1951 : DO j = 1, n_local_grid
1952 : DO i = j + 1, n_local_grid
1953 : D_local(i, j) = D_local(i, j)**2
1954 : D_local(i, j) = D_local(i, j)*d_vec_local(i)*d_vec_local(j)
1955 : D_local(j, i) = D_local(i, j)
1956 : END DO
1957 : END DO
1958 : !$OMP END PARALLEL DO
1959 :
1960 50 : CALL timestop(handle)
1961 :
1962 100 : END SUBROUTINE build_gram_jacobi_blas
1963 :
1964 : ! **************************************************************************************************
1965 : !> \brief Computes the Jacobi conditioning vector directly from phi for the distributed
1966 : !> (ScaLAPACK) Z_lP solve: d_l = 1 / Σ_μ Φ_μ(r_l)² = 1/sqrt(D_ll), without forming
1967 : !> the Gram matrix.
1968 : !> \param phi_local ...
1969 : !> \param n_local_grid ...
1970 : !> \param n_ao_used ...
1971 : !> \param d_vec_local ...
1972 : ! **************************************************************************************************
1973 0 : SUBROUTINE build_jacobi_diag_from_phi(phi_local, n_local_grid, n_ao_used, d_vec_local)
1974 :
1975 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_local
1976 : INTEGER, INTENT(IN) :: n_local_grid, n_ao_used
1977 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d_vec_local
1978 :
1979 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_jacobi_diag_from_phi'
1980 :
1981 : INTEGER :: handle, i, j
1982 :
1983 0 : CALL timeset(routineN, handle)
1984 :
1985 : !$OMP PARALLEL DO DEFAULT(NONE) &
1986 : !$OMP SHARED(n_local_grid, n_ao_used, phi_local, d_vec_local) &
1987 : !$OMP PRIVATE(i, j) &
1988 0 : !$OMP SCHEDULE(STATIC)
1989 : DO i = 1, n_local_grid
1990 : d_vec_local(i) = 0.0_dp
1991 : DO j = 1, n_ao_used
1992 : d_vec_local(i) = d_vec_local(i) + phi_local(i, j)*phi_local(i, j)
1993 : END DO
1994 : d_vec_local(i) = 1.0_dp/MAX(d_vec_local(i), 1.0E-16_dp)
1995 : END DO
1996 : !$OMP END PARALLEL DO
1997 :
1998 0 : CALL timestop(handle)
1999 :
2000 0 : END SUBROUTINE build_jacobi_diag_from_phi
2001 :
2002 : ! **************************************************************************************************
2003 : !> \brief Scales every row of a matrix by the corresponding diagonal entry,
2004 : !> A(l, :) <- d_l * A(l, :). Used in Z_lP solve: to pre-scale the RHS
2005 : !> (D'_lP = d_l D_lP) and to undo the conditioning of the solution (Z_lP = d_l Z'_lP).
2006 : !> \param d_lp_local ...
2007 : !> \param d_vec_local ...
2008 : !> \param n_local_grid ...
2009 : !> \param n_loc_ri ...
2010 : ! **************************************************************************************************
2011 100 : SUBROUTINE scale_rows_by_diag(d_lp_local, d_vec_local, n_local_grid, n_loc_ri)
2012 :
2013 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: d_lp_local
2014 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: d_vec_local
2015 : INTEGER, INTENT(IN) :: n_local_grid, n_loc_ri
2016 :
2017 : CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_rows_by_diag'
2018 :
2019 : INTEGER :: handle, i, j_ri
2020 :
2021 100 : CALL timeset(routineN, handle)
2022 :
2023 : !$OMP PARALLEL DO DEFAULT(NONE) &
2024 : !$OMP SHARED(n_loc_ri, n_local_grid, d_lp_local, d_vec_local) &
2025 : !$OMP PRIVATE(j_ri, i) &
2026 100 : !$OMP SCHEDULE(STATIC)
2027 : DO j_ri = 1, n_loc_ri
2028 : DO i = 1, n_local_grid
2029 : d_lp_local(i, j_ri) = d_lp_local(i, j_ri)*d_vec_local(i)
2030 : END DO
2031 : END DO
2032 : !$OMP END PARALLEL DO
2033 :
2034 100 : CALL timestop(handle)
2035 :
2036 100 : END SUBROUTINE scale_rows_by_diag
2037 :
2038 : ! **************************************************************************************************
2039 : !> \brief Scatters the solved Z columns of one atom P from the dense sphere-local solution back
2040 : !> into the global sparse mat_Z_lP
2041 : !> \param mat_Z_lP ...
2042 : !> \param d_lp_local ...
2043 : !> \param local_grid_idx ...
2044 : !> \param n_local_grid ...
2045 : !> \param n_loc_ri ...
2046 : !> \param atom_P ...
2047 : !> \param r_blk_sizes ...
2048 : !> \param row_offset ...
2049 : !> \param eps_filter ...
2050 : ! **************************************************************************************************
2051 50 : SUBROUTINE scatter_z_columns(mat_Z_lP, d_lp_local, local_grid_idx, n_local_grid, n_loc_ri, &
2052 50 : atom_P, r_blk_sizes, row_offset, eps_filter)
2053 :
2054 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_Z_lP
2055 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: d_lp_local
2056 : INTEGER, DIMENSION(:), INTENT(IN) :: local_grid_idx
2057 : INTEGER, INTENT(IN) :: n_local_grid, n_loc_ri, atom_P
2058 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes, row_offset
2059 : REAL(KIND=dp), INTENT(IN) :: eps_filter
2060 :
2061 : CHARACTER(LEN=*), PARAMETER :: routineN = 'scatter_z_columns'
2062 :
2063 : INTEGER :: current_chunk_size, g_pt, handle, i_blk, &
2064 : loc_ptr, r_end, r_start
2065 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Z_blk
2066 :
2067 50 : CALL timeset(routineN, handle)
2068 :
2069 870 : ALLOCATE (Z_blk(MAXVAL(r_blk_sizes), n_loc_ri))
2070 50 : loc_ptr = 1
2071 :
2072 720 : DO i_blk = 1, SIZE(r_blk_sizes)
2073 670 : r_start = row_offset(i_blk) + 1
2074 670 : r_end = row_offset(i_blk) + r_blk_sizes(i_blk)
2075 670 : current_chunk_size = r_blk_sizes(i_blk)
2076 :
2077 670 : Z_blk = 0.0_dp
2078 :
2079 : ! Copy the sphere points whose global grid index falls inside this block
2080 18689 : DO WHILE (loc_ptr <= n_local_grid)
2081 18635 : g_pt = local_grid_idx(loc_ptr)
2082 18635 : IF (g_pt > r_end) EXIT
2083 254814 : Z_blk(g_pt - r_start + 1, 1:n_loc_ri) = d_lp_local(loc_ptr, 1:n_loc_ri)
2084 18635 : loc_ptr = loc_ptr + 1
2085 : END DO
2086 :
2087 250858 : IF (MAXVAL(ABS(Z_blk(1:current_chunk_size, 1:n_loc_ri))) > eps_filter) THEN
2088 : CALL dbcsr_put_block(mat_Z_lP, row=i_blk, col=atom_P, &
2089 658 : block=Z_blk(1:current_chunk_size, 1:n_loc_ri))
2090 : END IF
2091 : END DO
2092 :
2093 50 : DEALLOCATE (Z_blk)
2094 :
2095 50 : CALL timestop(handle)
2096 :
2097 50 : END SUBROUTINE scatter_z_columns
2098 :
2099 : ! **************************************************************************************************
2100 : !> \brief LPT (longest-processing-time) assignment of the Z_lP atoms to workers (MPI ranks in
2101 : !> Phase A, rank subgroups in Phase B): sort by estimated solve cost n_local_grid^3
2102 : !> (the per-atom Cholesky dominates; the n^2 assembly terms order the atoms the same
2103 : !> way) and greedily give each atom to the least-loaded worker.
2104 : !> \param atom_list ...
2105 : !> \param n_atoms ...
2106 : !> \param n_local_grid_atom ...
2107 : !> \param n_workers ...
2108 : !> \param my_worker ...
2109 : !> \param my_atoms ...
2110 : !> \param max_over_mean ...
2111 : ! **************************************************************************************************
2112 42 : SUBROUTINE lpt_assign_atoms(atom_list, n_atoms, n_local_grid_atom, n_workers, my_worker, &
2113 : my_atoms, max_over_mean)
2114 :
2115 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
2116 : INTEGER, INTENT(IN) :: n_atoms
2117 : INTEGER, DIMENSION(:), INTENT(IN) :: n_local_grid_atom
2118 : INTEGER, INTENT(IN) :: n_workers, my_worker
2119 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: my_atoms
2120 : REAL(KIND=dp), INTENT(OUT) :: max_over_mean
2121 :
2122 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lpt_assign_atoms'
2123 :
2124 : INTEGER :: handle, i, iw, n_mine, w_min
2125 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: mine_tmp, perm
2126 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cost, load
2127 :
2128 42 : CALL timeset(routineN, handle)
2129 :
2130 42 : max_over_mean = 1.0_dp
2131 42 : IF (n_atoms <= 0) THEN
2132 0 : ALLOCATE (my_atoms(0))
2133 0 : CALL timestop(handle)
2134 0 : RETURN
2135 : END IF
2136 :
2137 336 : ALLOCATE (cost(n_atoms), perm(n_atoms), mine_tmp(n_atoms), load(n_workers))
2138 142 : DO i = 1, n_atoms
2139 142 : cost(i) = REAL(n_local_grid_atom(atom_list(i)), dp)**3
2140 : END DO
2141 42 : CALL sort(cost, n_atoms, perm) ! ascending; walk backwards for largest-first
2142 :
2143 42 : load(:) = 0.0_dp
2144 42 : n_mine = 0
2145 142 : DO i = n_atoms, 1, -1
2146 : w_min = 1
2147 200 : DO iw = 2, n_workers
2148 200 : IF (load(iw) < load(w_min)) w_min = iw
2149 : END DO
2150 100 : load(w_min) = load(w_min) + cost(i)
2151 142 : IF (w_min - 1 == my_worker) THEN
2152 50 : n_mine = n_mine + 1
2153 50 : mine_tmp(n_mine) = atom_list(perm(i))
2154 : END IF
2155 : END DO
2156 :
2157 126 : ALLOCATE (my_atoms(n_mine))
2158 92 : my_atoms(:) = mine_tmp(1:n_mine)
2159 294 : IF (SUM(load) > 0.0_dp) max_over_mean = MAXVAL(load)*REAL(n_workers, dp)/SUM(load)
2160 :
2161 42 : CALL timestop(handle)
2162 :
2163 84 : END SUBROUTINE lpt_assign_atoms
2164 :
2165 : ! **************************************************************************************************
2166 : !> \brief Computes the dense localized RHS d_lp(l,P) = Σ_{μν} Φ_μ(r_l)·Φ_ν(r_l)·(μν|P) for one
2167 : !> RI atom P, OMP-threaded over (atom_j, atom_k) AO-pair blocks: per thread, build the 3c
2168 : !> block, then contract grid-chunked pair densities into a private d_lp partial; partials
2169 : !> are reduced into d_lp at the end.
2170 : !> Pair screening is handled inside build_3c_integral_block_ctx via the `screened` output.
2171 : !> \param bs_env ...
2172 : !> \param ctx ...
2173 : !> \param phi_val ...
2174 : !> \param ao_col_map ...
2175 : !> \param d_lp ...
2176 : !> \param n_grid_total ...
2177 : !> \param n_loc_ri ...
2178 : !> \param atom_P ...
2179 : !> \param max_ao_size ...
2180 : !> \param atom_j_mepos ...
2181 : !> \param atom_j_stride ...
2182 : ! **************************************************************************************************
2183 50 : SUBROUTINE compute_d_lp(bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid_total, n_loc_ri, atom_P, &
2184 : max_ao_size, atom_j_mepos, atom_j_stride)
2185 :
2186 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2187 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
2188 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val
2189 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
2190 : INTEGER, INTENT(IN) :: n_grid_total, n_loc_ri
2191 : REAL(KIND=dp), INTENT(INOUT) :: d_lp(n_grid_total, n_loc_ri)
2192 : INTEGER, INTENT(IN) :: atom_P, max_ao_size, atom_j_mepos, &
2193 : atom_j_stride
2194 :
2195 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp'
2196 : INTEGER, PARAMETER :: grid_chunk = 1024
2197 :
2198 : INTEGER :: atom_j, atom_k, c, handle, handle_dgemm, &
2199 : j, jk_idx, jsize, jstart, k, ksize, &
2200 : kstart, l, l0, n_grid_pair, point, ri
2201 50 : INTEGER, ALLOCATABLE :: grid_index(:)
2202 : LOGICAL :: screened
2203 50 : LOGICAL, ALLOCATABLE :: skip_grid_point(:, :)
2204 : REAL(KIND=dp) :: pair_factor
2205 50 : REAL(KIND=dp), ALLOCATABLE :: grid_result(:, :)
2206 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_lp_prv, int_2d_prv, rho_chunk
2207 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c_prv
2208 50 : TYPE(gw_3c_ws_type) :: ws
2209 :
2210 50 : CALL timeset(routineN, handle)
2211 :
2212 : !$OMP PARALLEL DEFAULT(NONE) &
2213 : !$OMP SHARED(skip_grid_point, bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid_total, n_loc_ri, atom_P, max_ao_size, &
2214 : !$OMP atom_j_mepos, atom_j_stride) &
2215 : !$OMP PRIVATE(grid_index, grid_result, n_grid_pair, point, pair_factor, &
2216 : !$OMP atom_j, atom_k, c, handle_dgemm, j, jk_idx, jsize, jstart, k, ksize, kstart, &
2217 50 : !$OMP l, l0, ri, screened, d_lp_prv, int_2d_prv, rho_chunk, int_3c_prv, ws)
2218 :
2219 : CALL gw_3c_ws_create(ws, ctx)
2220 : ALLOCATE (int_3c_prv(max_ao_size, max_ao_size, n_loc_ri))
2221 : ALLOCATE (int_2d_prv(max_ao_size*max_ao_size, n_loc_ri))
2222 : ALLOCATE (rho_chunk(grid_chunk, max_ao_size*max_ao_size))
2223 : ALLOCATE (d_lp_prv(n_grid_total, n_loc_ri))
2224 : ALLOCATE (grid_index(n_grid_total), grid_result(grid_chunk, n_loc_ri))
2225 : d_lp_prv(:, :) = 0.0_dp
2226 :
2227 : !$OMP SINGLE
2228 : ALLOCATE (skip_grid_point(n_grid_total, bs_env%n_atom))
2229 : CALL compute_skip_grid_point(bs_env, phi_val, ao_col_map, skip_grid_point)
2230 : !$OMP END SINGLE
2231 :
2232 : ! MPI assigns each unordered pair through its first atom; OpenMP divides those atoms.
2233 : ! Skip only exact-zero pair-grid support, without a new screening threshold.
2234 : !$OMP DO SCHEDULE(DYNAMIC)
2235 : DO atom_j = atom_j_mepos + 1, bs_env%n_atom, atom_j_stride
2236 : DO atom_k = atom_j, bs_env%n_atom
2237 : jstart = ao_col_map(bs_env%i_ao_start_from_atom(atom_j))
2238 : kstart = ao_col_map(bs_env%i_ao_start_from_atom(atom_k))
2239 : IF (jstart == 0 .OR. kstart == 0) CYCLE
2240 : jsize = bs_env%i_ao_end_from_atom(atom_j) - bs_env%i_ao_start_from_atom(atom_j) + 1
2241 : ksize = bs_env%i_ao_end_from_atom(atom_k) - bs_env%i_ao_start_from_atom(atom_k) + 1
2242 :
2243 : n_grid_pair = 0
2244 : DO point = 1, n_grid_total
2245 : IF (skip_grid_point(point, atom_j) .OR. skip_grid_point(point, atom_k)) CYCLE
2246 : n_grid_pair = n_grid_pair + 1
2247 : grid_index(n_grid_pair) = point
2248 : END DO
2249 : IF (n_grid_pair == 0) CYCLE
2250 : ! (μν|P) = (νμ|P): distinct atom pairs contribute twice.
2251 : pair_factor = 1.0_dp
2252 : IF (atom_j /= atom_k) pair_factor = 2.0_dp
2253 :
2254 : int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp
2255 :
2256 : ! Compute B_{μν,P} = (μν|P); ctx-internal triangle-inequality screening on
2257 : ! kind_radius sets `screened=.TRUE.` for negligible triples.
2258 : CALL build_3c_integral_block_ctx(int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri), &
2259 : ctx, ws, atom_j=atom_j, atom_k=atom_k, atom_i=atom_P, &
2260 : screened=screened)
2261 :
2262 : IF (screened) CYCLE
2263 :
2264 : ! Flatten 3D B_{μν, P} tensor to 2D B_{(μν), P} matrix for BLAS
2265 : DO ri = 1, n_loc_ri
2266 : DO k = 1, ksize
2267 : DO j = 1, jsize
2268 : jk_idx = (k - 1)*jsize + j
2269 : int_2d_prv(jk_idx, ri) = int_3c_prv(j, k, ri)
2270 : END DO
2271 : END DO
2272 : END DO
2273 :
2274 : ! Pair density ρ(l, μν) = Φ_μ(r_l) Φ_ν(r_l) in grid chunks, contracted on the fly:
2275 : ! d_{l,P} += ρ(l, μν) B_{(μν),P} (dgemm runs serially inside the parallel region)
2276 : DO l0 = 1, n_grid_pair, grid_chunk
2277 : c = MIN(grid_chunk, n_grid_pair - l0 + 1)
2278 : DO k = 1, ksize
2279 : DO j = 1, jsize
2280 : jk_idx = (k - 1)*jsize + j
2281 : DO l = 1, c
2282 : point = grid_index(l0 + l - 1)
2283 : rho_chunk(l, jk_idx) = phi_val(point, jstart + j - 1)* &
2284 : phi_val(point, kstart + k - 1)
2285 : END DO
2286 : END DO
2287 : END DO
2288 : CALL timeset(routineN//"_dgemm", handle_dgemm)
2289 : CALL dgemm("N", "N", c, n_loc_ri, jsize*ksize, &
2290 : pair_factor, rho_chunk, grid_chunk, &
2291 : int_2d_prv, max_ao_size*max_ao_size, &
2292 : 0.0_dp, grid_result, grid_chunk)
2293 : DO ri = 1, n_loc_ri
2294 : DO l = 1, c
2295 : point = grid_index(l0 + l - 1)
2296 : d_lp_prv(point, ri) = d_lp_prv(point, ri) + grid_result(l, ri)
2297 : END DO
2298 : END DO
2299 : CALL timestop(handle_dgemm)
2300 : END DO
2301 : END DO
2302 : END DO
2303 : !$OMP END DO
2304 :
2305 : !$OMP CRITICAL (compute_d_lp_reduce)
2306 : d_lp(1:n_grid_total, 1:n_loc_ri) = d_lp(1:n_grid_total, 1:n_loc_ri) + &
2307 : d_lp_prv(1:n_grid_total, 1:n_loc_ri)
2308 : !$OMP END CRITICAL (compute_d_lp_reduce)
2309 :
2310 : DEALLOCATE (int_3c_prv, int_2d_prv, rho_chunk, d_lp_prv, grid_index, grid_result)
2311 : CALL gw_3c_ws_release(ws)
2312 :
2313 : !$OMP SINGLE
2314 : DEALLOCATE (skip_grid_point)
2315 : !$OMP END SINGLE
2316 :
2317 : !$OMP END PARALLEL
2318 :
2319 50 : CALL timestop(handle)
2320 :
2321 50 : END SUBROUTINE compute_d_lp
2322 :
2323 : ! **************************************************************************************************
2324 : !> \brief Marks grid points where all stored AO values of an atom are exactly zero.
2325 : !> \param bs_env ...
2326 : !> \param phi_val ...
2327 : !> \param ao_col_map ...
2328 : !> \param skip_grid_point ...
2329 : ! **************************************************************************************************
2330 50 : SUBROUTINE compute_skip_grid_point(bs_env, phi_val, ao_col_map, skip_grid_point)
2331 :
2332 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2333 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val
2334 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
2335 : LOGICAL, DIMENSION(:, :), INTENT(OUT) :: skip_grid_point
2336 :
2337 : INTEGER :: atom, first_ao, number_of_aos
2338 :
2339 45311 : skip_grid_point(:, :) = .TRUE.
2340 174 : DO atom = 1, bs_env%n_atom
2341 124 : first_ao = ao_col_map(bs_env%i_ao_start_from_atom(atom))
2342 124 : IF (first_ao == 0) CYCLE
2343 124 : number_of_aos = bs_env%i_ao_end_from_atom(atom) - bs_env%i_ao_start_from_atom(atom) + 1
2344 166003 : skip_grid_point(:, atom) = ALL(phi_val(:, first_ao:first_ao + number_of_aos - 1) == 0.0_dp, DIM=2)
2345 : END DO
2346 :
2347 50 : END SUBROUTINE compute_skip_grid_point
2348 :
2349 : ! **************************************************************************************************
2350 : !> \brief Distributed pdpotrf/pdpotrs solve of D x = b for one atom of the
2351 : !> RI-RS Z_lP build (Phase B, "big" atoms), called with a subgroup of
2352 : !> cooperating ranks and an associated BLACS context.
2353 : !> Each rank in the subgroup holds the (replicated) phi_local and the
2354 : !> (replicated) RHS d_lp; it builds its own block-cyclic slice of the
2355 : !> squared+Jacobi-scaled Gram matrix D via tiled DGEMM, factorizes via
2356 : !> cp_fm_cholesky_decompose (UPLO='U'), and solves with
2357 : !> cp_fm_cholesky_solve. The replicated d_lp is updated in place via
2358 : !> cp_fm_get_submatrix.
2359 : !> \param phi_local ...
2360 : !> \param d_vec ...
2361 : !> \param d_lp ...
2362 : !> \param n_loc ...
2363 : !> \param n_ao ...
2364 : !> \param n_rhs ...
2365 : !> \param tikhonov ...
2366 : !> \param para_env_sub ...
2367 : !> \param blacs_env_sub ...
2368 : !> \param fm_struct_D ...
2369 : !> \param fm_struct_b ...
2370 : !> \param fm_D ...
2371 : !> \param fm_b ...
2372 : !> \param info ...
2373 : ! **************************************************************************************************
2374 0 : SUBROUTINE solve_D_lp_distributed(phi_local, d_vec, d_lp, n_loc, n_ao, n_rhs, &
2375 : tikhonov, para_env_sub, blacs_env_sub, &
2376 : fm_struct_D, fm_struct_b, fm_D, fm_b, info)
2377 :
2378 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_local
2379 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: d_vec
2380 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: d_lp
2381 : INTEGER, INTENT(IN) :: n_loc, n_ao, n_rhs
2382 : REAL(KIND=dp), INTENT(IN) :: tikhonov
2383 : TYPE(mp_para_env_type), POINTER :: para_env_sub
2384 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
2385 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_D, fm_struct_b
2386 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_D, fm_b
2387 : INTEGER, INTENT(OUT) :: info
2388 :
2389 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_D_lp_distributed'
2390 :
2391 : INTEGER :: handle, i_loc, ig, j_loc, jg, &
2392 : ncol_local, nrow_local
2393 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2394 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2395 0 : POINTER :: local_data
2396 :
2397 0 : CALL timeset(routineN, handle)
2398 0 : info = 0
2399 :
2400 0 : NULLIFY (fm_struct_D, fm_struct_b)
2401 : CALL cp_fm_struct_create(fm_struct_D, para_env=para_env_sub, &
2402 : context=blacs_env_sub, &
2403 0 : nrow_global=n_loc, ncol_global=n_loc)
2404 : CALL cp_fm_struct_create(fm_struct_b, para_env=para_env_sub, &
2405 : context=blacs_env_sub, &
2406 0 : nrow_global=n_loc, ncol_global=n_rhs)
2407 0 : CALL cp_fm_create(fm_D, fm_struct_D)
2408 0 : CALL cp_fm_create(fm_b, fm_struct_b)
2409 :
2410 : ! ---- Build the local block-cyclic slice of fm_D --------------------
2411 : ! Tiled DGEMM build: for each (row_tile x col_tile) sub-block we
2412 : ! gather small phi strips, DGEMM into a small gram tile, then write
2413 : ! the squared & d_vec-scaled result directly into fm_D%local_data.
2414 : ! This avoids materialising the full (nrow_local x n_ao),
2415 : ! (n_ao x ncol_local), and (nrow_local x ncol_local) buffers.
2416 : CALL cp_fm_get_info(fm_D, nrow_local=nrow_local, ncol_local=ncol_local, &
2417 : row_indices=row_indices, col_indices=col_indices, &
2418 0 : local_data=local_data)
2419 :
2420 0 : IF (nrow_local > 0 .AND. ncol_local > 0) THEN
2421 0 : BLOCK
2422 : INTEGER, PARAMETER :: ntile = 1024
2423 : INTEGER :: ib, ie, jb, je, mb, kb, ti, tj, handle_dgemm
2424 0 : REAL(KIND=dp), ALLOCATABLE :: gram_t(:, :), phi_cols_t(:, :), phi_rows_t(:, :)
2425 0 : ALLOCATE (phi_rows_t(ntile, n_ao), phi_cols_t(n_ao, ntile), gram_t(ntile, ntile))
2426 0 : DO ib = 1, nrow_local, ntile
2427 0 : ie = MIN(ib + ntile - 1, nrow_local)
2428 0 : mb = ie - ib + 1
2429 : !$OMP PARALLEL DO DEFAULT(NONE) &
2430 : !$OMP SHARED(mb, n_ao, phi_rows_t, phi_local, row_indices, ib) &
2431 0 : !$OMP PRIVATE(ti, j_loc) SCHEDULE(STATIC)
2432 : DO j_loc = 1, n_ao
2433 : DO ti = 1, mb
2434 : phi_rows_t(ti, j_loc) = phi_local(row_indices(ib + ti - 1), j_loc)
2435 : END DO
2436 : END DO
2437 : !$OMP END PARALLEL DO
2438 0 : DO jb = 1, ncol_local, ntile
2439 0 : je = MIN(jb + ntile - 1, ncol_local)
2440 0 : kb = je - jb + 1
2441 : !$OMP PARALLEL DO DEFAULT(NONE) &
2442 : !$OMP SHARED(kb, n_ao, phi_cols_t, phi_local, col_indices, jb) &
2443 0 : !$OMP PRIVATE(tj, i_loc) SCHEDULE(STATIC)
2444 : DO tj = 1, kb
2445 : DO i_loc = 1, n_ao
2446 : phi_cols_t(i_loc, tj) = phi_local(col_indices(jb + tj - 1), i_loc)
2447 : END DO
2448 : END DO
2449 : !$OMP END PARALLEL DO
2450 0 : CALL timeset(routineN//"_dgemm", handle_dgemm)
2451 : CALL dgemm('N', 'N', mb, kb, n_ao, &
2452 : 1.0_dp, phi_rows_t, ntile, phi_cols_t, n_ao, &
2453 0 : 0.0_dp, gram_t, ntile)
2454 0 : CALL timestop(handle_dgemm)
2455 : !$OMP PARALLEL DO DEFAULT(NONE) &
2456 : !$OMP SHARED(mb, kb, gram_t, d_vec, row_indices, col_indices, ib, jb) &
2457 : !$OMP SHARED(local_data, tikhonov) &
2458 0 : !$OMP PRIVATE(ti, tj, ig, jg) SCHEDULE(STATIC)
2459 : DO tj = 1, kb
2460 : jg = col_indices(jb + tj - 1)
2461 : DO ti = 1, mb
2462 : ig = row_indices(ib + ti - 1)
2463 : local_data(ib + ti - 1, jb + tj - 1) = &
2464 : gram_t(ti, tj)*gram_t(ti, tj)*d_vec(ig)*d_vec(jg)
2465 : IF (ig == jg) THEN
2466 : local_data(ib + ti - 1, jb + tj - 1) = &
2467 : local_data(ib + ti - 1, jb + tj - 1) + tikhonov
2468 : END IF
2469 : END DO
2470 : END DO
2471 : !$OMP END PARALLEL DO
2472 : END DO
2473 : END DO
2474 0 : DEALLOCATE (phi_rows_t, phi_cols_t, gram_t)
2475 : END BLOCK
2476 : END IF
2477 :
2478 : ! Load the replicated d_lp into the block-cyclic fm_b
2479 0 : CALL cp_fm_set_submatrix(fm_b, d_lp)
2480 :
2481 : ! pdpotrf (Cholesky factorisation; cp_fm_cholesky_decompose
2482 : ! factors with UPLO='U', so pdpotrs must match)
2483 0 : CALL cp_fm_cholesky_decompose(fm_D, n=n_loc, info_out=info)
2484 0 : IF (info /= 0) THEN
2485 0 : CPABORT("pdpotrf failed in solve_D_lp_distributed")
2486 : END IF
2487 :
2488 : ! pdpotrs/dpotrs (solve in place on fm_b)
2489 0 : CALL cp_fm_cholesky_solve(fm_D, fm_b, n=n_loc, info_out=info)
2490 0 : IF (info /= 0) THEN
2491 0 : CPABORT("pdpotrs failed in solve_D_lp_distributed")
2492 : END IF
2493 :
2494 : ! Gather distributed solution back into the replicated d_lp
2495 0 : CALL cp_fm_get_submatrix(fm_b, d_lp)
2496 :
2497 0 : CALL cp_fm_release(fm_D)
2498 0 : CALL cp_fm_release(fm_b)
2499 0 : CALL cp_fm_struct_release(fm_struct_D)
2500 0 : CALL cp_fm_struct_release(fm_struct_b)
2501 :
2502 0 : CALL timestop(handle)
2503 :
2504 0 : END SUBROUTINE solve_D_lp_distributed
2505 :
2506 : ! **************************************************************************************************
2507 : !> \brief Computes the polarizability matrix in the RI basis for every
2508 : !> imaginary-time point:
2509 : !> G^occ/vir_μν(i|τ|) = Σ_n C_μn e^(-|(ϵ_n-ϵ_F)τ|) C_νn (AO x AO, build_G_ao)
2510 : !> χ_ll'(iτ) = [Σ_μν Φ_μ(r_l) G^occ_μν Φ_ν(r_l')] ∘ [Σ_μν Φ_μ(r_l) G^vir_μν Φ_ν(r_l')]
2511 : !> χ_PQ(iτ) = g_s Σ_ll' Z_lP χ_ll'(iτ) Z_l'Q (g_s = spin degeneracy)
2512 : !> The grid (l) index is streamed in panels (contract_grid_panels).
2513 : !> \param bs_env ...
2514 : !> \param mat_chi_Gamma_tau ...
2515 : !> \param mat_phi_mu_l ...
2516 : !> \param mat_Z_lP ...
2517 : ! **************************************************************************************************
2518 36 : SUBROUTINE get_mat_chi_Gamma_tau(bs_env, mat_chi_Gamma_tau, mat_phi_mu_l, mat_Z_lP)
2519 :
2520 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2521 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
2522 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
2523 :
2524 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
2525 :
2526 : INTEGER :: handle, i_t, ispin, n_panels
2527 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
2528 : REAL(KIND=dp) :: grid_occ, t1, tau
2529 : TYPE(dbcsr_type) :: matrix_G_occ_ao, matrix_G_vir_ao
2530 :
2531 36 : CALL timeset(routineN, handle)
2532 :
2533 : ! Panel boundaries for the grid-streaming contraction.
2534 : ! The panels are identical for χ, Σ^x and Σ^c, so the count is reported once here for all three stages.
2535 :
2536 36 : CALL resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
2537 :
2538 36 : n_panels = SIZE(pan_first)
2539 36 : IF (bs_env%unit_nr > 0) THEN
2540 18 : WRITE (bs_env%unit_nr, '(T2,A,I44)') 'Number of batches for χ, Σ matrices', n_panels
2541 18 : CALL m_flush(bs_env%unit_nr)
2542 : END IF
2543 :
2544 : ! =========================================================================
2545 : ! IMAGINARY TIME LOOP
2546 : ! χ_PQ(iτ) = Σ_s g_s · Z^T ( (φ G^occ_s φ^T) ∘ (φ G^vir_s φ^T) ) Z
2547 : ! (g_s = spin degeneracy)
2548 : ! =========================================================================
2549 546 : DO i_t = 1, bs_env%num_time_freq_points
2550 510 : t1 = m_walltime()
2551 510 : tau = bs_env%imag_time_points(i_t)
2552 :
2553 1100 : DO ispin = 1, bs_env%n_spin
2554 :
2555 : ! AO-space Green's functions G^occ_µν, G^vir_µν (dense AO x AO, small)
2556 590 : CALL build_G_ao(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_G_occ_ao)
2557 590 : CALL build_G_ao(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, matrix_G_vir_ao)
2558 :
2559 : ! χ_PQ += g_s · Z^T ( (φ G^occ φ^T) ∘ (φ G^vir φ^T) ) Z
2560 : CALL contract_grid_panels(L_A=mat_phi_mu_l, M_A=matrix_G_occ_ao, &
2561 : L_B=mat_phi_mu_l, M_B=matrix_G_vir_ao, &
2562 : L_out=mat_Z_lP, mat_out=mat_chi_Gamma_tau(i_t)%matrix, &
2563 : scale=bs_env%spin_degeneracy, eps=bs_env%eps_filter, &
2564 : para_env=bs_env%para_env, &
2565 : pan_first=pan_first, pan_last=pan_last, &
2566 : lb_eq_la=.TRUE., lout_eq_la=.FALSE., &
2567 : zero_out=(ispin == 1), &
2568 : keep_sparsity=bs_env%ri_rs%keep_sparsity_rirs, &
2569 : centroids=bs_env%ri_rs%chunk_centroids, &
2570 : cutoff=bs_env%ri_rs%cutoff_radius_v_w, &
2571 590 : grid_occupation=grid_occ)
2572 :
2573 590 : CALL dbcsr_release(matrix_G_occ_ao)
2574 1100 : CALL dbcsr_release(matrix_G_vir_ao)
2575 :
2576 : END DO ! ispin
2577 :
2578 : ! Sparsity reports
2579 510 : IF (i_t == 1 .AND. bs_env%unit_nr > 0) THEN
2580 18 : WRITE (bs_env%unit_nr, '(A)') ' '
2581 : WRITE (bs_env%unit_nr, '(T2,A,F17.2,A)') &
2582 18 : 'Occupation of non-zero elements of G(l,l''), χ(l,l''), W(l,l'')', &
2583 36 : grid_occ*100.0_dp, ' %'
2584 18 : WRITE (bs_env%unit_nr, '(A)') ' '
2585 18 : CALL m_flush(bs_env%unit_nr)
2586 : END IF
2587 510 : IF (i_t == 1) THEN
2588 36 : CALL print_matrix_occupation(mat_chi_Gamma_tau(i_t)%matrix, 'χ(P,Q)', bs_env)
2589 : END IF
2590 :
2591 546 : IF (bs_env%unit_nr > 0) THEN
2592 : WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
2593 255 : 'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
2594 510 : ', Execution time', m_walltime() - t1, ' s'
2595 : END IF
2596 :
2597 : END DO ! i_t
2598 :
2599 36 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
2600 :
2601 36 : CALL timestop(handle)
2602 :
2603 72 : END SUBROUTINE get_mat_chi_Gamma_tau
2604 :
2605 : ! **************************************************************************************************
2606 : !> \brief Prints the non-zero occupation percentage of a DBCSR matrix on one line.
2607 : !> \param matrix ...
2608 : !> \param label ...
2609 : !> \param bs_env ...
2610 : ! **************************************************************************************************
2611 122 : SUBROUTINE print_matrix_occupation(matrix, label, bs_env)
2612 :
2613 : TYPE(dbcsr_type), INTENT(IN) :: matrix
2614 : CHARACTER(LEN=*), INTENT(IN) :: label
2615 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2616 :
2617 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_matrix_occupation'
2618 :
2619 : CHARACTER(LEN=max_line_length) :: msg
2620 : INTEGER :: handle
2621 : REAL(KIND=dp) :: frac_2p31, max_loc, occ
2622 :
2623 122 : CALL timeset(routineN, handle)
2624 :
2625 122 : occ = dbcsr_get_occupation(matrix)
2626 122 : max_loc = REAL(dbcsr_get_data_size(matrix), dp)
2627 122 : CALL bs_env%para_env%max(max_loc)
2628 :
2629 122 : IF (bs_env%unit_nr > 0) THEN
2630 61 : frac_2p31 = max_loc/REAL(bs_env%dbcsr_msg_elem_limit, dp)
2631 61 : WRITE (bs_env%unit_nr, '(A)') ' '
2632 : WRITE (bs_env%unit_nr, '(T2,A,F36.2,A)') &
2633 61 : 'Occupation of non-zero elements of '//TRIM(label), occ*100.0_dp, ' %'
2634 61 : IF (frac_2p31 > 0.5_dp) THEN
2635 : WRITE (msg, '(3A,F0.2,A)') &
2636 0 : "The largest per-rank message of ", TRIM(label), " reaches ", frac_2p31, &
2637 : " of the 32-bit limit that DBCSR uses for its message length. Beyond it the "// &
2638 : "length overflows and multiply_cannon fails. Reduce the per-rank block size, "// &
2639 0 : "for instance with more MPI ranks or a larger N_PANELS."
2640 0 : CPWARN(TRIM(msg))
2641 : END IF
2642 61 : CALL m_flush(bs_env%unit_nr)
2643 : END IF
2644 :
2645 122 : CALL timestop(handle)
2646 :
2647 122 : END SUBROUTINE print_matrix_occupation
2648 :
2649 : ! **************************************************************************************************
2650 : !> \brief Marks the grid blocks whose centroid lies within cutoff of the bounding box of the
2651 : !> panel [blk0, blk1]'s chunk centroids.
2652 : !> \param centroids ...
2653 : !> \param blk0 ...
2654 : !> \param blk1 ...
2655 : !> \param cutoff ...
2656 : !> \param used ...
2657 : ! **************************************************************************************************
2658 92 : SUBROUTINE mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, used)
2659 :
2660 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
2661 : INTEGER, INTENT(IN) :: blk0, blk1
2662 : REAL(KIND=dp), INTENT(IN) :: cutoff
2663 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: used
2664 :
2665 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mask_grid_blocks_near_panel'
2666 :
2667 : INTEGER :: c, handle, k
2668 : REAL(KIND=dp) :: cutoff2, d2, dx
2669 : REAL(KIND=dp), DIMENSION(3) :: hi, lo
2670 :
2671 92 : CALL timeset(routineN, handle)
2672 :
2673 92 : cutoff2 = cutoff**2
2674 92 : lo(:) = MINVAL(centroids(:, blk0:blk1), DIM=2)
2675 92 : hi(:) = MAXVAL(centroids(:, blk0:blk1), DIM=2)
2676 :
2677 276 : ALLOCATE (used(SIZE(centroids, 2)))
2678 2024 : DO c = 1, SIZE(centroids, 2)
2679 : d2 = 0.0_dp
2680 7728 : DO k = 1, 3
2681 5796 : dx = MAX(0.0_dp, lo(k) - centroids(k, c), centroids(k, c) - hi(k))
2682 7728 : d2 = d2 + dx*dx
2683 : END DO
2684 2024 : used(c) = (d2 <= cutoff2)
2685 : END DO
2686 :
2687 92 : CALL timestop(handle)
2688 :
2689 92 : END SUBROUTINE mask_grid_blocks_near_panel
2690 :
2691 : ! **************************************************************************************************
2692 : !> \brief Exact allocated-element count of the geo template of panel [blk0, blk1]: the very same
2693 : !> per-block-pair centroid test as build_geo_template_panel, so this is the true DBCSR
2694 : !> data size of A_pan/B_pan/C_pan (DBCSR stores whole blocks).
2695 : !> \param r_blk_sizes ...
2696 : !> \param centroids ...
2697 : !> \param used ...
2698 : !> \param blk0 ...
2699 : !> \param blk1 ...
2700 : !> \param cutoff ...
2701 : !> \param nze_tmpl ...
2702 : ! **************************************************************************************************
2703 8 : SUBROUTINE panel_template_elems(r_blk_sizes, centroids, used, blk0, blk1, cutoff, nze_tmpl)
2704 :
2705 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes
2706 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
2707 : LOGICAL, DIMENSION(:), INTENT(IN) :: used
2708 : INTEGER, INTENT(IN) :: blk0, blk1
2709 : REAL(KIND=dp), INTENT(IN) :: cutoff
2710 : INTEGER(KIND=int_8), INTENT(OUT) :: nze_tmpl
2711 :
2712 : CHARACTER(LEN=*), PARAMETER :: routineN = 'panel_template_elems'
2713 :
2714 : INTEGER :: c, handle, ib, n_used
2715 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: used_idx
2716 : REAL(KIND=dp) :: cutoff2
2717 :
2718 8 : CALL timeset(routineN, handle)
2719 :
2720 : ! Compress the near mask once so the pair loop only visits candidate columns.
2721 176 : n_used = COUNT(used)
2722 24 : ALLOCATE (used_idx(n_used))
2723 8 : n_used = 0
2724 176 : DO c = 1, SIZE(used)
2725 176 : IF (used(c)) THEN
2726 168 : n_used = n_used + 1
2727 168 : used_idx(n_used) = c
2728 : END IF
2729 : END DO
2730 :
2731 8 : cutoff2 = cutoff**2
2732 8 : nze_tmpl = 0_int_8
2733 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(blk0, blk1, n_used, used_idx, centroids, cutoff2, &
2734 8 : !$OMP r_blk_sizes) PRIVATE(ib, c) REDUCTION(+:nze_tmpl)
2735 : DO ib = blk0, blk1
2736 : DO c = 1, n_used
2737 : IF (SUM((centroids(:, ib) - centroids(:, used_idx(c)))**2) <= cutoff2) THEN
2738 : nze_tmpl = nze_tmpl + INT(r_blk_sizes(ib), int_8)*INT(r_blk_sizes(used_idx(c)), int_8)
2739 : END IF
2740 : END DO
2741 : END DO
2742 : !$OMP END PARALLEL DO
2743 :
2744 8 : CALL timestop(handle)
2745 :
2746 16 : END SUBROUTINE panel_template_elems
2747 :
2748 : ! **************************************************************************************************
2749 : !> \brief Per-rank peak memory (GB) of one panel step of the neighborhood-restricted
2750 : !> contractions: three grid x grid panels of the template size (A_pan, B_pan, C_pan)
2751 : !> plus the grid x RI intermediates (tmp2 and the accumulation operand) and the
2752 : !> grid x AO intermediate (tmpA), whose column support is the panel's geometric
2753 : !> neighborhood fraction f_near = width/n_grid. Shared by the panel planner and
2754 : !> \param nze_tmpl ...
2755 : !> \param pan_rows ...
2756 : !> \param width ...
2757 : !> \param n_grid_total ...
2758 : !> \param n_RI ...
2759 : !> \param n_ao ...
2760 : !> \param n_procs ...
2761 : !> \param mem_GB ...
2762 : ! **************************************************************************************************
2763 8 : SUBROUTINE panel_mem_estimate_GB(nze_tmpl, pan_rows, width, n_grid_total, n_RI, n_ao, &
2764 : n_procs, mem_GB)
2765 :
2766 : INTEGER(KIND=int_8), INTENT(IN) :: nze_tmpl
2767 : INTEGER, INTENT(IN) :: pan_rows, width, n_grid_total, n_RI, &
2768 : n_ao, n_procs
2769 : REAL(KIND=dp), INTENT(OUT) :: mem_GB
2770 :
2771 : REAL(KIND=dp) :: f_near
2772 :
2773 8 : f_near = REAL(width, dp)/REAL(MAX(n_grid_total, 1), dp)
2774 : mem_GB = (3.0_dp*REAL(nze_tmpl, dp) + &
2775 : REAL(pan_rows, dp)*f_near*(2.0_dp*REAL(n_RI, dp) + REAL(n_ao, dp)))* &
2776 8 : 8.0_dp/REAL(MAX(n_procs, 1), dp)*1.0E-9_dp
2777 :
2778 8 : END SUBROUTINE panel_mem_estimate_GB
2779 :
2780 : ! **************************************************************************************************
2781 : !> \brief Plans the panel boundaries for the streaming contractions. Panels grow by whole grid
2782 : !> row-blocks towards ~panel_size rows. When the neighborhood restriction is active
2783 : !> (centroids+cutoff), each candidate panel is additionally checked against
2784 : !> (a) the 32-bit message bound with the panel's TRUE occupancy
2785 : !> (b) the per-rank memory budget: panel_mem_estimate_GB <= mem_budget_GB.
2786 : !> \param bs_env ...
2787 : !> \param r_blk_sizes ...
2788 : !> \param panel_size ...
2789 : !> \param min_dim ...
2790 : !> \param pan_first ...
2791 : !> \param pan_last ...
2792 : !> \param centroids ...
2793 : !> \param cutoff ...
2794 : !> \param n_RI ...
2795 : !> \param n_ao ...
2796 : !> \param n_procs ...
2797 : !> \param mem_budget_GB ...
2798 : !> \param honor_exact ...
2799 : !> \param unsafe ...
2800 : ! **************************************************************************************************
2801 36 : SUBROUTINE plan_grid_panels(bs_env, r_blk_sizes, panel_size, min_dim, pan_first, pan_last, &
2802 36 : centroids, cutoff, n_RI, n_ao, n_procs, mem_budget_GB, &
2803 : honor_exact, unsafe)
2804 :
2805 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2806 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes
2807 : INTEGER, INTENT(IN) :: panel_size, min_dim
2808 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: pan_first, pan_last
2809 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
2810 : OPTIONAL :: centroids
2811 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
2812 : INTEGER, INTENT(IN), OPTIONAL :: n_RI, n_ao, n_procs
2813 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: mem_budget_GB
2814 : LOGICAL, INTENT(IN), OPTIONAL :: honor_exact
2815 : LOGICAL, INTENT(OUT), OPTIONAL :: unsafe
2816 :
2817 : CHARACTER(LEN=*), PARAMETER :: routineN = 'plan_grid_panels'
2818 :
2819 : INTEGER :: blk0, blk1, handle, ib, n_grid_blocks, &
2820 : n_grid_total, n_panels, rows_acc, &
2821 : TARGET, width
2822 : INTEGER(KIND=int_8) :: msg, nze_tmpl, side
2823 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_first, tmp_last
2824 : LOGICAL :: fits, my_honor_exact, my_unsafe, &
2825 : use_cutoff
2826 36 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
2827 : REAL(KIND=dp) :: f_near, mem_GB
2828 :
2829 36 : CALL timeset(routineN, handle)
2830 :
2831 36 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
2832 36 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
2833 2 : IF (use_cutoff) THEN
2834 2 : CPASSERT(PRESENT(n_RI) .AND. PRESENT(n_ao) .AND. PRESENT(n_procs))
2835 : END IF
2836 :
2837 : ! honor_exact: use N_PANELS as requested -- do NOT split a panel further even if it trips
2838 : ! the message-overflow / memory-budget check; instead flag `unsafe` so the caller can warn.
2839 36 : my_honor_exact = .FALSE.
2840 36 : IF (PRESENT(honor_exact)) my_honor_exact = honor_exact
2841 36 : my_unsafe = .FALSE.
2842 :
2843 36 : n_grid_blocks = SIZE(r_blk_sizes)
2844 546 : n_grid_total = SUM(r_blk_sizes)
2845 144 : ALLOCATE (tmp_first(n_grid_blocks), tmp_last(n_grid_blocks))
2846 :
2847 36 : n_panels = 0
2848 36 : blk0 = 1
2849 74 : DO WHILE (blk0 <= n_grid_blocks)
2850 38 : TARGET = panel_size
2851 0 : DO
2852 38 : rows_acc = 0
2853 38 : blk1 = blk0
2854 512 : DO ib = blk0, n_grid_blocks
2855 510 : rows_acc = rows_acc + r_blk_sizes(ib)
2856 510 : blk1 = ib
2857 512 : IF (rows_acc >= TARGET) EXIT
2858 : END DO
2859 38 : IF (.NOT. use_cutoff .OR. blk1 == blk0) EXIT
2860 4 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, used)
2861 88 : width = SUM(r_blk_sizes, MASK=used)
2862 4 : CALL panel_template_elems(r_blk_sizes, centroids, used, blk0, blk1, cutoff, nze_tmpl)
2863 4 : f_near = REAL(width, dp)/REAL(MAX(n_grid_total, 1), dp)
2864 4 : side = INT(REAL(rows_acc, dp)*f_near*REAL(MAX(n_RI, n_ao), dp), int_8)
2865 4 : msg = MAX(nze_tmpl, side)/INT(MAX(min_dim, 1), int_8)
2866 4 : fits = (msg <= bs_env%dbcsr_msg_elem_limit/4)
2867 4 : IF (fits .AND. PRESENT(mem_budget_GB)) THEN
2868 4 : IF (mem_budget_GB > 0.0_dp) THEN
2869 : CALL panel_mem_estimate_GB(nze_tmpl, rows_acc, width, n_grid_total, &
2870 4 : n_RI, n_ao, n_procs, mem_GB)
2871 4 : fits = (mem_GB <= mem_budget_GB)
2872 : END IF
2873 : END IF
2874 : ! Panel size is bounded only by the message-overflow and memory checks above; there is
2875 : ! no neighborhood-width (f_near) cap. mp_waitall is dominated by the NUMBER of panel
2876 : ! multiplies, so fewer/larger panels are cheaper here -- panel count is driven DOWN by
2877 : ! the N_PANELS keyword (panel_size), not split up by a width heuristic.
2878 4 : IF (my_honor_exact) THEN
2879 : ! Keep exactly the requested grouping; just record if it exceeds a safety limit.
2880 4 : IF (.NOT. fits) my_unsafe = .TRUE.
2881 : EXIT
2882 : END IF
2883 0 : IF (fits) EXIT
2884 34 : TARGET = MAX(1, MIN(TARGET, rows_acc)/2)
2885 : END DO
2886 38 : n_panels = n_panels + 1
2887 38 : tmp_first(n_panels) = blk0
2888 38 : tmp_last(n_panels) = blk1
2889 38 : blk0 = blk1 + 1
2890 : END DO
2891 :
2892 144 : ALLOCATE (pan_first(n_panels), pan_last(n_panels))
2893 74 : pan_first(:) = tmp_first(1:n_panels)
2894 74 : pan_last(:) = tmp_last(1:n_panels)
2895 36 : DEALLOCATE (tmp_first, tmp_last)
2896 :
2897 36 : IF (PRESENT(unsafe)) unsafe = my_unsafe
2898 :
2899 36 : CALL timestop(handle)
2900 :
2901 72 : END SUBROUTINE plan_grid_panels
2902 :
2903 : ! **************************************************************************************************
2904 : !> \brief Resolves the panel boundaries for the streaming contractions from the bs_env settings:
2905 : !> \param bs_env ...
2906 : !> \param mat_phi_mu_l ...
2907 : !> \param pan_first ...
2908 : !> \param pan_last ...
2909 : ! **************************************************************************************************
2910 158 : SUBROUTINE resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
2911 :
2912 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2913 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
2914 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: pan_first, pan_last
2915 :
2916 : CHARACTER(LEN=*), PARAMETER :: routineN = 'resolve_grid_panels'
2917 :
2918 : CHARACTER(LEN=max_line_length) :: msg
2919 : INTEGER :: handle, min_dim, n_grid_total, &
2920 : n_panels_req, npcols, nprows, &
2921 : panel_size, safe_max
2922 158 : INTEGER, DIMENSION(:), POINTER :: r_blk_sizes
2923 : LOGICAL :: honor_exact, panels_unsafe, use_cutoff
2924 : REAL(KIND=dp) :: mem_avail_GB, mem_budget_GB
2925 : TYPE(dbcsr_distribution_type) :: dist
2926 :
2927 158 : CALL timeset(routineN, handle)
2928 :
2929 158 : IF (ALLOCATED(bs_env%ri_rs%pan_first)) THEN
2930 494 : ALLOCATE (pan_first, SOURCE=bs_env%ri_rs%pan_first)
2931 494 : ALLOCATE (pan_last, SOURCE=bs_env%ri_rs%pan_last)
2932 122 : CALL timestop(handle)
2933 : RETURN
2934 : END IF
2935 :
2936 : use_cutoff = bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp .AND. &
2937 36 : ALLOCATED(bs_env%ri_rs%chunk_centroids)
2938 :
2939 : ! MIN(nprows, npcols) is the divisor that bounds the worst-rank Cannon message: a
2940 : ! P x n_grid panel is replicated into block row strips (P/nprows x n_grid) or column
2941 : ! strips (P x n_grid/npcols) during multiply_cannon, so the largest single-rank
2942 : ! message is ~ P*n_grid / MIN(nprows,npcols) elements.
2943 : CALL dbcsr_get_info(mat_phi_mu_l, nfullrows_total=n_grid_total, row_blk_size=r_blk_sizes, &
2944 36 : distribution=dist)
2945 36 : CALL dbcsr_distribution_get(dist, nprows=nprows, npcols=npcols)
2946 36 : min_dim = MAX(MIN(nprows, npcols), 1)
2947 :
2948 : ! Panel height such that NO per-rank DBCSR message can overflow the 32-bit length field
2949 : ! (see bs_env%dbcsr_msg_elem_limit): requiring the worst-rank message to stay under
2950 : ! 0.5 * HUGE(int_4) gives the safe height P_safe = 0.5 * HUGE(int_4) * min_dim / n_grid.
2951 36 : IF (use_cutoff) THEN
2952 2 : safe_max = n_grid_total
2953 : ELSE
2954 : safe_max = INT(0.5_dp*REAL(bs_env%dbcsr_msg_elem_limit, dp)*REAL(min_dim, dp)/ &
2955 34 : REAL(n_grid_total, dp))
2956 34 : safe_max = MAX(1, MIN(safe_max, n_grid_total))
2957 : END IF
2958 :
2959 : ! A user-set N_PANELS ( > 1 ) is honored EXACTLY: the planner produces that many panels
2960 : ! (up to grid-block granularity) and never force-splits them for the message/memory safety
2961 : ! limits -- if a limit is tripped it warns instead of silently changing the count.
2962 36 : n_panels_req = bs_env%ri_rs%n_panels
2963 36 : honor_exact = (n_panels_req > 1)
2964 36 : panels_unsafe = .FALSE.
2965 36 : IF (n_panels_req > 1) THEN
2966 : ! ceil(n_grid/n_panels_req) rows per panel => exactly n_panels_req panels. With the
2967 : ! cutoff active safe_max = n_grid_total (no clamp, honored exactly); without it, safe_max
2968 : ! is the int32-overflow ceiling and MUST still bound the panel (the non-cutoff planner
2969 : ! loop has no in-loop message-size check).
2970 2 : panel_size = MIN((n_grid_total + n_panels_req - 1)/n_panels_req, safe_max)
2971 : ELSE
2972 : ! Default (<= 1): a single whole-grid panel, clamped to the overflow-safe ceiling.
2973 34 : panel_size = safe_max
2974 : END IF
2975 36 : panel_size = MAX(1, panel_size)
2976 :
2977 36 : IF (use_cutoff) THEN
2978 : ! Half of the measured free memory as panel budget.
2979 2 : CALL ri_rs_mem_avail_per_proc_GB(bs_env, mem_avail_GB)
2980 2 : mem_budget_GB = 0.5_dp*mem_avail_GB
2981 : CALL plan_grid_panels(bs_env, r_blk_sizes, panel_size, min_dim, pan_first, pan_last, &
2982 : centroids=bs_env%ri_rs%chunk_centroids, &
2983 : cutoff=bs_env%ri_rs%cutoff_radius_v_w, &
2984 : n_RI=bs_env%n_RI, n_ao=bs_env%n_ao, &
2985 : n_procs=bs_env%para_env%num_pe, mem_budget_GB=mem_budget_GB, &
2986 2 : honor_exact=honor_exact, unsafe=panels_unsafe)
2987 : ELSE
2988 34 : CALL plan_grid_panels(bs_env, r_blk_sizes, panel_size, min_dim, pan_first, pan_last)
2989 : END IF
2990 :
2991 36 : IF (honor_exact .AND. panels_unsafe) THEN
2992 : WRITE (msg, '(A,I0,A)') &
2993 0 : "N_PANELS = ", n_panels_req, " is used as requested, but one or more panels "// &
2994 : "exceed the DBCSR 32-bit message length or the memory budget. The run may abort "// &
2995 0 : "or swap; increase N_PANELS if it does."
2996 0 : CPWARN(TRIM(msg))
2997 : END IF
2998 :
2999 146 : ALLOCATE (bs_env%ri_rs%pan_first, SOURCE=pan_first)
3000 146 : ALLOCATE (bs_env%ri_rs%pan_last, SOURCE=pan_last)
3001 :
3002 36 : CALL timestop(handle)
3003 :
3004 194 : END SUBROUTINE resolve_grid_panels
3005 :
3006 : ! **************************************************************************************************
3007 : !> \brief Available memory per MPI process (GB). /proc/meminfo reports node-wide memory, so
3008 : !> every rank on a node reads the SAME MemLikelyFree; the per-process share is
3009 : !> node_free / ranks_per_node (ranks grouped by a hostname hash exchanged via allgather).
3010 : !> Returns the MIN across all ranks (most-constrained node)
3011 : !> \param bs_env ...
3012 : !> \param mem_avail_GB ...
3013 : ! **************************************************************************************************
3014 80 : SUBROUTINE ri_rs_mem_avail_per_proc_GB(bs_env, mem_avail_GB)
3015 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3016 : REAL(KIND=dp), INTENT(OUT) :: mem_avail_GB
3017 :
3018 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ri_rs_mem_avail_per_proc_GB'
3019 :
3020 : CHARACTER(LEN=default_string_length) :: hostname
3021 : INTEGER :: handle, host_hash, ic, n_procs, &
3022 : ranks_per_node
3023 : INTEGER(KIND=int_8) :: h8, mem_buffers, mem_cached, mem_free, &
3024 : mem_likely_free, mem_slab, &
3025 : mem_sreclaimable, mem_total
3026 : INTEGER, ALLOCATABLE, DIMENSION(:) :: all_host_hashes
3027 :
3028 80 : CALL timeset(routineN, handle)
3029 :
3030 80 : n_procs = bs_env%para_env%num_pe
3031 80 : CALL m_hostnm(hostname)
3032 80 : h8 = 0_int_8
3033 1040 : DO ic = 1, LEN_TRIM(hostname)
3034 1040 : h8 = MOD(h8*127_int_8 + INT(ICHAR(hostname(ic:ic)), int_8), 2147483647_int_8)
3035 : END DO
3036 80 : host_hash = INT(h8)
3037 240 : ALLOCATE (all_host_hashes(n_procs))
3038 80 : CALL bs_env%para_env%allgather(host_hash, all_host_hashes)
3039 240 : ranks_per_node = MAX(COUNT(all_host_hashes == host_hash), 1)
3040 80 : DEALLOCATE (all_host_hashes)
3041 :
3042 : CALL m_memory_details(MemTotal=mem_total, MemFree=mem_free, Buffers=mem_buffers, &
3043 : Cached=mem_cached, Slab=mem_slab, SReclaimable=mem_sreclaimable, &
3044 80 : MemLikelyFree=mem_likely_free)
3045 80 : mem_avail_GB = REAL(mem_likely_free, dp)*1.0E-9_dp/REAL(ranks_per_node, dp)
3046 80 : CALL bs_env%para_env%min(mem_avail_GB)
3047 :
3048 80 : CALL timestop(handle)
3049 :
3050 80 : END SUBROUTINE ri_rs_mem_avail_per_proc_GB
3051 :
3052 : ! **************************************************************************************************
3053 : !> \brief Estimates and prints per-process memory requirements for the RI-RS GW calculation.
3054 : !> \param qs_env ...
3055 : !> \param bs_env ...
3056 : ! **************************************************************************************************
3057 36 : SUBROUTINE print_ri_rs_memory_estimate(qs_env, bs_env)
3058 :
3059 : !$ USE OMP_LIB, ONLY: omp_get_max_threads
3060 :
3061 : TYPE(qs_environment_type), POINTER :: qs_env
3062 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3063 :
3064 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_ri_rs_memory_estimate'
3065 :
3066 : CHARACTER(LEN=max_line_length) :: msg
3067 : INTEGER :: handle, iatom, ipan, l, &
3068 : max_n_ao_used, max_n_local_grid, &
3069 : n_ao_used_atom, n_grid_total, &
3070 : n_local_grid, n_loc_ri_max, n_procs, &
3071 : n_procs_per_atom, n_RI, n_threads, &
3072 : natom, pan_rows, pan_width
3073 : INTEGER(KIND=int_8) :: nze_tmpl
3074 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
3075 36 : INTEGER, DIMENSION(:), POINTER :: r_blk_sizes
3076 : LOGICAL :: use_cutoff
3077 36 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used
3078 : REAL(KIND=dp) :: cutoff_ri, mem_avail_GB, mem_D_local_GB, &
3079 : mem_dlp_GB, mem_pan_GB, mem_panels_GB, &
3080 : mem_phi_local_GB, mem_Z_lP_GB, &
3081 : mem_Zlp_peak_GB, pos_P(3)
3082 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3083 :
3084 36 : CALL timeset(routineN, handle)
3085 :
3086 : CALL dbcsr_get_info(bs_env%ri_rs%mat_phi_mu_l, nfullrows_total=n_grid_total, &
3087 36 : row_blk_size=r_blk_sizes)
3088 36 : n_RI = bs_env%n_RI
3089 36 : n_procs = bs_env%para_env%num_pe
3090 :
3091 : ! Z_lP upper bound: dense n_grid × n_RI, distributed evenly across all ranks.
3092 : ! The actual sparse Z_lP is smaller due to the per-atom locality cutoff.
3093 : mem_Z_lP_GB = REAL(n_grid_total, dp)*REAL(n_RI, dp)*8.0_dp/ &
3094 36 : REAL(n_procs, dp)*1.0E-9_dp
3095 :
3096 : ! Peak panel memory during Σ^c: two G panels (A_occ, A_vir) + one W panel plus the
3097 : ! grid × RI / grid × AO intermediates. With the CUTOFF_RADIUS_RL_W restriction the panel
3098 : ! matrices only allocate the geo-template blocks, so use the same nze-aware model as the
3099 : ! panel planner (panel_mem_estimate_GB); without the cutoff, dense panel_rows × n_grid.
3100 : ! Plus the n_RI × n_RI W_aux matrix. Distributed over n_procs ranks.
3101 : use_cutoff = bs_env%ri_rs%cutoff_radius_v_w > 0.0_dp .AND. &
3102 36 : ALLOCATED(bs_env%ri_rs%chunk_centroids)
3103 36 : CALL resolve_grid_panels(bs_env, bs_env%ri_rs%mat_phi_mu_l, pan_first, pan_last)
3104 36 : mem_panels_GB = 0.0_dp
3105 74 : DO ipan = 1, SIZE(pan_first)
3106 548 : pan_rows = SUM(r_blk_sizes(pan_first(ipan):pan_last(ipan)))
3107 38 : IF (use_cutoff) THEN
3108 : CALL mask_grid_blocks_near_panel(bs_env%ri_rs%chunk_centroids, pan_first(ipan), &
3109 : pan_last(ipan), bs_env%ri_rs%cutoff_radius_v_w, &
3110 4 : grid_used)
3111 88 : pan_width = SUM(r_blk_sizes, MASK=grid_used)
3112 : CALL panel_template_elems(r_blk_sizes, bs_env%ri_rs%chunk_centroids, &
3113 : grid_used, pan_first(ipan), pan_last(ipan), &
3114 4 : bs_env%ri_rs%cutoff_radius_v_w, nze_tmpl)
3115 : CALL panel_mem_estimate_GB(nze_tmpl, pan_rows, pan_width, n_grid_total, &
3116 4 : n_RI, bs_env%n_ao, n_procs, mem_pan_GB)
3117 : ELSE
3118 34 : pan_width = n_grid_total
3119 : mem_pan_GB = (3.0_dp*REAL(pan_rows, dp)*REAL(pan_width, dp) + &
3120 : 2.0_dp*REAL(pan_rows, dp)*REAL(n_RI, dp))* &
3121 34 : 8.0_dp/REAL(n_procs, dp)*1.0E-9_dp
3122 : END IF
3123 74 : mem_panels_GB = MAX(mem_panels_GB, mem_pan_GB)
3124 : END DO
3125 : mem_panels_GB = mem_panels_GB + &
3126 36 : REAL(n_RI, dp)*REAL(n_RI, dp)*8.0_dp/REAL(n_procs, dp)*1.0E-9_dp
3127 :
3128 : ! Z_lP SOLVE peak (compute_coeff_Z_lP). For the atom P with the largest integration
3129 : ! sphere, one rank holds simultaneously:
3130 : ! D_local : n_local_grid x n_local_grid (dense Gram, BLAS path only; O(n_local_grid^2))
3131 : ! phi_local: n_local_grid x n_ao_used (AOs reaching into the sphere only)
3132 : ! d_lp : n_local_grid x n_loc_ri, replicated once + one private copy per OMP thread
3133 : ! n_local_grid = # grid points within cutoff_ri(P) = CUTOFF_RADIUS_RL_RI (if > 0) else
3134 : ! r_c(RI metric) + r_RI(P). This is NOT evenly distributed: n_local_grid depends on the
3135 : ! local density of atoms/grid, so the rank owning the densest atom peaks well above the
3136 : ! average. We report the worst-case (max over atoms) as a per-rank upper bound.
3137 36 : CALL get_qs_env(qs_env, particle_set=particle_set)
3138 36 : natom = SIZE(particle_set)
3139 :
3140 36 : max_n_local_grid = 0
3141 36 : n_loc_ri_max = 0
3142 36 : max_n_ao_used = 0
3143 126 : DO iatom = 1, natom
3144 90 : IF (bs_env%ri_rs%cutoff_radius_ri_rs > 0.0_dp) THEN
3145 6 : cutoff_ri = bs_env%ri_rs%cutoff_radius_ri_rs
3146 : ELSE
3147 84 : cutoff_ri = bs_env%ri_metric%cutoff_radius + bs_env%ri_rs%radius_ri_per_atom(iatom)
3148 : END IF
3149 360 : pos_P(:) = particle_set(iatom)%r(:)
3150 : n_local_grid = 0
3151 36562 : DO l = 1, n_grid_total
3152 145978 : IF (SUM((bs_env%ri_rs%grid_points(1:3, l) - pos_P(1:3))**2) <= cutoff_ri**2) THEN
3153 35086 : n_local_grid = n_local_grid + 1
3154 : END IF
3155 : END DO
3156 90 : max_n_local_grid = MAX(max_n_local_grid, n_local_grid)
3157 90 : CALL get_n_ao_in_sphere(bs_env, particle_set, iatom, cutoff_ri, n_ao_used_atom)
3158 90 : max_n_ao_used = MAX(max_n_ao_used, n_ao_used_atom)
3159 : n_loc_ri_max = MAX(n_loc_ri_max, &
3160 126 : bs_env%i_RI_end_from_atom(iatom) - bs_env%i_RI_start_from_atom(iatom) + 1)
3161 : END DO
3162 :
3163 36 : n_procs_per_atom = MIN(MAX(bs_env%ri_rs%n_procs_per_atom_z_lp, 1), n_procs)
3164 36 : n_threads = 1
3165 36 : !$ n_threads = omp_get_max_threads()
3166 :
3167 : ! D_local: dense on one rank for the BLAS path; block-cyclic over the subgroup (=> /G) for
3168 : ! the ScaLAPACK path (N_PROCS_PER_ATOM_Z_LP = G > 1). phi_local/d_lp stay per-rank either way.
3169 36 : IF (n_procs_per_atom > 1) THEN
3170 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
3171 : ELSE
3172 36 : mem_D_local_GB = REAL(max_n_local_grid, dp)**2*8.0_dp*1.0E-9_dp
3173 : END IF
3174 36 : mem_phi_local_GB = REAL(max_n_local_grid, dp)*REAL(max_n_ao_used, dp)*8.0_dp*1.0E-9_dp
3175 : mem_dlp_GB = REAL(max_n_local_grid, dp)*REAL(n_loc_ri_max, dp)*8.0_dp* &
3176 36 : REAL(1 + n_threads, dp)*1.0E-9_dp
3177 36 : mem_Zlp_peak_GB = mem_D_local_GB + mem_phi_local_GB + mem_dlp_GB
3178 :
3179 : ! Available memory per process = node MemLikelyFree / ranks-per-node, min across ranks
3180 : ! (0 on non-Linux => warnings suppressed below).
3181 36 : CALL ri_rs_mem_avail_per_proc_GB(bs_env, mem_avail_GB)
3182 :
3183 36 : IF (bs_env%unit_nr > 0) THEN
3184 18 : WRITE (bs_env%unit_nr, '(A)') ' '
3185 18 : WRITE (bs_env%unit_nr, '(T2,A)') 'RI-RS memory estimate per MPI process:'
3186 : WRITE (bs_env%unit_nr, '(T4,A,F37.2,A)') &
3187 18 : 'Available memory per process (system)', mem_avail_GB, ' GB'
3188 : WRITE (bs_env%unit_nr, '(T4,A,F18.2,A)') &
3189 18 : 'Required for Z_lP (dense upper bound; actual is sparser)', mem_Z_lP_GB, ' GB'
3190 : WRITE (bs_env%unit_nr, '(T4,A,F25.2,A)') &
3191 18 : 'Required for χ, W, Σ panels (peak per panel step)', mem_panels_GB, ' GB'
3192 : WRITE (bs_env%unit_nr, '(T4,A,F17.2,A)') &
3193 18 : 'Required for Z_lP solve peak (D_local+φ, worst-case atom)', mem_Zlp_peak_GB, ' GB'
3194 : WRITE (bs_env%unit_nr, '(T6,A,I21,A,F10.2,A)') &
3195 18 : 'worst-case n_local_grid', max_n_local_grid, ' points (D_local', mem_D_local_GB, ' GB)'
3196 18 : WRITE (bs_env%unit_nr, '(A)') ' '
3197 :
3198 : END IF
3199 :
3200 36 : IF (mem_avail_GB > 0.0_dp .AND. mem_Z_lP_GB > mem_avail_GB) THEN
3201 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
3202 0 : "The estimated memory for Z_lP, ", mem_Z_lP_GB, " GB per process, exceeds the "// &
3203 0 : "available ", mem_avail_GB, " GB. Z_lP (n_grid x n_RI) is distributed across all "// &
3204 : "MPI ranks, so add nodes, use fewer MPI ranks per node, or raise "// &
3205 : "N_PROCS_PER_ATOM_Z_LP to distribute each atom block via ScaLAPACK, which reduces "// &
3206 0 : "the per-rank memory roughly by the number of ranks per atom."
3207 0 : CPWARN(TRIM(msg))
3208 : END IF
3209 :
3210 36 : IF (mem_avail_GB > 0.0_dp .AND. mem_panels_GB > mem_avail_GB) THEN
3211 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
3212 0 : "The estimated peak memory of the chi/W/Sigma panels, ", mem_panels_GB, &
3213 0 : " GB per process, exceeds the available ", mem_avail_GB, " GB. Panel memory "// &
3214 : "scales roughly as 3*panel_size*n_grid/n_procs, so add nodes, use fewer MPI ranks "// &
3215 0 : "per node, or raise N_PANELS for more but smaller panels."
3216 0 : CPWARN(TRIM(msg))
3217 : END IF
3218 :
3219 36 : IF (mem_avail_GB > 0.0_dp .AND. mem_Zlp_peak_GB > mem_avail_GB) THEN
3220 : WRITE (msg, '(A,F0.2,A,F0.2,A)') &
3221 0 : "The estimated peak memory of the Z_lP solve, ", mem_Zlp_peak_GB, &
3222 0 : " GB per process, exceeds the available ", mem_avail_GB, &
3223 : " GB. The per-atom Gram matrix D_local(n_local_grid, n_local_grid) dominates and "// &
3224 : "scales as n_local_grid^2, and it is not balanced across ranks: the rank owning "// &
3225 : "the atom with the largest integration sphere peaks well above the average. "// &
3226 : "Either raise N_PROCS_PER_ATOM_Z_LP to distribute D_local block-cyclic via "// &
3227 : "ScaLAPACK, which reduces that term roughly by the number of ranks per atom at no "// &
3228 : "loss of accuracy, or lower CUTOFF_RADIUS_RL_RI, which shrinks D_local as "// &
3229 : "n_local_grid^2 but trades accuracy, or use fewer MPI ranks per node so that each "// &
3230 0 : "rank has more memory for the peak atom."
3231 0 : CPWARN(TRIM(msg))
3232 : END IF
3233 :
3234 36 : CALL timestop(handle)
3235 :
3236 108 : END SUBROUTINE print_ri_rs_memory_estimate
3237 :
3238 : ! **************************************************************************************************
3239 : !> \brief Creates an empty (panel_chunks x neighborhood_chunks) DBCSR matrix with zero blocks
3240 : !> pre-allocated only where |centroid(panel_row r) - centroid(column c)| <= cutoff.
3241 : !> Used with retain_sparsity=.TRUE. in the subsequent dbcsr_multiply so distant blocks
3242 : !> of the grid-basis panels (φ G φ^T, Z W Z^T, ...) are never computed at all.
3243 : !> \param L_pan ...
3244 : !> \param L_full ...
3245 : !> \param centroids ...
3246 : !> \param cutoff ...
3247 : !> \param blk0 ...
3248 : !> \param A_template ...
3249 : !> \param col_map ...
3250 : ! **************************************************************************************************
3251 124 : SUBROUTINE build_geo_template_panel(L_pan, L_full, centroids, cutoff, blk0, A_template, col_map)
3252 : TYPE(dbcsr_type), INTENT(IN) :: L_pan, L_full
3253 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
3254 : REAL(KIND=dp), INTENT(IN) :: cutoff
3255 : INTEGER, INTENT(IN) :: blk0
3256 : TYPE(dbcsr_type), INTENT(OUT) :: A_template
3257 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: col_map
3258 :
3259 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_geo_template_panel'
3260 :
3261 : INTEGER :: c, cg, cs, handle, my_pcol, my_prow, &
3262 : n_grid_blks, n_pan_blks, npcols, &
3263 : nprows, r, rs
3264 124 : INTEGER, DIMENSION(:), POINTER :: grid_blk_sizes, pan_blk_sizes
3265 : REAL(KIND=dp) :: cutoff2
3266 124 : REAL(KIND=dp), ALLOCATABLE :: zero_blk(:, :)
3267 : TYPE(dbcsr_distribution_type) :: dist
3268 :
3269 124 : CALL timeset(routineN, handle)
3270 :
3271 124 : cutoff2 = cutoff**2
3272 124 : CALL dbcsr_get_info(L_pan, nblkrows_total=n_pan_blks, row_blk_size=pan_blk_sizes)
3273 124 : CALL dbcsr_get_info(L_full, nblkrows_total=n_grid_blks, row_blk_size=grid_blk_sizes)
3274 :
3275 : ! create_product_matrix assigns row r to process MOD(r-1,nprows) and
3276 : ! col c to MOD(c-1,npcols), so we can determine local ownership analytically.
3277 124 : CALL create_product_matrix(L_pan, L_full, 'N', 'T', A_template)
3278 124 : CALL dbcsr_get_info(A_template, distribution=dist)
3279 : CALL dbcsr_distribution_get(dist, nprows=nprows, npcols=npcols, &
3280 124 : myprow=my_prow, mypcol=my_pcol)
3281 :
3282 : ALLOCATE (zero_blk(MAXVAL(pan_blk_sizes(1:n_pan_blks)), &
3283 4402 : MAXVAL(grid_blk_sizes(1:n_grid_blks))))
3284 124 : zero_blk(:, :) = 0.0_dp
3285 :
3286 1426 : DO r = 1, n_pan_blks
3287 1302 : IF (MOD(r - 1, nprows) /= my_prow) CYCLE
3288 651 : rs = pan_blk_sizes(r)
3289 14446 : DO c = 1, n_grid_blks
3290 13671 : IF (MOD(c - 1, npcols) /= my_pcol) CYCLE
3291 13671 : cg = c
3292 13671 : IF (PRESENT(col_map)) cg = col_map(c)
3293 13671 : IF ((centroids(1, blk0 + r - 1) - centroids(1, cg))**2 + &
3294 : (centroids(2, blk0 + r - 1) - centroids(2, cg))**2 + &
3295 1302 : (centroids(3, blk0 + r - 1) - centroids(3, cg))**2 <= cutoff2) THEN
3296 12741 : cs = grid_blk_sizes(c)
3297 12741 : CALL dbcsr_put_block(A_template, r, c, zero_blk(1:rs, 1:cs))
3298 : END IF
3299 : END DO
3300 : END DO
3301 124 : CALL dbcsr_finalize(A_template)
3302 :
3303 124 : DEALLOCATE (zero_blk)
3304 124 : CALL timestop(handle)
3305 :
3306 124 : END SUBROUTINE build_geo_template_panel
3307 :
3308 : ! **************************************************************************************************
3309 : !> \brief Slices a contiguous range of grid row-blocks [blk0, blk1] out of a (grid x n) DBCSR
3310 : !> matrix into a new (P x n) panel matrix: iterate the source's local blocks, put the
3311 : !> in-range ones into the panel with a remapped row-block index, then finalize. Row-block
3312 : !> index i of the panel corresponds to source row-block blk0+i-1.
3313 : !> \param mat_full ...
3314 : !> \param blk0 ...
3315 : !> \param blk1 ...
3316 : !> \param mat_panel ...
3317 : ! **************************************************************************************************
3318 2928 : SUBROUTINE extract_grid_panel(mat_full, blk0, blk1, mat_panel)
3319 :
3320 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_full
3321 : INTEGER, INTENT(IN) :: blk0, blk1
3322 : TYPE(dbcsr_type), INTENT(OUT) :: mat_panel
3323 :
3324 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extract_grid_panel'
3325 :
3326 : INTEGER :: handle, ib, jb, npb
3327 2928 : INTEGER, DIMENSION(:), POINTER :: col_blk_full, col_dist_full, &
3328 2928 : row_blk_full, row_blk_pan, &
3329 2928 : row_dist_full, row_dist_pan
3330 2928 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
3331 : TYPE(dbcsr_distribution_type) :: dist_full, dist_pan
3332 : TYPE(dbcsr_iterator_type) :: iter
3333 :
3334 2928 : CALL timeset(routineN, handle)
3335 :
3336 : CALL dbcsr_get_info(mat_full, distribution=dist_full, &
3337 2928 : row_blk_size=row_blk_full, col_blk_size=col_blk_full)
3338 2928 : CALL dbcsr_distribution_get(dist_full, row_dist=row_dist_full, col_dist=col_dist_full)
3339 :
3340 2928 : npb = blk1 - blk0 + 1
3341 11712 : ALLOCATE (row_dist_pan(npb), row_blk_pan(npb))
3342 90504 : row_dist_pan(:) = row_dist_full(blk0:blk1)
3343 90504 : row_blk_pan(:) = row_blk_full(blk0:blk1)
3344 :
3345 : CALL dbcsr_distribution_new(dist_pan, template=dist_full, &
3346 2928 : row_dist=row_dist_pan, col_dist=col_dist_full)
3347 : CALL dbcsr_create(mat_panel, name="grid_panel", dist=dist_pan, &
3348 : matrix_type=dbcsr_type_no_symmetry, &
3349 2928 : row_blk_size=row_blk_pan, col_blk_size=col_blk_full)
3350 :
3351 2928 : CALL dbcsr_iterator_start(iter, mat_full)
3352 54828 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3353 51900 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
3354 51900 : IF (ib < blk0 .OR. ib > blk1) CYCLE
3355 54828 : CALL dbcsr_put_block(mat_panel, ib - blk0 + 1, jb, blk)
3356 : END DO
3357 2928 : CALL dbcsr_iterator_stop(iter)
3358 2928 : CALL dbcsr_finalize(mat_panel)
3359 :
3360 2928 : CALL dbcsr_distribution_release(dist_pan)
3361 2928 : DEALLOCATE (row_dist_pan, row_blk_pan)
3362 :
3363 2928 : CALL timestop(handle)
3364 :
3365 8784 : END SUBROUTINE extract_grid_panel
3366 :
3367 : ! **************************************************************************************************
3368 : !> \brief Marks which column blocks of a DBCSR matrix carry at least one non-zero block anywhere
3369 : !> (global union). Used to restrict the inner index of the panel multiplies to the
3370 : !> AO/RI atoms that actually touch the panel (exact: dropped rows only meet zeros).
3371 : !> \param matrix ...
3372 : !> \param para_env ...
3373 : !> \param used ...
3374 : ! **************************************************************************************************
3375 2318 : SUBROUTINE collect_used_col_blocks(matrix, para_env, used)
3376 :
3377 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
3378 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
3379 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: used
3380 :
3381 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collect_used_col_blocks'
3382 :
3383 : INTEGER :: handle, ib, jb, nblkcols
3384 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iused
3385 2318 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
3386 : TYPE(dbcsr_iterator_type) :: iter
3387 :
3388 2318 : CALL timeset(routineN, handle)
3389 :
3390 2318 : CALL dbcsr_get_info(matrix, nblkcols_total=nblkcols)
3391 6954 : ALLOCATE (iused(nblkcols))
3392 2318 : iused(:) = 0
3393 :
3394 2318 : CALL dbcsr_iterator_start(iter, matrix)
3395 42591 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3396 40273 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
3397 40273 : iused(jb) = 1
3398 : END DO
3399 2318 : CALL dbcsr_iterator_stop(iter)
3400 :
3401 2318 : CALL para_env%sum(iused)
3402 :
3403 4636 : ALLOCATE (used(nblkcols))
3404 7784 : used(:) = (iused(:) > 0)
3405 2318 : DEALLOCATE (iused)
3406 :
3407 2318 : CALL timestop(handle)
3408 :
3409 4636 : END SUBROUTINE collect_used_col_blocks
3410 :
3411 : ! **************************************************************************************************
3412 : !> \brief Copies the flagged block rows (compress_rows=.TRUE.) or block columns (.FALSE.) of a
3413 : !> DBCSR matrix into a compressed matrix. The subset keeps the parent's process assignment
3414 : !> along the compressed dimension, so every block stays on its owning rank: the extraction
3415 : !> is purely local (zero communication), like extract_grid_panel.
3416 : !> \param mat_full ...
3417 : !> \param used ...
3418 : !> \param mat_out ...
3419 : !> \param compress_rows ...
3420 : !> \param blk_map ...
3421 : ! **************************************************************************************************
3422 6224 : SUBROUTINE extract_masked_blocks(mat_full, used, mat_out, compress_rows, blk_map)
3423 :
3424 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_full
3425 : LOGICAL, DIMENSION(:), INTENT(IN) :: used
3426 : TYPE(dbcsr_type), INTENT(OUT) :: mat_out
3427 : LOGICAL, INTENT(IN) :: compress_rows
3428 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT), &
3429 : OPTIONAL :: blk_map
3430 :
3431 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extract_masked_blocks'
3432 :
3433 : INTEGER :: handle, ib, jb, n_blk, n_sub, r
3434 6224 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inv_map
3435 6224 : INTEGER, DIMENSION(:), POINTER :: blk_full, blk_sub, col_blk_full, &
3436 6224 : col_dist_full, dist_full_1d, &
3437 6224 : dist_sub_1d, row_blk_full, &
3438 6224 : row_dist_full
3439 6224 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
3440 : TYPE(dbcsr_distribution_type) :: dist_full, dist_sub
3441 : TYPE(dbcsr_iterator_type) :: iter
3442 :
3443 6224 : CALL timeset(routineN, handle)
3444 :
3445 : CALL dbcsr_get_info(mat_full, distribution=dist_full, &
3446 6224 : row_blk_size=row_blk_full, col_blk_size=col_blk_full)
3447 6224 : CALL dbcsr_distribution_get(dist_full, row_dist=row_dist_full, col_dist=col_dist_full)
3448 :
3449 6224 : IF (compress_rows) THEN
3450 3906 : blk_full => row_blk_full
3451 3906 : dist_full_1d => row_dist_full
3452 : ELSE
3453 2318 : blk_full => col_blk_full
3454 2318 : dist_full_1d => col_dist_full
3455 : END IF
3456 6224 : n_blk = SIZE(blk_full)
3457 :
3458 24024 : n_sub = COUNT(used)
3459 6224 : CPASSERT(n_sub > 0)
3460 37344 : ALLOCATE (inv_map(n_blk), blk_sub(n_sub), dist_sub_1d(n_sub))
3461 6308 : IF (PRESENT(blk_map)) ALLOCATE (blk_map(n_sub))
3462 6224 : inv_map(:) = 0
3463 6224 : r = 0
3464 24024 : DO ib = 1, n_blk
3465 24024 : IF (used(ib)) THEN
3466 17800 : r = r + 1
3467 17800 : inv_map(ib) = r
3468 17800 : blk_sub(r) = blk_full(ib)
3469 17800 : dist_sub_1d(r) = dist_full_1d(ib)
3470 17800 : IF (PRESENT(blk_map)) blk_map(r) = ib
3471 : END IF
3472 : END DO
3473 :
3474 6224 : IF (compress_rows) THEN
3475 : CALL dbcsr_distribution_new(dist_sub, template=dist_full, &
3476 3906 : row_dist=dist_sub_1d, col_dist=col_dist_full)
3477 : CALL dbcsr_create(mat_out, name="row_subset", dist=dist_sub, &
3478 : matrix_type=dbcsr_type_no_symmetry, &
3479 3906 : row_blk_size=blk_sub, col_blk_size=col_blk_full)
3480 : ELSE
3481 : CALL dbcsr_distribution_new(dist_sub, template=dist_full, &
3482 2318 : row_dist=row_dist_full, col_dist=dist_sub_1d)
3483 : CALL dbcsr_create(mat_out, name="col_subset", dist=dist_sub, &
3484 : matrix_type=dbcsr_type_no_symmetry, &
3485 2318 : row_blk_size=row_blk_full, col_blk_size=blk_sub)
3486 : END IF
3487 :
3488 6224 : CALL dbcsr_iterator_start(iter, mat_full)
3489 61582 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3490 55358 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
3491 61582 : IF (compress_rows) THEN
3492 15085 : IF (inv_map(ib) > 0) CALL dbcsr_put_block(mat_out, inv_map(ib), jb, blk)
3493 : ELSE
3494 40273 : IF (inv_map(jb) > 0) CALL dbcsr_put_block(mat_out, ib, inv_map(jb), blk)
3495 : END IF
3496 : END DO
3497 6224 : CALL dbcsr_iterator_stop(iter)
3498 6224 : CALL dbcsr_finalize(mat_out)
3499 :
3500 6224 : CALL dbcsr_distribution_release(dist_sub)
3501 6224 : DEALLOCATE (inv_map, blk_sub, dist_sub_1d)
3502 :
3503 6224 : CALL timestop(handle)
3504 :
3505 12448 : END SUBROUTINE extract_masked_blocks
3506 :
3507 : ! **************************************************************************************************
3508 : !> \brief Pre-seeds a square blocked DBCSR matrix with zero blocks only for block pairs whose
3509 : !> centers lie within radius, for use with copy_fm_to_dbcsr(keep_sparsity=T) or
3510 : !> dbcsr_multiply(retain_sparsity=T). Consumers: the CUTOFF_RADIUS_G_W operator truncation
3511 : !> (atom-blocked, centers = atom_centers) and the RT-BSE CUTOFF_RADIUS_W0 truncation of the
3512 : !> grid-basis W^0 (grid-blocked, centers = chunk_centroids).
3513 : !> \param matrix ...
3514 : !> \param centers block positions, one column per block row/column of matrix
3515 : !> \param radius truncation radius, same units as centers (bohr)
3516 : ! **************************************************************************************************
3517 6 : SUBROUTINE reserve_blocks_within_radius(matrix, centers, radius)
3518 :
3519 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
3520 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centers
3521 : REAL(KIND=dp), INTENT(IN) :: radius
3522 :
3523 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reserve_blocks_within_radius'
3524 :
3525 : INTEGER :: handle, i, j, my_pcol, my_prow, &
3526 : nblkcols, nblkrows
3527 6 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist, row_blk, row_dist
3528 : REAL(KIND=dp) :: radius2
3529 6 : REAL(KIND=dp), ALLOCATABLE :: zero_blk(:, :)
3530 : TYPE(dbcsr_distribution_type) :: dist
3531 :
3532 6 : CALL timeset(routineN, handle)
3533 :
3534 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows, nblkcols_total=nblkcols, &
3535 6 : row_blk_size=row_blk, col_blk_size=col_blk, distribution=dist)
3536 : CALL dbcsr_distribution_get(dist, row_dist=row_dist, col_dist=col_dist, &
3537 6 : myprow=my_prow, mypcol=my_pcol)
3538 6 : CPASSERT(nblkrows == SIZE(centers, 2))
3539 6 : CPASSERT(nblkcols == SIZE(centers, 2))
3540 :
3541 6 : radius2 = radius**2
3542 144 : ALLOCATE (zero_blk(MAXVAL(row_blk(1:nblkrows)), MAXVAL(col_blk(1:nblkcols))))
3543 6 : zero_blk(:, :) = 0.0_dp
3544 :
3545 66 : DO i = 1, nblkrows
3546 60 : IF (row_dist(i) /= my_prow) CYCLE
3547 336 : DO j = 1, nblkcols
3548 300 : IF (col_dist(j) /= my_pcol) CYCLE
3549 300 : IF ((centers(1, i) - centers(1, j))**2 + (centers(2, i) - centers(2, j))**2 + &
3550 60 : (centers(3, i) - centers(3, j))**2 <= radius2) THEN
3551 168 : CALL dbcsr_put_block(matrix, i, j, zero_blk(1:row_blk(i), 1:col_blk(j)))
3552 : END IF
3553 : END DO
3554 : END DO
3555 6 : CALL dbcsr_finalize(matrix)
3556 :
3557 6 : DEALLOCATE (zero_blk)
3558 6 : CALL timestop(handle)
3559 :
3560 6 : END SUBROUTINE reserve_blocks_within_radius
3561 :
3562 : ! **************************************************************************************************
3563 : !> \brief Creates the (empty) result matrix of op(mat_left) * op(mat_right) with the correct block
3564 : !> structure and a distribution on the shared process grid, ready to be filled by
3565 : !> dbcsr_multiply. Row structure comes from op(left), column structure from op(right).
3566 : !> \param mat_left ...
3567 : !> \param mat_right ...
3568 : !> \param transa 'N' or 'T' applied to mat_left
3569 : !> \param transb 'N' or 'T' applied to mat_right
3570 : !> \param mat_out ...
3571 : ! **************************************************************************************************
3572 8286 : SUBROUTINE create_product_matrix(mat_left, mat_right, transa, transb, mat_out)
3573 :
3574 : TYPE(dbcsr_type), INTENT(IN) :: mat_left, mat_right
3575 : CHARACTER(LEN=1), INTENT(IN) :: transa, transb
3576 : TYPE(dbcsr_type), INTENT(OUT) :: mat_out
3577 :
3578 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_product_matrix'
3579 :
3580 : INTEGER :: handle, i, npcols, nprows
3581 8286 : INTEGER, DIMENSION(:), POINTER :: col_blk_l, col_blk_r, out_col_blk, &
3582 8286 : out_col_dist, out_row_blk, &
3583 8286 : out_row_dist, row_blk_l, row_blk_r
3584 : TYPE(dbcsr_distribution_type) :: dist_l, dist_out
3585 :
3586 8286 : CALL timeset(routineN, handle)
3587 :
3588 8286 : CALL dbcsr_get_info(mat_left, distribution=dist_l, row_blk_size=row_blk_l, col_blk_size=col_blk_l)
3589 8286 : CALL dbcsr_get_info(mat_right, row_blk_size=row_blk_r, col_blk_size=col_blk_r)
3590 8286 : CALL dbcsr_distribution_get(dist_l, nprows=nprows, npcols=npcols)
3591 :
3592 : ! block SIZES follow op(left)/op(right); DISTRIBUTIONS are freshly round-robined onto the
3593 : ! shared process grid (a transposed operand's row-dist is NOT a valid col-dist on a
3594 : ! non-square grid). dbcsr_multiply redistributes internally, so any valid mapping works.
3595 8286 : IF (transa == 'N') THEN
3596 8286 : out_row_blk => row_blk_l
3597 : ELSE
3598 0 : out_row_blk => col_blk_l
3599 : END IF
3600 8286 : IF (transb == 'N') THEN
3601 6012 : out_col_blk => col_blk_r
3602 : ELSE
3603 2274 : out_col_blk => row_blk_r
3604 : END IF
3605 :
3606 41430 : ALLOCATE (out_row_dist(SIZE(out_row_blk)), out_col_dist(SIZE(out_col_blk)))
3607 130224 : DO i = 1, SIZE(out_row_blk)
3608 130224 : out_row_dist(i) = MOD(i - 1, nprows)
3609 : END DO
3610 57344 : DO i = 1, SIZE(out_col_blk)
3611 57344 : out_col_dist(i) = MOD(i - 1, npcols)
3612 : END DO
3613 :
3614 : CALL dbcsr_distribution_new(dist_out, template=dist_l, &
3615 8286 : row_dist=out_row_dist, col_dist=out_col_dist)
3616 : CALL dbcsr_create(mat_out, name="panel_product", dist=dist_out, &
3617 : matrix_type=dbcsr_type_no_symmetry, &
3618 8286 : row_blk_size=out_row_blk, col_blk_size=out_col_blk)
3619 8286 : CALL dbcsr_distribution_release(dist_out)
3620 8286 : DEALLOCATE (out_row_dist, out_col_dist)
3621 :
3622 8286 : CALL timestop(handle)
3623 :
3624 24858 : END SUBROUTINE create_product_matrix
3625 :
3626 : ! **************************************************************************************************
3627 : !> \brief Builds the AO-space Green's function operator G^occ/vir_µν (AO x AO DBCSR)
3628 : !> \param bs_env ...
3629 : !> \param tau ...
3630 : !> \param ispin ...
3631 : !> \param occ ...
3632 : !> \param vir ...
3633 : !> \param template ...
3634 : !> \param matrix_G_ao ...
3635 : ! **************************************************************************************************
3636 5604 : SUBROUTINE build_G_ao(bs_env, tau, ispin, occ, vir, template, matrix_G_ao)
3637 :
3638 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
3639 : REAL(KIND=dp), INTENT(IN) :: tau
3640 : INTEGER, INTENT(IN) :: ispin
3641 : LOGICAL, INTENT(IN) :: occ, vir
3642 : TYPE(dbcsr_type), INTENT(INOUT) :: template
3643 : TYPE(dbcsr_type), INTENT(OUT) :: matrix_G_ao
3644 :
3645 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_G_ao'
3646 :
3647 : INTEGER :: handle
3648 2802 : INTEGER, DIMENSION(:), POINTER :: blk_ao, dist_row_ao
3649 : TYPE(cp_fm_type), POINTER :: fm_G
3650 : TYPE(dbcsr_distribution_type) :: dist_ao_ao
3651 :
3652 2802 : CALL timeset(routineN, handle)
3653 :
3654 2802 : IF (occ) THEN
3655 1422 : fm_G => bs_env%fm_Gocc
3656 : ELSE
3657 1380 : fm_G => bs_env%fm_Gvir
3658 : END IF
3659 :
3660 2802 : CALL G_occ_vir(bs_env, tau, fm_G, ispin, occ=occ, vir=vir)
3661 :
3662 2802 : CALL setup_square_topology(template, dist_ao_ao, blk_ao, dist_row_ao)
3663 : CALL dbcsr_create(matrix_G_ao, name="G_ao", dist=dist_ao_ao, &
3664 : matrix_type=dbcsr_type_no_symmetry, &
3665 2802 : row_blk_size=blk_ao, col_blk_size=blk_ao)
3666 :
3667 : ! Optional CUTOFF_RADIUS_G_W operator truncation: only atom-pair blocks within the radius
3668 : ! are reserved and filled.
3669 2802 : IF (bs_env%ri_rs%cutoff_radius_g_w > 0.0_dp .AND. &
3670 : ALLOCATED(bs_env%ri_rs%atom_centers)) THEN
3671 : CALL reserve_blocks_within_radius(matrix_G_ao, bs_env%ri_rs%atom_centers, &
3672 0 : bs_env%ri_rs%cutoff_radius_g_w)
3673 0 : CALL copy_fm_to_dbcsr(fm_G, matrix_G_ao, keep_sparsity=.TRUE.)
3674 : ELSE
3675 2802 : CALL copy_fm_to_dbcsr(fm_G, matrix_G_ao, keep_sparsity=.FALSE.)
3676 : END IF
3677 2802 : CALL dbcsr_filter(matrix_G_ao, bs_env%eps_filter)
3678 :
3679 : ! release only the topology; keep matrix_G_ao for the caller
3680 2802 : CALL release_square_topology(dist=dist_ao_ao, mapped_dist=dist_row_ao)
3681 :
3682 2802 : CALL timestop(handle)
3683 :
3684 2802 : END SUBROUTINE build_G_ao
3685 :
3686 : ! **************************************************************************************************
3687 : !> \brief Panel-streaming evaluation of out += scale * L_out^T (A_grid ∘ B_grid) L_out,
3688 : !> 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
3689 : !> grid x grid objects. The grid (row) index is processed in panels of ~panel_size rows; for
3690 : !> each panel only P x grid slabs are built, Hadamard-multiplied, and contracted into the
3691 : !> (small) output. Algebraically identical to L_out^T (A_grid ∘ B_grid) L_out summed over
3692 : !> grid rows, so the result matches the non-streamed path to eps_filter.
3693 : !>
3694 : !> Mapping (L in {phi (grid x AO), Z (grid x RI)}, M the AO/RI-space operator):
3695 : !> chi : L_A=L_B=phi, M_A=G_occ_ao, M_B=G_vir_ao, L_out=Z -> RI x RI
3696 : !> Sig : L_A=phi (M_A=D/G), L_B=Z (M_B=V/W), L_out=phi -> AO x AO
3697 : !> \param L_A ...
3698 : !> \param M_A ...
3699 : !> \param L_B ...
3700 : !> \param M_B ...
3701 : !> \param L_out ...
3702 : !> \param mat_out ...
3703 : !> \param scale ...
3704 : !> \param eps ...
3705 : !> \param para_env ...
3706 : !> \param pan_first ...
3707 : !> \param pan_last ...
3708 : !> \param lb_eq_la ...
3709 : !> \param lout_eq_la ...
3710 : !> \param zero_out ...
3711 : !> \param keep_sparsity ...
3712 : !> \param centroids ...
3713 : !> \param cutoff ...
3714 : !> \param grid_occupation ...
3715 : ! **************************************************************************************************
3716 632 : SUBROUTINE contract_grid_panels(L_A, M_A, L_B, M_B, L_out, mat_out, scale, eps, para_env, &
3717 632 : pan_first, pan_last, lb_eq_la, lout_eq_la, zero_out, &
3718 632 : keep_sparsity, centroids, cutoff, grid_occupation)
3719 :
3720 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_A
3721 : TYPE(dbcsr_type), INTENT(INOUT) :: M_A
3722 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_B
3723 : TYPE(dbcsr_type), INTENT(INOUT) :: M_B
3724 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_out
3725 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_out
3726 : REAL(KIND=dp), INTENT(IN) :: scale, eps
3727 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
3728 : INTEGER, DIMENSION(:), INTENT(IN) :: pan_first, pan_last
3729 : LOGICAL, INTENT(IN) :: lb_eq_la, lout_eq_la, zero_out
3730 : LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
3731 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
3732 : OPTIONAL :: centroids
3733 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
3734 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: grid_occupation
3735 :
3736 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_grid_panels'
3737 :
3738 : INTEGER :: blk0, blk1, handle, ipan, n_grid_total, &
3739 : ncols_pan, nrows_pan
3740 632 : INTEGER, ALLOCATABLE, DIMENSION(:) :: gmap
3741 : LOGICAL :: my_keep_sparsity, use_cutoff
3742 632 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used, usedA, usedB
3743 : TYPE(dbcsr_type) :: A_pan, B_pan, C_pan, LA_pan, LA_panC, &
3744 : LB_pan, LB_panC, Lout_pan, MA_sub, &
3745 : MB_sub, tmp2, tmpA, tmpB
3746 : TYPE(dbcsr_type), POINTER :: RB_A, RB_B, RB_out
3747 : TYPE(dbcsr_type), TARGET :: LA_near, LB_near, Lout_near
3748 :
3749 632 : CALL timeset(routineN, handle)
3750 :
3751 632 : my_keep_sparsity = .FALSE.
3752 632 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
3753 632 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
3754 632 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
3755 632 : IF (PRESENT(grid_occupation)) grid_occupation = 0.0_dp
3756 :
3757 632 : CALL dbcsr_get_info(L_A, nfullrows_total=n_grid_total)
3758 :
3759 632 : IF (zero_out) CALL dbcsr_set(mat_out, 0.0_dp)
3760 :
3761 1286 : DO ipan = 1, SIZE(pan_first)
3762 654 : blk0 = pan_first(ipan)
3763 654 : blk1 = pan_last(ipan)
3764 :
3765 : ! phi/Z panel slices (P x n)
3766 654 : CALL extract_grid_panel(L_A, blk0, blk1, LA_pan)
3767 654 : IF (.NOT. lb_eq_la) CALL extract_grid_panel(L_B, blk0, blk1, LB_pan)
3768 654 : IF (.NOT. lout_eq_la) CALL extract_grid_panel(L_out, blk0, blk1, Lout_pan)
3769 :
3770 : ! Which AO/RI atoms (column blocks) actually touch this panel: the inner index of
3771 : ! every multiply below is restricted to them, so only the matching rows of the
3772 : ! system-wide operators M_A/M_B ever enter Cannon (exact: dropped rows meet zeros).
3773 654 : CALL collect_used_col_blocks(LA_pan, para_env, usedA)
3774 654 : IF (.NOT. lb_eq_la) THEN
3775 44 : CALL collect_used_col_blocks(LB_pan, para_env, usedB)
3776 : ELSE
3777 610 : IF (ALLOCATED(usedB)) DEALLOCATE (usedB)
3778 3260 : ALLOCATE (usedB, SOURCE=usedA)
3779 : END IF
3780 654 : IF (.NOT. (ANY(usedA) .AND. ANY(usedB))) THEN
3781 : ! empty panel slice: its Hadamard contribution is exactly zero on all ranks
3782 0 : CALL dbcsr_release(LA_pan)
3783 0 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_pan)
3784 0 : IF (.NOT. lout_eq_la) CALL dbcsr_release(Lout_pan)
3785 : CYCLE
3786 : END IF
3787 :
3788 : ! Grid rows within reach of the panel: with the CUTOFF_RADIUS_RL_W truncation only
3789 : ! they can appear as columns of the panel products / inner rows of the L_out multiply,
3790 : ! so the system-wide phi/Z right operands are cut down to this neighborhood slice
3791 : ! (local extraction, zero communication; exact w.r.t. the geo template).
3792 654 : IF (use_cutoff) THEN
3793 44 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, grid_used)
3794 44 : CALL extract_masked_blocks(L_A, grid_used, LA_near, compress_rows=.TRUE., blk_map=gmap)
3795 44 : IF (.NOT. lb_eq_la) CALL extract_masked_blocks(L_B, grid_used, LB_near, compress_rows=.TRUE.)
3796 44 : IF (.NOT. lout_eq_la) CALL extract_masked_blocks(L_out, grid_used, Lout_near, compress_rows=.TRUE.)
3797 : RB_A => LA_near
3798 : ELSE
3799 : RB_A => L_A
3800 : END IF
3801 654 : IF (lb_eq_la) THEN
3802 : RB_B => RB_A
3803 44 : ELSE IF (use_cutoff) THEN
3804 : RB_B => LB_near
3805 : ELSE
3806 40 : RB_B => L_B
3807 : END IF
3808 654 : IF (lout_eq_la) THEN
3809 : RB_out => RB_A
3810 610 : ELSE IF (use_cutoff) THEN
3811 : RB_out => Lout_near
3812 : ELSE
3813 570 : RB_out => L_out
3814 : END IF
3815 :
3816 : ! A_pan = LA_pan * M_A * L_A^T (P x grid_near).
3817 : ! When cutoff is active, A_pan is pre-seeded with only nearby blocks via
3818 : ! build_geo_template_panel, and the multiply uses retain_sparsity to skip
3819 : ! computing distant blocks entirely (exact: they are zero by locality).
3820 654 : CALL extract_masked_blocks(LA_pan, usedA, LA_panC, compress_rows=.FALSE.)
3821 654 : CALL extract_masked_blocks(M_A, usedA, MA_sub, compress_rows=.TRUE.)
3822 654 : CALL create_product_matrix(LA_panC, MA_sub, 'N', 'N', tmpA)
3823 654 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LA_panC, MA_sub, 0.0_dp, tmpA, filter_eps=eps)
3824 654 : CALL dbcsr_release(MA_sub)
3825 654 : IF (use_cutoff) THEN
3826 : CALL build_geo_template_panel(LA_pan, LA_near, centroids, cutoff, blk0, A_pan, &
3827 44 : col_map=gmap)
3828 : ELSE
3829 610 : CALL create_product_matrix(tmpA, RB_A, 'N', 'T', A_pan)
3830 : END IF
3831 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpA, RB_A, 0.0_dp, A_pan, &
3832 654 : filter_eps=eps, retain_sparsity=use_cutoff)
3833 654 : CALL dbcsr_release(tmpA)
3834 :
3835 : ! Grid-basis occupation of A_pan = φ G φ^T, accumulated over ALL panels into the
3836 : ! occupation of the (never formed) full grid x grid object:
3837 : ! sum_panels nnz(A_pan) / n_grid^2, with nnz = occ * pan_rows * pan_cols.
3838 : ! Panel-independent by construction -- a single-panel sample would instead report the
3839 : ! local neighbor count of whichever region happens to land in that panel.
3840 654 : IF (PRESENT(grid_occupation)) THEN
3841 610 : CALL dbcsr_get_info(A_pan, nfullrows_total=nrows_pan, nfullcols_total=ncols_pan)
3842 : grid_occupation = grid_occupation + dbcsr_get_occupation(A_pan)* &
3843 : REAL(ncols_pan, dp)*REAL(nrows_pan, dp)/ &
3844 610 : (REAL(n_grid_total, dp)*REAL(n_grid_total, dp))
3845 : END IF
3846 :
3847 : ! B_pan = LB_pan * M_B * L_B^T (P x grid_near); reuse the L_A slices when L_B == L_A.
3848 : ! With keep_sparsity, B_pan is pre-populated with A_pan's block structure so that
3849 : ! retain_sparsity forces the final multiply to fill only those blocks (exact for ∘).
3850 654 : IF (lb_eq_la) THEN
3851 610 : CALL extract_masked_blocks(M_B, usedA, MB_sub, compress_rows=.TRUE.)
3852 610 : CALL create_product_matrix(LA_panC, MB_sub, 'N', 'N', tmpB)
3853 610 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LA_panC, MB_sub, 0.0_dp, tmpB, filter_eps=eps)
3854 : ELSE
3855 44 : CALL extract_masked_blocks(LB_pan, usedB, LB_panC, compress_rows=.FALSE.)
3856 44 : CALL extract_masked_blocks(M_B, usedB, MB_sub, compress_rows=.TRUE.)
3857 44 : CALL create_product_matrix(LB_panC, MB_sub, 'N', 'N', tmpB)
3858 44 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LB_panC, MB_sub, 0.0_dp, tmpB, filter_eps=eps)
3859 44 : CALL dbcsr_release(LB_panC)
3860 : END IF
3861 654 : CALL dbcsr_release(MB_sub)
3862 654 : IF (my_keep_sparsity) THEN
3863 654 : CALL dbcsr_create(B_pan, template=A_pan)
3864 654 : CALL dbcsr_copy(B_pan, A_pan)
3865 654 : CALL dbcsr_set(B_pan, 0.0_dp)
3866 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_B, 0.0_dp, B_pan, &
3867 654 : filter_eps=eps, retain_sparsity=.TRUE.)
3868 : ELSE
3869 0 : CALL create_product_matrix(tmpB, RB_B, 'N', 'T', B_pan)
3870 0 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_B, 0.0_dp, B_pan, filter_eps=eps)
3871 : END IF
3872 654 : CALL dbcsr_release(tmpB)
3873 654 : CALL dbcsr_release(LA_panC)
3874 :
3875 : ! C_pan = scale * (A_pan ∘ B_pan) (P x grid_near)
3876 654 : CALL dbcsr_create(C_pan, template=A_pan)
3877 654 : CALL hadamard_product(A_pan, B_pan, C_pan, scale)
3878 654 : CALL dbcsr_release(A_pan)
3879 654 : CALL dbcsr_release(B_pan)
3880 :
3881 : ! tmp2 = C_pan * L_out (P x n_out; inner index restricted to the neighborhood)
3882 654 : CALL create_product_matrix(C_pan, RB_out, 'N', 'N', tmp2)
3883 654 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_pan, RB_out, 0.0_dp, tmp2, filter_eps=eps)
3884 654 : CALL dbcsr_release(C_pan)
3885 :
3886 : ! mat_out += L_out_pan^T * tmp2 (accumulate: beta = 1)
3887 654 : IF (lout_eq_la) THEN
3888 44 : CALL dbcsr_multiply('T', 'N', 1.0_dp, LA_pan, tmp2, 1.0_dp, mat_out, filter_eps=eps)
3889 : ELSE
3890 610 : CALL dbcsr_multiply('T', 'N', 1.0_dp, Lout_pan, tmp2, 1.0_dp, mat_out, filter_eps=eps)
3891 610 : CALL dbcsr_release(Lout_pan)
3892 : END IF
3893 654 : CALL dbcsr_release(tmp2)
3894 654 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_pan)
3895 654 : CALL dbcsr_release(LA_pan)
3896 3248 : IF (use_cutoff) THEN
3897 44 : CALL dbcsr_release(LA_near)
3898 44 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_near)
3899 44 : IF (.NOT. lout_eq_la) CALL dbcsr_release(Lout_near)
3900 : END IF
3901 :
3902 : END DO
3903 :
3904 632 : CALL timestop(handle)
3905 :
3906 1264 : END SUBROUTINE contract_grid_panels
3907 :
3908 : ! **************************************************************************************************
3909 : !> \brief Σ^c-specific panel loop: computes both the occupied (neg) and virtual (pos) contributions
3910 : !> in a single pass over grid panels, forming W_pan = Z_panel × W_aux × Z^T only ONCE per
3911 : !> panel and reusing it for both the G^occ and G^vir Hadamard contractions.
3912 : !>
3913 : !> Computes:
3914 : !> mat_Sigma_neg = φ^T ( (φ G^occ φ^T) ∘ (Z W^MIC Z^T) ) φ
3915 : !> mat_Sigma_pos = φ^T ( (φ G^vir φ^T) ∘ (Z W^MIC Z^T) ) φ
3916 : !>
3917 : !> \param mat_phi ...
3918 : !> \param mat_Z ...
3919 : !> \param mat_G_occ_ao ...
3920 : !> \param mat_G_vir_ao ...
3921 : !> \param mat_W_aux ...
3922 : !> \param mat_Sigma_neg ...
3923 : !> \param mat_Sigma_pos ...
3924 : !> \param eps ...
3925 : !> \param para_env ...
3926 : !> \param pan_first ...
3927 : !> \param pan_last ...
3928 : !> \param keep_sparsity ...
3929 : !> \param centroids ...
3930 : !> \param cutoff ...
3931 : ! **************************************************************************************************
3932 790 : SUBROUTINE contract_grid_panels_sigma_c(mat_phi, mat_Z, mat_G_occ_ao, mat_G_vir_ao, &
3933 : mat_W_aux, mat_Sigma_neg, mat_Sigma_pos, eps, para_env, &
3934 790 : pan_first, pan_last, keep_sparsity, centroids, cutoff)
3935 :
3936 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: mat_phi, mat_Z
3937 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_G_occ_ao, mat_G_vir_ao, mat_W_aux, &
3938 : mat_Sigma_neg, mat_Sigma_pos
3939 : REAL(KIND=dp), INTENT(IN) :: eps
3940 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
3941 : INTEGER, DIMENSION(:), INTENT(IN) :: pan_first, pan_last
3942 : LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
3943 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
3944 : OPTIONAL :: centroids
3945 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
3946 :
3947 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_grid_panels_sigma_c'
3948 :
3949 : INTEGER :: blk0, blk1, handle, ipan
3950 790 : INTEGER, ALLOCATABLE, DIMENSION(:) :: gmap
3951 : LOGICAL :: my_keep_sparsity, use_cutoff
3952 790 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used, used_ao, used_ri
3953 : TYPE(dbcsr_type) :: A_occ, A_vir, C_pan, G_occ_sub, &
3954 : G_vir_sub, phi_pan, phi_panC, tmp2, &
3955 : tmpA, tmpB, W_pan, W_sub, Z_pan, Z_panC
3956 : TYPE(dbcsr_type), POINTER :: RB_phi, RB_Z
3957 : TYPE(dbcsr_type), TARGET :: phi_near, Z_near
3958 :
3959 790 : CALL timeset(routineN, handle)
3960 :
3961 790 : my_keep_sparsity = .FALSE.
3962 790 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
3963 790 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
3964 790 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
3965 :
3966 790 : CALL dbcsr_set(mat_Sigma_neg, 0.0_dp)
3967 790 : CALL dbcsr_set(mat_Sigma_pos, 0.0_dp)
3968 :
3969 1600 : DO ipan = 1, SIZE(pan_first)
3970 810 : blk0 = pan_first(ipan)
3971 810 : blk1 = pan_last(ipan)
3972 :
3973 810 : CALL extract_grid_panel(mat_phi, blk0, blk1, phi_pan)
3974 810 : CALL extract_grid_panel(mat_Z, blk0, blk1, Z_pan)
3975 :
3976 : ! AO/RI atoms touching this panel: only the matching rows of G_occ/G_vir/W ever
3977 : ! enter the multiplies below (exact: dropped rows meet zero columns of the panel).
3978 810 : CALL collect_used_col_blocks(phi_pan, para_env, used_ao)
3979 810 : CALL collect_used_col_blocks(Z_pan, para_env, used_ri)
3980 810 : IF (.NOT. (ANY(used_ao) .AND. ANY(used_ri))) THEN
3981 0 : CALL dbcsr_release(phi_pan)
3982 0 : CALL dbcsr_release(Z_pan)
3983 0 : CYCLE
3984 : END IF
3985 :
3986 : ! Neighborhood slices of phi/Z (grid rows within cutoff of the panel): they replace
3987 : ! the system-wide right operands in every multiply (local extraction, zero comm).
3988 810 : IF (use_cutoff) THEN
3989 40 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, grid_used)
3990 40 : CALL extract_masked_blocks(mat_phi, grid_used, phi_near, compress_rows=.TRUE., blk_map=gmap)
3991 40 : CALL extract_masked_blocks(mat_Z, grid_used, Z_near, compress_rows=.TRUE.)
3992 40 : RB_phi => phi_near
3993 40 : RB_Z => Z_near
3994 : ELSE
3995 : RB_phi => mat_phi
3996 : RB_Z => mat_Z
3997 : END IF
3998 :
3999 810 : CALL extract_masked_blocks(phi_pan, used_ao, phi_panC, compress_rows=.FALSE.)
4000 810 : CALL extract_masked_blocks(Z_pan, used_ri, Z_panC, compress_rows=.FALSE.)
4001 810 : CALL extract_masked_blocks(mat_G_occ_ao, used_ao, G_occ_sub, compress_rows=.TRUE.)
4002 810 : CALL extract_masked_blocks(mat_G_vir_ao, used_ao, G_vir_sub, compress_rows=.TRUE.)
4003 810 : CALL extract_masked_blocks(mat_W_aux, used_ri, W_sub, compress_rows=.TRUE.)
4004 :
4005 : ! A_occ = phi_pan × G_occ × phi^T (built first so W_pan can inherit its pattern).
4006 : ! With cutoff active, A_occ is pre-seeded with geo-local blocks so that the
4007 : ! phi^T multiply uses retain_sparsity and never computes distant blocks.
4008 810 : CALL create_product_matrix(phi_panC, G_occ_sub, 'N', 'N', tmpA)
4009 810 : CALL dbcsr_multiply('N', 'N', 1.0_dp, phi_panC, G_occ_sub, 0.0_dp, tmpA, filter_eps=eps)
4010 810 : IF (use_cutoff) THEN
4011 : CALL build_geo_template_panel(phi_pan, phi_near, centroids, cutoff, blk0, A_occ, &
4012 40 : col_map=gmap)
4013 : ELSE
4014 770 : CALL create_product_matrix(tmpA, RB_phi, 'N', 'T', A_occ)
4015 : END IF
4016 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpA, RB_phi, 0.0_dp, A_occ, &
4017 810 : filter_eps=eps, retain_sparsity=use_cutoff)
4018 810 : CALL dbcsr_release(tmpA)
4019 810 : CALL dbcsr_release(G_occ_sub)
4020 :
4021 : ! A_vir = phi_pan × G_vir × phi^T (same pre-screen as A_occ)
4022 810 : CALL create_product_matrix(phi_panC, G_vir_sub, 'N', 'N', tmpA)
4023 810 : CALL dbcsr_multiply('N', 'N', 1.0_dp, phi_panC, G_vir_sub, 0.0_dp, tmpA, filter_eps=eps)
4024 810 : IF (use_cutoff) THEN
4025 : CALL build_geo_template_panel(phi_pan, phi_near, centroids, cutoff, blk0, A_vir, &
4026 40 : col_map=gmap)
4027 : ELSE
4028 770 : CALL create_product_matrix(tmpA, RB_phi, 'N', 'T', A_vir)
4029 : END IF
4030 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpA, RB_phi, 0.0_dp, A_vir, &
4031 810 : filter_eps=eps, retain_sparsity=use_cutoff)
4032 810 : CALL dbcsr_release(tmpA)
4033 810 : CALL dbcsr_release(G_vir_sub)
4034 :
4035 : ! W_pan = Z_pan × W_aux × Z^T (computed once, reused for both Σ^c terms).
4036 : ! With keep_sparsity, W_pan is pre-seeded with the union of A_occ and A_vir block
4037 : ! patterns so that retain_sparsity forces the multiply to fill only those blocks:
4038 : ! exact since W outside G_occ∪G_vir is multiplied by zero in the Hadamard.
4039 810 : CALL create_product_matrix(Z_panC, W_sub, 'N', 'N', tmpB)
4040 810 : CALL dbcsr_multiply('N', 'N', 1.0_dp, Z_panC, W_sub, 0.0_dp, tmpB, filter_eps=eps)
4041 810 : IF (my_keep_sparsity) THEN
4042 810 : CALL dbcsr_create(W_pan, template=A_occ)
4043 810 : CALL dbcsr_copy(W_pan, A_occ)
4044 810 : CALL dbcsr_add(W_pan, A_vir, 1.0_dp, 1.0_dp)
4045 810 : CALL dbcsr_set(W_pan, 0.0_dp)
4046 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_Z, 0.0_dp, W_pan, &
4047 810 : filter_eps=eps, retain_sparsity=.TRUE.)
4048 : ELSE
4049 0 : CALL create_product_matrix(tmpB, RB_Z, 'N', 'T', W_pan)
4050 0 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_Z, 0.0_dp, W_pan, filter_eps=eps)
4051 : END IF
4052 810 : CALL dbcsr_release(tmpB)
4053 810 : CALL dbcsr_release(W_sub)
4054 810 : CALL dbcsr_release(phi_panC)
4055 810 : CALL dbcsr_release(Z_panC)
4056 :
4057 : ! Σ^c_neg: φ^T ( A_occ ∘ W_pan ) φ
4058 810 : CALL dbcsr_create(C_pan, template=A_occ)
4059 810 : CALL hadamard_product(A_occ, W_pan, C_pan, 1.0_dp)
4060 810 : CALL dbcsr_release(A_occ)
4061 810 : CALL create_product_matrix(C_pan, RB_phi, 'N', 'N', tmp2)
4062 810 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_pan, RB_phi, 0.0_dp, tmp2, filter_eps=eps)
4063 810 : CALL dbcsr_release(C_pan)
4064 810 : CALL dbcsr_multiply('T', 'N', 1.0_dp, phi_pan, tmp2, 1.0_dp, mat_Sigma_neg, filter_eps=eps)
4065 810 : CALL dbcsr_release(tmp2)
4066 :
4067 : ! Σ^c_pos: φ^T ( A_vir ∘ W_pan ) φ — W_pan reused
4068 810 : CALL dbcsr_create(C_pan, template=A_vir)
4069 810 : CALL hadamard_product(A_vir, W_pan, C_pan, 1.0_dp)
4070 810 : CALL dbcsr_release(A_vir)
4071 810 : CALL create_product_matrix(C_pan, RB_phi, 'N', 'N', tmp2)
4072 810 : CALL dbcsr_multiply('N', 'N', 1.0_dp, C_pan, RB_phi, 0.0_dp, tmp2, filter_eps=eps)
4073 810 : CALL dbcsr_release(C_pan)
4074 810 : CALL dbcsr_multiply('T', 'N', 1.0_dp, phi_pan, tmp2, 1.0_dp, mat_Sigma_pos, filter_eps=eps)
4075 810 : CALL dbcsr_release(tmp2)
4076 :
4077 810 : CALL dbcsr_release(W_pan)
4078 810 : CALL dbcsr_release(Z_pan)
4079 810 : CALL dbcsr_release(phi_pan)
4080 6460 : IF (use_cutoff) THEN
4081 40 : CALL dbcsr_release(phi_near)
4082 40 : CALL dbcsr_release(Z_near)
4083 : END IF
4084 :
4085 : END DO
4086 :
4087 790 : CALL timestop(handle)
4088 :
4089 1580 : END SUBROUTINE contract_grid_panels_sigma_c
4090 :
4091 : ! **************************************************************************************************
4092 : !> \brief Computes C = A ◦ B (Element-wise Hadamard product) for sparse DBCSR matrices.
4093 : !> \param matrix_A ...
4094 : !> \param matrix_B ...
4095 : !> \param matrix_C ...
4096 : !> \param fac (Scaling factor applied to the product)
4097 : ! **************************************************************************************************
4098 4548 : SUBROUTINE hadamard_product(matrix_A, matrix_B, matrix_C, fac)
4099 :
4100 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B, matrix_C
4101 : REAL(KIND=dp), INTENT(IN) :: fac
4102 :
4103 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hadamard_product'
4104 :
4105 : INTEGER :: col, handle, row
4106 : LOGICAL :: found
4107 2274 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk_B, blk_C
4108 : TYPE(dbcsr_iterator_type) :: iter
4109 :
4110 2274 : CALL timeset(routineN, handle)
4111 :
4112 2274 : CALL dbcsr_copy(matrix_C, matrix_A)
4113 :
4114 2274 : CALL dbcsr_iterator_start(iter, matrix_C)
4115 297836 : DO WHILE (dbcsr_iterator_blocks_left(iter))
4116 295562 : CALL dbcsr_iterator_next_block(iter, row, col, blk_C)
4117 :
4118 295562 : CALL dbcsr_get_block_p(matrix_B, row, col, blk_B, found)
4119 :
4120 297836 : IF (found) THEN
4121 527048742 : blk_C(:, :) = fac*blk_C(:, :)*blk_B(:, :)
4122 : ELSE
4123 : ! If B is sparse here, the product is zero
4124 0 : blk_C(:, :) = 0.0_dp
4125 : END IF
4126 : END DO
4127 2274 : CALL dbcsr_iterator_stop(iter)
4128 :
4129 2274 : CALL timestop(handle)
4130 :
4131 2274 : END SUBROUTINE hadamard_product
4132 :
4133 : ! **************************************************************************************************
4134 : !> \brief Computes the screened Coulomb interaction on the imaginary-time grid, entirely in the
4135 : !> RI auxiliary (PQ) basis:
4136 : !> χ_PQ(iω) = Σ_τ w(ω,τ) cos(ωτ) χ_PQ(iτ) (cosine transform)
4137 : !> ε(iω) = Id - V^0.5 M^-1 χ(iω) M^-1 V^0.5 (dielectric function)
4138 : !> W(iω) = V^0.5 ( ε^-1(iω) - Id ) V^0.5 (correlation part only)
4139 : !> W(iτ) = Σ_ω w̃(τ,ω) cos(ωτ) W(iω) (back transform)
4140 : !> W(iτ) <- M^-1 W(iτ) M^-1 (fold in the RI metric)
4141 : !> where V is the bare Coulomb matrix and M the RI metric.
4142 : !> \param bs_env ...
4143 : !> \param qs_env ...
4144 : !> \param mat_chi_Gamma_tau ...
4145 : !> \param fm_W_time ...
4146 : ! **************************************************************************************************
4147 144 : SUBROUTINE compute_W(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_time)
4148 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
4149 : TYPE(qs_environment_type), POINTER :: qs_env
4150 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
4151 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
4152 :
4153 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_W'
4154 :
4155 : INTEGER :: handle, i_t, j_w
4156 : REAL(KIND=dp) :: t1
4157 : TYPE(cp_fm_type) :: fm_M_inv_V_sqrt, fm_V, fm_V_sqrt
4158 :
4159 36 : CALL timeset(routineN, handle)
4160 :
4161 36 : t1 = m_walltime()
4162 :
4163 36 : CALL create_fm_W_MIC_time(bs_env, fm_W_time)
4164 :
4165 : ! 1. Allocate V and M matrices
4166 36 : CALL cp_fm_create(fm_V, bs_env%fm_RI_RI%matrix_struct)
4167 36 : CALL cp_fm_create(fm_V_sqrt, bs_env%fm_RI_RI%matrix_struct)
4168 36 : CALL cp_fm_create(fm_M_inv_V_sqrt, bs_env%fm_RI_RI%matrix_struct)
4169 :
4170 : ! Compute V and M^-1 * V^0.5
4171 36 : CALL compute_V_MinvVsqrt(bs_env, qs_env, fm_V, fm_V_sqrt, fm_M_inv_V_sqrt)
4172 :
4173 : ! 2. Loop over frequencies
4174 546 : DO j_w = 1, bs_env%num_time_freq_points
4175 : ! Fourier transformation of χ_PQ(iτ) to χ_PQ(iω_j)
4176 510 : CALL compute_fm_chi_Gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
4177 :
4178 : ! ε(iω_j) = Id - V^0.5*M^-1*χ(iω_j)*M^-1*V^0.5
4179 : ! W(iω_j) = V^0.5*(ε^-1(iω_j)-Id)*V^0.5
4180 : CALL compute_fm_W_freq(bs_env, bs_env%fm_chi_Gamma_freq, fm_V_sqrt, &
4181 510 : fm_M_inv_V_sqrt, bs_env%fm_W_MIC_freq)
4182 :
4183 : ! Fourier transform from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
4184 546 : CALL Fourier_transform_w_to_t(bs_env, fm_W_time, bs_env%fm_W_MIC_freq, j_w)
4185 : END DO
4186 :
4187 : ! M^-1*W^MIC(iτ)*M^-1
4188 36 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_W_time)
4189 :
4190 36 : IF (bs_env%unit_nr > 0) THEN
4191 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
4192 18 : 'Computed W(iτ),', ' Execution time', m_walltime() - t1, ' s'
4193 : END IF
4194 :
4195 36 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
4196 :
4197 : ! Cleanup
4198 36 : CALL cp_fm_release(fm_V)
4199 36 : CALL cp_fm_release(fm_V_sqrt)
4200 36 : CALL cp_fm_release(fm_M_inv_V_sqrt)
4201 :
4202 : ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
4203 36 : IF (bs_env%rtp_method == rtp_method_bse) THEN
4204 14 : t1 = m_walltime()
4205 14 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
4206 : ! Set to zero
4207 14 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
4208 : ! Sum over all times
4209 294 : DO i_t = 1, bs_env%num_time_freq_points
4210 : ! Add the relevant structure with correct weight
4211 : CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
4212 294 : bs_env%imag_time_weights_freq_zero(i_t), fm_W_time(i_t))
4213 : END DO
4214 : ! Done, save to file
4215 14 : CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
4216 : ! Report calculation
4217 14 : IF (bs_env%unit_nr > 0) THEN
4218 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
4219 7 : 'Computed W(0),', ' Execution time', m_walltime() - t1, ' s'
4220 : END IF
4221 : END IF
4222 :
4223 36 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
4224 :
4225 36 : CALL timestop(handle)
4226 :
4227 36 : END SUBROUTINE compute_W
4228 :
4229 : ! **************************************************************************************************
4230 : !> \brief Computes the static RI-basis Coulomb operators entering the dielectric function:
4231 : !> the bare Coulomb matrix V_PQ(k=0), its Cholesky/matrix square root V^0.5, and
4232 : !> M^-1 V^0.5 with M the RI metric (2c integrals of the RI_METRIC operator).
4233 : !> \param bs_env ...
4234 : !> \param qs_env ...
4235 : !> \param fm_V ...
4236 : !> \param fm_V_sqrt ...
4237 : !> \param fm_Minv_Vsqrt ...
4238 : ! **************************************************************************************************
4239 36 : SUBROUTINE compute_V_MinvVsqrt(bs_env, qs_env, fm_V, fm_V_sqrt, fm_Minv_Vsqrt)
4240 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
4241 : TYPE(qs_environment_type), POINTER :: qs_env
4242 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_V, fm_V_sqrt, fm_Minv_Vsqrt
4243 :
4244 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_MinvVsqrt'
4245 :
4246 : INTEGER :: handle, info, n_RI, ndep
4247 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
4248 : TYPE(cell_type), POINTER :: cell
4249 : TYPE(cp_fm_type) :: fm_work
4250 36 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_M
4251 36 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_kp
4252 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
4253 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
4254 :
4255 36 : CALL timeset(routineN, handle)
4256 :
4257 36 : n_RI = bs_env%n_RI
4258 36 : CALL cp_fm_create(fm_work, fm_V%matrix_struct)
4259 :
4260 : ! -----------------------------------------------------------------------
4261 : ! 1. Build Coulomb Matrix V(k=0) using the kp-routine but only for ikp=1
4262 : ! -----------------------------------------------------------------------
4263 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, &
4264 36 : qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
4265 :
4266 180 : ALLOCATE (mat_V_kp(1:1, 1:2))
4267 36 : NULLIFY (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
4268 36 : ALLOCATE (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
4269 :
4270 36 : CALL dbcsr_create(mat_V_kp(1, 1)%matrix, template=bs_env%mat_RI_RI%matrix)
4271 36 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 1)%matrix)
4272 36 : CALL dbcsr_set(mat_V_kp(1, 1)%matrix, 0.0_dp)
4273 :
4274 : ! Dummy imaginary part just to satisfy the routine
4275 36 : CALL dbcsr_create(mat_V_kp(1, 2)%matrix, template=bs_env%mat_RI_RI%matrix)
4276 36 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 2)%matrix)
4277 36 : CALL dbcsr_set(mat_V_kp(1, 2)%matrix, 0.0_dp)
4278 :
4279 144 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
4280 :
4281 : CALL build_2c_coulomb_matrix_kp(mat_V_kp, bs_env%kpoints_chi_eps_W, "RI_AUX", cell, &
4282 : particle_set, qs_kind_set, atomic_kind_set, &
4283 36 : bs_env%size_lattice_sum_V, operator_coulomb, 1, 1)
4284 :
4285 : ! Copy real part to fm_V
4286 36 : CALL copy_dbcsr_to_fm(mat_V_kp(1, 1)%matrix, fm_V)
4287 :
4288 36 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 1)%matrix)
4289 36 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 2)%matrix)
4290 36 : DEALLOCATE (mat_V_kp)
4291 :
4292 : ! -----------------------------------------------------------------------
4293 : ! 2. Get RI-Metric Matrix M(k=0)
4294 : ! -----------------------------------------------------------------------
4295 : CALL RI_2c_integral_mat(qs_env, fm_M, fm_V, n_RI, bs_env%ri_metric, &
4296 36 : do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI)
4297 :
4298 : ! -----------------------------------------------------------------------
4299 : ! 3. M -> M^-1
4300 : ! -----------------------------------------------------------------------
4301 36 : CALL cp_fm_cholesky_decompose(fm_M(1, 1), info_out=info)
4302 36 : IF (info == 0) THEN
4303 36 : CALL cp_fm_cholesky_invert(fm_M(1, 1))
4304 36 : CALL cp_fm_uplo_to_full(fm_M(1, 1), fm_work)
4305 : ELSE
4306 : ! Fallback if Cholesky fails due to conditioning
4307 0 : CALL cp_fm_power(fm_M(1, 1), fm_work, -1.0_dp, bs_env%eps_eigval_mat_RI, ndep)
4308 0 : CALL cp_fm_to_fm(fm_work, fm_M(1, 1))
4309 : END IF
4310 :
4311 : ! -----------------------------------------------------------------------
4312 : ! 4. V -> V^0.5
4313 : ! -----------------------------------------------------------------------
4314 36 : CALL cp_fm_to_fm(fm_V, fm_V_sqrt)
4315 36 : CALL cp_fm_cholesky_decompose(fm_V_sqrt, info_out=info)
4316 36 : IF (info == 0) THEN
4317 36 : CALL clean_lower_part(fm_V_sqrt)
4318 : ELSE
4319 0 : CALL cp_fm_power(fm_V, fm_V_sqrt, 0.5_dp, bs_env%eps_eigval_mat_RI, ndep)
4320 : END IF
4321 :
4322 : ! -----------------------------------------------------------------------
4323 : ! 5. M^-1 * V^0.5
4324 : ! -----------------------------------------------------------------------
4325 : CALL parallel_gemm("N", "T", n_RI, n_RI, n_RI, 1.0_dp, fm_M(1, 1), fm_V_sqrt, &
4326 36 : 0.0_dp, fm_Minv_Vsqrt)
4327 :
4328 36 : CALL cp_fm_release(fm_M)
4329 36 : CALL cp_fm_release(fm_work)
4330 :
4331 36 : CALL timestop(handle)
4332 :
4333 144 : END SUBROUTINE compute_V_MinvVsqrt
4334 :
4335 : ! **************************************************************************************************
4336 : !> \brief Computes the screened interaction at one imaginary frequency:
4337 : !> ε(iω_j) = Id - (M^-1 V^0.5)^T χ(iω_j) (M^-1 V^0.5)
4338 : !> W(iω_j) = V^0.5^T ( ε^-1(iω_j) - Id ) V^0.5
4339 : !> ε is inverted via Cholesky; if that fails due to conditioning, via
4340 : !> eigendecomposition (cp_fm_power) with eigenvalue filtering.
4341 : !> \param bs_env ...
4342 : !> \param fm_chi_freq_j ...
4343 : !> \param fm_V_sqrt ...
4344 : !> \param fm_Minv_Vsqrt ...
4345 : !> \param fm_W_freq_j ...
4346 : ! **************************************************************************************************
4347 2040 : SUBROUTINE compute_fm_W_freq(bs_env, fm_chi_freq_j, fm_V_sqrt, fm_Minv_Vsqrt, fm_W_freq_j)
4348 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
4349 : TYPE(cp_fm_type), INTENT(IN) :: fm_chi_freq_j, fm_V_sqrt, fm_Minv_Vsqrt
4350 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_W_freq_j
4351 :
4352 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_freq'
4353 :
4354 : INTEGER :: handle, info, n_RI, ndep
4355 : TYPE(cp_fm_type) :: fm_eps_freq_j, fm_work
4356 :
4357 510 : CALL timeset(routineN, handle)
4358 :
4359 510 : n_RI = bs_env%n_RI
4360 :
4361 510 : CALL cp_fm_create(fm_eps_freq_j, fm_chi_freq_j%matrix_struct)
4362 510 : CALL cp_fm_create(fm_work, fm_chi_freq_j%matrix_struct)
4363 :
4364 : ! -----------------------------------------------------------------------
4365 : ! 1. ε(iω_j) = Id - (M^-1 * V^0.5)^T * χ(iω_j) * (M^-1 * V^0.5)
4366 : ! -----------------------------------------------------------------------
4367 : ! work = χ(iω_j) * (M^-1 * V^0.5)
4368 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
4369 510 : fm_chi_freq_j, fm_Minv_Vsqrt, 0.0_dp, fm_work)
4370 :
4371 : ! eps_work = (M^-1 * V^0.5)^T * work
4372 : CALL parallel_gemm('T', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
4373 510 : fm_Minv_Vsqrt, fm_work, 0.0_dp, fm_eps_freq_j)
4374 :
4375 : ! ε(iω_j) = Id - eps_work --> -eps_work + Id
4376 510 : CALL fm_add_on_diag(fm_eps_freq_j, 1.0_dp)
4377 :
4378 : ! Force perfect symmetry before Cholesky to avoid info != 0 due to GEMM noise
4379 510 : CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
4380 :
4381 : ! -----------------------------------------------------------------------
4382 : ! 2. W(iω_j) = V^0.5^T * (ε^-1(iω_j) - Id) * V^0.5
4383 : ! -----------------------------------------------------------------------
4384 :
4385 : ! a) Cholesky decomposition of ε(iω_j)
4386 510 : CALL cp_fm_cholesky_decompose(fm_eps_freq_j, info_out=info)
4387 :
4388 : ! b) Inversion
4389 510 : IF (info == 0) THEN
4390 510 : CALL cp_fm_cholesky_invert(fm_eps_freq_j)
4391 510 : CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
4392 : ELSE
4393 : ! Fallback to expensive diagonalization if Cholesky fails
4394 0 : CALL cp_fm_power(fm_eps_freq_j, fm_work, -1.0_dp, bs_env%eps_eigval_mat_RI, ndep)
4395 0 : CALL cp_fm_to_fm(fm_work, fm_eps_freq_j)
4396 : END IF
4397 :
4398 : ! c) ε^-1(iω_j) - Id
4399 510 : CALL fm_add_on_diag(fm_eps_freq_j, -1.0_dp)
4400 :
4401 : ! d) work = (ε^-1(iω_j) - Id) * V^0.5
4402 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_eps_freq_j, fm_V_sqrt, &
4403 510 : 0.0_dp, fm_work)
4404 :
4405 : ! e) W(iw) = V^0.5^T * work
4406 : CALL parallel_gemm('T', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_V_sqrt, fm_work, &
4407 510 : 0.0_dp, fm_W_freq_j)
4408 :
4409 : ! Cleanup
4410 510 : CALL cp_fm_release(fm_work)
4411 510 : CALL cp_fm_release(fm_eps_freq_j)
4412 :
4413 510 : CALL timestop(handle)
4414 :
4415 510 : END SUBROUTINE compute_fm_W_freq
4416 :
4417 : ! **************************************************************************************************
4418 : !> \brief Adds a real scalar value to the diagonal of a real full matrix
4419 : !> \param fm ...
4420 : !> \param alpha ...
4421 : ! **************************************************************************************************
4422 2040 : SUBROUTINE fm_add_on_diag(fm, alpha)
4423 : TYPE(cp_fm_type), INTENT(INOUT) :: fm
4424 : REAL(KIND=dp), INTENT(IN) :: alpha
4425 :
4426 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_add_on_diag'
4427 :
4428 : INTEGER :: handle, i_global, i_row, j_col, &
4429 : j_global, ncol_local, nrow_local
4430 1020 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
4431 :
4432 1020 : CALL timeset(routineN, handle)
4433 :
4434 : CALL cp_fm_get_info(matrix=fm, &
4435 : nrow_local=nrow_local, &
4436 : ncol_local=ncol_local, &
4437 : row_indices=row_indices, &
4438 1020 : col_indices=col_indices)
4439 :
4440 36520 : DO j_col = 1, ncol_local
4441 35500 : j_global = col_indices(j_col)
4442 732150 : DO i_row = 1, nrow_local
4443 695630 : i_global = row_indices(i_row)
4444 731130 : IF (j_global == i_global) THEN
4445 17750 : fm%local_data(i_row, j_col) = fm%local_data(i_row, j_col) + alpha
4446 : END IF
4447 : END DO
4448 : END DO
4449 :
4450 1020 : CALL timestop(handle)
4451 :
4452 1020 : END SUBROUTINE fm_add_on_diag
4453 :
4454 : ! **************************************************************************************************
4455 : !> \brief Zeroes out the strictly lower triangular part of a real matrix
4456 : !> \param fm_mat ...
4457 : ! **************************************************************************************************
4458 72 : SUBROUTINE clean_lower_part(fm_mat)
4459 : TYPE(cp_fm_type) :: fm_mat
4460 :
4461 : CHARACTER(LEN=*), PARAMETER :: routineN = 'clean_lower_part'
4462 :
4463 : INTEGER :: handle, i_row, j_col, j_global, &
4464 : ncol_local, nrow_local
4465 36 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
4466 :
4467 36 : CALL timeset(routineN, handle)
4468 :
4469 : CALL cp_fm_get_info(matrix=fm_mat, &
4470 : nrow_local=nrow_local, ncol_local=ncol_local, &
4471 36 : row_indices=row_indices, col_indices=col_indices)
4472 :
4473 1146 : DO j_col = 1, ncol_local
4474 1110 : j_global = col_indices(j_col)
4475 20895 : DO i_row = 1, nrow_local
4476 20859 : IF (j_global < row_indices(i_row)) fm_mat%local_data(i_row, j_col) = 0.0_dp
4477 : END DO
4478 : END DO
4479 :
4480 36 : CALL timestop(handle)
4481 :
4482 36 : END SUBROUTINE clean_lower_part
4483 :
4484 : ! **************************************************************************************************
4485 : !> \brief Computes the exact-exchange part of the GW self-energy:
4486 : !> D_μν = Σ_n^occ C_μn C_νn (density matrix = G^occ at τ=0)
4487 : !> V^tr_PQ = M^-1 (P|Q)_trunc M^-1 (truncated Coulomb, RI basis)
4488 : !> Σ^x_λσ(k=0) = -Σ_ll' Φ_λ(r_l) [ (φ D φ^T)_ll' ∘ (Z V^tr Z^T)_ll' ] Φ_σ(r_l')
4489 : !> \param bs_env ...
4490 : !> \param qs_env ...
4491 : !> \param mat_phi_mu_l ...
4492 : !> \param mat_Z_lP ...
4493 : !> \param fm_Sigma_x_Gamma ...
4494 : ! **************************************************************************************************
4495 36 : SUBROUTINE compute_Sigma_x(bs_env, qs_env, mat_phi_mu_l, mat_Z_lP, fm_Sigma_x_Gamma)
4496 :
4497 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
4498 : TYPE(qs_environment_type), POINTER :: qs_env
4499 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
4500 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
4501 :
4502 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_x'
4503 :
4504 : INTEGER :: handle, ispin
4505 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
4506 36 : INTEGER, DIMENSION(:), POINTER :: blk_aux, dist_row_aux
4507 : REAL(KIND=dp) :: t1
4508 36 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
4509 : TYPE(dbcsr_distribution_type) :: dist_aux_aux
4510 : TYPE(dbcsr_type) :: mat_Sigma_x_Gamma, matrix_D_ao, &
4511 : matrix_V_aux
4512 :
4513 36 : CALL timeset(routineN, handle)
4514 :
4515 36 : t1 = m_walltime()
4516 :
4517 150 : ALLOCATE (fm_Sigma_x_Gamma(bs_env%n_spin))
4518 78 : DO ispin = 1, bs_env%n_spin
4519 78 : CALL cp_fm_create(fm_Sigma_x_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
4520 : END DO
4521 :
4522 36 : CALL dbcsr_create(mat_Sigma_x_Gamma, template=bs_env%mat_ao_ao%matrix)
4523 :
4524 36 : CALL resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
4525 :
4526 : ! =========================================================================
4527 : ! 1. COMPUTE V^tr_PQ (RI x RI)
4528 : ! =========================================================================
4529 36 : CALL setup_square_topology(mat_Z_lP, dist_aux_aux, blk_aux, dist_row_aux)
4530 :
4531 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
4532 36 : bs_env%trunc_coulomb, do_kpoints=.FALSE.)
4533 :
4534 : ! M^-1 * V^tr * M^-1 directly modifies fm_Vtr_Gamma(:, 1)
4535 36 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_Vtr_Gamma(:, 1))
4536 :
4537 36 : CALL dbcsr_create(matrix_V_aux, "V_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
4538 : ! Optional CUTOFF_RADIUS_G_W operator truncation + filter (see build_G_ao)
4539 36 : IF (bs_env%ri_rs%cutoff_radius_g_w > 0.0_dp .AND. &
4540 : ALLOCATED(bs_env%ri_rs%atom_centers)) THEN
4541 : CALL reserve_blocks_within_radius(matrix_V_aux, bs_env%ri_rs%atom_centers, &
4542 0 : bs_env%ri_rs%cutoff_radius_g_w)
4543 0 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.TRUE.)
4544 : ELSE
4545 36 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.FALSE.)
4546 : END IF
4547 36 : CALL dbcsr_filter(matrix_V_aux, bs_env%eps_filter)
4548 :
4549 : ! =========================================================================
4550 : ! 2. SPIN LOOP FOR EXACT EXCHANGE
4551 : ! Σ^x_λσ = -Σ_ll' Φ_λ(r_l) ( D_ll' V^tr_ll' ) Φ_σ(r_l')
4552 : ! = -φ^T ( (φ D φ^T) ∘ (Z V^tr Z^T) ) φ
4553 : ! =========================================================================
4554 78 : DO ispin = 1, bs_env%n_spin
4555 :
4556 : ! AO-space density matrix D_µν = G^occ at τ = 0
4557 42 : CALL build_G_ao(bs_env, 0.0_dp, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_D_ao)
4558 :
4559 : CALL contract_grid_panels(L_A=mat_phi_mu_l, M_A=matrix_D_ao, &
4560 : L_B=mat_Z_lP, M_B=matrix_V_aux, &
4561 : L_out=mat_phi_mu_l, mat_out=mat_Sigma_x_Gamma, &
4562 : scale=1.0_dp, eps=bs_env%eps_filter, &
4563 : para_env=bs_env%para_env, &
4564 : pan_first=pan_first, pan_last=pan_last, &
4565 : lb_eq_la=.FALSE., lout_eq_la=.TRUE., zero_out=.TRUE., &
4566 : keep_sparsity=bs_env%ri_rs%keep_sparsity_rirs, &
4567 : centroids=bs_env%ri_rs%chunk_centroids, &
4568 42 : cutoff=bs_env%ri_rs%cutoff_radius_v_w)
4569 42 : CALL dbcsr_scale(mat_Sigma_x_Gamma, -1.0_dp)
4570 :
4571 42 : CALL dbcsr_release(matrix_D_ao)
4572 :
4573 : ! Data I/O and Export to CP2K Full Matrices
4574 78 : CALL copy_dbcsr_to_fm(mat_Sigma_x_Gamma, fm_Sigma_x_Gamma(ispin))
4575 :
4576 : END DO ! ispin
4577 :
4578 36 : IF (bs_env%unit_nr > 0) THEN
4579 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
4580 18 : 'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
4581 18 : WRITE (bs_env%unit_nr, '(A)') ' '
4582 : END IF
4583 :
4584 : ! =========================================================================
4585 : ! 3. CLEANUP
4586 : ! =========================================================================
4587 36 : CALL dbcsr_release(matrix_V_aux)
4588 36 : CALL dbcsr_release(mat_Sigma_x_Gamma)
4589 36 : CALL release_square_topology(dist=dist_aux_aux, mapped_dist=dist_row_aux)
4590 :
4591 36 : CALL cp_fm_release(fm_Vtr_Gamma)
4592 :
4593 36 : CALL timestop(handle)
4594 :
4595 108 : END SUBROUTINE compute_Sigma_x
4596 :
4597 : ! **************************************************************************************************
4598 : !> \brief Computes the correlation part of the GW self-energy on the imaginary-time grid:
4599 : !> Σ^c_λσ(iτ<0) = -Σ_ll' Φ_λ(r_l) [ (φ G^occ φ^T)_ll' ∘ (Z W^MIC Z^T)_ll' ] Φ_σ(r_l')
4600 : !> Σ^c_λσ(iτ>0) = +Σ_ll' Φ_λ(r_l) [ (φ G^vir φ^T)_ll' ∘ (Z W^MIC Z^T)_ll' ] Φ_σ(r_l')
4601 : !> \param bs_env ...
4602 : !> \param fm_W_time ...
4603 : !> \param mat_phi_mu_l ...
4604 : !> \param mat_Z_lP ...
4605 : !> \param fm_Sigma_c_Gamma_time ...
4606 : ! **************************************************************************************************
4607 50 : SUBROUTINE compute_Sigma_c(bs_env, fm_W_time, mat_phi_mu_l, mat_Z_lP, fm_Sigma_c_Gamma_time)
4608 :
4609 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
4610 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
4611 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
4612 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
4613 :
4614 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c'
4615 :
4616 : INTEGER :: handle, i_t, ispin
4617 50 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pan_first, pan_last
4618 50 : INTEGER, DIMENSION(:), POINTER :: blk_aux, dist_row_aux
4619 : REAL(KIND=dp) :: t1, tau
4620 : TYPE(dbcsr_distribution_type) :: dist_aux_aux
4621 50 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
4622 : TYPE(dbcsr_type) :: matrix_G_occ_ao, matrix_G_vir_ao, &
4623 : matrix_W_aux
4624 :
4625 50 : CALL timeset(routineN, handle)
4626 :
4627 : ! =========================================================================
4628 : ! 1. SETUP AUXILIARY TOPOLOGY AND PRE-ALLOCATE OUTPUT ARRAYS
4629 : ! =========================================================================
4630 50 : CALL setup_square_topology(mat_Z_lP, dist_aux_aux, blk_aux, dist_row_aux)
4631 :
4632 50 : CALL resolve_grid_panels(bs_env, mat_phi_mu_l, pan_first, pan_last)
4633 :
4634 : ! Pre-allocate local DBCSR matrices to act as targets for final output
4635 50 : NULLIFY (mat_Sigma_neg_tau, mat_Sigma_pos_tau)
4636 1052 : ALLOCATE (mat_Sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
4637 1052 : ALLOCATE (mat_Sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
4638 :
4639 700 : DO i_t = 1, bs_env%num_time_freq_points
4640 1490 : DO ispin = 1, bs_env%n_spin
4641 790 : ALLOCATE (mat_Sigma_neg_tau(i_t, ispin)%matrix)
4642 790 : ALLOCATE (mat_Sigma_pos_tau(i_t, ispin)%matrix)
4643 790 : CALL dbcsr_create(mat_Sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
4644 1440 : CALL dbcsr_create(mat_Sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
4645 : END DO
4646 : END DO
4647 :
4648 : ! =========================================================================
4649 : ! 2. IMAGINARY TIME LOOP
4650 : ! Σ^c_neg_λσ(iτ) = -φ^T ( (φ G^occ φ^T) ∘ (Z W^MIC Z^T) ) φ
4651 : ! Σ^c_pos_λσ(iτ) = φ^T ( (φ G^vir φ^T) ∘ (Z W^MIC Z^T) ) φ
4652 : ! =========================================================================
4653 700 : DO i_t = 1, bs_env%num_time_freq_points
4654 650 : tau = bs_env%imag_time_points(i_t)
4655 :
4656 : CALL dbcsr_create(matrix_W_aux, "W_aux", dist_aux_aux, dbcsr_type_no_symmetry, &
4657 650 : blk_aux, blk_aux)
4658 650 : IF (bs_env%ri_rs%cutoff_radius_g_w > 0.0_dp .AND. &
4659 : ALLOCATED(bs_env%ri_rs%atom_centers)) THEN
4660 : CALL reserve_blocks_within_radius(matrix_W_aux, bs_env%ri_rs%atom_centers, &
4661 0 : bs_env%ri_rs%cutoff_radius_g_w)
4662 0 : CALL copy_fm_to_dbcsr(fm_W_time(i_t), matrix_W_aux, keep_sparsity=.TRUE.)
4663 : ELSE
4664 650 : CALL copy_fm_to_dbcsr(fm_W_time(i_t), matrix_W_aux, keep_sparsity=.FALSE.)
4665 : END IF
4666 650 : CALL dbcsr_filter(matrix_W_aux, bs_env%eps_filter)
4667 :
4668 1440 : DO ispin = 1, bs_env%n_spin
4669 790 : t1 = m_walltime()
4670 :
4671 : ! AO-space Green's functions G^occ_µν, G^vir_µν (dense AO x AO, small)
4672 790 : CALL build_G_ao(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_G_occ_ao)
4673 790 : CALL build_G_ao(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, matrix_G_vir_ao)
4674 :
4675 : ! Σ^c_neg and Σ^c_pos in a single panel loop: W_pan = Z_panel × W × Z^T built once
4676 : CALL contract_grid_panels_sigma_c(mat_phi=mat_phi_mu_l, mat_Z=mat_Z_lP, &
4677 : mat_G_occ_ao=matrix_G_occ_ao, &
4678 : mat_G_vir_ao=matrix_G_vir_ao, &
4679 : mat_W_aux=matrix_W_aux, &
4680 : mat_Sigma_neg=mat_Sigma_neg_tau(i_t, ispin)%matrix, &
4681 : mat_Sigma_pos=mat_Sigma_pos_tau(i_t, ispin)%matrix, &
4682 : eps=bs_env%eps_filter, &
4683 : para_env=bs_env%para_env, &
4684 : pan_first=pan_first, pan_last=pan_last, &
4685 : keep_sparsity=bs_env%ri_rs%keep_sparsity_rirs, &
4686 : centroids=bs_env%ri_rs%chunk_centroids, &
4687 790 : cutoff=bs_env%ri_rs%cutoff_radius_v_w)
4688 790 : CALL dbcsr_scale(mat_Sigma_neg_tau(i_t, ispin)%matrix, -1.0_dp)
4689 :
4690 790 : CALL dbcsr_release(matrix_G_occ_ao)
4691 790 : CALL dbcsr_release(matrix_G_vir_ao)
4692 :
4693 1440 : IF (bs_env%unit_nr > 0) THEN
4694 : WRITE (bs_env%unit_nr, '(T2,A,I15,A,I3,A,F7.1,A)') &
4695 395 : 'Computed Σ^c(iτ) for time point', i_t, ' /', bs_env%num_time_freq_points, &
4696 790 : ', Execution time', m_walltime() - t1, ' s'
4697 : END IF
4698 :
4699 : END DO ! ispin
4700 :
4701 700 : CALL dbcsr_release(matrix_W_aux)
4702 :
4703 : END DO ! i_t
4704 :
4705 50 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
4706 :
4707 : ! -------------------------------------------------------------------------
4708 : ! 3. FINALIZE AND CLEANUP
4709 : ! -------------------------------------------------------------------------
4710 : CALL fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
4711 50 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
4712 :
4713 : ! fm_W_time and the scratch files are released by the caller: in an evGW0 cycle this
4714 : ! routine is entered once per iteration and both have to survive until it is done.
4715 50 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_neg_tau)
4716 50 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_pos_tau)
4717 :
4718 50 : CALL release_square_topology(dist=dist_aux_aux, mapped_dist=dist_row_aux)
4719 :
4720 50 : CALL timestop(handle)
4721 :
4722 100 : END SUBROUTINE compute_Sigma_c
4723 :
4724 : ! **************************************************************************************************
4725 : !> \brief Builds the DBCSR distribution.
4726 : !> \param matrix_template ...
4727 : !> \param square_dist ...
4728 : !> \param blk_sizes ...
4729 : !> \param mapped_dist ...
4730 : ! **************************************************************************************************
4731 2888 : SUBROUTINE setup_square_topology(matrix_template, square_dist, blk_sizes, mapped_dist)
4732 :
4733 : TYPE(dbcsr_type), INTENT(IN) :: matrix_template
4734 : TYPE(dbcsr_distribution_type), INTENT(OUT) :: square_dist
4735 : INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: blk_sizes, mapped_dist
4736 :
4737 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_square_topology'
4738 :
4739 : INTEGER :: handle, i, nprows
4740 2888 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist
4741 : TYPE(dbcsr_distribution_type) :: dist_template
4742 :
4743 2888 : CALL timeset(routineN, handle)
4744 :
4745 2888 : CALL dbcsr_get_info(matrix_template, distribution=dist_template, col_blk_size=col_blk)
4746 2888 : CALL dbcsr_distribution_get(dist_template, col_dist=col_dist, nprows=nprows)
4747 :
4748 2888 : blk_sizes => col_blk
4749 8664 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
4750 9646 : DO i = 1, SIZE(blk_sizes)
4751 9646 : mapped_dist(i) = MOD(i - 1, nprows)
4752 : END DO
4753 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
4754 2888 : row_dist=mapped_dist, col_dist=col_dist)
4755 :
4756 2888 : CALL timestop(handle)
4757 :
4758 5776 : END SUBROUTINE setup_square_topology
4759 :
4760 : ! **************************************************************************************************
4761 : !> \brief Releases a distribution created by setup_square_topology.
4762 : !> \param dist ...
4763 : !> \param mapped_dist ...
4764 : ! **************************************************************************************************
4765 2888 : SUBROUTINE release_square_topology(dist, mapped_dist)
4766 :
4767 : TYPE(dbcsr_distribution_type), INTENT(INOUT) :: dist
4768 : INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: mapped_dist
4769 :
4770 2888 : CALL dbcsr_distribution_release(dist)
4771 2888 : IF (ASSOCIATED(mapped_dist)) THEN
4772 2888 : DEALLOCATE (mapped_dist)
4773 : NULLIFY (mapped_dist)
4774 : END IF
4775 :
4776 2888 : END SUBROUTINE release_square_topology
4777 :
4778 : ! **************************************************************************************************
4779 : !> \brief Σ^c_λσ(iτ) -> Σ^c_nn(ϵ) and the quasi-particle levels of the non-periodic RI-RS path,
4780 : !> ϵ_n^GW = ϵ_n^DFT + Σ^c_nn(ϵ_n^GW) + Σ^x_nn - v^xc_nn.
4781 : !> \param bs_env ...
4782 : !> \param fm_Sigma_x_Gamma ...
4783 : !> \param fm_Sigma_c_Gamma_time ...
4784 : ! **************************************************************************************************
4785 50 : SUBROUTINE compute_QP_energies(bs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
4786 :
4787 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
4788 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
4789 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
4790 :
4791 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
4792 :
4793 : INTEGER :: handle, ispin, j_t
4794 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Sigma_x_n, V_xc_n
4795 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Sigma_c_n_freq, Sigma_c_n_time
4796 : TYPE(cp_fm_type) :: fm_ks, fm_mos, fm_s, fm_work
4797 :
4798 50 : CALL timeset(routineN, handle)
4799 :
4800 50 : CALL cp_fm_create(fm_ks, bs_env%fm_s_Gamma%matrix_struct)
4801 50 : CALL cp_fm_create(fm_s, bs_env%fm_s_Gamma%matrix_struct)
4802 50 : CALL cp_fm_create(fm_mos, bs_env%fm_s_Gamma%matrix_struct)
4803 50 : CALL cp_fm_create(fm_work, bs_env%fm_s_Gamma%matrix_struct)
4804 :
4805 200 : ALLOCATE (V_xc_n(bs_env%n_ao), Sigma_x_n(bs_env%n_ao))
4806 250 : ALLOCATE (Sigma_c_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
4807 150 : ALLOCATE (Sigma_c_n_freq(bs_env%n_ao, bs_env%num_time_freq_points, 2))
4808 :
4809 112 : DO ispin = 1, bs_env%n_spin
4810 :
4811 : ! 1. Roothaan-Hall H^KS_µν C_νn = S_µν C_νn ϵ_n
4812 62 : CALL cp_fm_to_fm(bs_env%fm_ks_Gamma(ispin), fm_ks)
4813 62 : CALL cp_fm_to_fm(bs_env%fm_s_Gamma, fm_s)
4814 62 : CALL cp_fm_geeig(fm_ks, fm_s, fm_mos, bs_env%eigenval_scf(:, 1, ispin), fm_work)
4815 :
4816 : ! 2. v^xc_µν -> v^xc_nn and Σ^x_µν -> Σ^x_nn
4817 62 : CALL to_Gamma_and_mo_real(V_xc_n, bs_env%fm_V_xc_Gamma(ispin), bs_env, fm_mos)
4818 62 : CALL to_Gamma_and_mo_real(Sigma_x_n, fm_Sigma_x_Gamma(ispin), bs_env, fm_mos)
4819 :
4820 : ! 3. Σ^c_µν(+/-i|τ_j|) -> Σ^c_nn(+/-i|τ_j|)
4821 852 : DO j_t = 1, bs_env%num_time_freq_points
4822 : CALL to_Gamma_and_mo_real(Sigma_c_n_time(:, j_t, 1), &
4823 790 : fm_Sigma_c_Gamma_time(j_t, 1, ispin), bs_env, fm_mos)
4824 : CALL to_Gamma_and_mo_real(Sigma_c_n_time(:, j_t, 2), &
4825 852 : fm_Sigma_c_Gamma_time(j_t, 2, ispin), bs_env, fm_mos)
4826 : END DO
4827 :
4828 : ! 4. Σ^c_nn(iτ) -> Σ^c_nn(iω)
4829 62 : CALL time_to_freq(bs_env, Sigma_c_n_time, Sigma_c_n_freq, ispin)
4830 :
4831 : ! 5. Analytic continuation Σ^c_nn(iω) -> Σ^c_nn(ϵ) and the QP levels
4832 : CALL analyt_conti_and_print(bs_env, Sigma_c_n_freq, Sigma_x_n, V_xc_n, &
4833 112 : bs_env%eigenval_scf(:, 1, ispin), 1, ispin)
4834 :
4835 : END DO ! ispin
4836 :
4837 50 : CALL get_all_VBM_CBM_bandgaps(bs_env)
4838 :
4839 50 : IF (bs_env%gw_flavour == G0W0) CALL cp_fm_release(fm_Sigma_x_Gamma)
4840 50 : CALL cp_fm_release(fm_Sigma_c_Gamma_time)
4841 :
4842 50 : CALL cp_fm_release(fm_ks)
4843 50 : CALL cp_fm_release(fm_s)
4844 50 : CALL cp_fm_release(fm_mos)
4845 50 : CALL cp_fm_release(fm_work)
4846 :
4847 50 : CALL timestop(handle)
4848 :
4849 100 : END SUBROUTINE compute_QP_energies
4850 :
4851 : ! **************************************************************************************************
4852 : !> \brief AO -> MO transform of a Γ-point matrix
4853 : !> \param array_n ...
4854 : !> \param fm_Gamma ...
4855 : !> \param bs_env ...
4856 : !> \param fm_mos ...
4857 : ! **************************************************************************************************
4858 1704 : SUBROUTINE to_Gamma_and_mo_real(array_n, fm_Gamma, bs_env, fm_mos)
4859 :
4860 : REAL(KIND=dp), DIMENSION(:) :: array_n
4861 : TYPE(cp_fm_type) :: fm_Gamma
4862 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
4863 : TYPE(cp_fm_type) :: fm_mos
4864 :
4865 : CHARACTER(LEN=*), PARAMETER :: routineN = 'to_Gamma_and_mo_real'
4866 :
4867 : INTEGER :: handle, nmo
4868 : TYPE(cp_fm_type) :: fm_mo, fm_tmp
4869 :
4870 1704 : CALL timeset(routineN, handle)
4871 :
4872 1704 : CALL cp_fm_create(fm_tmp, fm_Gamma%matrix_struct)
4873 1704 : CALL cp_fm_create(fm_mo, fm_Gamma%matrix_struct)
4874 :
4875 1704 : nmo = bs_env%n_ao
4876 1704 : CALL parallel_gemm('N', 'N', nmo, nmo, nmo, 1.0_dp, fm_Gamma, fm_mos, 0.0_dp, fm_tmp)
4877 1704 : CALL parallel_gemm('T', 'N', nmo, nmo, nmo, 1.0_dp, fm_mos, fm_tmp, 0.0_dp, fm_mo)
4878 :
4879 1704 : CALL cp_fm_get_diag(fm_mo, array_n)
4880 :
4881 1704 : CALL cp_fm_release(fm_mo)
4882 1704 : CALL cp_fm_release(fm_tmp)
4883 :
4884 1704 : CALL timestop(handle)
4885 :
4886 1704 : END SUBROUTINE to_Gamma_and_mo_real
4887 :
4888 : END MODULE gw_non_periodic_ri_rs
|
