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