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