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 Computes the RI-RS fitting matrix Z_lP.
10 : !> \par History
11 : !> 09.2026 created Jan Wilhelm
12 : !> 09.2026 moved code by Ritaj Tyagi from gw_ri_rs_non_periodic.F
13 : !> 09.2026 added routine to compute Z_lP for automatic RI optimization
14 : ! **************************************************************************************************
15 : MODULE gw_ri_rs_compute_Z_lP
16 : USE ai_contraction_sphi, ONLY: abc_contract_xsmm
17 : USE atomic_kind_types, ONLY: atomic_kind_type
18 : USE basis_set_types, ONLY: get_gto_basis_set
19 : USE cell_types, ONLY: cell_type,&
20 : pbc
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_binary_read, dbcsr_binary_write, dbcsr_create, dbcsr_distribution_get, &
26 : dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, &
27 : dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_multiply, dbcsr_put_block, &
28 : dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
29 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
30 : USE cp_fm_struct, ONLY: cp_fm_struct_type
31 : USE cp_fm_types, ONLY: cp_fm_type
32 : USE cp_log_handling, ONLY: cp_get_default_logger,&
33 : cp_logger_type
34 : USE cp_output_handling, ONLY: cp_p_file,&
35 : cp_print_key_should_output
36 : USE gw_auto_ri_types, ONLY: auto_ri_type
37 : USE gw_ri_rs_compute_Z_lP_utils, ONLY: build_gram_jacobi_blas,&
38 : build_jacobi_diag_from_phi,&
39 : scale_rows_by_diag,&
40 : solve_D_lp_distributed,&
41 : store_Z_lP_columns
42 : USE gw_ri_rs_utils, ONLY: evaluate_ao_on_points
43 : USE gw_utils_compute_integrals, ONLY: build_3c_integral_block_ctx,&
44 : gw_3c_ctx_create,&
45 : gw_3c_ctx_release,&
46 : gw_3c_ctx_type,&
47 : gw_3c_ws_create,&
48 : gw_3c_ws_release,&
49 : gw_3c_ws_type
50 : USE input_section_types, ONLY: section_vals_type
51 : USE kinds, ONLY: dp
52 : USE libint_2c_3c, ONLY: eri_3center
53 : USE machine, ONLY: m_flush,&
54 : m_walltime
55 : USE message_passing, ONLY: mp_mem_avail_per_rank_GB,&
56 : mp_para_env_type
57 : USE orbital_pointers, ONLY: ncoset
58 : USE particle_types, ONLY: particle_type
59 : USE physcon, ONLY: angstrom
60 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type
63 : USE qs_kind_types, ONLY: qs_kind_type
64 : USE util, ONLY: sort
65 : #include "./base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 : PRIVATE
69 :
70 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_ri_rs_compute_Z_lP'
71 :
72 : PUBLIC :: compute_Z_lP
73 :
74 : CONTAINS
75 :
76 : ! **************************************************************************************************
77 : !> \brief Computes Z_lP for (1) a tabulated RI basis set or (2) an on-the-fly generated RI basis set.
78 : !>
79 : !> In both cases, the Z_lP coefficients are computed by solving the linear system
80 : !>
81 : !> Σ_l' D_ll' Z_l'P = d_lP,
82 : !>
83 : !> with D_ll' and d_lP given by
84 : !>
85 : !> D_ll' = [Σ_μ ϕ_μ(r_l) ϕ_μ(r_l')]²,
86 : !>
87 : !> d_lP = Σ_μν ϕ_μ(r_l) ϕ_ν(r_l) (μν|P).
88 : !>
89 : !> For an on-the-fly generated RI basis set, an RI function φ_P may be a contraction of
90 : !> Gaussians on neighboring atoms, which requires special computation of d_lP.
91 : !>
92 : !> \param qs_env ...
93 : !> \param bs_env Band-structure environment containing GW parameters.
94 : !> \param ri_rs_grid_points ...
95 : !> \param mat_phi_mu_l ...
96 : !> \param mat_Z_lP ...
97 : ! **************************************************************************************************
98 48 : SUBROUTINE compute_Z_lP(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
99 : TYPE(qs_environment_type), POINTER :: qs_env
100 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
101 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
102 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
103 : TYPE(dbcsr_type), INTENT(OUT) :: mat_Z_lP
104 :
105 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Z_lP'
106 :
107 : INTEGER :: handle
108 :
109 48 : CALL timeset(routineN, handle)
110 :
111 48 : IF (bs_env%auto_ri%enabled) THEN
112 4 : CALL compute_Z_lP_auto_ri(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
113 : ELSE
114 44 : CALL compute_Z_lP_standard(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
115 : END IF
116 :
117 48 : CALL timestop(handle)
118 :
119 48 : END SUBROUTINE compute_Z_lP
120 :
121 : ! **************************************************************************************************
122 : !> \brief Computes the RI-RS fitting coefficients Z_lP by solving, independently for every RI
123 : !> atom P, a Jacobi-conditioned, Tikhonov-regularized linear system restricted to the
124 : !> grid points r_l inside P's integration sphere |r_l - R_P| <= cutoff_ri(P):
125 : !>
126 : !> 1. D_ll' = [ Σ_μ ϕ_μ(r_l) ϕ_μ(r_l') ]²
127 : !> 2. d_l = 1 / sqrt(D_ll) (Jacobi conditioning vector)
128 : !> 3. D'_ll' = d_l D_ll' d_l' + λ δ_ll' (λ = TIKHONOV regularization)
129 : !> 4. d_lP = Σ_μν ϕ_μ(r_l) ϕ_ν(r_l) (μν|P)
130 : !> 5. Σ_l' D'_ll' Z'_l'P = d_l d_lP (solve linear system for Z'_lP)
131 : !> 6. Z_lP = d_l Z'_l'P (undo the conditioning)
132 : !>
133 : !> Work is distributed over atoms in two phases (planned by classify_z_lp_atoms and
134 : !> lpt_assign_atoms): Phase A solves "small" atoms with single-rank LAPACK
135 : !> (dpotrf/dpotrs); Phase B solves "big" atoms, whose dense matrix D'_ll' would exceed one
136 : !> rank's memory, with ScaLAPACK (pdpotrf/pdpotrs) over rank subgroups of size G.
137 : !> The solved Z columns are scattered into the sparse global mat_Z_lP.
138 : !> If a Z_lP restart file exists, it is read instead and the solve is skipped entirely.
139 : !> \param qs_env ...
140 : !> \param bs_env ...
141 : !> \param ri_rs_grid_points ...
142 : !> \param mat_phi_mu_l ...
143 : !> \param mat_Z_lP ...
144 : ! **************************************************************************************************
145 88 : SUBROUTINE compute_Z_lP_standard(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
146 :
147 : TYPE(qs_environment_type), POINTER :: qs_env
148 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
149 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
150 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
151 : TYPE(dbcsr_type), INTENT(OUT) :: mat_Z_lP
152 :
153 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
154 : routineN = 'compute_Z_lP_standard'
155 :
156 : INTEGER :: atom_j_mepos, atom_j_stride, atom_P, G, handle, handle_dpotrf, handle_dpotrs, &
157 : i_blk, iatom, idx, info, iphase, j, max_ao_size, my_group, n_ao_total, n_ao_used, n_big, &
158 : n_done, n_groups, n_loc_ri, n_local_grid, n_my_atoms, n_small, natom, npcol_phi, &
159 : num_grid_chunks, phase_hi
160 88 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_col_map, big_list, local_grid_idx, &
161 44 : my_atoms_A, my_atoms_B, &
162 44 : n_local_grid_atom, row_offset, &
163 44 : small_list
164 88 : INTEGER, DIMENSION(:), POINTER :: col_dist_ri, r_blk_sizes, ri_blk_sizes, &
165 44 : row_dist_grid
166 : LOGICAL :: do_scatter, use_dist
167 : REAL(KIND=dp) :: balance_A, balance_B, cutoff_ri, &
168 : item_start_time, r_c, t1
169 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cutoff_ri_per_atom, d_vec_local
170 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: D_local, d_lp_local, phi_local
171 : TYPE(cell_type), POINTER :: cell
172 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
173 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_b, fm_struct_D
174 : TYPE(cp_fm_type) :: fm_b, fm_D
175 : TYPE(cp_logger_type), POINTER :: logger
176 : TYPE(dbcsr_distribution_type) :: dist_phi, dist_Z
177 572 : TYPE(gw_3c_ctx_type) :: ctx_3c
178 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
179 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
180 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
181 : TYPE(section_vals_type), POINTER :: input
182 :
183 44 : CALL timeset(routineN, handle)
184 :
185 44 : t1 = m_walltime()
186 :
187 : CALL get_qs_env(qs_env, para_env=para_env, particle_set=particle_set, input=input, &
188 44 : qs_kind_set=qs_kind_set, cell=cell)
189 :
190 44 : NULLIFY (para_env_sub, blacs_env_sub)
191 :
192 44 : natom = bs_env%n_atom
193 44 : n_ao_total = bs_env%i_ao_end_from_atom(natom)
194 :
195 : ! =========================================================================
196 : ! 1. SETUP DBCSR TOPOLOGY & EXACT OFFSETS
197 : ! mat_Z_lP inherits the grid row blocking (and row distribution) of
198 : ! mat_phi_mu_l; its columns are one block per RI atom.
199 : ! =========================================================================
200 44 : CALL dbcsr_get_info(mat_phi_mu_l, row_blk_size=r_blk_sizes, distribution=dist_phi)
201 44 : CALL dbcsr_distribution_get(dist_phi, row_dist=row_dist_grid, npcols=npcol_phi)
202 :
203 44 : num_grid_chunks = SIZE(r_blk_sizes)
204 :
205 132 : ALLOCATE (row_offset(num_grid_chunks))
206 44 : row_offset(1) = 0
207 602 : DO i_blk = 2, num_grid_chunks
208 602 : row_offset(i_blk) = row_offset(i_blk - 1) + r_blk_sizes(i_blk - 1)
209 : END DO
210 :
211 176 : ALLOCATE (ri_blk_sizes(natom), col_dist_ri(natom))
212 150 : DO iatom = 1, natom
213 : ri_blk_sizes(iatom) = &
214 106 : bs_env%i_RI_end_from_atom(iatom) - bs_env%i_RI_start_from_atom(iatom) + 1
215 150 : col_dist_ri(iatom) = MOD(iatom - 1, npcol_phi)
216 : END DO
217 :
218 : CALL dbcsr_distribution_new(dist_Z, template=dist_phi, &
219 44 : row_dist=row_dist_grid, col_dist=col_dist_ri)
220 :
221 44 : IF (bs_env%ri_rs%Z_lP_exists) THEN
222 : CALL dbcsr_binary_read(filepath=TRIM(bs_env%prefix)//"Z_lP.matrix", &
223 : distribution=dist_Z, &
224 2 : matrix_new=mat_Z_lP)
225 2 : IF (bs_env%unit_nr > 0) THEN
226 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
227 1 : 'Read Z_lP from file ', ' Execution time', m_walltime() - t1, ' s'
228 : ! The grid rows use Morton order (spatial_atom_order). A Z_lP.matrix
229 : ! written with another grid ordering would be silently read into
230 : ! the current row order. Delete stale Z_lP.matrix files and recompute if in doubt.
231 : WRITE (bs_env%unit_nr, '(T2,A)') &
232 1 : '*** NOTE: Z_lP restart must match the current (spatial) grid row ordering ***'
233 1 : WRITE (bs_env%unit_nr, '(A)') ' '
234 : END IF
235 : ELSE
236 :
237 42 : IF (bs_env%unit_nr > 0) THEN
238 21 : WRITE (bs_env%unit_nr, '(A)') ' '
239 21 : WRITE (bs_env%unit_nr, '(T2,A)') 'Started computing Z_lP'
240 21 : CALL m_flush(bs_env%unit_nr)
241 : END IF
242 : CALL dbcsr_create(mat_Z_lP, name="mat_Z_lP", dist=dist_Z, &
243 : matrix_type=dbcsr_type_no_symmetry, &
244 42 : row_blk_size=r_blk_sizes, col_blk_size=ri_blk_sizes)
245 :
246 : ! Largest per-atom AO block, needed to size the 3c-integral work buffers.
247 42 : max_ao_size = 0
248 142 : DO j = 1, bs_env%n_atom
249 : max_ao_size = MAX(max_ao_size, &
250 : bs_env%i_ao_end_from_atom(j) - &
251 142 : bs_env%i_ao_start_from_atom(j) + 1)
252 : END DO
253 :
254 : ! Per-atom RI-RS integration sphere:
255 : ! cutoff_ri(P) = r_c + r_RI(P)
256 : ! where r_c is the truncated-Coulomb cutoff of the RI metric and r_RI the radius of
257 : ! the most diffuse RI auxiliary Gaussian on P.
258 126 : ALLOCATE (cutoff_ri_per_atom(natom))
259 :
260 42 : r_c = bs_env%ri_metric%cutoff_radius
261 142 : DO iatom = 1, natom
262 142 : cutoff_ri_per_atom(iatom) = r_c + bs_env%ri_rs%radius_ri_per_atom(iatom)
263 : END DO
264 :
265 42 : CALL print_sphere_cutoff_table(bs_env, cutoff_ri_per_atom)
266 :
267 : ! =========================================================================
268 : ! 2. PER-ATOM SOLVER CLASSIFICATION
269 : ! Split the atoms into "small" (single-rank LAPACK, Phase A) and "big"
270 : ! (distributed ScaLAPACK over subgroups of G ranks, Phase B) by comparing
271 : ! each atom's estimated solve peak memory against the measured budget.
272 : ! =========================================================================
273 : CALL classify_z_lp_atoms(bs_env, ri_rs_grid_points, cutoff_ri_per_atom, ri_blk_sizes, &
274 42 : n_local_grid_atom, small_list, n_small, big_list, n_big, G)
275 :
276 : ! LPT scheduling: sort the atoms of each phase by estimated solve cost
277 : ! (n_local_grid^3, Cholesky-dominated) and greedily assign to the least-loaded
278 : ! rank (Phase A) / subgroup (Phase B).
279 : CALL lpt_assign_atoms(small_list, n_small, n_local_grid_atom, para_env%num_pe, &
280 42 : para_env%mepos, my_atoms_A, balance_A)
281 42 : IF (n_big > 0) THEN
282 0 : n_groups = para_env%num_pe/G
283 0 : my_group = MIN(para_env%mepos/G, n_groups - 1)
284 : CALL lpt_assign_atoms(big_list, n_big, n_local_grid_atom, n_groups, my_group, &
285 0 : my_atoms_B, balance_B)
286 : ELSE
287 42 : ALLOCATE (my_atoms_B(0))
288 42 : balance_B = 1.0_dp
289 : END IF
290 :
291 : ! Atoms this rank will process across both phases for rank-0 progress
292 42 : n_my_atoms = SIZE(my_atoms_A) + SIZE(my_atoms_B)
293 42 : n_done = 0
294 :
295 : ! Shared context for the three-center integrals (μν|P) of the RHS build
296 : CALL gw_3c_ctx_create(ctx_3c, bs_env, bs_env%ri_metric, &
297 : basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, &
298 42 : basis_i=bs_env%basis_set_RI)
299 :
300 : ! =========================================================================
301 : ! 3. TWO-PHASE LOOP OVER ATOMS
302 : ! Phase A processes the "small" atoms with the single-rank BLAS path
303 : ! Phase B processes the "big" atoms with the distributed ScaLAPACK path
304 : ! over rank subgroups of size G. phi_local for each atom's cutoff sphere
305 : ! is built on the fly to avoid replicating a global grid x AO matrix.
306 : ! =========================================================================
307 126 : DO iphase = 1, 2
308 84 : IF (iphase == 1) THEN
309 42 : use_dist = .FALSE.
310 42 : atom_j_mepos = 0
311 42 : atom_j_stride = 1
312 42 : phase_hi = SIZE(my_atoms_A)
313 : ELSE
314 42 : IF (n_big == 0) CYCLE
315 0 : use_dist = .TRUE.
316 0 : n_groups = para_env%num_pe/G
317 0 : my_group = MIN(para_env%mepos/G, n_groups - 1)
318 0 : ALLOCATE (para_env_sub)
319 0 : CALL para_env_sub%from_split(para_env, my_group)
320 0 : CALL cp_blacs_env_create(blacs_env=blacs_env_sub, para_env=para_env_sub)
321 0 : atom_j_mepos = para_env_sub%mepos
322 0 : atom_j_stride = para_env_sub%num_pe
323 : ! All ranks of a subgroup share my_group, hence the identical my_atoms_B list
324 : ! (the per-atom ScaLAPACK solve is collective over the subgroup).
325 0 : phase_hi = SIZE(my_atoms_B)
326 : END IF
327 :
328 92 : DO idx = 1, phase_hi
329 50 : item_start_time = m_walltime()
330 50 : IF (iphase == 1) THEN
331 50 : atom_P = my_atoms_A(idx)
332 : ELSE
333 0 : atom_P = my_atoms_B(idx)
334 : END IF
335 :
336 50 : n_loc_ri = ri_blk_sizes(atom_P)
337 50 : cutoff_ri = cutoff_ri_per_atom(atom_P)
338 :
339 : ! ---------------------------------------------------------------------
340 : ! A. Sphere-local AO matrix ϕ_μ(r_l): select the grid points with
341 : ! |r_l - R_P| <= cutoff_ri(P), evaluate every AO on them, and drop
342 : ! points whose largest AO amplitude is below EPS_FILTER.
343 : ! ---------------------------------------------------------------------
344 : CALL build_phi_on_sphere(bs_env, qs_kind_set, &
345 : ri_rs_grid_points, atom_P, cutoff_ri, n_ao_total, &
346 : local_grid_idx, n_local_grid, phi_local, &
347 50 : ao_col_map, n_ao_used)
348 :
349 : ! ---------------------------------------------------------------------
350 : ! B. Right-hand side D_lP = Σ_μν ϕ_μ(r_l) ϕ_ν(r_l) (μν|P)
351 : ! ---------------------------------------------------------------------
352 200 : ALLOCATE (d_lp_local(n_local_grid, n_loc_ri))
353 50 : d_lp_local = 0.0_dp
354 :
355 : CALL compute_d_lp(bs_env, ctx_3c, phi_local, ao_col_map, d_lp_local, n_local_grid, &
356 50 : n_loc_ri, atom_P, max_ao_size, atom_j_mepos, atom_j_stride)
357 :
358 : ! Reduce per-subgroup-rank partials into the replicated d_lp_local.
359 : ! Skipped for BLAS path: each rank has the full sum locally.
360 50 : IF (use_dist) THEN
361 0 : CALL para_env_sub%sum(d_lp_local)
362 : END IF
363 :
364 : ! ---------------------------------------------------------------------
365 : ! C. Jacobi conditioning vector d_l = 1/sqrt(D_ll) and, on the BLAS path,
366 : ! the dense conditioned matrix D'_ll' = d_l D_ll' d_l' + λδ_ll'.
367 : ! ---------------------------------------------------------------------
368 150 : ALLOCATE (d_vec_local(n_local_grid))
369 :
370 50 : IF (.NOT. use_dist) THEN
371 : CALL build_gram_jacobi_blas(phi_local, n_local_grid, n_ao_used, &
372 50 : bs_env%ri_rs%tikhonov, D_local, d_vec_local)
373 : ELSE
374 : ! ScaLAPACK path: only d_vec is needed here (= 1/||phi(r_l)||^2);
375 : ! solve_D_lp_distributed builds its block-cyclic slice of D' internally
376 : ! with the squared+scaled values, so no dense D_local on this rank.
377 : CALL build_jacobi_diag_from_phi(phi_local, n_local_grid, n_ao_used, &
378 0 : d_vec_local)
379 : END IF
380 :
381 : ! ---------------------------------------------------------------------
382 : ! D. Pre-scale the RHS: D'_lP = d_l * D_lP
383 : ! ---------------------------------------------------------------------
384 50 : CALL scale_rows_by_diag(d_lp_local, d_vec_local, n_local_grid, n_loc_ri)
385 :
386 : ! ---------------------------------------------------------------------
387 : ! E. Cholesky solve Σ_l' D'_ll' Z'_l'P = D'_lP
388 : ! (BLAS dpotrf/dpotrs or ScaLAPACK pdpotrf/pdpotrs)
389 : ! ---------------------------------------------------------------------
390 50 : IF (.NOT. use_dist) THEN
391 50 : CALL timeset(routineN//"_dpotrf", handle_dpotrf)
392 50 : CALL dpotrf('L', n_local_grid, D_local, n_local_grid, info)
393 50 : CALL timestop(handle_dpotrf)
394 50 : IF (info /= 0) CPABORT("RI-RS Cholesky factorization failed")
395 50 : CALL timeset(routineN//"_dpotrs", handle_dpotrs)
396 : CALL dpotrs('L', n_local_grid, n_loc_ri, D_local, n_local_grid, &
397 50 : d_lp_local, n_local_grid, info)
398 50 : CALL timestop(handle_dpotrs)
399 50 : IF (info /= 0) CPABORT("RI-RS Cholesky solve failed")
400 50 : DEALLOCATE (D_local)
401 : ELSE
402 : CALL solve_D_lp_distributed(phi_local, d_vec_local, d_lp_local, &
403 : n_local_grid, n_ao_used, n_loc_ri, &
404 : bs_env%ri_rs%tikhonov, &
405 : para_env_sub, blacs_env_sub, &
406 0 : fm_struct_D, fm_struct_b, fm_D, fm_b, info)
407 0 : IF (info /= 0) CPABORT("Distributed RI-RS Cholesky solve failed")
408 : END IF
409 :
410 : ! ---------------------------------------------------------------------
411 : ! F. Undo the conditioning: Z_lP = d_l * Z'_lP
412 : ! ---------------------------------------------------------------------
413 50 : CALL scale_rows_by_diag(d_lp_local, d_vec_local, n_local_grid, n_loc_ri)
414 :
415 : ! ---------------------------------------------------------------------
416 : ! G. Scatter the solved Z columns back into the global sparse mat_Z_lP.
417 : ! ---------------------------------------------------------------------
418 50 : do_scatter = .TRUE.
419 50 : IF (use_dist) do_scatter = (para_env_sub%mepos == 0)
420 0 : IF (do_scatter) THEN
421 : CALL store_Z_lP_columns(mat_Z_lP, d_lp_local, local_grid_idx, n_local_grid, &
422 : n_loc_ri, atom_P, r_blk_sizes, row_offset, &
423 50 : bs_env%eps_filter)
424 : END IF
425 :
426 50 : DEALLOCATE (d_vec_local, d_lp_local)
427 50 : DEALLOCATE (local_grid_idx, phi_local, ao_col_map)
428 :
429 : ! Report each atom completed by the printing rank. No inter-rank communication is
430 : ! needed; the final message below is printed only after the global synchronization.
431 50 : n_done = n_done + 1
432 : CALL print_Z_lP_progress(bs_env, n_done, n_my_atoms, &
433 92 : m_walltime() - item_start_time)
434 : END DO ! idx: atoms of this phase owned by this rank / subgroup
435 :
436 : ! Tear down the Phase-B subgroup (all ranks created it collectively).
437 84 : IF (iphase == 2) THEN
438 0 : CALL cp_blacs_env_release(blacs_env_sub)
439 0 : CALL para_env_sub%free()
440 0 : DEALLOCATE (para_env_sub)
441 : END IF
442 : END DO ! iphase
443 :
444 42 : DEALLOCATE (cutoff_ri_per_atom)
445 42 : DEALLOCATE (small_list, big_list)
446 :
447 42 : CALL gw_3c_ctx_release(ctx_3c)
448 :
449 42 : CALL dbcsr_finalize(mat_Z_lP)
450 :
451 42 : CALL para_env%sync()
452 : CALL print_Z_lP_progress(bs_env, natom, natom, m_walltime() - t1, &
453 42 : all_mpi_ranks=.TRUE.)
454 :
455 42 : logger => cp_get_default_logger()
456 :
457 42 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
458 4 : CALL dbcsr_binary_write(matrix=mat_Z_lP, filepath=TRIM(bs_env%prefix)//"Z_lP.matrix")
459 : END IF
460 :
461 : END IF
462 :
463 44 : DEALLOCATE (row_offset, ri_blk_sizes, col_dist_ri)
464 44 : CALL dbcsr_distribution_release(dist_Z)
465 :
466 44 : DEALLOCATE (ri_rs_grid_points)
467 :
468 44 : CALL timestop(handle)
469 :
470 220 : END SUBROUTINE compute_Z_lP_standard
471 :
472 : ! **************************************************************************************************
473 : !> \brief Prints the per-kind maximum RI-RS integration-sphere cutoff table.
474 : !> \param bs_env ...
475 : !> \param cutoff_ri_per_atom ...
476 : ! **************************************************************************************************
477 42 : SUBROUTINE print_sphere_cutoff_table(bs_env, cutoff_ri_per_atom)
478 :
479 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
480 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutoff_ri_per_atom
481 :
482 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_sphere_cutoff_table'
483 :
484 : INTEGER :: handle, iatom, ikind, nkind
485 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cutoff_ri_per_kind
486 42 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
487 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
488 :
489 42 : CALL timeset(routineN, handle)
490 :
491 42 : atomic_kind_set => bs_env%ri_rs%atomic_kind_set
492 42 : particle_set => bs_env%ri_rs%particle_set
493 :
494 42 : IF (bs_env%unit_nr <= 0) THEN
495 21 : CALL timestop(handle)
496 : RETURN
497 : END IF
498 :
499 21 : nkind = SIZE(atomic_kind_set)
500 63 : ALLOCATE (cutoff_ri_per_kind(nkind))
501 21 : cutoff_ri_per_kind(:) = 0.0_dp
502 :
503 71 : DO iatom = 1, bs_env%n_atom
504 50 : ikind = particle_set(iatom)%atomic_kind%kind_number
505 71 : cutoff_ri_per_kind(ikind) = MAX(cutoff_ri_per_kind(ikind), cutoff_ri_per_atom(iatom))
506 : END DO
507 :
508 21 : WRITE (bs_env%unit_nr, '(T2,A)') 'Per-kind maximum RI-RS sphere cutoff (Å):'
509 21 : WRITE (bs_env%unit_nr, '(T4,A4,A14)') 'Kind', 'cutoff (Å)'
510 50 : DO ikind = 1, nkind
511 : WRITE (bs_env%unit_nr, '(T4,A4,F14.4)') &
512 29 : atomic_kind_set(ikind)%element_symbol, &
513 79 : cutoff_ri_per_kind(ikind)*angstrom
514 : END DO
515 21 : WRITE (bs_env%unit_nr, '(A)') ' '
516 :
517 21 : DEALLOCATE (cutoff_ri_per_kind)
518 :
519 21 : CALL timestop(handle)
520 :
521 42 : END SUBROUTINE print_sphere_cutoff_table
522 :
523 : ! **************************************************************************************************
524 : !> \brief LPT (longest-processing-time) assignment of the Z_lP atoms to workers (MPI ranks in
525 : !> Phase A, rank subgroups in Phase B): sort by estimated solve cost n_local_grid^3
526 : !> (the per-atom Cholesky dominates; the n^2 assembly terms order the atoms the same
527 : !> way) and greedily give each atom to the least-loaded worker.
528 : !> \param atom_list ...
529 : !> \param n_atoms ...
530 : !> \param n_local_grid_atom ...
531 : !> \param n_workers ...
532 : !> \param my_worker ...
533 : !> \param my_atoms ...
534 : !> \param max_over_mean ...
535 : ! **************************************************************************************************
536 42 : SUBROUTINE lpt_assign_atoms(atom_list, n_atoms, n_local_grid_atom, n_workers, my_worker, &
537 : my_atoms, max_over_mean)
538 :
539 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
540 : INTEGER, INTENT(IN) :: n_atoms
541 : INTEGER, DIMENSION(:), INTENT(IN) :: n_local_grid_atom
542 : INTEGER, INTENT(IN) :: n_workers, my_worker
543 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: my_atoms
544 : REAL(KIND=dp), INTENT(OUT) :: max_over_mean
545 :
546 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lpt_assign_atoms'
547 :
548 : INTEGER :: handle, i, iw, n_mine, w_min
549 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: mine_tmp, perm
550 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cost, load
551 :
552 42 : CALL timeset(routineN, handle)
553 :
554 42 : max_over_mean = 1.0_dp
555 42 : IF (n_atoms <= 0) THEN
556 0 : ALLOCATE (my_atoms(0))
557 0 : CALL timestop(handle)
558 0 : RETURN
559 : END IF
560 :
561 336 : ALLOCATE (cost(n_atoms), perm(n_atoms), mine_tmp(n_atoms), load(n_workers))
562 142 : DO i = 1, n_atoms
563 142 : cost(i) = REAL(n_local_grid_atom(atom_list(i)), dp)**3
564 : END DO
565 42 : CALL sort(cost, n_atoms, perm) ! ascending; walk backwards for largest-first
566 :
567 42 : load(:) = 0.0_dp
568 42 : n_mine = 0
569 142 : DO i = n_atoms, 1, -1
570 : w_min = 1
571 200 : DO iw = 2, n_workers
572 200 : IF (load(iw) < load(w_min)) w_min = iw
573 : END DO
574 100 : load(w_min) = load(w_min) + cost(i)
575 142 : IF (w_min - 1 == my_worker) THEN
576 50 : n_mine = n_mine + 1
577 50 : mine_tmp(n_mine) = atom_list(perm(i))
578 : END IF
579 : END DO
580 :
581 126 : ALLOCATE (my_atoms(n_mine))
582 92 : my_atoms(:) = mine_tmp(1:n_mine)
583 294 : IF (SUM(load) > 0.0_dp) max_over_mean = MAXVAL(load)*REAL(n_workers, dp)/SUM(load)
584 :
585 42 : CALL timestop(handle)
586 :
587 84 : END SUBROUTINE lpt_assign_atoms
588 :
589 : ! **************************************************************************************************
590 : !> \brief Computes the dense localized right-hand side for one RI atom P,
591 : !>
592 : !> d_lP = Σ_μν ϕ_μ(r_l) ϕ_ν(r_l) (μν|P).
593 : !>
594 : !> RI atom P, OMP-threaded over (atom_j, atom_k) AO-pair blocks: per thread, build the 3c
595 : !> block, then contract grid-chunked pair densities into a private d_lp partial; partials
596 : !> are reduced into d_lp at the end.
597 : !> Pair screening is handled inside build_3c_integral_block_ctx via the `screened` output.
598 : !> \param bs_env ...
599 : !> \param ctx ...
600 : !> \param phi_val ...
601 : !> \param ao_col_map ...
602 : !> \param d_lp ...
603 : !> \param n_grid_total ...
604 : !> \param n_loc_ri ...
605 : !> \param atom_P ...
606 : !> \param max_ao_size ...
607 : !> \param atom_j_mepos ...
608 : !> \param atom_j_stride ...
609 : ! **************************************************************************************************
610 50 : SUBROUTINE compute_d_lp(bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid_total, n_loc_ri, atom_P, &
611 : max_ao_size, atom_j_mepos, atom_j_stride)
612 :
613 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
614 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
615 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val
616 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
617 : INTEGER, INTENT(IN) :: n_grid_total, n_loc_ri
618 : REAL(KIND=dp), INTENT(INOUT) :: d_lp(n_grid_total, n_loc_ri)
619 : INTEGER, INTENT(IN) :: atom_P, max_ao_size, atom_j_mepos, &
620 : atom_j_stride
621 :
622 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp'
623 : INTEGER, PARAMETER :: grid_chunk = 1024
624 :
625 : INTEGER :: atom_j, atom_k, c, handle, handle_dgemm, &
626 : j, jk_idx, jsize, jstart, k, ksize, &
627 : kstart, l, l0, n_grid_pair, point, ri
628 50 : INTEGER, ALLOCATABLE :: grid_index(:)
629 : LOGICAL :: screened
630 50 : LOGICAL, ALLOCATABLE :: skip_grid_point(:, :)
631 : REAL(KIND=dp) :: pair_factor
632 50 : REAL(KIND=dp), ALLOCATABLE :: grid_result(:, :)
633 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_lp_prv, int_2d_prv, rho_chunk
634 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c_prv
635 50 : TYPE(gw_3c_ws_type) :: ws
636 :
637 50 : CALL timeset(routineN, handle)
638 :
639 : !$OMP PARALLEL DEFAULT(NONE) &
640 : !$OMP SHARED(skip_grid_point, bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid_total, n_loc_ri, atom_P, max_ao_size, &
641 : !$OMP atom_j_mepos, atom_j_stride) &
642 : !$OMP PRIVATE(grid_index, grid_result, n_grid_pair, point, pair_factor, &
643 : !$OMP atom_j, atom_k, c, handle_dgemm, j, jk_idx, jsize, jstart, k, ksize, kstart, &
644 50 : !$OMP l, l0, ri, screened, d_lp_prv, int_2d_prv, rho_chunk, int_3c_prv, ws)
645 :
646 : CALL gw_3c_ws_create(ws, ctx)
647 : ALLOCATE (int_3c_prv(max_ao_size, max_ao_size, n_loc_ri))
648 : ALLOCATE (int_2d_prv(max_ao_size*max_ao_size, n_loc_ri))
649 : ALLOCATE (rho_chunk(grid_chunk, max_ao_size*max_ao_size))
650 : ALLOCATE (d_lp_prv(n_grid_total, n_loc_ri))
651 : ALLOCATE (grid_index(n_grid_total), grid_result(grid_chunk, n_loc_ri))
652 : d_lp_prv(:, :) = 0.0_dp
653 :
654 : !$OMP SINGLE
655 : ALLOCATE (skip_grid_point(n_grid_total, bs_env%n_atom))
656 : CALL compute_skip_grid_point(bs_env, phi_val, ao_col_map, skip_grid_point)
657 : !$OMP END SINGLE
658 :
659 : ! MPI assigns each unordered pair through its first atom; OpenMP divides those atoms.
660 : ! Skip only exact-zero pair-grid support, without a new screening threshold.
661 : !$OMP DO SCHEDULE(DYNAMIC)
662 : DO atom_j = atom_j_mepos + 1, bs_env%n_atom, atom_j_stride
663 : DO atom_k = atom_j, bs_env%n_atom
664 : jstart = ao_col_map(bs_env%i_ao_start_from_atom(atom_j))
665 : kstart = ao_col_map(bs_env%i_ao_start_from_atom(atom_k))
666 : IF (jstart == 0 .OR. kstart == 0) CYCLE
667 : jsize = bs_env%i_ao_end_from_atom(atom_j) - bs_env%i_ao_start_from_atom(atom_j) + 1
668 : ksize = bs_env%i_ao_end_from_atom(atom_k) - bs_env%i_ao_start_from_atom(atom_k) + 1
669 :
670 : n_grid_pair = 0
671 : DO point = 1, n_grid_total
672 : IF (skip_grid_point(point, atom_j) .OR. skip_grid_point(point, atom_k)) CYCLE
673 : n_grid_pair = n_grid_pair + 1
674 : grid_index(n_grid_pair) = point
675 : END DO
676 : IF (n_grid_pair == 0) CYCLE
677 : ! (μν|P) = (νμ|P): distinct atom pairs contribute twice.
678 : pair_factor = 1.0_dp
679 : IF (atom_j /= atom_k) pair_factor = 2.0_dp
680 :
681 : int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp
682 :
683 : ! Compute B_{μν,P} = (μν|P); ctx-internal triangle-inequality screening on
684 : ! kind_radius sets `screened=.TRUE.` for negligible triples.
685 : CALL build_3c_integral_block_ctx(int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri), &
686 : ctx, ws, atom_j=atom_j, atom_k=atom_k, atom_i=atom_P, &
687 : screened=screened)
688 :
689 : IF (screened) CYCLE
690 :
691 : ! Flatten 3D B_{μν, P} tensor to 2D B_{(μν), P} matrix for BLAS
692 : DO ri = 1, n_loc_ri
693 : DO k = 1, ksize
694 : DO j = 1, jsize
695 : jk_idx = (k - 1)*jsize + j
696 : int_2d_prv(jk_idx, ri) = int_3c_prv(j, k, ri)
697 : END DO
698 : END DO
699 : END DO
700 :
701 : ! Pair density ρ(l, μν) = ϕ_μ(r_l) ϕ_ν(r_l) in grid chunks, contracted on the fly:
702 : ! d_{l,P} += ρ(l, μν) B_{(μν),P} (dgemm runs serially inside the parallel region)
703 : DO l0 = 1, n_grid_pair, grid_chunk
704 : c = MIN(grid_chunk, n_grid_pair - l0 + 1)
705 : DO k = 1, ksize
706 : DO j = 1, jsize
707 : jk_idx = (k - 1)*jsize + j
708 : DO l = 1, c
709 : point = grid_index(l0 + l - 1)
710 : rho_chunk(l, jk_idx) = phi_val(point, jstart + j - 1)* &
711 : phi_val(point, kstart + k - 1)
712 : END DO
713 : END DO
714 : END DO
715 : CALL timeset(routineN//"_dgemm", handle_dgemm)
716 : CALL dgemm("N", "N", c, n_loc_ri, jsize*ksize, &
717 : pair_factor, rho_chunk, grid_chunk, &
718 : int_2d_prv, max_ao_size*max_ao_size, &
719 : 0.0_dp, grid_result, grid_chunk)
720 : DO ri = 1, n_loc_ri
721 : DO l = 1, c
722 : point = grid_index(l0 + l - 1)
723 : d_lp_prv(point, ri) = d_lp_prv(point, ri) + grid_result(l, ri)
724 : END DO
725 : END DO
726 : CALL timestop(handle_dgemm)
727 : END DO
728 : END DO
729 : END DO
730 : !$OMP END DO
731 :
732 : !$OMP CRITICAL (compute_d_lp_reduce)
733 : d_lp(1:n_grid_total, 1:n_loc_ri) = d_lp(1:n_grid_total, 1:n_loc_ri) + &
734 : d_lp_prv(1:n_grid_total, 1:n_loc_ri)
735 : !$OMP END CRITICAL (compute_d_lp_reduce)
736 :
737 : DEALLOCATE (int_3c_prv, int_2d_prv, rho_chunk, d_lp_prv, grid_index, grid_result)
738 : CALL gw_3c_ws_release(ws)
739 :
740 : !$OMP SINGLE
741 : DEALLOCATE (skip_grid_point)
742 : !$OMP END SINGLE
743 :
744 : !$OMP END PARALLEL
745 :
746 50 : CALL timestop(handle)
747 :
748 50 : END SUBROUTINE compute_d_lp
749 :
750 : ! **************************************************************************************************
751 : !> \brief Marks grid points where all stored AO values of an atom are exactly zero.
752 : !> \param bs_env ...
753 : !> \param phi_val ...
754 : !> \param ao_col_map ...
755 : !> \param skip_grid_point ...
756 : ! **************************************************************************************************
757 50 : SUBROUTINE compute_skip_grid_point(bs_env, phi_val, ao_col_map, skip_grid_point)
758 :
759 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
760 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val
761 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
762 : LOGICAL, DIMENSION(:, :), INTENT(OUT) :: skip_grid_point
763 :
764 : INTEGER :: atom, first_ao, number_of_aos
765 :
766 46697 : skip_grid_point(:, :) = .TRUE.
767 174 : DO atom = 1, bs_env%n_atom
768 124 : first_ao = ao_col_map(bs_env%i_ao_start_from_atom(atom))
769 124 : IF (first_ao == 0) CYCLE
770 124 : number_of_aos = bs_env%i_ao_end_from_atom(atom) - bs_env%i_ao_start_from_atom(atom) + 1
771 169237 : skip_grid_point(:, atom) = ALL(phi_val(:, first_ao:first_ao + number_of_aos - 1) == 0.0_dp, DIM=2)
772 : END DO
773 :
774 50 : END SUBROUTINE compute_skip_grid_point
775 :
776 : ! **************************************************************************************************
777 : !> \brief Computes the RI-RS matrix Z_lP for an automatically optimized RI basis.
778 : !>
779 : !> The fitting coefficients Z_lP satisfy
780 : !>
781 : !> Σ_l' D_ll' Z_l'P = d_lP,
782 : !> D_ll' = [Σ_μ ϕ_μ(r_l) ϕ_μ(r_l')]²,
783 : !> d_lP = Σ_μν ϕ_μ(r_l) ϕ_ν(r_l) (μν|P).
784 : !>
785 : !> For the standard RI-RS fit, P is an atom-centered RI function. Here an optimized
786 : !> function may contain reference functions on neighboring atoms,
787 : !>
788 : !> φ_p(r) = Σ_A Σ_{P∈A} U_Pp^A φ_P^A(r),
789 : !>
790 : !> so all atomic contributions to d_lp must be accumulated before solving for Z_lp.
791 : !> If every optimized column uses the complete RI-RS grid, all columns are solved in one
792 : !> system. Otherwise, each atom-blocked column set is fitted on its own integration sphere;
793 : !> grid points outside that sphere are excluded from its fitting equations.
794 : !> \param qs_env ...
795 : !> \param bs_env ...
796 : !> \param ri_rs_grid_points ...
797 : !> \param mat_phi_mu_l ...
798 : !> \param mat_Z_lP ...
799 : ! **************************************************************************************************
800 4 : SUBROUTINE compute_Z_lP_auto_ri(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
801 : TYPE(qs_environment_type), POINTER :: qs_env
802 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
803 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ri_rs_grid_points
804 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
805 : TYPE(dbcsr_type), INTENT(OUT) :: mat_Z_lP
806 :
807 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
808 : routineN = 'compute_Z_lP_auto_ri'
809 :
810 : INTEGER :: AB_block, atom_A, atom_B, block_size_b, column_first, first_p_AB, fit_atom, &
811 : handle, handle_dpotrf, handle_dpotrs, info, max_ao_size, max_nRI_ref, mypcol, myprow, &
812 : n_ao_total, n_ao_used, n_done, n_my_atoms, n_to_a, natom, ncol, ngrid, nri, nRI_ref_A, &
813 : nRI_ref_B, output_offset, ri_atom
814 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_col_map, local_grid_idx, row_offset
815 8 : INTEGER, DIMENSION(:), POINTER :: col_dist_ri, ri_blk_sizes, &
816 4 : row_dist_grid, row_size_grid
817 : LOGICAL :: AB_block_local, common_grid_available, &
818 : have_fitted_columns, &
819 : reuse_atomic_integrals
820 4 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: active_atom
821 : REAL(KIND=dp) :: cutoff_ri, item_start_time, r_c, t1
822 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: d_vec
823 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: D_local, d_lp_all, d_lp_local, &
824 4 : phi_local, U_Pp_AB
825 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: U_Pp_by_atom
826 : REAL(KIND=dp), DIMENSION(3) :: center
827 : TYPE(cell_type), POINTER :: cell
828 : TYPE(cp_logger_type), POINTER :: logger
829 : TYPE(dbcsr_distribution_type) :: dist_Z
830 : TYPE(dbcsr_type) :: mat_rhs
831 52 : TYPE(gw_3c_ctx_type) :: ctx_3c
832 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_col
833 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
834 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
835 : TYPE(section_vals_type), POINTER :: input
836 :
837 4 : CALL timeset(routineN, handle)
838 4 : t1 = m_walltime()
839 4 : IF (bs_env%unit_nr > 0) THEN
840 2 : WRITE (bs_env%unit_nr, '(A)') ' '
841 2 : WRITE (bs_env%unit_nr, '(T2,A)') 'Started computing Z_lP'
842 2 : CALL m_flush(bs_env%unit_nr)
843 : END IF
844 : CALL get_qs_env(qs_env, para_env=para_env, particle_set=particle_set, input=input, &
845 4 : qs_kind_set=qs_kind_set, cell=cell)
846 4 : NULLIFY (para_env_col)
847 :
848 4 : natom = bs_env%n_atom
849 4 : n_ao_total = bs_env%i_ao_end_from_atom(natom)
850 4 : CPASSERT(bs_env%auto_ri%AB_block_count > 0)
851 4 : CPASSERT(SIZE(particle_set) == natom)
852 :
853 : CALL prepare_Z_lP_auto_ri(bs_env, mat_phi_mu_l, mat_Z_lP, dist_Z, row_size_grid, &
854 : row_dist_grid, col_dist_ri, ri_blk_sizes, row_offset, &
855 4 : myprow, mypcol, max_ao_size)
856 :
857 : CALL gw_3c_ctx_create(ctx_3c, bs_env, bs_env%ri_metric, &
858 : basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, &
859 4 : basis_i=bs_env%basis_set_RI)
860 :
861 : ! Check whether one grid can represent d_lp for every optimized atomic column block.
862 : CALL common_Z_lP_grid_available(bs_env, ri_rs_grid_points, &
863 4 : common_grid_available, have_fitted_columns)
864 :
865 4 : IF (common_grid_available .AND. have_fitted_columns) THEN
866 : ! Evaluate ϕ_μ(r_l) once on the common grid.
867 : CALL build_phi_on_complete_grid(bs_env, qs_kind_set, &
868 : ri_rs_grid_points, n_ao_total, local_grid_idx, ngrid, &
869 4 : phi_local, ao_col_map, n_ao_used)
870 16 : nri = SUM(ri_blk_sizes)
871 16 : ALLOCATE (d_lp_all(ngrid, nri), source=0.0_dp)
872 : ! d_lp = Σ_A Σ_{P∈A} d_lP U_Pp for every optimized RI function q.
873 : CALL compute_d_lp_auto_ri_batch(bs_env, ctx_3c, phi_local, ao_col_map, d_lp_all, ngrid, &
874 4 : max_ao_size, para_env%mepos, para_env%num_pe)
875 4 : CALL para_env%sum(d_lp_all)
876 12 : ALLOCATE (d_vec(ngrid))
877 : ! D_ll' = [Σ_μ ϕ_μ(r_l) ϕ_μ(r_l')]², followed by D Z = d.
878 : CALL build_gram_jacobi_blas(phi_local, ngrid, n_ao_used, bs_env%ri_rs%tikhonov, &
879 4 : D_local, d_vec)
880 4 : CALL scale_rows_by_diag(d_lp_all, d_vec, ngrid, nri)
881 4 : CALL timeset(routineN//'_dpotrf', handle_dpotrf)
882 4 : CALL dpotrf('L', ngrid, D_local, ngrid, info)
883 4 : CALL timestop(handle_dpotrf)
884 4 : CPASSERT(info == 0)
885 4 : CALL timeset(routineN//'_dpotrs', handle_dpotrs)
886 4 : CALL dpotrs('L', ngrid, nri, D_local, ngrid, d_lp_all, ngrid, info)
887 4 : CALL timestop(handle_dpotrs)
888 4 : CPASSERT(info == 0)
889 4 : CALL scale_rows_by_diag(d_lp_all, d_vec, ngrid, nri)
890 :
891 22 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
892 18 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
893 18 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
894 18 : AB_block_local = col_dist_ri(atom_A) == mypcol
895 18 : IF (atom_B /= atom_A) AB_block_local = AB_block_local .OR. col_dist_ri(atom_B) == mypcol
896 12 : IF (.NOT. AB_block_local) CYCLE
897 18 : ncol = bs_env%auto_ri%AB_size_opt_RI(AB_block)
898 18 : n_to_a = bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
899 18 : IF (n_to_a > 0 .AND. col_dist_ri(atom_A) == mypcol) THEN
900 : output_offset = SUM(ri_blk_sizes(:atom_A - 1)) + &
901 28 : bs_env%auto_ri%AB_first_p_A(AB_block) - 1
902 : CALL add_Z_lP_columns(mat_Z_lP, &
903 : d_lp_all(:, output_offset + 1:output_offset + n_to_a), &
904 : local_grid_idx, ngrid, atom_A, &
905 : bs_env%auto_ri%AB_first_p_A(AB_block), &
906 : ri_blk_sizes(atom_A), row_size_grid, row_offset, &
907 : row_dist_grid, &
908 14 : myprow, bs_env%eps_filter)
909 : END IF
910 22 : IF (n_to_a < ncol) THEN
911 4 : CPASSERT(atom_B /= atom_A)
912 4 : IF (col_dist_ri(atom_B) == mypcol) THEN
913 4 : block_size_b = ri_blk_sizes(atom_B)
914 : output_offset = SUM(ri_blk_sizes(:atom_B - 1)) + &
915 10 : bs_env%auto_ri%AB_first_p_B(AB_block) - 1
916 : CALL add_Z_lP_columns(mat_Z_lP, &
917 : d_lp_all(:, output_offset + 1: &
918 : output_offset + ncol - n_to_a), &
919 : local_grid_idx, ngrid, atom_B, &
920 : bs_env%auto_ri%AB_first_p_B(AB_block), block_size_b, &
921 : row_size_grid, row_offset, row_dist_grid, myprow, &
922 4 : bs_env%eps_filter)
923 : END IF
924 : END IF
925 : END DO
926 12 : DEALLOCATE (D_local, d_vec, d_lp_all, local_grid_idx, phi_local, ao_col_map)
927 : ELSE
928 : reuse_atomic_integrals = &
929 : bs_env%ri_rs%cutoff_radius_ri_ao <= 0.0_dp .OR. &
930 : bs_env%ri_rs%cutoff_radius_ri_ao <= &
931 0 : MINVAL(bs_env%ri_rs%radius_ao_per_atom)
932 : IF (reuse_atomic_integrals) THEN
933 : ! Form d_lp = Σ_A Σ_{P∈A} d_lP U_Pp before fitting each atomic column block.
934 : CALL compute_auto_ri_d_lp(qs_env, bs_env, ctx_3c, &
935 : ri_rs_grid_points, mat_phi_mu_l, mat_rhs, &
936 0 : max_ao_size)
937 0 : CALL dbcsr_release(mat_Z_lP)
938 : CALL dbcsr_create(mat_Z_lP, name='mat_Z_lP localized AA/AB', dist=dist_Z, &
939 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_size_grid, &
940 0 : col_blk_size=ri_blk_sizes)
941 : ! Solve Σ_l' D_ll' Z_l'q = d_lp on each atomic fitting grid.
942 0 : CALL fit_auto_ri_z_lp(qs_env, bs_env, ri_rs_grid_points, mat_rhs, mat_Z_lP)
943 0 : CALL dbcsr_release(mat_rhs)
944 : ELSE
945 0 : ALLOCATE (para_env_col)
946 0 : CALL para_env_col%from_split(para_env, mypcol)
947 0 : n_my_atoms = 0
948 0 : DO fit_atom = 1, natom
949 0 : IF (col_dist_ri(fit_atom) /= mypcol) CYCLE
950 0 : IF (ri_blk_sizes(fit_atom) == 0) CYCLE
951 0 : n_my_atoms = n_my_atoms + 1
952 : END DO
953 0 : n_done = 0
954 0 : max_nRI_ref = 0
955 0 : DO ri_atom = 1, natom
956 0 : max_nRI_ref = MAX(max_nRI_ref, get_ref_RI_size(bs_env, ri_atom))
957 : END DO
958 0 : DO fit_atom = 1, natom
959 0 : IF (col_dist_ri(fit_atom) /= mypcol) CYCLE
960 0 : ncol = ri_blk_sizes(fit_atom)
961 0 : IF (ncol == 0) CYCLE
962 0 : item_start_time = m_walltime()
963 0 : ALLOCATE (U_Pp_by_atom(max_nRI_ref, ncol, natom), source=0.0_dp)
964 0 : ALLOCATE (active_atom(natom), source=.FALSE.)
965 :
966 0 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
967 0 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
968 0 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
969 0 : n_to_a = bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
970 0 : IF (fit_atom == atom_A .AND. n_to_a > 0) THEN
971 0 : column_first = bs_env%auto_ri%AB_first_p_A(AB_block)
972 0 : first_p_AB = 1
973 0 : ncol = n_to_a
974 0 : ELSE IF (fit_atom == atom_B .AND. &
975 : n_to_a < bs_env%auto_ri%AB_size_opt_RI(AB_block)) THEN
976 0 : column_first = bs_env%auto_ri%AB_first_p_B(AB_block)
977 0 : first_p_AB = n_to_a + 1
978 0 : ncol = bs_env%auto_ri%AB_size_opt_RI(AB_block) - n_to_a
979 : ELSE
980 : CYCLE
981 : END IF
982 0 : CALL get_U_Pp_AB(bs_env%auto_ri, AB_block, U_Pp_AB)
983 :
984 0 : nRI_ref_A = get_ref_RI_size(bs_env, atom_A)
985 : U_Pp_by_atom(1:nRI_ref_A, &
986 : column_first:column_first + ncol - 1, atom_A) = &
987 : U_Pp_AB( &
988 0 : 1:nRI_ref_A, first_p_AB:first_p_AB + ncol - 1)
989 0 : active_atom(atom_A) = .TRUE.
990 0 : IF (atom_B /= atom_A) THEN
991 0 : nRI_ref_B = get_ref_RI_size(bs_env, atom_B)
992 : U_Pp_by_atom(1:nRI_ref_B, &
993 : column_first:column_first + ncol - 1, atom_B) = &
994 : U_Pp_AB( &
995 : nRI_ref_A + 1:nRI_ref_A + nRI_ref_B, &
996 0 : first_p_AB:first_p_AB + ncol - 1)
997 0 : active_atom(atom_B) = .TRUE.
998 : END IF
999 0 : DEALLOCATE (U_Pp_AB)
1000 : END DO
1001 :
1002 0 : center = particle_set(fit_atom)%r
1003 0 : r_c = bs_env%ri_metric%cutoff_radius
1004 0 : cutoff_ri = 0.0_dp
1005 0 : DO ri_atom = 1, natom
1006 0 : IF (.NOT. active_atom(ri_atom)) CYCLE
1007 : cutoff_ri = MAX(cutoff_ri, &
1008 : r_c + bs_env%ri_rs%radius_ri_per_atom(ri_atom) + &
1009 0 : NORM2(center - particle_set(ri_atom)%r))
1010 : END DO
1011 : CALL build_phi_on_sphere(bs_env, qs_kind_set, &
1012 : ri_rs_grid_points, fit_atom, cutoff_ri, n_ao_total, &
1013 : local_grid_idx, ngrid, phi_local, ao_col_map, n_ao_used, &
1014 0 : center=center)
1015 0 : ncol = ri_blk_sizes(fit_atom)
1016 0 : ALLOCATE (d_lp_local(ngrid, ncol), source=0.0_dp)
1017 : ! d_lp = Σ_A Σ_{P∈A} d_lP U_Pp on the grid of this atomic column block.
1018 : CALL compute_d_lp_auto_ri_atoms(bs_env, ctx_3c, phi_local, ao_col_map, &
1019 : d_lp_local, ngrid, U_Pp_by_atom, &
1020 : active_atom, max_ao_size, &
1021 0 : para_env_col%mepos, para_env_col%num_pe)
1022 0 : CALL para_env_col%sum(d_lp_local)
1023 0 : ALLOCATE (d_vec(ngrid))
1024 : ! D_ll' = [Σ_μ ϕ_μ(r_l) ϕ_μ(r_l')]², followed by D Z = d.
1025 : CALL build_gram_jacobi_blas(phi_local, ngrid, n_ao_used, bs_env%ri_rs%tikhonov, &
1026 0 : D_local, d_vec)
1027 0 : CALL scale_rows_by_diag(d_lp_local, d_vec, ngrid, ncol)
1028 0 : CALL timeset(routineN//'_dpotrf', handle_dpotrf)
1029 0 : CALL dpotrf('L', ngrid, D_local, ngrid, info)
1030 0 : CALL timestop(handle_dpotrf)
1031 0 : CPASSERT(info == 0)
1032 0 : CALL timeset(routineN//'_dpotrs', handle_dpotrs)
1033 0 : CALL dpotrs('L', ngrid, ncol, D_local, ngrid, d_lp_local, ngrid, info)
1034 0 : CALL timestop(handle_dpotrs)
1035 0 : CPASSERT(info == 0)
1036 0 : CALL scale_rows_by_diag(d_lp_local, d_vec, ngrid, ncol)
1037 : CALL add_Z_lP_columns(mat_Z_lP, d_lp_local, local_grid_idx, ngrid, fit_atom, 1, &
1038 : ri_blk_sizes(fit_atom), row_size_grid, row_offset, &
1039 : row_dist_grid, &
1040 0 : myprow, bs_env%eps_filter)
1041 0 : DEALLOCATE (D_local, d_vec, d_lp_local, local_grid_idx, phi_local, ao_col_map)
1042 0 : DEALLOCATE (U_Pp_by_atom, active_atom)
1043 0 : n_done = n_done + 1
1044 : CALL print_Z_lP_progress(bs_env, n_done, n_my_atoms, &
1045 0 : m_walltime() - item_start_time)
1046 : END DO
1047 0 : CALL para_env_col%free()
1048 0 : DEALLOCATE (para_env_col)
1049 : END IF
1050 : END IF
1051 4 : CALL gw_3c_ctx_release(ctx_3c)
1052 :
1053 4 : CALL dbcsr_filter(mat_Z_lP, bs_env%eps_filter)
1054 4 : CALL dbcsr_finalize(mat_Z_lP)
1055 4 : CALL para_env%sync()
1056 : CALL print_Z_lP_progress(bs_env, natom, natom, m_walltime() - t1, &
1057 4 : all_mpi_ranks=.TRUE.)
1058 4 : logger => cp_get_default_logger()
1059 4 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
1060 0 : CALL dbcsr_binary_write(matrix=mat_Z_lP, filepath=TRIM(bs_env%prefix)//'Z_lP.matrix')
1061 : END IF
1062 :
1063 4 : DEALLOCATE (row_offset, ri_blk_sizes, col_dist_ri)
1064 4 : CALL dbcsr_distribution_release(dist_Z)
1065 4 : CALL timestop(handle)
1066 :
1067 12 : END SUBROUTINE compute_Z_lP_auto_ri
1068 :
1069 : ! **************************************************************************************************
1070 : !> \brief Prints completed Z_lP atoms and execution time from the printing rank.
1071 : !> \param bs_env ...
1072 : !> \param n_done Number of atoms completed by the printing rank.
1073 : !> \param n_total Total number of atoms assigned to the printing rank.
1074 : !> \param execution_time Wall-clock time used for the reported work.
1075 : !> \param all_mpi_ranks Whether all MPI ranks have completed the calculation.
1076 : ! **************************************************************************************************
1077 102 : SUBROUTINE print_Z_lP_progress(bs_env, n_done, n_total, execution_time, all_mpi_ranks)
1078 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1079 : INTEGER, INTENT(IN) :: n_done, n_total
1080 : REAL(KIND=dp), INTENT(IN) :: execution_time
1081 : LOGICAL, INTENT(IN), OPTIONAL :: all_mpi_ranks
1082 :
1083 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_Z_lP_progress'
1084 :
1085 : INTEGER :: handle
1086 : LOGICAL :: completed
1087 :
1088 102 : CALL timeset(routineN, handle)
1089 :
1090 102 : IF (bs_env%unit_nr > 0) THEN
1091 56 : completed = .FALSE.
1092 56 : IF (PRESENT(all_mpi_ranks)) completed = all_mpi_ranks
1093 23 : IF (completed) THEN
1094 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A,/)') &
1095 23 : 'Computed Z_lP (all MPI ranks) for all atoms,', &
1096 46 : 'Execution time', execution_time, ' s'
1097 : ELSE
1098 : WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F7.1,A)') &
1099 33 : 'Computed Z_lP (MPI rank 0) for atom', n_done, ' /', n_total, &
1100 66 : ', Execution time', execution_time, ' s'
1101 : END IF
1102 56 : CALL m_flush(bs_env%unit_nr)
1103 : END IF
1104 :
1105 102 : CALL timestop(handle)
1106 :
1107 102 : END SUBROUTINE print_Z_lP_progress
1108 :
1109 : ! **************************************************************************************************
1110 : !> \brief Prepares the distributed Z_lP matrix and its atom-blocked column layout. Column block A
1111 : !> contains all optimized functions assigned to atom A. Every local block is reserved once
1112 : !> because several AA/AB contraction blocks may contribute to the same block.
1113 : !> \param bs_env ...
1114 : !> \param mat_phi_mu_l ...
1115 : !> \param mat_Z_lP ...
1116 : !> \param dist_Z ...
1117 : !> \param row_size_grid ...
1118 : !> \param row_dist_grid ...
1119 : !> \param col_dist_ri ...
1120 : !> \param ri_blk_sizes ...
1121 : !> \param row_offset ...
1122 : !> \param myprow ...
1123 : !> \param mypcol ...
1124 : !> \param max_ao_size ...
1125 : ! **************************************************************************************************
1126 4 : SUBROUTINE prepare_Z_lP_auto_ri(bs_env, mat_phi_mu_l, mat_Z_lP, dist_Z, row_size_grid, &
1127 : row_dist_grid, col_dist_ri, ri_blk_sizes, row_offset, &
1128 : myprow, mypcol, max_ao_size)
1129 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1130 : TYPE(dbcsr_type), INTENT(IN) :: mat_phi_mu_l
1131 : TYPE(dbcsr_type), INTENT(OUT) :: mat_Z_lP
1132 : TYPE(dbcsr_distribution_type), INTENT(OUT) :: dist_Z
1133 : INTEGER, DIMENSION(:), POINTER :: row_size_grid, row_dist_grid, &
1134 : col_dist_ri, ri_blk_sizes
1135 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: row_offset
1136 : INTEGER, INTENT(OUT) :: myprow, mypcol, max_ao_size
1137 :
1138 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_Z_lP_auto_ri'
1139 :
1140 : INTEGER :: handle, i_blk, iatom, npcol
1141 : TYPE(dbcsr_distribution_type) :: dist_phi
1142 :
1143 4 : CALL timeset(routineN, handle)
1144 :
1145 4 : CALL dbcsr_get_info(mat_phi_mu_l, row_blk_size=row_size_grid, distribution=dist_phi)
1146 : CALL dbcsr_distribution_get(dist_phi, row_dist=row_dist_grid, npcols=npcol, &
1147 4 : myprow=myprow, mypcol=mypcol)
1148 12 : ALLOCATE (row_offset(SIZE(row_size_grid)))
1149 4 : row_offset(1) = 0
1150 12 : DO i_blk = 2, SIZE(row_size_grid)
1151 12 : row_offset(i_blk) = row_offset(i_blk - 1) + row_size_grid(i_blk - 1)
1152 : END DO
1153 :
1154 20 : ALLOCATE (ri_blk_sizes(bs_env%n_atom), col_dist_ri(bs_env%n_atom))
1155 32 : ri_blk_sizes = bs_env%auto_ri%sizes_opt_RI
1156 16 : DO iatom = 1, bs_env%n_atom
1157 16 : col_dist_ri(iatom) = MOD(iatom - 1, npcol)
1158 : END DO
1159 : CALL dbcsr_distribution_new(dist_Z, template=dist_phi, row_dist=row_dist_grid, &
1160 4 : col_dist=col_dist_ri)
1161 : CALL dbcsr_create(mat_Z_lP, name='mat_Z_lP localized AA/AB', dist=dist_Z, &
1162 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_size_grid, &
1163 4 : col_blk_size=ri_blk_sizes)
1164 4 : CALL dbcsr_reserve_all_blocks(mat_Z_lP)
1165 4 : CALL dbcsr_set(mat_Z_lP, 0.0_dp)
1166 :
1167 4 : max_ao_size = 0
1168 16 : DO iatom = 1, bs_env%n_atom
1169 : max_ao_size = MAX(max_ao_size, bs_env%i_ao_end_from_atom(iatom) - &
1170 16 : bs_env%i_ao_start_from_atom(iatom) + 1)
1171 : END DO
1172 :
1173 4 : CALL timestop(handle)
1174 :
1175 8 : END SUBROUTINE prepare_Z_lP_auto_ri
1176 :
1177 : ! **************************************************************************************************
1178 : !> \brief Tests whether every fitted atomic block can use the complete RI-RS grid.
1179 : !>
1180 : !> For a block centered on A, every grid point must satisfy
1181 : !>
1182 : !> |r_l - R_A| <= R_A^fit,
1183 : !>
1184 : !> and every atom B whose AOs can contribute must satisfy
1185 : !>
1186 : !> |R_B - R_A| <= R_B^AO + R_A^fit.
1187 : !>
1188 : !> Pair-midpoint grids do not satisfy this atom-centered criterion in general.
1189 : !> \param bs_env ...
1190 : !> \param ri_rs_grid_points ...
1191 : !> \param common_grid_available ...
1192 : !> \param have_fitted_columns ...
1193 : ! **************************************************************************************************
1194 4 : SUBROUTINE common_Z_lP_grid_available(bs_env, ri_rs_grid_points, &
1195 : common_grid_available, have_fitted_columns)
1196 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1197 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ri_rs_grid_points
1198 : LOGICAL, INTENT(OUT) :: common_grid_available, &
1199 : have_fitted_columns
1200 :
1201 : CHARACTER(LEN=*), PARAMETER :: routineN = 'common_Z_lP_grid_available'
1202 :
1203 : INTEGER :: AB_block, atom_A, atom_B, fit_atom, &
1204 : handle, iatom, igrid, n_to_a, natom
1205 4 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: active_atom
1206 : REAL(KIND=dp) :: cutoff_ri
1207 : REAL(KIND=dp), DIMENSION(3) :: center
1208 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1209 :
1210 4 : CALL timeset(routineN, handle)
1211 :
1212 4 : particle_set => bs_env%ri_rs%particle_set
1213 4 : natom = bs_env%n_atom
1214 4 : common_grid_available = .TRUE.
1215 4 : have_fitted_columns = .FALSE.
1216 12 : ALLOCATE (active_atom(natom))
1217 16 : DO fit_atom = 1, natom
1218 12 : IF (bs_env%auto_ri%sizes_opt_RI(fit_atom) == 0) CYCLE
1219 12 : have_fitted_columns = .TRUE.
1220 12 : active_atom = .FALSE.
1221 66 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
1222 54 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
1223 54 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
1224 54 : n_to_a = bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
1225 54 : IF (.NOT. (fit_atom == atom_A .AND. n_to_a > 0) .AND. &
1226 : .NOT. (fit_atom == atom_B .AND. &
1227 : n_to_a < bs_env%auto_ri%AB_size_opt_RI(AB_block))) CYCLE
1228 18 : active_atom(atom_A) = .TRUE.
1229 30 : IF (atom_B /= atom_A) active_atom(atom_B) = .TRUE.
1230 : END DO
1231 :
1232 48 : center = particle_set(fit_atom)%r
1233 : cutoff_ri = 0.0_dp
1234 48 : DO iatom = 1, natom
1235 36 : IF (.NOT. active_atom(iatom)) CYCLE
1236 : cutoff_ri = MAX(cutoff_ri, bs_env%ri_metric%cutoff_radius + &
1237 : bs_env%ri_rs%radius_ri_per_atom(iatom) + &
1238 102 : NORM2(center - particle_set(iatom)%r))
1239 : END DO
1240 :
1241 390 : DO igrid = 1, bs_env%ri_rs%n_grid_points
1242 1524 : IF (NORM2(ri_rs_grid_points(1:3, igrid) - center) > cutoff_ri) THEN
1243 0 : common_grid_available = .FALSE.
1244 0 : EXIT
1245 : END IF
1246 : END DO
1247 12 : IF (.NOT. common_grid_available) EXIT
1248 :
1249 48 : DO iatom = 1, natom
1250 144 : IF (NORM2(particle_set(iatom)%r - center) > &
1251 12 : bs_env%ri_rs%radius_ao_per_atom(iatom) + cutoff_ri) THEN
1252 0 : common_grid_available = .FALSE.
1253 0 : EXIT
1254 : END IF
1255 : END DO
1256 16 : IF (.NOT. common_grid_available) EXIT
1257 : END DO
1258 4 : DEALLOCATE (active_atom)
1259 :
1260 4 : CALL timestop(handle)
1261 :
1262 4 : END SUBROUTINE common_Z_lP_grid_available
1263 :
1264 : ! **************************************************************************************************
1265 : !> \brief Evaluates the AO collocation matrix on the complete RI-RS grid,
1266 : !>
1267 : !> ϕ_lμ = ϕ_μ(r_l), l = 1, ..., N_grid.
1268 : !>
1269 : !> "Complete common grid" means that the same full set of grid points is valid for every
1270 : !> optimized RI column block. The enclosing sphere is only an implementation device passed
1271 : !> to build_phi_on_sphere; it contains every r_l and every AO center, so its chosen center
1272 : !> does not select an AA or AB contraction.
1273 : !> \param bs_env ...
1274 : !> \param qs_kind_set ...
1275 : !> \param ri_rs_grid_points ...
1276 : !> \param n_ao_total ...
1277 : !> \param local_grid_idx ...
1278 : !> \param ngrid ...
1279 : !> \param phi_local ...
1280 : !> \param ao_col_map ...
1281 : !> \param n_ao_used ...
1282 : ! **************************************************************************************************
1283 8 : SUBROUTINE build_phi_on_complete_grid(bs_env, qs_kind_set, &
1284 4 : ri_rs_grid_points, n_ao_total, local_grid_idx, ngrid, &
1285 : phi_local, ao_col_map, n_ao_used)
1286 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1287 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1288 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ri_rs_grid_points
1289 : INTEGER, INTENT(IN) :: n_ao_total
1290 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: local_grid_idx
1291 : INTEGER, INTENT(OUT) :: ngrid
1292 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1293 : INTENT(OUT) :: phi_local
1294 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_col_map
1295 : INTEGER, INTENT(OUT) :: n_ao_used
1296 :
1297 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_phi_on_complete_grid'
1298 :
1299 : INTEGER :: handle, iatom, igrid, reference_atom
1300 : REAL(KIND=dp) :: cutoff_ri
1301 : REAL(KIND=dp), DIMENSION(3) :: center
1302 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1303 :
1304 4 : CALL timeset(routineN, handle)
1305 :
1306 4 : particle_set => bs_env%ri_rs%particle_set
1307 4 : reference_atom = 1
1308 16 : center = particle_set(reference_atom)%r
1309 4 : cutoff_ri = 0.0_dp
1310 130 : DO igrid = 1, bs_env%ri_rs%n_grid_points
1311 508 : cutoff_ri = MAX(cutoff_ri, NORM2(ri_rs_grid_points(1:3, igrid) - center))
1312 : END DO
1313 16 : DO iatom = 1, bs_env%n_atom
1314 52 : cutoff_ri = MAX(cutoff_ri, NORM2(particle_set(iatom)%r - center))
1315 : END DO
1316 4 : cutoff_ri = cutoff_ri + 1.0_dp
1317 :
1318 : CALL build_phi_on_sphere(bs_env, qs_kind_set, &
1319 : ri_rs_grid_points, reference_atom, cutoff_ri, n_ao_total, &
1320 : local_grid_idx, ngrid, phi_local, ao_col_map, n_ao_used, &
1321 4 : center=center)
1322 :
1323 4 : CALL timestop(handle)
1324 :
1325 4 : END SUBROUTINE build_phi_on_complete_grid
1326 :
1327 : ! **************************************************************************************************
1328 : !> \brief Writes Z_l,p0+p += z_lp into the process-owned grid-row blocks of one atomic column
1329 : !> block. The routine performs no explicit MPI communication; DBCSR summation combines
1330 : !> contributions when a block receives columns from more than one contraction group.
1331 : !> \param mat_Z_lP ...
1332 : !> \param z_block ...
1333 : !> \param local_grid_idx ...
1334 : !> \param n_local_grid ...
1335 : !> \param atom_index ...
1336 : !> \param first_column ...
1337 : !> \param atom_block_size ...
1338 : !> \param r_blk_sizes ...
1339 : !> \param row_offset ...
1340 : !> \param row_dist ...
1341 : !> \param myprow ...
1342 : !> \param eps_filter ...
1343 : ! **************************************************************************************************
1344 18 : SUBROUTINE add_Z_lP_columns(mat_Z_lP, z_block, local_grid_idx, n_local_grid, atom_index, &
1345 18 : first_column, atom_block_size, r_blk_sizes, row_offset, row_dist, &
1346 : myprow, eps_filter)
1347 :
1348 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_Z_lP
1349 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: z_block
1350 : INTEGER, DIMENSION(:), INTENT(IN) :: local_grid_idx
1351 : INTEGER, INTENT(IN) :: n_local_grid, atom_index, first_column, &
1352 : atom_block_size
1353 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes, row_offset, row_dist
1354 : INTEGER, INTENT(IN) :: myprow
1355 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1356 :
1357 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_Z_lP_columns'
1358 :
1359 : INTEGER :: current_chunk_size, g_pt, handle, i_blk, &
1360 : loc_ptr, ncolumn, r_end, r_start
1361 : LOGICAL :: row_owned
1362 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Z_blk
1363 :
1364 18 : CALL timeset(routineN, handle)
1365 :
1366 18 : ncolumn = SIZE(z_block, 2)
1367 18 : CPASSERT(first_column > 0)
1368 18 : CPASSERT(first_column + ncolumn - 1 <= atom_block_size)
1369 126 : ALLOCATE (Z_blk(MAXVAL(r_blk_sizes), atom_block_size), source=0.0_dp)
1370 18 : loc_ptr = 1
1371 72 : DO i_blk = 1, SIZE(r_blk_sizes)
1372 54 : r_start = row_offset(i_blk) + 1
1373 54 : r_end = row_offset(i_blk) + r_blk_sizes(i_blk)
1374 54 : current_chunk_size = r_blk_sizes(i_blk)
1375 54 : row_owned = row_dist(i_blk) == myprow
1376 54 : Z_blk = 0.0_dp
1377 684 : DO WHILE (loc_ptr <= n_local_grid)
1378 666 : g_pt = local_grid_idx(loc_ptr)
1379 666 : IF (g_pt > r_end) EXIT
1380 630 : IF (row_owned) THEN
1381 : Z_blk(g_pt - r_start + 1, first_column:first_column + ncolumn - 1) = &
1382 1302 : z_block(loc_ptr, 1:ncolumn)
1383 : END IF
1384 666 : loc_ptr = loc_ptr + 1
1385 : END DO
1386 3750 : IF (row_owned .AND. MAXVAL(ABS(Z_blk(1:current_chunk_size, :))) > eps_filter) THEN
1387 : CALL dbcsr_put_block(mat_Z_lP, row=i_blk, col=atom_index, &
1388 27 : block=Z_blk(1:current_chunk_size, :), summation=.TRUE.)
1389 : END IF
1390 : END DO
1391 18 : DEALLOCATE (Z_blk)
1392 :
1393 18 : CALL timestop(handle)
1394 :
1395 18 : END SUBROUTINE add_Z_lP_columns
1396 :
1397 : ! **************************************************************************************************
1398 : !> \brief Computes m_l^A = OR_{μ∈A}[ϕ_μ(r_l) /= 0]. The product
1399 : !>
1400 : !> ϕ_μ(r_l) ϕ_ν(r_l), μ∈A, ν∈B,
1401 : !>
1402 : !> is evaluated only where m_l^A AND m_l^B is true. This avoids products and contractions
1403 : !> whenever either atom has no nonzero AO on r_l; no additional numerical threshold is used.
1404 : !> \param bs_env ...
1405 : !> \param phi_val ...
1406 : !> \param ao_col_map ...
1407 : !> \param nonzero_ao ...
1408 : ! **************************************************************************************************
1409 6 : SUBROUTINE compute_nonzero_AO_grid_mask(bs_env, phi_val, ao_col_map, nonzero_ao)
1410 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1411 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val
1412 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
1413 : LOGICAL, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: nonzero_ao
1414 :
1415 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_nonzero_AO_grid_mask'
1416 :
1417 : INTEGER :: col, first, handle, iatom, last, natom, &
1418 : ngrid, point
1419 :
1420 6 : CALL timeset(routineN, handle)
1421 :
1422 6 : ngrid = SIZE(phi_val, 1)
1423 6 : natom = bs_env%n_atom
1424 24 : ALLOCATE (nonzero_ao(ngrid, natom), source=.FALSE.)
1425 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
1426 : !$OMP SHARED(bs_env, phi_val, ao_col_map, nonzero_ao, ngrid, natom) &
1427 6 : !$OMP PRIVATE(iatom, first, last, col, point)
1428 : DO iatom = 1, natom
1429 : first = ao_col_map(bs_env%i_ao_start_from_atom(iatom))
1430 : IF (first == 0) CYCLE
1431 : last = ao_col_map(bs_env%i_ao_end_from_atom(iatom))
1432 : DO col = first, last
1433 : DO point = 1, ngrid
1434 : nonzero_ao(point, iatom) = &
1435 : nonzero_ao(point, iatom) .OR. phi_val(point, col) /= 0.0_dp
1436 : END DO
1437 : END DO
1438 : END DO
1439 : !$OMP END PARALLEL DO
1440 :
1441 6 : CALL timestop(handle)
1442 6 : END SUBROUTINE compute_nonzero_AO_grid_mask
1443 :
1444 : ! **************************************************************************************************
1445 : !> \brief Computes d_lp = Σ_{μνP} ϕ_μ(r_l)ϕ_ν(r_l)(μν|P)U_Pp for all optimized
1446 : !> columns p that contain reference RI functions P on atom I.
1447 : !> \param bs_env ...
1448 : !> \param ctx ...
1449 : !> \param phi_val ...
1450 : !> \param ao_col_map ...
1451 : !> \param d_lp ...
1452 : !> \param n_grid ...
1453 : !> \param iatom ...
1454 : !> \param U_Pp ...
1455 : !> \param max_ao_size ...
1456 : ! **************************************************************************************************
1457 6 : SUBROUTINE compute_d_lp_auto_ri_atom(bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid, iatom, &
1458 6 : U_Pp, max_ao_size)
1459 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
1460 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1461 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
1462 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1463 : INTENT(IN) :: phi_val
1464 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
1465 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: d_lp
1466 : INTEGER, INTENT(IN) :: n_grid, iatom
1467 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1468 : INTENT(IN) :: U_Pp
1469 : INTEGER, INTENT(IN) :: max_ao_size
1470 :
1471 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp_auto_ri_atom'
1472 : INTEGER, PARAMETER :: grid_chunk = 1024
1473 :
1474 : INTEGER :: jatom, katom, c, handle, i_thread, j, jk_idx, jsize, jstart, k, ksize, &
1475 : kstart, l, l0, ncol, nRI_ref, nthreads, ri, thread_id
1476 : LOGICAL :: screened
1477 : REAL(KIND=dp) :: pair_factor
1478 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: int_2d_prv, rho_chunk
1479 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: d_lp_threads, int_3c_prv
1480 6 : TYPE(gw_3c_ws_type) :: ws
1481 :
1482 6 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: nonzero_ao
1483 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: grid_index
1484 : INTEGER :: n_grid_pair, grid_l, point
1485 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: grid_result
1486 :
1487 6 : CALL timeset(routineN, handle)
1488 6 : ncol = SIZE(U_Pp, 2)
1489 6 : nRI_ref = get_ref_RI_size(bs_env, iatom)
1490 6 : nthreads = 1
1491 6 : !$ nthreads = omp_get_max_threads()
1492 6 : CPASSERT(SIZE(d_lp, 1) == n_grid)
1493 6 : CPASSERT(SIZE(d_lp, 2) == ncol)
1494 6 : CPASSERT(SIZE(U_Pp, 1) == nRI_ref)
1495 30 : ALLOCATE (d_lp_threads(n_grid, ncol, nthreads), source=0.0_dp)
1496 :
1497 6 : CALL compute_nonzero_AO_grid_mask(bs_env, phi_val, ao_col_map, nonzero_ao)
1498 :
1499 : !$OMP PARALLEL DEFAULT(NONE) &
1500 : !$OMP SHARED(nonzero_ao, bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid, iatom, &
1501 : !$OMP U_Pp, &
1502 : !$OMP max_ao_size, ncol, nRI_ref, d_lp_threads, nthreads) &
1503 : !$OMP PRIVATE(grid_index, n_grid_pair, grid_l, point, grid_result, &
1504 : !$OMP jatom, katom, c, &
1505 : !$OMP i_thread, j, jk_idx, jsize, jstart, k, ksize, &
1506 : !$OMP kstart, l, l0, ri, screened, pair_factor, int_2d_prv, rho_chunk, &
1507 6 : !$OMP int_3c_prv, ws, thread_id)
1508 :
1509 : thread_id = 1
1510 : !$ thread_id = omp_get_thread_num() + 1
1511 :
1512 : CALL gw_3c_ws_create(ws, ctx)
1513 : ALLOCATE (int_3c_prv(max_ao_size, max_ao_size, ncol))
1514 : ALLOCATE (int_2d_prv(max_ao_size*max_ao_size, ncol))
1515 : ALLOCATE (rho_chunk(grid_chunk, max_ao_size*max_ao_size))
1516 : ALLOCATE (grid_index(n_grid), grid_result(grid_chunk, ncol))
1517 :
1518 : !$OMP DO SCHEDULE(DYNAMIC)
1519 : DO jatom = 1, bs_env%n_atom
1520 : DO katom = jatom, bs_env%n_atom
1521 : jstart = ao_col_map(bs_env%i_ao_start_from_atom(jatom))
1522 : kstart = ao_col_map(bs_env%i_ao_start_from_atom(katom))
1523 : IF (jstart == 0 .OR. kstart == 0) CYCLE
1524 : jsize = bs_env%i_ao_end_from_atom(jatom) - bs_env%i_ao_start_from_atom(jatom) + 1
1525 : ksize = bs_env%i_ao_end_from_atom(katom) - bs_env%i_ao_start_from_atom(katom) + 1
1526 : n_grid_pair = 0
1527 : DO grid_l = 1, n_grid
1528 : IF (.NOT. (nonzero_ao(grid_l, jatom) .AND. nonzero_ao(grid_l, katom))) CYCLE
1529 : n_grid_pair = n_grid_pair + 1
1530 : grid_index(n_grid_pair) = grid_l
1531 : END DO
1532 : IF (n_grid_pair == 0) CYCLE
1533 : int_3c_prv(1:jsize, 1:ksize, 1:ncol) = 0.0_dp
1534 : CALL build_3c_integral_block_auto_ri_ctx( &
1535 : int_3c_prv(1:jsize, 1:ksize, 1:ncol), ctx, ws, &
1536 : atom_j=jatom, atom_k=katom, atom_i=iatom, &
1537 : transform=U_Pp, transform_row=1, screened=screened)
1538 : IF (screened) CYCLE
1539 : DO ri = 1, ncol
1540 : DO k = 1, ksize
1541 : DO j = 1, jsize
1542 : jk_idx = (k - 1)*jsize + j
1543 : int_2d_prv(jk_idx, ri) = int_3c_prv(j, k, ri)
1544 : END DO
1545 : END DO
1546 : END DO
1547 : pair_factor = 1.0_dp
1548 : IF (jatom /= katom) pair_factor = 2.0_dp
1549 : DO l0 = 1, n_grid_pair, grid_chunk
1550 : c = MIN(grid_chunk, n_grid_pair - l0 + 1)
1551 : DO k = 1, ksize
1552 : DO j = 1, jsize
1553 : jk_idx = (k - 1)*jsize + j
1554 : DO l = 1, c
1555 : point = grid_index(l0 + l - 1)
1556 : rho_chunk(l, jk_idx) = phi_val(point, jstart + j - 1)* &
1557 : phi_val(point, kstart + k - 1)
1558 : END DO
1559 : END DO
1560 : END DO
1561 : CALL dgemm('N', 'N', c, ncol, jsize*ksize, pair_factor, rho_chunk, grid_chunk, &
1562 : int_2d_prv, max_ao_size*max_ao_size, 0.0_dp, &
1563 : grid_result, grid_chunk)
1564 : DO ri = 1, ncol
1565 : DO l = 1, c
1566 : point = grid_index(l0 + l - 1)
1567 : d_lp_threads(point, ri, thread_id) = &
1568 : d_lp_threads(point, ri, thread_id) + grid_result(l, ri)
1569 : END DO
1570 : END DO
1571 : END DO
1572 : END DO
1573 : END DO
1574 : !$OMP END DO
1575 :
1576 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1577 : DO ri = 1, ncol
1578 : DO l = 1, n_grid
1579 : DO i_thread = 1, nthreads
1580 : d_lp(l, ri) = d_lp(l, ri) + d_lp_threads(l, ri, i_thread)
1581 : END DO
1582 : END DO
1583 : END DO
1584 : !$OMP END DO
1585 : DEALLOCATE (int_3c_prv, int_2d_prv, rho_chunk)
1586 : DEALLOCATE (grid_index, grid_result)
1587 : CALL gw_3c_ws_release(ws)
1588 : !$OMP END PARALLEL
1589 :
1590 6 : DEALLOCATE (d_lp_threads)
1591 :
1592 6 : DEALLOCATE (nonzero_ao)
1593 :
1594 6 : CALL timestop(handle)
1595 :
1596 12 : END SUBROUTINE compute_d_lp_auto_ri_atom
1597 :
1598 : ! **************************************************************************************************
1599 : !> \brief Computes d_lp = Σ_{μνP} ϕ_μ(r_l)ϕ_ν(r_l)(μν|P)U_Pp for all AA and AB
1600 : !> columns p of one RI-RS matrix block. P can belong to either atom of an AB contraction.
1601 : !> \param bs_env ...
1602 : !> \param ctx ...
1603 : !> \param phi_val ...
1604 : !> \param ao_col_map ...
1605 : !> \param d_lp ...
1606 : !> \param n_grid ...
1607 : !> \param ri_coefficients ...
1608 : !> \param active_atom ...
1609 : !> \param max_ao_size ...
1610 : !> \param atom_j_mepos ...
1611 : !> \param atom_j_stride ...
1612 : ! **************************************************************************************************
1613 0 : SUBROUTINE compute_d_lp_auto_ri_atoms(bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid, &
1614 0 : ri_coefficients, active_atom, max_ao_size, atom_j_mepos, &
1615 : atom_j_stride)
1616 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
1617 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1618 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
1619 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1620 : INTENT(IN) :: phi_val
1621 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
1622 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: d_lp
1623 : INTEGER, INTENT(IN) :: n_grid
1624 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: ri_coefficients
1625 : LOGICAL, DIMENSION(:), INTENT(IN) :: active_atom
1626 : INTEGER, INTENT(IN) :: max_ao_size, atom_j_mepos, atom_j_stride
1627 :
1628 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp_auto_ri_atoms'
1629 : INTEGER, PARAMETER :: grid_chunk = 1024
1630 :
1631 : INTEGER :: active_column, iatom, jatom, katom, c, handle, i_thread, j, jk_idx, jsize, &
1632 : jstart, k, ksize, kstart, l, l0, &
1633 : max_active, nactive, ncol, nRI_ref, nthreads, ri, thread_id
1634 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_ncol
1635 0 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: atom_column
1636 : LOGICAL :: any_integral, screened
1637 : REAL(KIND=dp) :: pair_factor
1638 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: int_2d_prv, rho_chunk
1639 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coefficient_compact, d_lp_threads, &
1640 0 : int_3c_prv, int_3c_atom
1641 0 : TYPE(gw_3c_ws_type) :: ws
1642 :
1643 0 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: nonzero_ao
1644 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: grid_index
1645 : INTEGER :: n_grid_pair, grid_l, point
1646 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: grid_result
1647 :
1648 0 : CALL timeset(routineN, handle)
1649 0 : ncol = SIZE(d_lp, 2)
1650 0 : nthreads = 1
1651 0 : !$ nthreads = omp_get_max_threads()
1652 0 : CPASSERT(SIZE(d_lp, 1) == n_grid)
1653 0 : CPASSERT(SIZE(ri_coefficients, 2) == ncol)
1654 0 : CPASSERT(SIZE(ri_coefficients, 3) == SIZE(active_atom))
1655 0 : ALLOCATE (atom_ncol(SIZE(active_atom)), atom_column(ncol, SIZE(active_atom)))
1656 0 : atom_ncol = 0
1657 0 : atom_column = 0
1658 0 : DO iatom = 1, SIZE(active_atom)
1659 0 : IF (.NOT. active_atom(iatom)) CYCLE
1660 0 : nRI_ref = get_ref_RI_size(bs_env, iatom)
1661 0 : DO ri = 1, ncol
1662 0 : IF (.NOT. ANY(ri_coefficients(1:nRI_ref, ri, iatom) /= 0.0_dp)) CYCLE
1663 0 : atom_ncol(iatom) = atom_ncol(iatom) + 1
1664 0 : atom_column(atom_ncol(iatom), iatom) = ri
1665 : END DO
1666 : END DO
1667 0 : max_active = MAXVAL(atom_ncol)
1668 0 : CPASSERT(max_active > 0)
1669 0 : ALLOCATE (coefficient_compact(SIZE(ri_coefficients, 1), max_active, &
1670 0 : SIZE(active_atom)), source=0.0_dp)
1671 0 : DO iatom = 1, SIZE(active_atom)
1672 0 : nRI_ref = get_ref_RI_size(bs_env, iatom)
1673 0 : DO active_column = 1, atom_ncol(iatom)
1674 0 : ri = atom_column(active_column, iatom)
1675 : coefficient_compact(1:nRI_ref, active_column, iatom) = &
1676 0 : ri_coefficients(1:nRI_ref, ri, iatom)
1677 : END DO
1678 : END DO
1679 0 : ALLOCATE (d_lp_threads(n_grid, ncol, nthreads), source=0.0_dp)
1680 :
1681 0 : CALL compute_nonzero_AO_grid_mask(bs_env, phi_val, ao_col_map, nonzero_ao)
1682 :
1683 : !$OMP PARALLEL DEFAULT(NONE) &
1684 : !$OMP SHARED(nonzero_ao, bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid, ri_coefficients, &
1685 : !$OMP active_atom, max_ao_size, atom_j_mepos, atom_j_stride, ncol, &
1686 : !$OMP d_lp_threads, nthreads, atom_ncol, &
1687 : !$OMP atom_column, coefficient_compact, max_active) &
1688 : !$OMP PRIVATE(grid_index, n_grid_pair, grid_l, point, grid_result, &
1689 : !$OMP active_column, iatom, jatom, katom, c, i_thread, j, jk_idx, jsize, jstart, &
1690 : !$OMP k, ksize, kstart, l, l0, &
1691 : !$OMP nactive, nRI_ref, ri, any_integral, screened, pair_factor, &
1692 : !$OMP int_2d_prv, rho_chunk, &
1693 0 : !$OMP int_3c_prv, int_3c_atom, ws, thread_id)
1694 :
1695 : thread_id = 1
1696 : !$ thread_id = omp_get_thread_num() + 1
1697 :
1698 : CALL gw_3c_ws_create(ws, ctx)
1699 : ALLOCATE (int_3c_prv(max_ao_size, max_ao_size, ncol))
1700 : ALLOCATE (int_3c_atom(max_ao_size, max_ao_size, max_active))
1701 : ALLOCATE (int_2d_prv(max_ao_size*max_ao_size, ncol))
1702 : ALLOCATE (rho_chunk(grid_chunk, max_ao_size*max_ao_size))
1703 : ALLOCATE (grid_index(n_grid), grid_result(grid_chunk, ncol))
1704 :
1705 : !$OMP DO SCHEDULE(DYNAMIC)
1706 : DO jatom = atom_j_mepos + 1, bs_env%n_atom, atom_j_stride
1707 : DO katom = jatom, bs_env%n_atom
1708 : jstart = ao_col_map(bs_env%i_ao_start_from_atom(jatom))
1709 : kstart = ao_col_map(bs_env%i_ao_start_from_atom(katom))
1710 : IF (jstart == 0 .OR. kstart == 0) CYCLE
1711 : jsize = bs_env%i_ao_end_from_atom(jatom) - bs_env%i_ao_start_from_atom(jatom) + 1
1712 : ksize = bs_env%i_ao_end_from_atom(katom) - bs_env%i_ao_start_from_atom(katom) + 1
1713 : n_grid_pair = 0
1714 : DO grid_l = 1, n_grid
1715 : IF (.NOT. (nonzero_ao(grid_l, jatom) .AND. nonzero_ao(grid_l, katom))) CYCLE
1716 : n_grid_pair = n_grid_pair + 1
1717 : grid_index(n_grid_pair) = grid_l
1718 : END DO
1719 : IF (n_grid_pair == 0) CYCLE
1720 : int_3c_prv(1:jsize, 1:ksize, 1:ncol) = 0.0_dp
1721 : any_integral = .FALSE.
1722 : DO iatom = 1, SIZE(active_atom)
1723 : IF (.NOT. active_atom(iatom)) CYCLE
1724 : nRI_ref = get_ref_RI_size(bs_env, iatom)
1725 : nactive = atom_ncol(iatom)
1726 : int_3c_atom(1:jsize, 1:ksize, 1:nactive) = 0.0_dp
1727 : CALL build_3c_integral_block_auto_ri_ctx( &
1728 : int_3c_atom(1:jsize, 1:ksize, 1:nactive), ctx, ws, &
1729 : atom_j=jatom, atom_k=katom, atom_i=iatom, &
1730 : transform=coefficient_compact(1:nRI_ref, 1:nactive, iatom), &
1731 : transform_row=1, screened=screened)
1732 : IF (.NOT. screened) THEN
1733 : any_integral = .TRUE.
1734 : DO active_column = 1, nactive
1735 : ri = atom_column(active_column, iatom)
1736 : int_3c_prv(1:jsize, 1:ksize, ri) = &
1737 : int_3c_prv(1:jsize, 1:ksize, ri) + &
1738 : int_3c_atom(1:jsize, 1:ksize, active_column)
1739 : END DO
1740 : END IF
1741 : END DO
1742 : IF (.NOT. any_integral) CYCLE
1743 :
1744 : DO ri = 1, ncol
1745 : DO k = 1, ksize
1746 : DO j = 1, jsize
1747 : jk_idx = (k - 1)*jsize + j
1748 : int_2d_prv(jk_idx, ri) = int_3c_prv(j, k, ri)
1749 : END DO
1750 : END DO
1751 : END DO
1752 :
1753 : pair_factor = 1.0_dp
1754 : IF (jatom /= katom) pair_factor = 2.0_dp
1755 : DO l0 = 1, n_grid_pair, grid_chunk
1756 : c = MIN(grid_chunk, n_grid_pair - l0 + 1)
1757 : DO k = 1, ksize
1758 : DO j = 1, jsize
1759 : jk_idx = (k - 1)*jsize + j
1760 : DO l = 1, c
1761 : point = grid_index(l0 + l - 1)
1762 : rho_chunk(l, jk_idx) = phi_val(point, jstart + j - 1)* &
1763 : phi_val(point, kstart + k - 1)
1764 : END DO
1765 : END DO
1766 : END DO
1767 : CALL dgemm('N', 'N', c, ncol, jsize*ksize, pair_factor, rho_chunk, grid_chunk, &
1768 : int_2d_prv, max_ao_size*max_ao_size, 0.0_dp, &
1769 : grid_result, grid_chunk)
1770 : DO ri = 1, ncol
1771 : DO l = 1, c
1772 : point = grid_index(l0 + l - 1)
1773 : d_lp_threads(point, ri, thread_id) = &
1774 : d_lp_threads(point, ri, thread_id) + grid_result(l, ri)
1775 : END DO
1776 : END DO
1777 : END DO
1778 : END DO
1779 : END DO
1780 : !$OMP END DO
1781 :
1782 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1783 : DO ri = 1, ncol
1784 : DO l = 1, n_grid
1785 : DO i_thread = 1, nthreads
1786 : d_lp(l, ri) = d_lp(l, ri) + d_lp_threads(l, ri, i_thread)
1787 : END DO
1788 : END DO
1789 : END DO
1790 : !$OMP END DO
1791 : DEALLOCATE (int_3c_prv, int_3c_atom, int_2d_prv, rho_chunk)
1792 : DEALLOCATE (grid_index, grid_result)
1793 : CALL gw_3c_ws_release(ws)
1794 : !$OMP END PARALLEL
1795 :
1796 0 : DEALLOCATE (coefficient_compact, d_lp_threads, atom_column, atom_ncol)
1797 :
1798 0 : DEALLOCATE (nonzero_ao)
1799 :
1800 0 : CALL timestop(handle)
1801 :
1802 0 : END SUBROUTINE compute_d_lp_auto_ri_atoms
1803 :
1804 : ! **************************************************************************************************
1805 : !> \brief Computes d_lp = Σ_A Σ_{P∈A} d_lP U_Pp in rank-sized batches on a common grid.
1806 : !> \param bs_env ...
1807 : !> \param ctx ...
1808 : !> \param phi_val ...
1809 : !> \param ao_col_map ...
1810 : !> \param d_lp ...
1811 : !> \param n_grid ...
1812 : !> \param max_ao_size ...
1813 : !> \param mepos ...
1814 : !> \param num_pe ...
1815 : ! **************************************************************************************************
1816 4 : SUBROUTINE compute_d_lp_auto_ri_batch(bs_env, ctx, phi_val, ao_col_map, d_lp, n_grid, &
1817 : max_ao_size, mepos, num_pe)
1818 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1819 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
1820 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1821 : INTENT(IN) :: phi_val
1822 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_col_map
1823 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: d_lp
1824 : INTEGER, INTENT(IN) :: n_grid, max_ao_size, mepos, num_pe
1825 :
1826 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp_auto_ri_batch'
1827 :
1828 : INTEGER :: column, handle, iatom, n_done, n_total, &
1829 : ncol_batch
1830 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: global_map
1831 : REAL(KIND=dp) :: item_start_time
1832 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_lp_batch, U_Pp
1833 :
1834 4 : CALL timeset(routineN, handle)
1835 :
1836 4 : n_total = 0
1837 4 : DO iatom = mepos + 1, bs_env%n_atom, num_pe
1838 6 : n_total = n_total + 1
1839 : END DO
1840 4 : n_done = 0
1841 :
1842 10 : DO iatom = mepos + 1, bs_env%n_atom, num_pe
1843 6 : item_start_time = m_walltime()
1844 6 : CALL collect_auto_ri_columns_for_atom(bs_env, iatom, U_Pp, global_map)
1845 6 : ncol_batch = SIZE(global_map)
1846 6 : IF (ncol_batch == 0) THEN
1847 0 : DEALLOCATE (U_Pp, global_map)
1848 0 : n_done = n_done + 1
1849 : CALL print_Z_lP_progress(bs_env, n_done, n_total, &
1850 0 : m_walltime() - item_start_time)
1851 0 : CYCLE
1852 : END IF
1853 :
1854 24 : ALLOCATE (d_lp_batch(n_grid, ncol_batch), source=0.0_dp)
1855 : CALL compute_d_lp_auto_ri_atom(bs_env, ctx, phi_val, ao_col_map, d_lp_batch, &
1856 6 : n_grid, iatom, U_Pp, max_ao_size)
1857 44 : DO column = 1, ncol_batch
1858 1409 : d_lp(:, global_map(column)) = d_lp(:, global_map(column)) + d_lp_batch(:, column)
1859 : END DO
1860 6 : DEALLOCATE (U_Pp, global_map, d_lp_batch)
1861 6 : n_done = n_done + 1
1862 : CALL print_Z_lP_progress(bs_env, n_done, n_total, &
1863 10 : m_walltime() - item_start_time)
1864 : END DO
1865 4 : CALL timestop(handle)
1866 :
1867 4 : END SUBROUTINE compute_d_lp_auto_ri_batch
1868 :
1869 : ! **************************************************************************************************
1870 : !> \brief Collects every optimized column p containing reference functions P on atom A. It returns
1871 : !>
1872 : !> U_Pp^A
1873 : !>
1874 : !> together with the global optimized-column index of q.
1875 : !> \param bs_env ...
1876 : !> \param iatom ...
1877 : !> \param U_Pp ...
1878 : !> \param global_map ...
1879 : ! **************************************************************************************************
1880 6 : SUBROUTINE collect_auto_ri_columns_for_atom(bs_env, iatom, U_Pp, global_map)
1881 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1882 : INTEGER, INTENT(IN) :: iatom
1883 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: U_Pp(:, :)
1884 : INTEGER, ALLOCATABLE, INTENT(OUT) :: global_map(:)
1885 :
1886 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collect_auto_ri_columns_for_atom'
1887 :
1888 : INTEGER :: AB_block, atom_A, atom_B, column, global_column, handle, local_column, ncol, &
1889 : ncol_batch, nRI_ref, nRI_ref_A, output_offset, row_first
1890 6 : REAL(KIND=dp), ALLOCATABLE :: U_Pp_AB(:, :)
1891 :
1892 6 : CALL timeset(routineN, handle)
1893 :
1894 6 : ncol_batch = 0
1895 33 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
1896 27 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
1897 27 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
1898 33 : IF (iatom == atom_A .OR. iatom == atom_B) THEN
1899 12 : ncol_batch = ncol_batch + bs_env%auto_ri%AB_size_opt_RI(AB_block)
1900 : END IF
1901 : END DO
1902 :
1903 6 : nRI_ref = get_ref_RI_size(bs_env, iatom)
1904 24 : ALLOCATE (U_Pp(nRI_ref, ncol_batch), source=0.0_dp)
1905 18 : ALLOCATE (global_map(ncol_batch))
1906 6 : column = 0
1907 33 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
1908 27 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
1909 27 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
1910 27 : IF (iatom /= atom_A .AND. iatom /= atom_B) CYCLE
1911 12 : ncol = bs_env%auto_ri%AB_size_opt_RI(AB_block)
1912 12 : CALL get_U_Pp_AB(bs_env%auto_ri, AB_block, U_Pp_AB)
1913 12 : row_first = 1
1914 12 : IF (iatom == atom_B .AND. atom_B /= atom_A) THEN
1915 3 : nRI_ref_A = get_ref_RI_size(bs_env, atom_A)
1916 3 : row_first = 1 + nRI_ref_A
1917 : END IF
1918 : U_Pp(:, column + 1:column + ncol) = &
1919 : U_Pp_AB( &
1920 3892 : row_first:row_first + nRI_ref - 1, 1:ncol)
1921 :
1922 50 : DO local_column = 1, ncol
1923 38 : IF (local_column <= bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)) THEN
1924 : output_offset = SUM(bs_env%auto_ri%sizes_opt_RI(:atom_A - 1)) + &
1925 30 : bs_env%auto_ri%AB_first_p_A(AB_block) - 1
1926 22 : global_column = output_offset + local_column
1927 : ELSE
1928 16 : CPASSERT(atom_B /= atom_A)
1929 : output_offset = SUM(bs_env%auto_ri%sizes_opt_RI(:atom_B - 1)) + &
1930 40 : bs_env%auto_ri%AB_first_p_B(AB_block) - 1
1931 : global_column = output_offset + local_column - &
1932 16 : bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
1933 : END IF
1934 50 : global_map(column + local_column) = global_column
1935 : END DO
1936 12 : column = column + ncol
1937 33 : DEALLOCATE (U_Pp_AB)
1938 : END DO
1939 6 : CPASSERT(column == ncol_batch)
1940 :
1941 6 : CALL timestop(handle)
1942 :
1943 12 : END SUBROUTINE collect_auto_ri_columns_for_atom
1944 :
1945 : ! **************************************************************************************************
1946 : !> \brief Unpacks one AB contraction matrix U_Pp; B=A denotes an AA block.
1947 : !> \param auto_ri ...
1948 : !> \param AB_block ...
1949 : !> \param U_Pp ...
1950 : ! **************************************************************************************************
1951 12 : SUBROUTINE get_U_Pp_AB(auto_ri, AB_block, U_Pp)
1952 : TYPE(auto_ri_type), INTENT(IN) :: auto_ri
1953 : INTEGER, INTENT(IN) :: AB_block
1954 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: U_Pp(:, :)
1955 :
1956 : INTEGER :: first, last, ncolumn, nrow
1957 :
1958 12 : nrow = auto_ri%AB_size_ref_RI(AB_block)
1959 12 : ncolumn = auto_ri%AB_size_opt_RI(AB_block)
1960 12 : first = auto_ri%U_Pp_AB_offset(AB_block)
1961 12 : last = first + nrow*ncolumn - 1
1962 48 : ALLOCATE (U_Pp(nrow, ncolumn))
1963 36 : U_Pp(:, :) = RESHAPE(auto_ri%U_Pp_AB(first:last), [nrow, ncolumn])
1964 :
1965 12 : END SUBROUTINE get_U_Pp_AB
1966 :
1967 : ! **************************************************************************************************
1968 : !> \brief Computes the fitting radius for every optimized atomic column block. If q assigned to A
1969 : !> contains reference functions on B, the required radius is
1970 : !>
1971 : !> R_A^fit = max_B [r_c + R_B^RI + |R_A - R_B|].
1972 : !> \param bs_env ...
1973 : !> \param radius ...
1974 : ! **************************************************************************************************
1975 0 : SUBROUTINE compute_auto_ri_grid_radii(bs_env, radius)
1976 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1977 : REAL(KIND=dp), INTENT(OUT) :: radius(:)
1978 :
1979 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_auto_ri_grid_radii'
1980 :
1981 : INTEGER :: a, AB_block, b, handle, n_to_a, ncol
1982 : REAL(KIND=dp) :: cutoff, distance
1983 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1984 :
1985 0 : CALL timeset(routineN, handle)
1986 :
1987 0 : particle_set => bs_env%ri_rs%particle_set
1988 0 : radius = 0.0_dp
1989 0 : cutoff = bs_env%ri_metric%cutoff_radius
1990 0 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
1991 0 : a = bs_env%auto_ri%AB_atom_A(AB_block)
1992 0 : b = bs_env%auto_ri%AB_atom_B(AB_block)
1993 0 : n_to_a = bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
1994 0 : ncol = bs_env%auto_ri%AB_size_opt_RI(AB_block)
1995 0 : IF (n_to_a > 0) THEN
1996 0 : radius(a) = MAX(radius(a), cutoff + bs_env%ri_rs%radius_ri_per_atom(a))
1997 0 : IF (b /= a) THEN
1998 0 : distance = NORM2(particle_set(a)%r - particle_set(b)%r)
1999 0 : radius(a) = MAX(radius(a), cutoff + bs_env%ri_rs%radius_ri_per_atom(b) + distance)
2000 : END IF
2001 : END IF
2002 0 : IF (n_to_a < ncol) THEN
2003 0 : CPASSERT(b /= a)
2004 0 : distance = NORM2(particle_set(b)%r - particle_set(a)%r)
2005 : radius(b) = MAX(radius(b), cutoff + bs_env%ri_rs%radius_ri_per_atom(b), &
2006 0 : cutoff + bs_env%ri_rs%radius_ri_per_atom(a) + distance)
2007 : END IF
2008 : END DO
2009 :
2010 0 : CALL timestop(handle)
2011 0 : END SUBROUTINE compute_auto_ri_grid_radii
2012 :
2013 : ! **************************************************************************************************
2014 : !> \brief Computes
2015 : !>
2016 : !> d_lp = Σ_A Σ_{P∈A} d_lP U_Pp^A,
2017 : !> d_lP = Σ_μν ϕ_μ(r_l) ϕ_ν(r_l) (μν|P),
2018 : !>
2019 : !> on the union of the atomic RI-RS grids needed by q. Each atom's three-center integrals are
2020 : !> evaluated once and transformed with U_Pp^A.
2021 : !> \param qs_env ...
2022 : !> \param bs_env ...
2023 : !> \param ctx ...
2024 : !> \param grid ...
2025 : !> \param mat_phi ...
2026 : !> \param mat_rhs ...
2027 : !> \param max_ao_size ...
2028 : ! **************************************************************************************************
2029 0 : SUBROUTINE compute_auto_ri_d_lp(qs_env, bs_env, ctx, grid, mat_phi, mat_rhs, max_ao_size)
2030 : TYPE(qs_environment_type), POINTER :: qs_env
2031 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2032 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
2033 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: grid
2034 : TYPE(dbcsr_type), INTENT(IN) :: mat_phi
2035 : TYPE(dbcsr_type), INTENT(OUT) :: mat_rhs
2036 : INTEGER, INTENT(IN) :: max_ao_size
2037 :
2038 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_auto_ri_d_lp'
2039 :
2040 : INTEGER :: AB_block, atom_A, atom_B, first, fit_atom, handle, handle_project, handle_rhs, l, &
2041 : last, n_done, n_first_p_ABs, n_to_a, n_total, n_union, nao, natom, ncol, ngrid, npcol, &
2042 : nprow, ri_atom
2043 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_map, global_map, local_index, &
2044 0 : row_offset, union_index
2045 0 : INTEGER, DIMENSION(:), POINTER :: AB_row_dist, col_dist, &
2046 0 : first_p_ABs_per_atom, retained_size, &
2047 0 : row_dist, row_size
2048 0 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: needed_fit_atom, union_mask
2049 : REAL(KIND=dp) :: item_start_time, radius
2050 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fit_radius
2051 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: column_map, phi, rhs, U_Pp, union_grid
2052 : TYPE(cell_type), POINTER :: cell
2053 : TYPE(dbcsr_distribution_type) :: dist_AB, dist_phi, dist_t
2054 : TYPE(dbcsr_type) :: AB_d_lp, transform
2055 : TYPE(mp_para_env_type), POINTER :: para_env
2056 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2057 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2058 :
2059 0 : CALL timeset(routineN, handle)
2060 : CALL get_qs_env(qs_env, para_env=para_env, particle_set=particle_set, &
2061 0 : qs_kind_set=qs_kind_set, cell=cell)
2062 0 : natom = bs_env%n_atom
2063 0 : CALL dbcsr_get_info(mat_phi, row_blk_size=row_size, distribution=dist_phi)
2064 0 : CALL dbcsr_distribution_get(dist_phi, row_dist=row_dist, nprows=nprow, npcols=npcol)
2065 : ALLOCATE (first_p_ABs_per_atom(natom), col_dist(natom), AB_row_dist(natom), &
2066 0 : retained_size(natom), row_offset(SIZE(row_size)))
2067 0 : retained_size(:) = bs_env%auto_ri%sizes_opt_RI
2068 0 : DO ri_atom = 1, natom
2069 0 : first_p_ABs_per_atom(ri_atom) = 0
2070 0 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
2071 0 : IF (ri_atom == bs_env%auto_ri%AB_atom_A(AB_block) .OR. &
2072 0 : ri_atom == bs_env%auto_ri%AB_atom_B(AB_block)) THEN
2073 : first_p_ABs_per_atom(ri_atom) = &
2074 : first_p_ABs_per_atom(ri_atom) + &
2075 0 : bs_env%auto_ri%AB_size_opt_RI(AB_block)
2076 : END IF
2077 : END DO
2078 0 : col_dist(ri_atom) = MOD(ri_atom - 1, npcol)
2079 0 : AB_row_dist(ri_atom) = MOD(ri_atom - 1, nprow)
2080 : END DO
2081 0 : row_offset(1) = 0
2082 0 : DO l = 2, SIZE(row_size)
2083 0 : row_offset(l) = row_offset(l - 1) + row_size(l - 1)
2084 : END DO
2085 : CALL dbcsr_distribution_new(dist_AB, template=dist_phi, &
2086 0 : row_dist=row_dist, col_dist=col_dist)
2087 : CALL dbcsr_create(AB_d_lp, name='AUTO_RI AA/AB d_lp', dist=dist_AB, &
2088 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_size, &
2089 0 : col_blk_size=first_p_ABs_per_atom)
2090 : CALL dbcsr_distribution_new(dist_t, template=dist_phi, &
2091 0 : row_dist=AB_row_dist, col_dist=col_dist)
2092 : CALL dbcsr_create(transform, name='AUTO_RI U_Pp', dist=dist_t, &
2093 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=first_p_ABs_per_atom, &
2094 0 : col_blk_size=retained_size)
2095 : CALL dbcsr_create(mat_rhs, name='AUTO_RI optimized d_lp', dist=dist_AB, &
2096 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_size, &
2097 0 : col_blk_size=retained_size)
2098 0 : ALLOCATE (fit_radius(natom))
2099 0 : CALL compute_auto_ri_grid_radii(bs_env, fit_radius)
2100 0 : ALLOCATE (needed_fit_atom(natom), union_mask(SIZE(grid, 2)))
2101 0 : n_total = 0
2102 0 : DO ri_atom = para_env%mepos + 1, natom, para_env%num_pe
2103 0 : n_total = n_total + 1
2104 : END DO
2105 0 : n_done = 0
2106 0 : CALL timeset(routineN//'_AB_d_lp', handle_rhs)
2107 0 : DO ri_atom = para_env%mepos + 1, natom, para_env%num_pe
2108 0 : item_start_time = m_walltime()
2109 0 : needed_fit_atom = .FALSE.
2110 0 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
2111 0 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
2112 0 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
2113 0 : IF (ri_atom /= atom_A .AND. ri_atom /= atom_B) CYCLE
2114 0 : n_to_a = bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
2115 0 : IF (n_to_a > 0) needed_fit_atom(atom_A) = .TRUE.
2116 0 : IF (n_to_a < bs_env%auto_ri%AB_size_opt_RI(AB_block)) THEN
2117 0 : needed_fit_atom(atom_B) = .TRUE.
2118 : END IF
2119 : END DO
2120 0 : IF (.NOT. ANY(needed_fit_atom)) THEN
2121 0 : n_done = n_done + 1
2122 : CALL print_Z_lP_progress(bs_env, n_done, n_total, &
2123 0 : m_walltime() - item_start_time)
2124 0 : CYCLE
2125 : END IF
2126 0 : union_mask = .FALSE.
2127 0 : radius = 0.0_dp
2128 0 : DO fit_atom = 1, natom
2129 0 : IF (.NOT. needed_fit_atom(fit_atom)) CYCLE
2130 : radius = MAX(radius, fit_radius(fit_atom) + &
2131 0 : NORM2(particle_set(fit_atom)%r - particle_set(ri_atom)%r))
2132 0 : DO l = 1, SIZE(grid, 2)
2133 0 : IF (NORM2(grid(:, l) - particle_set(fit_atom)%r) <= &
2134 0 : fit_radius(fit_atom)) union_mask(l) = .TRUE.
2135 : END DO
2136 : END DO
2137 0 : n_union = COUNT(union_mask)
2138 0 : ALLOCATE (union_index(n_union), union_grid(3, n_union))
2139 0 : n_union = 0
2140 0 : DO l = 1, SIZE(grid, 2)
2141 0 : IF (.NOT. union_mask(l)) CYCLE
2142 0 : n_union = n_union + 1
2143 0 : union_index(n_union) = l
2144 0 : union_grid(:, n_union) = grid(:, l)
2145 : END DO
2146 : CALL build_phi_on_sphere(bs_env, qs_kind_set, union_grid, ri_atom, &
2147 : radius + 1.0_dp, bs_env%i_ao_end_from_atom(natom), local_index, &
2148 0 : ngrid, phi, ao_map, nao)
2149 0 : local_index(1:ngrid) = union_index(local_index(1:ngrid))
2150 0 : n_first_p_ABs = first_p_ABs_per_atom(ri_atom)
2151 0 : CALL collect_auto_ri_columns_for_atom(bs_env, ri_atom, U_Pp, global_map)
2152 0 : CPASSERT(SIZE(global_map) == n_first_p_ABs)
2153 0 : first = 1
2154 0 : DO fit_atom = 1, natom
2155 0 : ncol = retained_size(fit_atom)
2156 0 : last = first + ncol - 1
2157 0 : IF (ANY(global_map >= first .AND. global_map <= last)) THEN
2158 0 : ALLOCATE (column_map(n_first_p_ABs, ncol), source=0.0_dp)
2159 0 : DO l = 1, n_first_p_ABs
2160 0 : IF (global_map(l) < first .OR. global_map(l) > last) CYCLE
2161 0 : column_map(l, global_map(l) - first + 1) = 1.0_dp
2162 : END DO
2163 0 : CALL dbcsr_put_block(transform, ri_atom, fit_atom, column_map)
2164 0 : DEALLOCATE (column_map)
2165 : END IF
2166 0 : first = last + 1
2167 : END DO
2168 0 : ALLOCATE (rhs(ngrid, n_first_p_ABs), source=0.0_dp)
2169 : CALL compute_d_lp_auto_ri_atom(bs_env, ctx, phi, ao_map, rhs, ngrid, &
2170 0 : ri_atom, U_Pp, max_ao_size)
2171 : CALL store_Z_lP_columns(AB_d_lp, rhs, local_index, ngrid, n_first_p_ABs, ri_atom, &
2172 0 : row_size, row_offset, 0.0_dp)
2173 0 : DEALLOCATE (rhs, U_Pp, global_map, phi, ao_map, local_index, &
2174 0 : union_grid, union_index)
2175 0 : n_done = n_done + 1
2176 : CALL print_Z_lP_progress(bs_env, n_done, n_total, &
2177 0 : m_walltime() - item_start_time)
2178 : END DO
2179 0 : CALL dbcsr_finalize(transform)
2180 0 : CALL dbcsr_finalize(AB_d_lp)
2181 0 : CALL timestop(handle_rhs)
2182 0 : CALL timeset(routineN//'_transform', handle_project)
2183 : CALL dbcsr_multiply('N', 'N', 1.0_dp, AB_d_lp, transform, 0.0_dp, &
2184 0 : mat_rhs, filter_eps=0.0_dp)
2185 0 : CALL timestop(handle_project)
2186 0 : CALL dbcsr_release(AB_d_lp)
2187 0 : CALL dbcsr_release(transform)
2188 0 : CALL dbcsr_distribution_release(dist_AB)
2189 0 : CALL dbcsr_distribution_release(dist_t)
2190 0 : DEALLOCATE (first_p_ABs_per_atom, retained_size, col_dist, AB_row_dist, row_offset, &
2191 0 : fit_radius, needed_fit_atom, union_mask)
2192 0 : CALL timestop(handle)
2193 0 : END SUBROUTINE compute_auto_ri_d_lp
2194 :
2195 : ! **************************************************************************************************
2196 : !> \brief Extracts rhs(i,q) = d_{local_index(i),q} from one optimized atomic DBCSR column block.
2197 : !> \param mat_rhs ...
2198 : !> \param atom_index ...
2199 : !> \param local_index ...
2200 : !> \param row_offset ...
2201 : !> \param rhs ...
2202 : ! **************************************************************************************************
2203 0 : SUBROUTINE extract_atom_d_lp(mat_rhs, atom_index, local_index, row_offset, rhs)
2204 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_rhs
2205 : INTEGER, INTENT(IN) :: atom_index
2206 : INTEGER, DIMENSION(:), INTENT(IN) :: local_index, row_offset
2207 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rhs
2208 :
2209 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extract_atom_d_lp'
2210 :
2211 : INTEGER :: handle, l, next_row, row
2212 : LOGICAL :: found
2213 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
2214 :
2215 0 : CALL timeset(routineN, handle)
2216 :
2217 0 : rhs = 0.0_dp
2218 0 : row = 0
2219 0 : next_row = 1
2220 0 : NULLIFY (block)
2221 0 : DO l = 1, SIZE(rhs, 1)
2222 0 : DO WHILE (next_row <= SIZE(row_offset))
2223 0 : IF (row_offset(next_row) >= local_index(l)) EXIT
2224 0 : row = next_row
2225 0 : next_row = next_row + 1
2226 0 : CALL dbcsr_get_block_p(mat_rhs, row, atom_index, block, found)
2227 0 : IF (.NOT. found) NULLIFY (block)
2228 : END DO
2229 0 : IF (ASSOCIATED(block)) rhs(l, :) = block(local_index(l) - row_offset(row), :)
2230 : END DO
2231 :
2232 0 : CALL timestop(handle)
2233 0 : END SUBROUTINE extract_atom_d_lp
2234 :
2235 : ! **************************************************************************************************
2236 : !> \brief Solves Σ_l' D_ll' Z_l'q = d_lp for each optimized atomic column block. Small grids are
2237 : !> gathered within one process column and solved by a local Cholesky factorization.
2238 : !> \param qs_env ...
2239 : !> \param bs_env ...
2240 : !> \param grid ...
2241 : !> \param mat_rhs ...
2242 : !> \param mat_z ...
2243 : ! **************************************************************************************************
2244 0 : SUBROUTINE fit_auto_ri_z_lp(qs_env, bs_env, grid, mat_rhs, mat_z)
2245 : TYPE(qs_environment_type), POINTER :: qs_env
2246 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2247 : REAL(KIND=dp), INTENT(IN) :: grid(:, :)
2248 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_rhs, mat_z
2249 :
2250 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_auto_ri_z_lp'
2251 :
2252 : INTEGER :: atom_index, base, group_size, handle, &
2253 : handle_gather, i, info, mypcol, n_big, &
2254 : n_small, nao, natom, ncol, ngrid, &
2255 : npcol, slot
2256 0 : INTEGER, ALLOCATABLE :: all_index(:), ao_map(:), big_list(:), &
2257 0 : grid_size(:), local_index(:), &
2258 0 : row_offset(:), small_list(:)
2259 0 : INTEGER, POINTER :: row_size(:)
2260 0 : LOGICAL, ALLOCATABLE :: single_rank(:)
2261 0 : REAL(KIND=dp), ALLOCATABLE :: atom_rhs(:, :), buffer(:, :), &
2262 0 : diagonal(:), gram(:, :), &
2263 0 : local_rhs(:, :), phi(:, :), radius(:)
2264 : TYPE(dbcsr_distribution_type) :: distribution
2265 : TYPE(mp_para_env_type), POINTER :: column_env, para_env
2266 0 : TYPE(qs_kind_type), POINTER :: qs_kinds(:)
2267 :
2268 0 : CALL timeset(routineN, handle)
2269 0 : CALL get_qs_env(qs_env, para_env=para_env, qs_kind_set=qs_kinds)
2270 0 : CALL dbcsr_get_info(mat_rhs, distribution=distribution, row_blk_size=row_size)
2271 0 : CALL dbcsr_distribution_get(distribution, mypcol=mypcol, npcols=npcol)
2272 0 : ALLOCATE (column_env)
2273 0 : CALL column_env%from_split(para_env, mypcol)
2274 0 : natom = bs_env%n_atom
2275 0 : ALLOCATE (radius(natom), all_index(SIZE(grid, 2)), row_offset(SIZE(row_size)))
2276 0 : CALL compute_auto_ri_grid_radii(bs_env, radius)
2277 : CALL classify_z_lp_atoms(bs_env, grid, radius, &
2278 : bs_env%auto_ri%sizes_opt_RI, grid_size, &
2279 0 : small_list, n_small, big_list, n_big, group_size)
2280 0 : ALLOCATE (single_rank(natom), source=.FALSE.)
2281 0 : single_rank(small_list(:n_small)) = .TRUE.
2282 0 : DO i = 1, SIZE(all_index)
2283 0 : all_index(i) = i
2284 : END DO
2285 0 : row_offset(1) = 0
2286 0 : DO i = 2, SIZE(row_size)
2287 0 : row_offset(i) = row_offset(i - 1) + row_size(i - 1)
2288 : END DO
2289 0 : DO base = mypcol + 1, natom, npcol*column_env%num_pe
2290 0 : CALL timeset(routineN//'_gather_rhs', handle_gather)
2291 0 : DO slot = 0, column_env%num_pe - 1
2292 0 : atom_index = base + slot*npcol
2293 0 : IF (atom_index > natom) EXIT
2294 0 : IF (.NOT. single_rank(atom_index)) CYCLE
2295 0 : ncol = bs_env%auto_ri%sizes_opt_RI(atom_index)
2296 0 : IF (ncol == 0) CYCLE
2297 0 : ALLOCATE (buffer(SIZE(grid, 2), ncol))
2298 0 : CALL extract_atom_d_lp(mat_rhs, atom_index, all_index, row_offset, buffer)
2299 0 : CALL column_env%sum(buffer, slot)
2300 0 : IF (column_env%mepos == slot) THEN
2301 0 : CALL MOVE_ALLOC(buffer, atom_rhs)
2302 : ELSE
2303 0 : DEALLOCATE (buffer)
2304 : END IF
2305 : END DO
2306 0 : CALL timestop(handle_gather)
2307 0 : atom_index = base + column_env%mepos*npcol
2308 0 : IF (atom_index > natom) CYCLE
2309 0 : IF (.NOT. single_rank(atom_index)) CYCLE
2310 0 : ncol = bs_env%auto_ri%sizes_opt_RI(atom_index)
2311 0 : IF (ncol == 0) CYCLE
2312 : CALL build_phi_on_sphere(bs_env, qs_kinds, grid, atom_index, &
2313 : radius(atom_index), bs_env%i_ao_end_from_atom(natom), &
2314 0 : local_index, ngrid, phi, ao_map, nao)
2315 0 : ALLOCATE (local_rhs(ngrid, ncol), diagonal(ngrid))
2316 0 : DO i = 1, ngrid
2317 0 : local_rhs(i, :) = atom_rhs(local_index(i), :)
2318 : END DO
2319 0 : DEALLOCATE (atom_rhs)
2320 0 : CALL build_gram_jacobi_blas(phi, ngrid, nao, bs_env%ri_rs%tikhonov, gram, diagonal)
2321 0 : CALL scale_rows_by_diag(local_rhs, diagonal, ngrid, ncol)
2322 0 : CALL dpotrf('L', ngrid, gram, ngrid, info)
2323 0 : CPASSERT(info == 0)
2324 0 : CALL dpotrs('L', ngrid, ncol, gram, ngrid, local_rhs, ngrid, info)
2325 0 : CPASSERT(info == 0)
2326 0 : CALL scale_rows_by_diag(local_rhs, diagonal, ngrid, ncol)
2327 : CALL store_Z_lP_columns(mat_z, local_rhs, local_index, ngrid, ncol, atom_index, &
2328 0 : row_size, row_offset, 0.0_dp)
2329 0 : DEALLOCATE (local_rhs, diagonal, gram, phi, ao_map, local_index)
2330 : END DO
2331 0 : IF (n_big > 0) THEN
2332 : CALL fit_distributed_auto_ri_z_lp(qs_env, bs_env, grid, mat_rhs, mat_z, radius, &
2333 : big_list(:n_big), group_size, row_size, row_offset, &
2334 0 : all_index)
2335 : END IF
2336 0 : DEALLOCATE (radius, all_index, row_offset, single_rank, grid_size, small_list, big_list)
2337 0 : CALL column_env%free()
2338 0 : DEALLOCATE (column_env)
2339 0 : CALL dbcsr_finalize(mat_z)
2340 0 : CALL timestop(handle)
2341 0 : END SUBROUTINE fit_auto_ri_z_lp
2342 :
2343 : ! **************************************************************************************************
2344 : !> \brief Solves Σ_l' D_ll' Z_l'q = d_lp for large atomic grids with distributed Cholesky
2345 : !> factorization. Each d_lp block is gathered once and distributed over its assigned rank
2346 : !> group.
2347 : !> \param qs_env ...
2348 : !> \param bs_env ...
2349 : !> \param grid ...
2350 : !> \param mat_rhs ...
2351 : !> \param mat_z ...
2352 : !> \param radius ...
2353 : !> \param atom_list ...
2354 : !> \param group_size ...
2355 : !> \param row_size ...
2356 : !> \param row_offset ...
2357 : !> \param all_index ...
2358 : ! **************************************************************************************************
2359 0 : SUBROUTINE fit_distributed_auto_ri_z_lp(qs_env, bs_env, grid, mat_rhs, mat_z, radius, &
2360 0 : atom_list, group_size, row_size, row_offset, &
2361 0 : all_index)
2362 : TYPE(qs_environment_type), POINTER :: qs_env
2363 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2364 : REAL(KIND=dp), INTENT(IN) :: grid(:, :)
2365 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_rhs, mat_z
2366 : REAL(KIND=dp), INTENT(IN) :: radius(:)
2367 : INTEGER, INTENT(IN) :: atom_list(:), group_size, row_size(:), &
2368 : row_offset(:), all_index(:)
2369 :
2370 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_distributed_auto_ri_z_lp'
2371 :
2372 : INTEGER :: atom_index, base, handle, i, info, &
2373 : my_group, nao, ncol, ngrid, ngroups, &
2374 : root, slot
2375 0 : INTEGER, ALLOCATABLE :: ao_map(:), local_index(:)
2376 0 : REAL(KIND=dp), ALLOCATABLE :: atom_rhs(:, :), buffer(:, :), &
2377 0 : diagonal(:), local_rhs(:, :), phi(:, :)
2378 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
2379 : TYPE(cp_fm_struct_type), POINTER :: gram_struct, rhs_struct
2380 : TYPE(cp_fm_type) :: gram, rhs
2381 : TYPE(mp_para_env_type), POINTER :: group_env, para_env
2382 0 : TYPE(qs_kind_type), POINTER :: qs_kinds(:)
2383 :
2384 0 : CALL timeset(routineN, handle)
2385 0 : CALL get_qs_env(qs_env, para_env=para_env, qs_kind_set=qs_kinds)
2386 0 : ngroups = MAX(1, para_env%num_pe/group_size)
2387 0 : my_group = MIN(para_env%mepos/group_size, ngroups - 1)
2388 0 : ALLOCATE (group_env)
2389 0 : CALL group_env%from_split(para_env, my_group)
2390 0 : NULLIFY (blacs_env, gram_struct, rhs_struct)
2391 0 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=group_env)
2392 0 : DO base = 1, SIZE(atom_list), ngroups
2393 0 : DO slot = 0, ngroups - 1
2394 0 : IF (base + slot > SIZE(atom_list)) EXIT
2395 0 : atom_index = atom_list(base + slot)
2396 0 : ncol = bs_env%auto_ri%sizes_opt_RI(atom_index)
2397 0 : IF (ncol == 0) CYCLE
2398 0 : root = slot*group_size
2399 0 : ALLOCATE (buffer(SIZE(grid, 2), ncol))
2400 0 : CALL extract_atom_d_lp(mat_rhs, atom_index, all_index, row_offset, buffer)
2401 0 : CALL para_env%sum(buffer, root)
2402 0 : IF (para_env%mepos == root) THEN
2403 0 : CALL MOVE_ALLOC(buffer, atom_rhs)
2404 : ELSE
2405 0 : DEALLOCATE (buffer)
2406 : END IF
2407 : END DO
2408 0 : IF (base + my_group > SIZE(atom_list)) CYCLE
2409 0 : atom_index = atom_list(base + my_group)
2410 0 : ncol = bs_env%auto_ri%sizes_opt_RI(atom_index)
2411 0 : IF (ncol == 0) CYCLE
2412 0 : IF (group_env%mepos /= 0) ALLOCATE (atom_rhs(SIZE(grid, 2), ncol))
2413 0 : CALL group_env%bcast(atom_rhs, 0)
2414 : CALL build_phi_on_sphere(bs_env, qs_kinds, grid, atom_index, &
2415 : radius(atom_index), bs_env%i_ao_end_from_atom(bs_env%n_atom), &
2416 0 : local_index, ngrid, phi, ao_map, nao)
2417 0 : ALLOCATE (local_rhs(ngrid, ncol), diagonal(ngrid))
2418 0 : DO i = 1, ngrid
2419 0 : local_rhs(i, :) = atom_rhs(local_index(i), :)
2420 : END DO
2421 0 : DEALLOCATE (atom_rhs)
2422 0 : CALL build_jacobi_diag_from_phi(phi, ngrid, nao, diagonal)
2423 0 : CALL scale_rows_by_diag(local_rhs, diagonal, ngrid, ncol)
2424 : CALL solve_D_lp_distributed(phi, diagonal, local_rhs, ngrid, nao, ncol, &
2425 : bs_env%ri_rs%tikhonov, group_env, blacs_env, &
2426 0 : gram_struct, rhs_struct, gram, rhs, info)
2427 0 : CPASSERT(info == 0)
2428 0 : CALL scale_rows_by_diag(local_rhs, diagonal, ngrid, ncol)
2429 0 : IF (group_env%mepos == 0) THEN
2430 : CALL store_Z_lP_columns(mat_z, local_rhs, local_index, ngrid, ncol, atom_index, &
2431 0 : row_size, row_offset, 0.0_dp)
2432 : END IF
2433 0 : DEALLOCATE (local_rhs, diagonal, phi, ao_map, local_index)
2434 : END DO
2435 0 : CALL cp_blacs_env_release(blacs_env)
2436 0 : CALL group_env%free()
2437 0 : DEALLOCATE (group_env)
2438 0 : CALL timestop(handle)
2439 0 : END SUBROUTINE fit_distributed_auto_ri_z_lp
2440 :
2441 : ! **************************************************************************************************
2442 : !> \brief Splits the atoms of the Z_lP solve into a single-rank list ("small", Phase A: LAPACK
2443 : !> dpotrf/dpotrs on one rank) and a distributed list ("big", Phase B: ScaLAPACK
2444 : !> pdpotrf/pdpotrs over a rank subgroup of size G), and sizes G.
2445 : !> Each atom's single-rank peak memory is estimated as
2446 : !> peak(P) = 8*n_local_grid(P)^2 (dense matrix D'_ll', stored in D_local)
2447 : !> + 8*n_local_grid(P)*n_ao_used(P) (phi_local)
2448 : !> + 8*n_local_grid(P)*n_RI(P)*(1+n_threads) (d_lp + OMP partials)
2449 : !> and send atoms whose peak exceeds mem_safety * available-memory-per-proc to the
2450 : !> distributed path; G is auto-sized so the biggest atom's distributed D_local (/G)
2451 : !> fits alongside the replicated phi_local + d_lp. G is floored by the ScaLAPACK 32-bit
2452 : !> index limit (a local block-cyclic slice of ~n_local_grid^2/G elements must stay below
2453 : !> 2^31 or pdpotrf segfaults).
2454 : !> \param bs_env ...
2455 : !> \param ri_rs_grid_points ...
2456 : !> \param cutoff_ri_per_atom ...
2457 : !> \param ri_blk_sizes per-atom ...
2458 : !> \param n_local_grid_atom ...
2459 : !> \param small_list ...
2460 : !> \param n_small ...
2461 : !> \param big_list ...
2462 : !> \param n_big ...
2463 : !> \param G ...
2464 : ! **************************************************************************************************
2465 42 : SUBROUTINE classify_z_lp_atoms(bs_env, ri_rs_grid_points, cutoff_ri_per_atom, ri_blk_sizes, &
2466 : n_local_grid_atom, small_list, n_small, big_list, n_big, G)
2467 :
2468 : !$ USE OMP_LIB, ONLY: omp_get_max_threads
2469 :
2470 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2471 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ri_rs_grid_points
2472 :
2473 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutoff_ri_per_atom
2474 : INTEGER, DIMENSION(:), INTENT(IN) :: ri_blk_sizes
2475 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: n_local_grid_atom, small_list, big_list
2476 : INTEGER, INTENT(OUT) :: n_small, n_big, G
2477 : CHARACTER(LEN=*), PARAMETER :: routineN = 'classify_z_lp_atoms'
2478 :
2479 : INTEGER :: handle
2480 :
2481 : ! Conservative fraction of measured available memory usable per rank for the Z_lP
2482 : REAL(KIND=dp), PARAMETER :: mem_safety = 0.8_dp
2483 :
2484 : ! ScaLAPACK/BLACS index the per-rank local block-cyclic slice (~n_local_grid^2/G
2485 : ! elements) with 32-bit integers; keep it safely below 2^31 or pdpotrf segfaults.
2486 : REAL(KIND=dp), PARAMETER :: scalapack_loc_limit = 2.0E9_dp
2487 :
2488 : INTEGER :: G_atom, G_int32, l, n_ao_used_atom, &
2489 : n_grid_total, n_local_grid, natom, &
2490 : nthreads_cls, P_loop_atom
2491 : REAL(KIND=dp) :: budget_bytes, cutoff_ri, dlp_bytes, &
2492 : mem_avail_GB, ng, nri, peak_bytes, &
2493 : phi_bytes
2494 : REAL(KIND=dp), DIMENSION(3) :: pos_P
2495 : TYPE(mp_para_env_type), POINTER :: para_env
2496 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2497 :
2498 42 : CALL timeset(routineN, handle)
2499 :
2500 42 : para_env => bs_env%para_env
2501 42 : particle_set => bs_env%ri_rs%particle_set
2502 42 : natom = bs_env%n_atom
2503 42 : n_grid_total = bs_env%ri_rs%n_grid_points
2504 42 : CPASSERT(SIZE(ri_rs_grid_points, 2) == n_grid_total)
2505 :
2506 : ! Per-atom sphere size: n_local_grid(P) = number of l with |r_l - R_P| <= cutoff_ri(P).
2507 : ! It sets both the memory footprint (D_local is n_local_grid^2) and the solve
2508 : ! cost (~n_local_grid^3), so it drives classification and the LPT load balancing.
2509 210 : ALLOCATE (n_local_grid_atom(natom), small_list(natom), big_list(natom))
2510 142 : DO P_loop_atom = 1, natom
2511 400 : pos_P(:) = particle_set(P_loop_atom)%r(:)
2512 100 : cutoff_ri = cutoff_ri_per_atom(P_loop_atom)
2513 100 : n_local_grid = 0
2514 38088 : DO l = 1, n_grid_total
2515 152052 : IF (SUM((ri_rs_grid_points(1:3, l) - pos_P(1:3))**2) <= cutoff_ri**2) THEN
2516 37850 : n_local_grid = n_local_grid + 1
2517 : END IF
2518 : END DO
2519 142 : n_local_grid_atom(P_loop_atom) = n_local_grid
2520 : END DO
2521 :
2522 42 : nthreads_cls = 1
2523 42 : !$ nthreads_cls = omp_get_max_threads()
2524 42 : CALL mp_mem_avail_per_rank_GB(bs_env%para_env, mem_avail_GB) ! collective over all ranks
2525 42 : budget_bytes = mem_safety*mem_avail_GB*1.0E9_dp
2526 :
2527 42 : n_small = 0
2528 42 : n_big = 0
2529 42 : G = 1
2530 42 : G_atom = 1 ! max G a big atom needs (memory + ScaLAPACK int32 floor)
2531 42 : IF (mem_avail_GB <= 0.0_dp) THEN
2532 : ! No /proc/meminfo => cannot size by memory.
2533 0 : IF (bs_env%unit_nr > 0) THEN
2534 0 : CPWARN("RI-RS Z_lP: no meminfo; single-rank solve for all atoms")
2535 : END IF
2536 0 : DO P_loop_atom = 1, natom
2537 0 : n_small = n_small + 1
2538 0 : small_list(n_small) = P_loop_atom
2539 : END DO
2540 : ELSE
2541 : ! Classify by memory: peak (D_local + phi_local + d_lp) vs budget. Small -> BLAS,
2542 : ! big -> distributed.
2543 142 : DO P_loop_atom = 1, natom
2544 100 : ng = REAL(n_local_grid_atom(P_loop_atom), dp)
2545 100 : nri = REAL(ri_blk_sizes(P_loop_atom), dp)
2546 : CALL get_n_ao_in_sphere(bs_env, P_loop_atom, &
2547 100 : cutoff_ri_per_atom(P_loop_atom), n_ao_used_atom)
2548 100 : phi_bytes = 8.0_dp*ng*REAL(n_ao_used_atom, dp)
2549 100 : dlp_bytes = 8.0_dp*ng*nri*REAL(1 + nthreads_cls, dp)
2550 100 : peak_bytes = 8.0_dp*ng*ng + phi_bytes + dlp_bytes
2551 142 : IF (peak_bytes <= budget_bytes) THEN
2552 100 : n_small = n_small + 1
2553 100 : small_list(n_small) = P_loop_atom
2554 : ELSE
2555 0 : n_big = n_big + 1
2556 0 : big_list(n_big) = P_loop_atom
2557 : ! G must satisfy BOTH: (a) memory — distributed D_local (/G) fits next to the
2558 : ! replicated phi_local + d_lp; (b) ScaLAPACK — local ~ng^2/G below the int32 limit.
2559 0 : G_int32 = CEILING(ng*ng/scalapack_loc_limit)
2560 : G_atom = MAX(G_atom, G_int32, &
2561 0 : CEILING(8.0_dp*ng*ng/MAX(budget_bytes - phi_bytes - dlp_bytes, 1.0_dp)))
2562 : END IF
2563 : END DO
2564 42 : IF (n_big > 0) THEN
2565 : ! Auto-size G from the most demanding big atom.
2566 0 : IF (G_atom > para_env%num_pe) THEN
2567 : CALL cp_abort(__LOCATION__, &
2568 : "RI-RS Z_lP: an atom is too large to fit even when "// &
2569 : "distributed over all ranks. Add nodes, use fewer MPI ranks "// &
2570 : "per node, lower CUTOFF_RADIUS_RI, or raise EPS_FILTER "// &
2571 0 : "for more grid screening.")
2572 : END IF
2573 0 : G = MIN(MAX(G_atom, 2), para_env%num_pe)
2574 : END IF
2575 : END IF
2576 :
2577 42 : CALL timestop(handle)
2578 :
2579 42 : END SUBROUTINE classify_z_lp_atoms
2580 :
2581 : ! **************************************************************************************************
2582 : !> \brief Number of AO basis functions that can be non-zero inside the RI-RS integration sphere
2583 : !> \param bs_env ...
2584 : !> \param atom_P ...
2585 : !> \param cutoff_ri ...
2586 : !> \param n_ao_used ...
2587 : ! **************************************************************************************************
2588 100 : SUBROUTINE get_n_ao_in_sphere(bs_env, atom_P, cutoff_ri, n_ao_used)
2589 :
2590 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2591 : INTEGER, INTENT(IN) :: atom_P
2592 : REAL(KIND=dp), INTENT(IN) :: cutoff_ri
2593 : INTEGER, INTENT(OUT) :: n_ao_used
2594 :
2595 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_n_ao_in_sphere'
2596 :
2597 : INTEGER :: handle, ri_atom
2598 100 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2599 :
2600 100 : CALL timeset(routineN, handle)
2601 :
2602 100 : particle_set => bs_env%ri_rs%particle_set
2603 100 : n_ao_used = 0
2604 348 : DO ri_atom = 1, bs_env%n_atom
2605 992 : IF (NORM2(particle_set(ri_atom)%r(:) - particle_set(atom_P)%r(:)) > &
2606 : bs_env%ri_rs%radius_ao_per_atom(ri_atom) + cutoff_ri) CYCLE
2607 : n_ao_used = n_ao_used + bs_env%i_ao_end_from_atom(ri_atom) - &
2608 348 : bs_env%i_ao_start_from_atom(ri_atom) + 1
2609 : END DO
2610 :
2611 100 : CALL timestop(handle)
2612 :
2613 100 : END SUBROUTINE get_n_ao_in_sphere
2614 :
2615 : ! **************************************************************************************************
2616 : !> \brief Builds the sphere-local AO matrix phi_local(l, μ) = ϕ_μ(r_l) for one RI atom P
2617 : !> \param bs_env ...
2618 : !> \param qs_kind_set ...
2619 : !> \param ri_rs_grid_points ...
2620 : !> \param atom_P ...
2621 : !> \param cutoff_ri ...
2622 : !> \param n_ao_total ...
2623 : !> \param local_grid_idx ...
2624 : !> \param n_local_grid ...
2625 : !> \param phi_local ...
2626 : !> \param ao_col_map ...
2627 : !> \param n_ao_used ...
2628 : !> \param center ...
2629 : ! **************************************************************************************************
2630 54 : SUBROUTINE build_phi_on_sphere(bs_env, qs_kind_set, ri_rs_grid_points, &
2631 : atom_P, cutoff_ri, n_ao_total, local_grid_idx, n_local_grid, &
2632 : phi_local, ao_col_map, n_ao_used, center)
2633 :
2634 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2635 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2636 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ri_rs_grid_points
2637 : INTEGER, INTENT(IN) :: atom_P
2638 : REAL(KIND=dp), INTENT(IN) :: cutoff_ri
2639 : INTEGER, INTENT(IN) :: n_ao_total
2640 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: local_grid_idx
2641 : INTEGER, INTENT(OUT) :: n_local_grid
2642 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
2643 : INTENT(OUT) :: phi_local
2644 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_col_map
2645 : INTEGER, INTENT(OUT) :: n_ao_used
2646 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: center
2647 :
2648 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_phi_on_sphere'
2649 :
2650 : INTEGER :: col_end, col_start, handle, j, k, l, &
2651 : loc_idx, n_grid_total, n_keep, ri_atom
2652 : REAL(KIND=dp) :: d_sP, dist, r2_threshold
2653 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w_pt
2654 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: phi_keep, sphere_grid
2655 : REAL(KIND=dp), DIMENSION(3) :: pos_P
2656 : TYPE(cell_type), POINTER :: cell
2657 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2658 :
2659 54 : CALL timeset(routineN, handle)
2660 :
2661 54 : cell => bs_env%ri_rs%cell
2662 54 : particle_set => bs_env%ri_rs%particle_set
2663 : ! AUTO_RI may pass only the global grid points needed for the current atom block.
2664 54 : n_grid_total = SIZE(ri_rs_grid_points, 2)
2665 54 : IF (PRESENT(center)) THEN
2666 4 : pos_P(:) = center(:)
2667 : ELSE
2668 200 : pos_P(:) = particle_set(atom_P)%r(:)
2669 : END IF
2670 :
2671 54 : n_local_grid = 0
2672 19174 : DO l = 1, n_grid_total
2673 76480 : dist = NORM2(ri_rs_grid_points(1:3, l) - pos_P(1:3))
2674 19174 : IF (dist <= cutoff_ri) n_local_grid = n_local_grid + 1
2675 : END DO
2676 :
2677 162 : ALLOCATE (local_grid_idx(n_local_grid))
2678 :
2679 54 : n_local_grid = 0
2680 19174 : DO l = 1, n_grid_total
2681 76480 : dist = NORM2(ri_rs_grid_points(1:3, l) - pos_P(1:3))
2682 19174 : IF (dist <= cutoff_ri) THEN
2683 19051 : n_local_grid = n_local_grid + 1
2684 19051 : local_grid_idx(n_local_grid) = l
2685 : END IF
2686 : END DO
2687 :
2688 162 : ALLOCATE (sphere_grid(3, n_local_grid))
2689 19105 : DO loc_idx = 1, n_local_grid
2690 76258 : sphere_grid(:, loc_idx) = ri_rs_grid_points(:, local_grid_idx(loc_idx))
2691 : END DO
2692 :
2693 : ! Only AOs on atoms that reach into the sphere can be non-zero here
2694 162 : ALLOCATE (ao_col_map(n_ao_total))
2695 54 : ao_col_map(:) = 0
2696 54 : n_ao_used = 0
2697 190 : DO ri_atom = 1, bs_env%n_atom
2698 544 : d_sP = NORM2(particle_set(ri_atom)%r(:) - pos_P(:))
2699 136 : IF (d_sP > bs_env%ri_rs%radius_ao_per_atom(ri_atom) + cutoff_ri) CYCLE
2700 :
2701 690 : DO j = bs_env%i_ao_start_from_atom(ri_atom), bs_env%i_ao_end_from_atom(ri_atom)
2702 500 : n_ao_used = n_ao_used + 1
2703 636 : ao_col_map(j) = n_ao_used
2704 : END DO
2705 : END DO
2706 :
2707 216 : ALLOCATE (phi_local(n_local_grid, n_ao_used))
2708 54 : phi_local = 0.0_dp
2709 :
2710 190 : DO ri_atom = 1, bs_env%n_atom
2711 544 : d_sP = NORM2(particle_set(ri_atom)%r(:) - pos_P(:))
2712 136 : IF (d_sP > bs_env%ri_rs%radius_ao_per_atom(ri_atom) + cutoff_ri) CYCLE
2713 :
2714 136 : col_start = ao_col_map(bs_env%i_ao_start_from_atom(ri_atom))
2715 136 : col_end = ao_col_map(bs_env%i_ao_end_from_atom(ri_atom))
2716 : ! A positive CUTOFF_RADIUS_RI_AO overrides the per-atom Gaussian radius
2717 : ! with a user-defined hard cutoff.
2718 136 : IF (bs_env%ri_rs%cutoff_radius_ri_ao > 0.0_dp) THEN
2719 9 : r2_threshold = bs_env%ri_rs%cutoff_radius_ri_ao**2
2720 : ELSE
2721 127 : r2_threshold = bs_env%ri_rs%radius_ao_per_atom(ri_atom)**2
2722 : END IF
2723 :
2724 : CALL evaluate_ao_on_points(phi_local(:, col_start:col_end), sphere_grid, &
2725 : ri_atom, particle_set, qs_kind_set, cell, &
2726 190 : cutoff_squared=r2_threshold)
2727 : END DO
2728 :
2729 54 : DEALLOCATE (sphere_grid)
2730 :
2731 54 : IF (n_local_grid > 0) THEN
2732 162 : ALLOCATE (w_pt(n_local_grid))
2733 : !$OMP PARALLEL DO DEFAULT(NONE) &
2734 : !$OMP SHARED(n_local_grid, n_ao_used, phi_local, w_pt) &
2735 54 : !$OMP PRIVATE(l, j) SCHEDULE(STATIC)
2736 : DO l = 1, n_local_grid
2737 : w_pt(l) = 0.0_dp
2738 : DO j = 1, n_ao_used
2739 : w_pt(l) = MAX(w_pt(l), ABS(phi_local(l, j)))
2740 : END DO
2741 : END DO
2742 : !$OMP END PARALLEL DO
2743 19105 : n_keep = COUNT(w_pt > bs_env%eps_filter)
2744 54 : IF (n_keep < n_local_grid) THEN
2745 24 : ALLOCATE (phi_keep(n_keep, n_ao_used))
2746 6 : k = 0
2747 2601 : DO l = 1, n_local_grid
2748 2601 : IF (w_pt(l) > bs_env%eps_filter) THEN
2749 2151 : k = k + 1
2750 17208 : phi_keep(k, :) = phi_local(l, :)
2751 2151 : local_grid_idx(k) = local_grid_idx(l)
2752 : END IF
2753 : END DO
2754 6 : CALL MOVE_ALLOC(phi_keep, phi_local)
2755 6 : n_local_grid = n_keep
2756 : END IF
2757 54 : DEALLOCATE (w_pt)
2758 : END IF
2759 :
2760 54 : CALL timestop(handle)
2761 :
2762 108 : END SUBROUTINE build_phi_on_sphere
2763 :
2764 : ! **************************************************************************************************
2765 : !> \brief Computes a three-center integral block directly in AUTO_RI contraction columns.
2766 : !>
2767 : !> For each angular momentum l, the primitive RI Gaussians on the atoms contributing
2768 : !> to one optimized function are collected before the contraction
2769 : !>
2770 : !> (μν|p) = Σ_P (μν|P) U_Pp .
2771 : !>
2772 : !> This avoids constructing and retaining the complete reference-basis block (μν|P).
2773 : !> \param int_3c ...
2774 : !> \param ctx ...
2775 : !> \param ws ...
2776 : !> \param atom_j ...
2777 : !> \param atom_k ...
2778 : !> \param atom_i ...
2779 : !> \param transform ...
2780 : !> \param transform_row ...
2781 : !> \param screened ...
2782 : ! **************************************************************************************************
2783 108 : SUBROUTINE build_3c_integral_block_auto_ri_ctx(int_3c, ctx, ws, atom_j, atom_k, atom_i, &
2784 36 : transform, transform_row, screened)
2785 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: int_3c
2786 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
2787 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
2788 : INTEGER, INTENT(IN) :: atom_j, atom_k, atom_i
2789 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: transform
2790 : INTEGER, INTENT(IN) :: transform_row
2791 : LOGICAL, INTENT(OUT), OPTIONAL :: screened
2792 :
2793 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integral_block_auto_ri_ctx'
2794 :
2795 : INTEGER :: handle_contract, handle_eri, ikind, iset, jkind, jset, kkind, kset, l, ncoi, &
2796 : ncoj, ncok, ncol, npgf_group, nseti, nsetj, nsetk, primitive_first, sgfi, sgfj, sgfk
2797 36 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
2798 36 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
2799 36 : nsgfj, nsgfk
2800 36 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
2801 : REAL(KIND=dp) :: dij, dik, djk, group_radius, &
2802 : kind_radius_i, kind_radius_j, &
2803 : kind_radius_k, sijk_ext
2804 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rpgf_group, zet_group
2805 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: spi_group
2806 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sijk, sijk_contr
2807 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
2808 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
2809 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, zeti, zetj, zetk
2810 :
2811 36 : IF (PRESENT(screened)) screened = .FALSE.
2812 36 : ri = pbc(ctx%particle_set(atom_i)%r(1:3), ctx%cell)
2813 36 : rj = pbc(ctx%particle_set(atom_j)%r(1:3), ctx%cell)
2814 36 : rk = pbc(ctx%particle_set(atom_k)%r(1:3), ctx%cell)
2815 144 : rjk = rk - rj
2816 144 : rij = rj - ri
2817 144 : rik = rk - ri
2818 144 : djk = NORM2(rjk)
2819 144 : dij = NORM2(rij)
2820 144 : dik = NORM2(rik)
2821 :
2822 36 : ikind = ctx%kind_of(atom_i)
2823 36 : jkind = ctx%kind_of(atom_j)
2824 36 : kkind = ctx%kind_of(atom_k)
2825 : CALL get_gto_basis_set(ctx%basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, &
2826 : lmax=lmax_i, lmin=lmin_i, npgf=npgfi, nset=nseti, &
2827 : nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
2828 36 : zet=zeti, kind_radius=kind_radius_i)
2829 : CALL get_gto_basis_set(ctx%basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, &
2830 : lmax=lmax_j, lmin=lmin_j, npgf=npgfj, nset=nsetj, &
2831 : nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
2832 36 : zet=zetj, kind_radius=kind_radius_j)
2833 : CALL get_gto_basis_set(ctx%basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, &
2834 : lmax=lmax_k, lmin=lmin_k, npgf=npgfk, nset=nsetk, &
2835 : nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
2836 36 : zet=zetk, kind_radius=kind_radius_k)
2837 :
2838 : IF (kind_radius_j + kind_radius_i + ctx%dr_ij < dij .OR. &
2839 36 : kind_radius_j + kind_radius_k + ctx%dr_jk < djk .OR. &
2840 : kind_radius_k + kind_radius_i + ctx%dr_ik < dik) THEN
2841 0 : IF (PRESENT(screened)) screened = .TRUE.
2842 0 : RETURN
2843 : END IF
2844 :
2845 36 : ncol = SIZE(transform, 2)
2846 36 : CPASSERT(SIZE(int_3c, 3) == ncol)
2847 252 : DO l = 0, ctx%maxli
2848 216 : npgf_group = 0
2849 216 : group_radius = 0.0_dp
2850 4392 : DO iset = 1, nseti
2851 4176 : IF (lmin_i(iset) /= l .OR. lmax_i(iset) /= l) CYCLE
2852 696 : npgf_group = npgf_group + npgfi(iset)
2853 4392 : group_radius = MAX(group_radius, set_radius_i(iset))
2854 : END DO
2855 216 : IF (npgf_group == 0) CYCLE
2856 192 : ncoi = npgf_group*ncoset(l)
2857 768 : ALLOCATE (zet_group(npgf_group), rpgf_group(npgf_group))
2858 768 : ALLOCATE (spi_group(ncoi, ncol), source=0.0_dp)
2859 192 : primitive_first = 1
2860 4008 : DO iset = 1, nseti
2861 3816 : IF (lmin_i(iset) /= l .OR. lmax_i(iset) /= l) CYCLE
2862 : zet_group(primitive_first:primitive_first + npgfi(iset) - 1) = &
2863 1392 : zeti(1:npgfi(iset), iset)
2864 : rpgf_group(primitive_first:primitive_first + npgfi(iset) - 1) = &
2865 1392 : rpgf_i(1:npgfi(iset), iset)
2866 696 : sgfi = first_sgf_i(1, iset)
2867 : spi_group((primitive_first - 1)*ncoset(l) + 1: &
2868 : (primitive_first + npgfi(iset) - 1)*ncoset(l), :) = &
2869 696 : MATMUL(ctx%spi(iset, ikind)%array, &
2870 : transform(transform_row + sgfi - 1: &
2871 1437180 : transform_row + sgfi + nsgfi(iset) - 2, :))
2872 4008 : primitive_first = primitive_first + npgfi(iset)
2873 : END DO
2874 :
2875 480 : DO jset = 1, nsetj
2876 288 : IF (set_radius_j(jset) + group_radius + ctx%dr_ij < dij) CYCLE
2877 832 : DO kset = 1, nsetk
2878 352 : IF (set_radius_j(jset) + set_radius_k(kset) + ctx%dr_jk < djk) CYCLE
2879 352 : IF (set_radius_k(kset) + group_radius + ctx%dr_ik < dik) CYCLE
2880 352 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
2881 352 : ncok = npgfk(kset)*ncoset(lmax_k(kset))
2882 352 : sgfj = first_sgf_j(1, jset)
2883 352 : sgfk = first_sgf_k(1, kset)
2884 352 : IF (ncoj*ncok*ncoi <= 0) CYCLE
2885 1760 : ALLOCATE (sijk(ncoj, ncok, ncoi), source=0.0_dp)
2886 352 : CALL timeset(routineN//'_eri', handle_eri)
2887 : CALL eri_3center(sijk, &
2888 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
2889 : rpgf_j(:, jset), rj, &
2890 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), &
2891 : rpgf_k(:, kset), rk, l, l, npgf_group, zet_group, &
2892 : rpgf_group, ri, djk, dij, dik, ws%lib, ctx%potential_parameter, &
2893 352 : int_abc_ext=sijk_ext)
2894 352 : CALL timestop(handle_eri)
2895 1760 : ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), ncol))
2896 352 : CALL timeset(routineN//'_contract', handle_contract)
2897 : CALL abc_contract_xsmm(sijk_contr, sijk, ctx%tspj(jset, jkind)%array, &
2898 : ctx%spk(kset, kkind)%array, spi_group, ncoj, ncok, ncoi, &
2899 352 : nsgfj(jset), nsgfk(kset), ncol, ws%cpp_buffer, ws%ccp_buffer)
2900 352 : CALL timestop(handle_contract)
2901 352 : DEALLOCATE (sijk)
2902 : int_3c(sgfj:sgfj + nsgfj(jset) - 1, sgfk:sgfk + nsgfk(kset) - 1, :) = &
2903 : int_3c(sgfj:sgfj + nsgfj(jset) - 1, sgfk:sgfk + nsgfk(kset) - 1, :) + &
2904 14476 : sijk_contr
2905 1344 : DEALLOCATE (sijk_contr)
2906 : END DO
2907 : END DO
2908 252 : DEALLOCATE (zet_group, rpgf_group, spi_group)
2909 : END DO
2910 :
2911 72 : END SUBROUTINE build_3c_integral_block_auto_ri_ctx
2912 :
2913 : ! **************************************************************************************************
2914 : !> \brief Returns the reference RI basis size of one atom.
2915 : !> \param bs_env ...
2916 : !> \param iatom ...
2917 : !> \return ...
2918 : ! **************************************************************************************************
2919 15 : INTEGER FUNCTION get_ref_RI_size(bs_env, iatom) RESULT(n)
2920 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2921 : INTEGER, INTENT(IN) :: iatom
2922 :
2923 : INTEGER :: ikind
2924 :
2925 15 : ikind = bs_env%ri_rs%particle_set(iatom)%atomic_kind%kind_number
2926 15 : n = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
2927 :
2928 15 : END FUNCTION get_ref_RI_size
2929 :
2930 : END MODULE gw_ri_rs_compute_Z_lP
|
