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