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 Automatic RI basis set optimization for molecular GW.
10 : !> \par History
11 : !> 09.2026 created [Jan Wilhelm]
12 : ! **************************************************************************************************
13 : MODULE gw_auto_ri_generate_RI_basis
14 : USE ai_contraction_sphi, ONLY: ab_contract,&
15 : abc_contract
16 : USE atomic_kind_types, ONLY: get_atomic_kind_set
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE constants_operator, ONLY: operator_coulomb,&
20 : operator_truncated
21 : USE coulomb_integral_interface, ONLY: compute_coulomb_2c,&
22 : compute_coulomb_3c,&
23 : coulomb_integral_cleanup,&
24 : coulomb_integral_context_type,&
25 : coulomb_integral_init,&
26 : coulomb_operator_type
27 : USE cp_blacs_env, ONLY: cp_blacs_env_type
28 : USE cp_dbcsr_api, ONLY: &
29 : dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, dbcsr_distribution_release, &
30 : dbcsr_distribution_type, dbcsr_finalize, dbcsr_put_block, dbcsr_release, dbcsr_type, &
31 : dbcsr_type_symmetric
32 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
33 : USE cp_files, ONLY: close_file,&
34 : discover_file,&
35 : open_file
36 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
37 : cp_fm_syrk,&
38 : cp_fm_uplo_to_full
39 : USE cp_fm_diag, ONLY: cp_fm_geeig
40 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
41 : cp_fm_struct_release,&
42 : cp_fm_struct_type
43 : USE cp_fm_types, ONLY: cp_fm_create,&
44 : cp_fm_get_info,&
45 : cp_fm_get_submatrix,&
46 : cp_fm_release,&
47 : cp_fm_type
48 : USE generic_os_integrals, ONLY: int_operators_r12_ab_os
49 : USE gw_auto_ri_types, ONLY: auto_ri_type
50 : USE gw_utils_compute_integrals, ONLY: build_3c_integral_block_ctx,&
51 : gw_3c_ctx_create,&
52 : gw_3c_ctx_release,&
53 : gw_3c_ctx_type,&
54 : gw_3c_ws_create,&
55 : gw_3c_ws_release,&
56 : gw_3c_ws_type
57 : USE ieee_arithmetic, ONLY: ieee_is_finite
58 : USE input_constants, ONLY: do_potential_coulomb,&
59 : do_potential_truncated,&
60 : tensor_small_cell_full_kp
61 : USE integral_library_types, ONLY: active_integral_library,&
62 : library_libint
63 : USE kinds, ONLY: default_path_length,&
64 : dp
65 : USE libint_2c_3c, ONLY: libint_potential_type
66 : USE local_gemm_api, ONLY: LOCAL_GEMM_PU_HOST,&
67 : local_gemm_ctxt_type
68 : USE machine, ONLY: m_walltime
69 : USE mathlib, ONLY: diamat_all
70 : USE message_passing, ONLY: mp_para_env_type
71 : USE orbital_pointers, ONLY: init_orbital_pointers,&
72 : ncoset
73 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
74 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
75 : USE qs_environment_types, ONLY: get_qs_env,&
76 : qs_environment_type
77 : USE t_c_g0, ONLY: get_lmax_init,&
78 : init_t_c_g0 => init
79 : USE util, ONLY: sort
80 : #include "base/base_uses.f90"
81 :
82 : IMPLICIT NONE
83 : PRIVATE
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_auto_ri_generate_RI_basis'
86 :
87 : ! Numerical cutoffs for the reference M_PQ eigenspace and final M_pq pseudoinverse.
88 : REAL(KIND=dp), PARAMETER :: reference_M_eigenval_threshold = 1.0E-10_dp
89 : REAL(KIND=dp), PARAMETER :: optimized_M_eigenval_threshold = 1.0E-6_dp
90 : ! Discard null target directions and reject non-positive-semidefinite M_PQ or B_PQ spectra.
91 : REAL(KIND=dp), PARAMETER :: relative_target_eigenval_threshold = 1.0E-12_dp
92 : REAL(KIND=dp), PARAMETER :: AB_absolute_target_eigenval_threshold = 1.0E-12_dp
93 : REAL(KIND=dp), PARAMETER :: negative_eigenval_tolerance = 1.0E-10_dp
94 : REAL(KIND=dp), PARAMETER :: conditioned_negative_eigenval_tolerance = SQRT(EPSILON(1.0_dp))
95 :
96 : TYPE AB_block_list_type
97 : INTEGER, ALLOCATABLE :: block(:)
98 : INTEGER, ALLOCATABLE :: first_p_AB(:)
99 : INTEGER, ALLOCATABLE :: first_p_atom(:)
100 : INTEGER, ALLOCATABLE :: size_opt_RI(:)
101 : END TYPE AB_block_list_type
102 :
103 : TYPE optimized_atom_U_type
104 : INTEGER, ALLOCATABLE :: ref_atom(:)
105 : INTEGER, ALLOCATABLE :: first_P(:)
106 : REAL(KIND=dp), ALLOCATABLE :: U_Pp(:, :)
107 : END TYPE optimized_atom_U_type
108 :
109 : PUBLIC :: generate_auto_ri_basis
110 :
111 : CONTAINS
112 :
113 : !**************************************************************************************************!
114 : ! The full derivation of the AUTO_RI algorithm will be published separately. !
115 : ! !
116 : ! The equations below define the notation used throughout this automatic RI basis optimizer. !
117 : ! !
118 : ! The RI scalar product uses the truncated Coulomb operator: !
119 : ! !
120 : ! (f|g) = ∫dr∫dr' f(r) V_rc(r-r') g(r'). (1) !
121 : ! !
122 : ! For the RI basis-set optimization, we consider the orbital-product density formed by !
123 : ! AOs on atoms A and B, !
124 : ! !
125 : ! ρ_mn^AB(r) = Σ_{μ∈A} Σ_{ν∈B} C_μm C_νn ϕ_μ(r) ϕ_ν(r). (2) !
126 : ! !
127 : ! The atom-local correlation target B_PQ^{A,C} mimics contributions to the density response χ, !
128 : ! and the exchange target B_PQ^{A,X} mimics contributions to the exchange self-energy: !
129 : ! !
130 : ! B_PQ^{A,C} = (1/N_C^A) Σ_ia^oe w_ia (φ_P^A|ρ_ia^AA)(ρ_ia^AA|φ_Q^A). (3) !
131 : ! !
132 : ! B_PQ^{A,X} = (1/N_X^A) Σ_{n∈W_occ} Σ_i^occ (φ_P^A|ρ_ni^AA)(ρ_ni^AA|φ_Q^A). (4) !
133 : ! !
134 : ! W_occ = [ε_HOMO-E_win, ε_LUMO+E_win], where E_win is OCC_EMPTY_FRONTIER_ORBITAL_WINDOW. !
135 : ! !
136 : ! Here w_ia gives larger weight to frontier-orbital products close to the HOMO and LUMO: !
137 : ! !
138 : ! w_ia = 1/(|ε_i-E_F| |ε_a-E_F|), E_F = (ε_HOMO+ε_LUMO)/2. (5) !
139 : ! !
140 : ! Diagonalizing B^A = B^{A,C}+B^{A,X} identifies the most important atom-local RI functions !
141 : ! that we use in the optimized RI basis set. The contraction coefficients U_Pp^A [see Eq. (7)] !
142 : ! used for optimized RI basis functions are obtained from the eigenvalue problem !
143 : ! !
144 : ! B^A U^A = M^A U^A λ^A, Σ_Q B_PQ^A U_Qp^A = λ_p^A Σ_Q M_PQ^A U_Qp^A. (6) !
145 : ! !
146 : ! The eigenvectors with the largest eigenvalues are the contracted RI functions φ_p^A(r) !
147 : ! used for the optimized RI basis set: !
148 : ! !
149 : ! φ_p^A(r) = Σ_P U_Pp^A φ_P^A(r). (7) !
150 : ! !
151 : ! We also use superpositions of Gaussians on neighboring atoms in the RI basis. Before optimizing !
152 : ! the expansion coefficients, the AA from Eq. (7) used in AB contractions are projected out: !
153 : ! !
154 : ! span{φ_P^{AB,⊥}} = span{φ_P^A,φ_Q^B} \ span{φ_p^A,φ_q^B}. (8) !
155 : ! !
156 : ! For obtaining the AB contraction coefficients, we use the symmetric off-diagonal density: !
157 : ! !
158 : ! ρ̅_mn^AB(r) = ρ_mn^AB(r)+ρ_mn^BA(r). (9) !
159 : ! !
160 : ! To identify the most important functions for a neighboring pair A-B, define !
161 : ! correlation and exchange targets in this projected basis: !
162 : ! !
163 : ! B_PQ^{AB,C} = (1/N_C^{AB}) Σ_ia^oe w_ia (φ_P^{AB,⊥}|ρ̅_ia^AB)(ρ̅_ia^AB|φ_Q^{AB,⊥}). (10) !
164 : ! !
165 : ! B_PQ^{AB,X} = (1/N_X^{AB}) Σ_{n∈W_occ} Σ_i^occ (φ_P^{AB,⊥}|ρ̅_ni^AB)(ρ̅_ni^AB|φ_Q^{AB,⊥}). (11) !
166 : ! !
167 : ! Diagonalizing their sum gives the AB-neighbor-pair contraction coefficients: !
168 : ! !
169 : ! B^AB U^AB = U^AB λ^AB, B^AB = B^{AB,C}+B^{AB,X}. (12) !
170 : ! !
171 : ! The eigenvectors with the largest eigenvalues are used as AB contractions in the optimized !
172 : ! RI basis set: !
173 : ! !
174 : ! φ_p^AB(r) = Σ_P U_Pp^AB φ_P^{AB,⊥}(r). (13) !
175 : ! !
176 : ! Eq. (13) means that the RI basis is not an atom-centered basis any more, but can be delocalized !
177 : ! over two neighboring atoms AB. !
178 : ! !
179 : ! AUTO_RI uses the same number of AA and AB contractions in the final optimized RI basis: !
180 : ! !
181 : ! Σ_A N_A = Σ_{A>B} N_AB = N_RI/2. (14) !
182 : ! !
183 : ! If no atom pairs are available, as for an isolated atom or when NEIGHBOR_RADIUS is below the !
184 : ! minimum bond length, there is no AB space, so the full requested basis is assigned to AA: !
185 : ! !
186 : ! Σ_A N_A^target = N_RI^target. (15) !
187 : ! !
188 : ! The retained size can be smaller only when B^A or M^A has fewer numerically non-null modes. !
189 : ! !
190 : !**************************************************************************************************!
191 :
192 : ! **************************************************************************************************
193 : !> \brief Executes the AUTO_RI algorithm defined by Eqs. (1)-(15).
194 : !> \param qs_env ...
195 : !> \param bs_env Band-structure environment containing GW parameters.
196 : ! **************************************************************************************************
197 4 : SUBROUTINE generate_auto_ri_basis(qs_env, bs_env)
198 : TYPE(qs_environment_type), POINTER :: qs_env
199 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
200 :
201 : CHARACTER(LEN=*), PARAMETER :: routineN = 'generate_auto_ri_basis'
202 :
203 : INTEGER :: handle
204 : REAL(KIND=dp) :: time_start
205 4 : REAL(KIND=dp), ALLOCATABLE :: U_Pp_A(:, :, :)
206 :
207 4 : CALL timeset(routineN, handle)
208 :
209 4 : time_start = m_walltime()
210 :
211 : ! Compute AA target matrices B^{AA,C}, B^{AA,X} and contraction U_Pp^A in Eqs. (3)-(7).
212 4 : CALL compute_AA_optimized_RI_basis(bs_env, U_Pp_A)
213 :
214 : ! Compute AB target matrices B^{AB,C}, B^{AB,X} and contraction U_Pp^AB in Eqs. (8)-(13).
215 4 : CALL compute_AB_optimized_RI_basis(bs_env, U_Pp_A)
216 4 : DEALLOCATE (U_Pp_A)
217 :
218 : ! Compute and invert M_pq = (φ_p|φ_q) of the optimized RI basis in Eqs. (7) and (13).
219 4 : CALL compute_M_pq_inv_and_V_pq(qs_env, bs_env)
220 :
221 : ! Print the optimized RI basis dimension and total optimization time.
222 4 : CALL print_auto_ri_summary(bs_env, time_start)
223 :
224 4 : CALL timestop(handle)
225 :
226 4 : END SUBROUTINE generate_auto_ri_basis
227 :
228 : ! **************************************************************************************************
229 : !> \brief Computes the atom-local contraction coefficients U^A from Eqs. (3)-(7):
230 : !>
231 : !> B_PQ^{A,C} = (1/N_C^A) Σ_ia^oe w_ia (φ_P^A|ρ_ia^AA)(ρ_ia^AA|φ_Q^A), (3)
232 : !>
233 : !> B_PQ^{A,X} = (1/N_X^A) Σ_{n∈W_occ} Σ_i^occ (φ_P^A|ρ_ni^AA)(ρ_ni^AA|φ_Q^A), (4)
234 : !>
235 : !> w_ia = 1/(|ε_i-E_F| |ε_a-E_F|), E_F = (ε_HOMO+ε_LUMO)/2, (5)
236 : !>
237 : !> B^A U^A = M^A U^A λ^A, (6)
238 : !>
239 : !> Contracted RI basis functions for optimized RI basis:
240 : !>
241 : !> φ_p^A(r) = Σ_P U_Pp^A φ_P^A(r). (7)
242 : !>
243 : !> \param bs_env ...
244 : !> \param U_Pp_A temporary atom-local contraction coefficients used to construct the AB space
245 : ! **************************************************************************************************
246 4 : SUBROUTINE compute_AA_optimized_RI_basis(bs_env, U_Pp_A)
247 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
248 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: U_Pp_A(:, :, :)
249 :
250 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_AA_optimized_RI_basis'
251 :
252 : INTEGER :: handle, iatom, ikind, max_nAO_kind, max_nRI_kind, nAO_kind, natom, nkind, &
253 : nRI_kind, RI_size_opt_AA_requested, root_rank
254 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of_atom, sizes_opt_RI_AA
255 : REAL(KIND=dp) :: target_norm
256 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
257 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: B_PQ, eigenval_AA, M_PQ, U_Pp_A_atom
258 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: M_PQ_by_kind
259 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: Int_3c
260 :
261 4 : CALL timeset(routineN, handle)
262 :
263 : ! Determine reference dimensions and the number of AA RI functions (i.e. number of
264 : ! eigenvectors U_Pp^A used for RI-basis contractions).
265 : CALL initialize_AA_optimization(bs_env, kind_of_atom, natom, nkind, &
266 4 : max_nRI_kind, RI_size_opt_AA_requested)
267 16 : ALLOCATE (eigenval_AA(max_nRI_kind, natom), source=0.0_dp)
268 12 : ALLOCATE (sizes_opt_RI_AA(natom))
269 :
270 4 : CALL initialize_auto_ri_integrals(bs_env)
271 16 : max_nAO_kind = MAXVAL(bs_env%sizes_AO)
272 0 : ALLOCATE (U_Pp_A(max_nRI_kind, max_nRI_kind, natom), &
273 36 : M_PQ_by_kind(max_nRI_kind, max_nRI_kind, nkind), source=0.0_dp)
274 24 : ALLOCATE (Int_3c(max_nAO_kind, max_nAO_kind, max_nRI_kind, nkind), source=0.0_dp)
275 :
276 : ! The few atom-kind integral blocks are inexpensive; compute them on every rank.
277 12 : DO ikind = 1, nkind
278 8 : nAO_kind = bs_env%basis_set_AO(ikind)%gto_basis_set%nsgf
279 8 : nRI_kind = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
280 : CALL compute_int_3c( &
281 : ao_basis=bs_env%basis_set_AO(ikind)%gto_basis_set, &
282 : ri_basis=bs_env%basis_set_RI(ikind)%gto_basis_set, &
283 : potential=bs_env%ri_metric, M_PQ=M_PQ, &
284 8 : Int_3c=Int_3c(1:nAO_kind, 1:nAO_kind, 1:nRI_kind, ikind))
285 : CALL add_RI_metric_regularization(bs_env%basis_set_RI(ikind)%gto_basis_set, &
286 8 : bs_env%regularization_RI, M_PQ)
287 78376 : M_PQ_by_kind(1:nRI_kind, 1:nRI_kind, ikind) = M_PQ
288 12 : DEALLOCATE (M_PQ)
289 : END DO
290 :
291 16 : bs_env%auto_ri%sizes_opt_RI = 0
292 16 : DO iatom = 1, natom
293 12 : root_rank = MODULO(iatom - 1, bs_env%para_env%num_pe)
294 12 : ikind = kind_of_atom(iatom)
295 12 : nRI_kind = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
296 12 : nAO_kind = bs_env%sizes_AO(iatom)
297 48 : ALLOCATE (M_PQ(nRI_kind, nRI_kind))
298 90700 : M_PQ(:, :) = M_PQ_by_kind(1:nRI_kind, 1:nRI_kind, ikind)
299 : ! Eqs. (3)-(5): form B^{A,C}, B^{A,X}, and their weighted sum B^A.
300 : CALL compute_B_PQ_AA(bs_env, iatom, root_rank, &
301 : Int_3c(1:nAO_kind, 1:nAO_kind, 1:nRI_kind, ikind), &
302 17204 : M_PQ, B_PQ, target_norm)
303 12 : DEALLOCATE (M_PQ)
304 12 : IF (root_rank /= bs_env%para_env%mepos) THEN
305 6 : DEALLOCATE (B_PQ)
306 6 : CYCLE
307 : END IF
308 : ! Eq. (6): solve B^A U^A = M^A U^A λ^A.
309 : CALL solve_AA_eigenproblem(M_PQ_by_kind(1:nRI_kind, 1:nRI_kind, ikind), &
310 6 : B_PQ, target_norm, U_Pp_A_atom, eigenval)
311 40 : eigenval_AA(1:SIZE(U_Pp_A_atom, 2), iatom) = eigenval(1:SIZE(U_Pp_A_atom, 2))
312 6 : bs_env%auto_ri%sizes_opt_RI(iatom) = SIZE(U_Pp_A_atom, 2)
313 4100 : U_Pp_A(1:SIZE(U_Pp_A_atom, 1), 1:SIZE(U_Pp_A_atom, 2), iatom) = U_Pp_A_atom
314 10 : DEALLOCATE (U_Pp_A_atom, eigenval, B_PQ)
315 : END DO
316 4 : CALL bs_env%para_env%sum(bs_env%auto_ri%sizes_opt_RI)
317 4 : CALL bs_env%para_env%sum(U_Pp_A)
318 4 : CALL bs_env%para_env%sum(eigenval_AA)
319 4 : CALL select_RI_sizes(eigenval_AA, RI_size_opt_AA_requested, 1, sizes_opt_RI_AA)
320 16 : bs_env%auto_ri%sizes_opt_RI(:) = sizes_opt_RI_AA
321 16 : DO iatom = 1, natom
322 193000 : U_Pp_A(:, sizes_opt_RI_AA(iatom) + 1:, iatom) = 0.0_dp
323 : END DO
324 4 : DEALLOCATE (Int_3c, M_PQ_by_kind)
325 :
326 4 : CALL timestop(handle)
327 :
328 8 : END SUBROUTINE compute_AA_optimized_RI_basis
329 :
330 : ! **************************************************************************************************
331 : !> \brief Initializes dimensions for the AA optimization.
332 : !>
333 : !> This routine maps every atom to its reference RI basis, allocates the persistent
334 : !> per-atom arrays, and determines the AA part of the requested basis size from
335 : !>
336 : !> Σ_A N_A = N_RI/2. (14)
337 : !>
338 : !> If no atom pair lies within NEIGHBOR_RADIUS, so that no AB subspace
339 : !> exists, it instead assigns
340 : !>
341 : !> Σ_A N_A^target = N_RI^target. (15)
342 : !>
343 : !> \param bs_env ...
344 : !> \param kind_of_atom ...
345 : !> \param natom ...
346 : !> \param nkind ...
347 : !> \param max_nRI_kind ...
348 : !> \param RI_size_opt_AA_requested ...
349 : ! **************************************************************************************************
350 4 : SUBROUTINE initialize_AA_optimization(bs_env, kind_of_atom, &
351 : natom, nkind, max_nRI_kind, &
352 : RI_size_opt_AA_requested)
353 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
354 : INTEGER, ALLOCATABLE, INTENT(OUT) :: kind_of_atom(:)
355 : INTEGER, INTENT(OUT) :: natom, nkind, max_nRI_kind, &
356 : RI_size_opt_AA_requested
357 :
358 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_AA_optimization'
359 :
360 : INTEGER :: AB_pair_count, handle, ikind, &
361 : RI_size_requested
362 4 : INTEGER, ALLOCATABLE :: sizes_ref_RI_by_kind(:)
363 :
364 4 : CALL timeset(routineN, handle)
365 :
366 4 : natom = bs_env%n_atom
367 4 : CPASSERT(ASSOCIATED(bs_env%ri_rs%atomic_kind_set))
368 4 : CPASSERT(ASSOCIATED(bs_env%ri_rs%particle_set))
369 4 : CPASSERT(ASSOCIATED(bs_env%para_env) .AND. natom > 0)
370 4 : CPASSERT(ALLOCATED(bs_env%basis_set_AO))
371 4 : CPASSERT(ALLOCATED(bs_env%basis_set_RI))
372 4 : nkind = SIZE(bs_env%basis_set_RI)
373 4 : CPASSERT(nkind > 0)
374 4 : CPASSERT(SIZE(bs_env%basis_set_AO) == nkind)
375 4 : CPASSERT(bs_env%do_gw_ri_rs)
376 4 : CPASSERT(.NOT. bs_env%do_periodic)
377 4 : CPASSERT(bs_env%gw_implementation /= tensor_small_cell_full_kp)
378 :
379 : CALL get_atomic_kind_set(atomic_kind_set=bs_env%ri_rs%atomic_kind_set, &
380 4 : kind_of=kind_of_atom)
381 12 : ALLOCATE (bs_env%auto_ri%sizes_opt_RI(natom))
382 12 : ALLOCATE (sizes_ref_RI_by_kind(nkind))
383 12 : DO ikind = 1, nkind
384 8 : CPASSERT(ASSOCIATED(bs_env%basis_set_AO(ikind)%gto_basis_set))
385 8 : CPASSERT(ASSOCIATED(bs_env%basis_set_RI(ikind)%gto_basis_set))
386 12 : sizes_ref_RI_by_kind(ikind) = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
387 : END DO
388 12 : max_nRI_kind = MAXVAL(sizes_ref_RI_by_kind)
389 : RI_size_requested = NINT( &
390 : MIN(bs_env%auto_ri%ri_ao_ratio, &
391 : REAL(SUM(sizes_ref_RI_by_kind(kind_of_atom)), dp)/ &
392 28 : REAL(SUM(bs_env%sizes_AO), dp))*REAL(SUM(bs_env%sizes_AO), dp))
393 4 : CALL count_AB_neighbor_pairs(bs_env, AB_pair_count)
394 4 : IF (AB_pair_count == 0) THEN
395 2 : RI_size_opt_AA_requested = RI_size_requested
396 : ELSE
397 2 : RI_size_opt_AA_requested = MAX(natom, NINT(0.5_dp*REAL(RI_size_requested, dp)))
398 : END IF
399 4 : DEALLOCATE (sizes_ref_RI_by_kind)
400 :
401 4 : CALL timestop(handle)
402 :
403 8 : END SUBROUTINE initialize_AA_optimization
404 :
405 : ! **************************************************************************************************
406 : !> \brief Counts atom pairs eligible for the AB contractions in Eqs. (8)-(13).
407 : !> \param bs_env ...
408 : !> \param AB_pair_count number of pairs separated by at most NEIGHBOR_RADIUS
409 : ! **************************************************************************************************
410 8 : SUBROUTINE count_AB_neighbor_pairs(bs_env, AB_pair_count)
411 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
412 : INTEGER, INTENT(OUT) :: AB_pair_count
413 :
414 : CHARACTER(LEN=*), PARAMETER :: routineN = 'count_AB_neighbor_pairs'
415 :
416 : INTEGER :: handle, iatom, jatom
417 :
418 8 : CALL timeset(routineN, handle)
419 :
420 8 : AB_pair_count = 0
421 24 : DO iatom = 1, bs_env%n_atom - 1
422 48 : DO jatom = iatom + 1, bs_env%n_atom
423 96 : IF (NORM2(bs_env%ri_rs%particle_set(iatom)%r - bs_env%ri_rs%particle_set(jatom)%r) <= &
424 16 : bs_env%auto_ri%neighbor_radius) THEN
425 12 : AB_pair_count = AB_pair_count + 1
426 : END IF
427 : END DO
428 : END DO
429 :
430 8 : CALL timestop(handle)
431 :
432 8 : END SUBROUTINE count_AB_neighbor_pairs
433 :
434 : ! **************************************************************************************************
435 : !> \brief Sets the RI-RS block sizes to the optimized RI basis dimensions.
436 : !> \param bs_env ...
437 : ! **************************************************************************************************
438 4 : SUBROUTINE set_optimized_RI_basis_sizes(bs_env)
439 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
440 :
441 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_optimized_RI_basis_sizes'
442 :
443 : INTEGER :: handle
444 :
445 4 : CALL timeset(routineN, handle)
446 :
447 16 : bs_env%sizes_RI(:) = bs_env%auto_ri%sizes_opt_RI
448 4 : bs_env%auto_ri%ready = .TRUE.
449 :
450 4 : CALL timestop(handle)
451 :
452 4 : END SUBROUTINE set_optimized_RI_basis_sizes
453 :
454 : ! **************************************************************************************************
455 : !> \brief Prints the optimized RI basis size and total optimization time.
456 : !> \param bs_env ...
457 : !> \param time_start ...
458 : ! **************************************************************************************************
459 4 : SUBROUTINE print_auto_ri_summary(bs_env, time_start)
460 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
461 : REAL(KIND=dp), INTENT(IN) :: time_start
462 :
463 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_auto_ri_summary'
464 :
465 : INTEGER :: handle, RI_size_optimized, unit_nr
466 : REAL(KIND=dp) :: elapsed
467 :
468 4 : CALL timeset(routineN, handle)
469 :
470 16 : RI_size_optimized = SUM(bs_env%auto_ri%sizes_opt_RI)
471 4 : elapsed = m_walltime() - time_start
472 4 : CALL bs_env%para_env%max(elapsed)
473 4 : unit_nr = bs_env%unit_nr
474 4 : IF (unit_nr > 0) THEN
475 2 : WRITE (unit_nr, '(A)') ' '
476 : WRITE (unit_nr, '(T2,A,T74,I9)') &
477 2 : 'AUTO_RI| Number of automatic RI functions for χ, ε, W:', RI_size_optimized
478 : WRITE (unit_nr, '(T2,A,T71,F8.1,A)') &
479 2 : 'AUTO_RI| Total RI basis optimization time:', elapsed, ' s'
480 : END IF
481 :
482 4 : CALL timestop(handle)
483 :
484 4 : END SUBROUTINE print_auto_ri_summary
485 :
486 : ! **************************************************************************************************
487 : !> \brief Adds REGULARIZATION_RI to a two-center RI metric matrix (P|Q).
488 : !> \param basis ...
489 : !> \param regularization ...
490 : !> \param M_PQ reference RI metric M_PQ
491 : !> \param offset ...
492 : ! **************************************************************************************************
493 14 : SUBROUTINE add_RI_metric_regularization(basis, regularization, M_PQ, offset)
494 : USE basis_set_types, ONLY: gto_basis_set_type
495 : TYPE(gto_basis_set_type), INTENT(IN) :: basis
496 : REAL(KIND=dp), INTENT(IN) :: regularization
497 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: M_PQ
498 : INTEGER, INTENT(IN), OPTIONAL :: offset
499 :
500 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_RI_metric_regularization'
501 :
502 : INTEGER :: first, handle, i, iset, nset, start
503 14 : INTEGER, DIMENSION(:), POINTER :: npgf, nsgf_set
504 14 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
505 : REAL(KIND=dp) :: shift
506 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
507 :
508 14 : CALL timeset(routineN, handle)
509 :
510 14 : start = 0
511 14 : IF (PRESENT(offset)) start = offset
512 : CALL get_gto_basis_set(basis, first_sgf=first_sgf, npgf=npgf, nset=nset, &
513 14 : nsgf_set=nsgf_set, zet=zet)
514 14 : CPASSERT(start >= 0 .AND. start + basis%nsgf <= SIZE(M_PQ, 1))
515 14 : CPASSERT(start + basis%nsgf <= SIZE(M_PQ, 2))
516 302 : DO iset = 1, nset
517 576 : shift = regularization*MAX(1.0_dp, 1.0_dp/MINVAL(zet(1:npgf(iset), iset)))
518 288 : first = start + first_sgf(1, iset)
519 1510 : DO i = first, first + nsgf_set(iset) - 1
520 1496 : M_PQ(i, i) = M_PQ(i, i) + shift
521 : END DO
522 : END DO
523 :
524 14 : CALL timestop(handle)
525 :
526 14 : END SUBROUTINE add_RI_metric_regularization
527 :
528 : ! **************************************************************************************************
529 : !> \brief Copies atom-local diagonal blocks M^AA from a distributed metric M.
530 : !> \param M_pq_local locally stored entries of the optimized RI metric M_pq
531 : !> \param M_pq_AA_local locally stored entries of the atom-block-diagonal metric
532 : !> \param row_indices global indices of the locally owned rows
533 : !> \param col_indices global indices of the locally owned columns
534 : !> \param sizes_RI optimized RI functions in each atomic matrix block
535 : ! **************************************************************************************************
536 4 : SUBROUTINE copy_M_pq_AA_blocks(M_pq_local, M_pq_AA_local, row_indices, &
537 4 : col_indices, sizes_RI)
538 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: M_pq_local
539 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: M_pq_AA_local
540 : INTEGER, DIMENSION(:), INTENT(IN) :: row_indices, col_indices, sizes_RI
541 :
542 : CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_M_pq_AA_blocks'
543 :
544 : INTEGER :: handle, i, j, last_index
545 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: block_ends, col_atom, row_atom
546 :
547 4 : CALL timeset(routineN, handle)
548 :
549 4 : CPASSERT(SIZE(M_pq_local, 1) == SIZE(row_indices))
550 4 : CPASSERT(SIZE(M_pq_local, 2) == SIZE(col_indices))
551 12 : CPASSERT(ALL(SHAPE(M_pq_AA_local) == SHAPE(M_pq_local)))
552 16 : CPASSERT(SIZE(sizes_RI) > 0 .AND. ALL(sizes_RI > 0))
553 :
554 12 : ALLOCATE (block_ends(SIZE(sizes_RI)))
555 4 : block_ends(1) = sizes_RI(1)
556 12 : DO i = 2, SIZE(sizes_RI)
557 12 : block_ends(i) = block_ends(i - 1) + sizes_RI(i)
558 : END DO
559 4 : last_index = block_ends(SIZE(block_ends))
560 33 : CPASSERT(ALL(row_indices >= 1 .AND. row_indices <= last_index))
561 62 : CPASSERT(ALL(col_indices >= 1 .AND. col_indices <= last_index))
562 :
563 20 : ALLOCATE (row_atom(SIZE(row_indices)), col_atom(SIZE(col_indices)))
564 33 : DO i = 1, SIZE(row_indices)
565 33 : row_atom(i) = M_pq_index_atom(row_indices(i), block_ends)
566 : END DO
567 62 : DO j = 1, SIZE(col_indices)
568 62 : col_atom(j) = M_pq_index_atom(col_indices(j), block_ends)
569 : END DO
570 :
571 507 : M_pq_AA_local = 0.0_dp
572 62 : DO j = 1, SIZE(col_indices)
573 507 : DO i = 1, SIZE(row_indices)
574 503 : IF (row_atom(i) == col_atom(j)) THEN
575 193 : M_pq_AA_local(i, j) = M_pq_local(i, j)
576 : END IF
577 : END DO
578 : END DO
579 :
580 4 : CALL timestop(handle)
581 :
582 8 : END SUBROUTINE copy_M_pq_AA_blocks
583 :
584 : ! **************************************************************************************************
585 : !> \brief Forms f_i=1/sqrt(λ_i) for eigenvalues of the RI metric matrix M_PQ.
586 : !> \param eigenval generalized metric eigenvalues
587 : !> \param factors inverse-square-root weights, or zero for discarded modes
588 : ! **************************************************************************************************
589 4 : SUBROUTINE compute_optimized_M_pq_inv_factors(eigenval, factors)
590 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenval
591 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
592 : INTENT(OUT) :: factors
593 :
594 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_optimized_M_pq_inv_factors'
595 :
596 : INTEGER :: handle, i, n_dependent
597 : REAL(KIND=dp) :: largest, threshold
598 :
599 4 : CALL timeset(routineN, handle)
600 :
601 4 : CPASSERT(SIZE(eigenval) > 0)
602 :
603 62 : largest = MAXVAL(eigenval)
604 4 : CPASSERT(largest > 0.0_dp)
605 4 : threshold = optimized_M_eigenval_threshold*largest
606 12 : ALLOCATE (factors(SIZE(eigenval)))
607 4 : factors = 0.0_dp
608 4 : n_dependent = 0
609 62 : DO i = 1, SIZE(eigenval)
610 62 : IF (eigenval(i) > 0.0_dp .AND. eigenval(i) >= threshold) THEN
611 58 : factors(i) = 1.0_dp/SQRT(eigenval(i))
612 : ELSE
613 0 : n_dependent = n_dependent + 1
614 : END IF
615 : END DO
616 4 : CPASSERT(n_dependent < SIZE(eigenval))
617 :
618 4 : CALL timestop(handle)
619 :
620 4 : END SUBROUTINE compute_optimized_M_pq_inv_factors
621 :
622 : ! **************************************************************************************************
623 : !> \brief Finds the atom owning a global RI metric index.
624 : !> \param global_index global metric index
625 : !> \param block_ends inclusive end index of each consecutive atom block
626 : !> \return atom index
627 : ! **************************************************************************************************
628 87 : PURE FUNCTION M_pq_index_atom(global_index, block_ends) RESULT(atom_index)
629 : INTEGER, INTENT(IN) :: global_index
630 : INTEGER, DIMENSION(:), INTENT(IN) :: block_ends
631 : INTEGER :: atom_index
632 :
633 : INTEGER :: high, low, middle
634 :
635 87 : low = 1
636 87 : high = SIZE(block_ends)
637 243 : DO WHILE (low < high)
638 156 : middle = (low + high)/2
639 243 : IF (global_index <= block_ends(middle)) THEN
640 : high = middle
641 : ELSE
642 39 : low = middle + 1
643 : END IF
644 : END DO
645 87 : atom_index = low
646 :
647 87 : END FUNCTION M_pq_index_atom
648 :
649 : ! **************************************************************************************************
650 : !> \brief Computes X=U diag(1/sqrt(λ)) with X^T M X=I.
651 : !>
652 : !> The columns of X span the non-null eigenspace of the symmetric
653 : !> two-center RI metric matrix M=U diag(λ) U^T.
654 : !> \param M_PQ ...
655 : !> \param X_Pi ...
656 : !> \param rank ...
657 : ! **************************************************************************************************
658 21 : SUBROUTINE compute_reference_M_PQ_inv_factor(M_PQ, X_Pi, rank)
659 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: M_PQ
660 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
661 : INTENT(OUT) :: X_Pi
662 : INTEGER, INTENT(OUT) :: rank
663 :
664 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_reference_M_PQ_inv_factor'
665 :
666 : INTEGER :: column, handle, i, n
667 : REAL(KIND=dp) :: largest
668 21 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
669 21 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvec
670 :
671 21 : CALL timeset(routineN, handle)
672 :
673 21 : n = SIZE(M_PQ, 1)
674 21 : CPASSERT(n > 0 .AND. SIZE(M_PQ, 2) == n)
675 499976 : ALLOCATE (eigenvec, source=0.5_dp*(M_PQ + TRANSPOSE(M_PQ)))
676 63 : ALLOCATE (eigenval(n))
677 21 : CALL diamat_all(eigenvec, eigenval, dac=.TRUE.)
678 1943 : largest = MAXVAL(eigenval)
679 21 : IF (largest <= 0.0_dp) CPABORT("AUTO_RI metric has no positive non-null range")
680 1943 : IF (MINVAL(eigenval) < -negative_eigenval_tolerance*largest) THEN
681 0 : CPABORT("AUTO_RI metric has no positive non-null range")
682 : END IF
683 1943 : rank = COUNT(eigenval > reference_M_eigenval_threshold*largest)
684 21 : IF (rank == 0) CPABORT("AUTO_RI metric has no positive non-null range")
685 84 : ALLOCATE (X_Pi(n, rank))
686 21 : column = 0
687 1943 : DO i = 1, n
688 1922 : IF (eigenval(i) <= reference_M_eigenval_threshold*largest) CYCLE
689 1922 : column = column + 1
690 249967 : X_Pi(:, column) = eigenvec(:, i)/SQRT(eigenval(i))
691 : END DO
692 :
693 21 : CALL timestop(handle)
694 :
695 42 : END SUBROUTINE compute_reference_M_PQ_inv_factor
696 :
697 : ! **************************************************************************************************
698 : !> \brief Solves eigenvalue problem B^A U^A=M^A U^A λ^A in Eq. (6).
699 : !>
700 : !> Columns are ordered from the largest to the smallest λ^A.
701 : !> Columns are retained when λ_i/λ_max >= eigenval_threshold.
702 : !> \param M_PQ metric of the large reference RI basis
703 : !> \param B_PQ normalized atom-local target matrix B^A
704 : !> \param target_normalization ...
705 : !> \param U_Pp_A selected AA contractions from the large RI basis P to optimized functions p
706 : !> \param eigenval normalized target eigenvalues in descending order
707 : ! **************************************************************************************************
708 6 : SUBROUTINE solve_AA_eigenproblem(M_PQ, B_PQ, target_normalization, U_Pp_A, eigenval)
709 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: M_PQ, B_PQ
710 : REAL(KIND=dp), INTENT(IN) :: target_normalization
711 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
712 : INTENT(OUT) :: U_Pp_A
713 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
714 : INTENT(OUT) :: eigenval
715 :
716 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_AA_eigenproblem'
717 :
718 : INTEGER :: handle, i, j, M_PQ_rank, n, nret
719 : REAL(KIND=dp) :: largest_score, relative_fitting_defect
720 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_B, eigenval_M
721 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: B_ij, M_pp_check, U_ij, U_Pp_A_rotated, &
722 6 : X_Pi
723 :
724 6 : CALL timeset(routineN, handle)
725 :
726 6 : n = SIZE(M_PQ, 1)
727 6 : CPASSERT(n > 0)
728 6 : CPASSERT(SIZE(M_PQ, 2) == n)
729 18 : CPASSERT(ALL(SHAPE(B_PQ) == [n, n]))
730 6 : CPASSERT(target_normalization > 0.0_dp)
731 6 : CALL compute_reference_M_PQ_inv_factor(M_PQ, X_Pi, M_PQ_rank)
732 :
733 24 : ALLOCATE (B_ij(M_PQ_rank, M_PQ_rank))
734 6 : B_ij(:, :) = MATMUL(TRANSPOSE(X_Pi), &
735 90712 : MATMUL(0.5_dp*(B_PQ + TRANSPOSE(B_PQ)), X_Pi))/target_normalization
736 45368 : ALLOCATE (U_ij, source=B_ij)
737 18 : ALLOCATE (eigenval_B(M_PQ_rank))
738 6 : CALL diamat_all(U_ij, eigenval_B, dac=.TRUE.)
739 482 : largest_score = MAXVAL(eigenval_B)
740 6 : IF (largest_score <= 0.0_dp) THEN
741 0 : CPABORT("AUTO_RI conditioned atom-local eigensolver failed")
742 : END IF
743 482 : IF (MINVAL(eigenval_B) < -conditioned_negative_eigenval_tolerance*largest_score) THEN
744 0 : CPABORT("AUTO_RI conditioned atom-local eigensolver failed")
745 : END IF
746 482 : eigenval_B(:) = MAX(eigenval_B, 0.0_dp)
747 482 : nret = COUNT(eigenval_B/largest_score >= relative_target_eigenval_threshold)
748 30 : ALLOCATE (U_Pp_A(n, nret), eigenval(M_PQ_rank))
749 482 : DO i = 1, M_PQ_rank
750 482 : eigenval(i) = eigenval_B(M_PQ_rank - i + 1)
751 : END DO
752 40 : DO i = 1, nret
753 34 : j = M_PQ_rank - i + 1
754 1519020 : U_Pp_A(:, i) = MATMUL(X_Pi, U_ij(:, j))
755 : END DO
756 :
757 : relative_fitting_defect = 1.0_dp
758 40 : DO i = 1, nret
759 34 : j = M_PQ_rank - i + 1
760 : relative_fitting_defect = relative_fitting_defect - &
761 511848 : DOT_PRODUCT(U_ij(:, j), MATMUL(B_ij, U_ij(:, j)))
762 : END DO
763 6 : IF (ABS(relative_fitting_defect) < 100.0_dp*EPSILON(1.0_dp)) THEN
764 : relative_fitting_defect = 0.0_dp
765 : END IF
766 : ! A projection cannot capture more than the exact product norm.
767 6 : IF (relative_fitting_defect < -1.0E-8_dp .OR. &
768 : relative_fitting_defect > 1.0_dp + 1.0E-8_dp) THEN
769 0 : CPABORT("AUTO_RI conditioned atom-local eigensolver failed")
770 : END IF
771 6 : IF (nret > 0) THEN
772 48 : ALLOCATE (M_pp_check(nret, nret), eigenval_M(nret), U_Pp_A_rotated(n, nret))
773 6 : M_pp_check(:, :) = MATMUL(TRANSPOSE(U_Pp_A), &
774 70682 : MATMUL(0.5_dp*(M_PQ + TRANSPOSE(M_PQ)), U_Pp_A))
775 6 : CALL diamat_all(M_pp_check, eigenval_M, dac=.TRUE.)
776 40 : IF (MINVAL(eigenval_M) <= reference_M_eigenval_threshold) THEN
777 0 : CPABORT("AUTO_RI conditioned atom-local eigensolver failed")
778 : END IF
779 894 : U_Pp_A_rotated(:, :) = MATMUL(U_Pp_A, M_pp_check)
780 40 : DO i = 1, nret
781 4100 : U_Pp_A_rotated(:, i) = U_Pp_A_rotated(:, i)/SQRT(eigenval_M(i))
782 : END DO
783 894 : U_Pp_A(:, :) = MATMUL(U_Pp_A_rotated, TRANSPOSE(M_pp_check))
784 6 : M_pp_check(:, :) = MATMUL(TRANSPOSE(U_Pp_A), &
785 71122 : MATMUL(0.5_dp*(M_PQ + TRANSPOSE(M_PQ)), U_Pp_A))
786 40 : DO i = 1, nret
787 40 : M_pp_check(i, i) = M_pp_check(i, i) - 1.0_dp
788 : END DO
789 494 : IF (MAXVAL(ABS(M_pp_check)) >= 1.0E-7_dp) THEN
790 0 : CPABORT("AUTO_RI conditioned atom-local eigensolver failed")
791 : END IF
792 : END IF
793 :
794 6 : CALL timestop(handle)
795 :
796 12 : END SUBROUTINE solve_AA_eigenproblem
797 :
798 : ! **************************************************************************************************
799 : !> \brief Solves eigenvalue problem B^AB U^AB=U^AB λ^AB of Eq. (8).
800 : !>
801 : !> Columns are ordered from the largest to the smallest λ^AB.
802 : !> \param M_PQ combined AB reference-basis metric
803 : !> \param B_PQ normalized two-center target matrix B^AB
804 : !> \param U_Pp_AA retained AA columns embedded in the combined AB reference space
805 : !> \param U_Pp_AB selected AB contractions from the large RI basis P to optimized functions p
806 : !> \param eigenval normalized scores in descending order
807 : ! **************************************************************************************************
808 3 : SUBROUTINE solve_AB_eigenproblem(M_PQ, B_PQ, U_Pp_AA, U_Pp_AB, eigenval)
809 : REAL(KIND=dp), INTENT(IN) :: M_PQ(:, :), B_PQ(:, :), U_Pp_AA(:, :)
810 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: U_Pp_AB(:, :), eigenval(:)
811 :
812 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_AB_eigenproblem'
813 :
814 : INTEGER :: handle, i, k, n, nkeep, rank
815 : REAL(KIND=dp) :: threshold
816 3 : REAL(KIND=dp), ALLOCATABLE :: B_ij(:, :), B_ij_residual(:, :), eigenval_residual(:), &
817 3 : M_AA(:, :), P_ij(:, :), Q_ia(:, :), U_Aa(:, :), U_ia(:, :), X_Pi(:, :)
818 :
819 3 : CALL timeset(routineN, handle)
820 :
821 3 : n = SIZE(M_PQ, 1)
822 3 : CALL compute_reference_M_PQ_inv_factor(M_PQ, X_Pi, rank)
823 :
824 : ! Xᵀ M X=I; Q spans the retained atomic columns in these orthonormal coordinates.
825 12 : ALLOCATE (Q_ia(rank, SIZE(U_Pp_AA, 2)))
826 121560 : Q_ia(:, :) = MATMUL(TRANSPOSE(X_Pi), MATMUL(M_PQ, U_Pp_AA))
827 12 : ALLOCATE (P_ij(rank, rank), source=0.0_dp)
828 479 : DO i = 1, rank
829 479 : P_ij(i, i) = 1.0_dp
830 : END DO
831 3 : IF (SIZE(U_Pp_AA, 2) > 0) THEN
832 12 : ALLOCATE (M_AA(SIZE(U_Pp_AA, 2), SIZE(U_Pp_AA, 2)))
833 24017 : M_AA(:, :) = MATMUL(TRANSPOSE(Q_ia), Q_ia)
834 3 : CALL compute_reference_M_PQ_inv_factor(M_AA, U_Aa, k)
835 12 : ALLOCATE (U_ia(rank, SIZE(U_Aa, 2)))
836 47881 : U_ia(:, :) = MATMUL(Q_ia, U_Aa)
837 643217 : P_ij(:, :) = P_ij - MATMUL(U_ia, TRANSPOSE(U_ia))
838 3 : DEALLOCATE (U_ia)
839 : END IF
840 9 : ALLOCATE (B_ij(rank, rank))
841 6 : B_ij(:, :) = MATMUL(TRANSPOSE(X_Pi), MATMUL(B_PQ, X_Pi))
842 9 : ALLOCATE (B_ij_residual(rank, rank))
843 6 : B_ij_residual(:, :) = MATMUL(P_ij, MATMUL(B_ij, P_ij))
844 9 : ALLOCATE (eigenval_residual(rank))
845 3 : CALL diamat_all(B_ij_residual, eigenval_residual, dac=.TRUE.)
846 : threshold = MAX(AB_absolute_target_eigenval_threshold, &
847 : relative_target_eigenval_threshold*MAXVAL(eigenval_residual), &
848 80036 : 100.0_dp*EPSILON(1.0_dp)*MAXVAL(ABS(B_ij)))
849 479 : nkeep = COUNT(eigenval_residual > threshold)
850 18 : ALLOCATE (U_Pp_AB(n, nkeep), eigenval(nkeep))
851 3 : k = 0
852 479 : DO i = rank, 1, -1
853 476 : IF (eigenval_residual(i) <= threshold) CYCLE
854 14 : k = k + 1
855 1349880 : U_Pp_AB(:, k) = MATMUL(X_Pi, MATMUL(P_ij, B_ij_residual(:, i)))
856 17 : eigenval(k) = eigenval_residual(i)
857 : END DO
858 :
859 3 : CALL timestop(handle)
860 :
861 6 : END SUBROUTINE solve_AB_eigenproblem
862 :
863 : ! **************************************************************************************************
864 : !> \brief Selects the requested RI sizes from the largest positive eigenvalues.
865 : !> \param eigenval_by_block descending eigenvalues; zero pads unavailable modes
866 : !> \param RI_size_requested requested total number of functions
867 : !> \param RI_size_minimum minimum retained size in every nonempty block
868 : !> \param sizes_RI_selected number retained from each block
869 : ! **************************************************************************************************
870 8 : SUBROUTINE select_RI_sizes(eigenval_by_block, RI_size_requested, RI_size_minimum, &
871 8 : sizes_RI_selected)
872 : REAL(KIND=dp), INTENT(IN) :: eigenval_by_block(:, :)
873 : INTEGER, INTENT(IN) :: RI_size_requested, RI_size_minimum
874 : INTEGER, INTENT(OUT) :: sizes_RI_selected(:)
875 :
876 : CHARACTER(LEN=*), PARAMETER :: routineN = 'select_RI_sizes'
877 :
878 : INTEGER :: block, handle, i, n, nmode, remaining
879 8 : INTEGER, ALLOCATABLE :: indices(:)
880 8 : REAL(KIND=dp), ALLOCATABLE :: values(:)
881 :
882 8 : CALL timeset(routineN, handle)
883 :
884 8 : nmode = SIZE(eigenval_by_block, 1)
885 24 : n = SIZE(eigenval_by_block)
886 8 : CPASSERT(SIZE(sizes_RI_selected) == SIZE(eigenval_by_block, 2))
887 8 : CPASSERT(RI_size_minimum >= 0)
888 26 : sizes_RI_selected(:) = 0
889 8 : IF (n > 0) THEN
890 30 : ALLOCATE (values(n), indices(n))
891 3084 : values(:) = -RESHAPE(eigenval_by_block, [n])
892 24 : DO block = 1, SIZE(sizes_RI_selected)
893 : sizes_RI_selected(block) = &
894 3090 : MIN(RI_size_minimum, COUNT(eigenval_by_block(:, block) > 0.0_dp))
895 : values((block - 1)*nmode + 1: &
896 36 : (block - 1)*nmode + sizes_RI_selected(block)) = 0.0_dp
897 : END DO
898 6 : CALL sort(values, n, indices)
899 3096 : remaining = MIN(MAX(0, RI_size_requested - SUM(sizes_RI_selected)), COUNT(values < 0.0_dp))
900 52 : DO i = 1, remaining
901 46 : block = (indices(i) - 1)/nmode + 1
902 52 : sizes_RI_selected(block) = sizes_RI_selected(block) + 1
903 : END DO
904 : END IF
905 :
906 8 : CALL timestop(handle)
907 :
908 8 : END SUBROUTINE select_RI_sizes
909 :
910 : ! **************************************************************************************************
911 : !> \brief Computes orbital factors for the occupied-empty weight in Eq. (5).
912 : !>
913 : !> weights(n) = 1/|ε_n-E_F|, so weights(i)*weights(a) = w_ia.
914 : !>
915 : !> \param eigenval_MO molecular-orbital eigenvalues, occupied first
916 : !> \param homo Last occupied state in this spin channel.
917 : !> \param weights factors whose occupied-empty products give w_ia
918 : ! **************************************************************************************************
919 15 : SUBROUTINE compute_orbital_weights(eigenval_MO, homo, weights)
920 : REAL(KIND=dp), INTENT(IN) :: eigenval_MO(:)
921 : INTEGER, INTENT(IN) :: homo
922 : REAL(KIND=dp), INTENT(OUT) :: weights(:)
923 :
924 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_orbital_weights'
925 :
926 : INTEGER :: handle, n
927 : REAL(KIND=dp) :: gap, middle
928 :
929 15 : CALL timeset(routineN, handle)
930 :
931 15 : CPASSERT(SIZE(weights) == SIZE(eigenval_MO))
932 15 : IF (homo < 1 .OR. homo >= SIZE(eigenval_MO)) THEN
933 0 : CPABORT("AUTO_RI orbital weights require occupied and empty states")
934 : END IF
935 120 : IF (.NOT. ALL(ieee_is_finite(eigenval_MO))) THEN
936 0 : CPABORT("AUTO_RI orbital weights require finite eigenvalues")
937 : END IF
938 15 : gap = eigenval_MO(homo + 1) - eigenval_MO(homo)
939 15 : IF (.NOT. ieee_is_finite(gap) .OR. gap <= 0.0_dp) THEN
940 0 : CPABORT("AUTO_RI orbital weights require a finite positive gap")
941 : END IF
942 15 : middle = 0.5_dp*(eigenval_MO(homo) + eigenval_MO(homo + 1))
943 120 : DO n = 1, SIZE(eigenval_MO)
944 120 : weights(n) = 1.0_dp/ABS(eigenval_MO(n) - middle)
945 : END DO
946 :
947 15 : CALL timestop(handle)
948 :
949 15 : END SUBROUTINE compute_orbital_weights
950 :
951 : ! **************************************************************************************************
952 : !> \brief Computes the occupied-empty contribution B_PQ^C in Eqs. (3) or (10).
953 : !>
954 : !> B_PQ^C = Σ_ia w_ia (P|ρ_ia)(ρ_ia|Q).
955 : !>
956 : !> \param Int_3c atom-local or atom-pair three-center integrals (μν|P)
957 : !> \param C_mu_n molecular-orbital coefficients restricted to the relevant AO rows
958 : !> \param weights Per-orbital weights for this spin channel.
959 : !> \param homo Last occupied state.
960 : !> \param B_PQ unnormalized transition matrix for this spin channel
961 : ! **************************************************************************************************
962 12 : SUBROUTINE compute_B_PQ_transition(Int_3c, C_mu_n, weights, homo, B_PQ)
963 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN) :: Int_3c(:, :, :), C_mu_n(:, :), weights(:)
964 : INTEGER, INTENT(IN) :: homo
965 : REAL(KIND=dp), INTENT(OUT) :: B_PQ(:, :)
966 :
967 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_B_PQ_transition'
968 :
969 : INTEGER :: handle, n, nAO, nRI, nstates
970 12 : REAL(KIND=dp), ALLOCATABLE :: empty(:, :), occupied(:, :), scaled(:, :)
971 : TYPE(local_gemm_ctxt_type) :: gemm_context
972 :
973 12 : CALL timeset(routineN, handle)
974 :
975 12 : CALL gemm_context%create(LOCAL_GEMM_PU_HOST)
976 12 : nAO = SIZE(C_mu_n, 1)
977 12 : nstates = SIZE(C_mu_n, 2)
978 12 : nRI = SIZE(Int_3c, 3)
979 12 : CPASSERT(SIZE(Int_3c, 1) == nAO .AND. SIZE(Int_3c, 2) == nAO)
980 12 : CPASSERT(SIZE(B_PQ, 1) == nRI .AND. SIZE(B_PQ, 2) == nRI)
981 12 : CPASSERT(SIZE(weights) == nstates .AND. homo > 0 .AND. homo < nstates)
982 108 : ALLOCATE (scaled(nAO, nstates), occupied(nAO, nAO), empty(nAO, nAO))
983 96 : DO n = 1, nstates
984 390 : scaled(:, n) = C_mu_n(:, n)*SQRT(weights(n))
985 : END DO
986 : CALL gemm_context%gemm( &
987 : 'N', 'T', nAO, nAO, homo, 1.0_dp, scaled, nAO, scaled, nAO, &
988 12 : 0.0_dp, occupied, nAO)
989 : CALL gemm_context%gemm('N', 'T', nAO, nAO, nstates - homo, 1.0_dp, &
990 : scaled(:, homo + 1:), nAO, scaled(:, homo + 1:), nAO, &
991 12 : 0.0_dp, empty, nAO)
992 12 : CALL contract_B_PQ(Int_3c, occupied, empty, B_PQ)
993 12 : CALL gemm_context%destroy()
994 :
995 12 : CALL timestop(handle)
996 :
997 24 : END SUBROUTINE compute_B_PQ_transition
998 :
999 : ! **************************************************************************************************
1000 : !> \brief Computes the frontier-exchange contribution B_PQ^X in Eqs. (4) or (11).
1001 : !>
1002 : !> B_PQ^X = Σ_{n∈W_occ} Σ_i^occ (P|ρ_ni)(ρ_ni|Q).
1003 : !>
1004 : !> \param Int_3c atom-local or atom-pair three-center integrals (μν|P)
1005 : !> \param C_mu_n molecular-orbital coefficients restricted to the relevant AO rows
1006 : !> \param eigenval_MO molecular-orbital eigenvalues, occupied first, in Hartree
1007 : !> \param homo Last occupied state in this spin channel.
1008 : !> \param occ_energy_window frontier-window extension in Hartree
1009 : !> \param B_PQ unnormalized exchange matrix
1010 : ! **************************************************************************************************
1011 12 : SUBROUTINE compute_B_PQ_exchange(Int_3c, C_mu_n, eigenval_MO, homo, occ_energy_window, B_PQ)
1012 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN) :: Int_3c(:, :, :), C_mu_n(:, :), &
1013 : eigenval_MO(:)
1014 : INTEGER, INTENT(IN) :: homo
1015 : REAL(KIND=dp), INTENT(IN) :: occ_energy_window
1016 : REAL(KIND=dp), INTENT(OUT) :: B_PQ(:, :)
1017 :
1018 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_B_PQ_exchange'
1019 :
1020 : INTEGER :: handle, n, nAO, nstates
1021 : REAL(KIND=dp), ALLOCATABLE :: frontier(:, :), occupied(:, :), &
1022 12 : selected(:, :)
1023 : TYPE(local_gemm_ctxt_type) :: gemm_context
1024 :
1025 12 : CALL timeset(routineN, handle)
1026 :
1027 204464 : B_PQ(:, :) = 0.0_dp
1028 12 : nstates = SIZE(eigenval_MO)
1029 12 : IF (homo < 1 .OR. homo >= nstates) THEN
1030 0 : CPABORT("AUTO_RI exchange target requires occupied and empty states")
1031 : END IF
1032 96 : IF (.NOT. ieee_is_finite(occ_energy_window) .OR. &
1033 : .NOT. ALL(ieee_is_finite(eigenval_MO))) THEN
1034 0 : CPABORT("AUTO_RI exchange target requires finite input")
1035 : END IF
1036 12 : IF (occ_energy_window < 0.0_dp .OR. eigenval_MO(homo + 1) <= eigenval_MO(homo)) THEN
1037 0 : CPABORT("AUTO_RI exchange target requires a valid frontier window")
1038 : END IF
1039 12 : nAO = SIZE(C_mu_n, 1)
1040 12 : CPASSERT(SIZE(C_mu_n, 2) == nstates)
1041 108 : ALLOCATE (occupied(nAO, nAO), frontier(nAO, nAO), selected(nAO, nstates), source=0.0_dp)
1042 96 : DO n = 1, nstates
1043 84 : IF (eigenval_MO(n) < eigenval_MO(homo) - occ_energy_window .OR. &
1044 : eigenval_MO(n) > eigenval_MO(homo + 1) + occ_energy_window) CYCLE
1045 240 : selected(:, n) = C_mu_n(:, n)
1046 : END DO
1047 12 : CALL gemm_context%create(LOCAL_GEMM_PU_HOST)
1048 : ! D_μν = Σ_i^occ C_μi C_νi.
1049 : ! F_μν = Σ_{n: ε_n ∈ [ε_H−E_win,ε_L+E_win]} C_μn C_νn.
1050 : CALL gemm_context%gemm('N', 'T', nAO, nAO, homo, 1.0_dp, C_mu_n, nAO, C_mu_n, nAO, &
1051 12 : 0.0_dp, occupied, nAO)
1052 : CALL gemm_context%gemm('N', 'T', nAO, nAO, nstates, 1.0_dp, selected, nAO, selected, nAO, &
1053 12 : 0.0_dp, frontier, nAO)
1054 12 : CALL contract_B_PQ(Int_3c, occupied, frontier, B_PQ)
1055 12 : CALL gemm_context%destroy()
1056 :
1057 12 : CALL timestop(handle)
1058 :
1059 24 : END SUBROUTINE compute_B_PQ_exchange
1060 :
1061 : ! **************************************************************************************************
1062 : !> \brief Contracts AO density matrices and (μν|P) to form B_PQ.
1063 : !>
1064 : !> B_PQ = Σ_μν (μν|P) [left (..|Q) right]_μν.
1065 : !>
1066 : !> \param Int_3c atom-local or atom-pair three-center integrals (μν|P)
1067 : !> \param left First density block, for occupied partners.
1068 : !> \param right Second density block, for empty or hard-window target states.
1069 : !> \param B_PQ symmetric RI-basis product matrix
1070 : ! **************************************************************************************************
1071 24 : SUBROUTINE contract_B_PQ(Int_3c, left, right, B_PQ)
1072 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN) :: Int_3c(:, :, :), left(:, :), right(:, :)
1073 : REAL(KIND=dp), INTENT(OUT) :: B_PQ(:, :)
1074 :
1075 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_B_PQ'
1076 :
1077 : INTEGER :: b, handle, mu, naux, nu, p, q
1078 24 : REAL(KIND=dp), ALLOCATABLE :: weighted(:, :), work(:, :)
1079 : TYPE(local_gemm_ctxt_type) :: gemm_context
1080 :
1081 24 : CALL timeset(routineN, handle)
1082 :
1083 24 : b = SIZE(Int_3c, 1)
1084 24 : naux = SIZE(Int_3c, 3)
1085 24 : CPASSERT(SIZE(Int_3c, 2) == b)
1086 24 : CPASSERT(SIZE(left, 1) == b .AND. SIZE(left, 2) == b)
1087 24 : CPASSERT(SIZE(right, 1) == b .AND. SIZE(right, 2) == b)
1088 24 : CPASSERT(SIZE(B_PQ, 1) == naux .AND. SIZE(B_PQ, 2) == naux)
1089 144 : ALLOCATE (work(b, b), weighted(b, b))
1090 24 : CALL gemm_context%create(LOCAL_GEMM_PU_HOST)
1091 2880 : DO q = 1, naux
1092 : ! B_PQ = Σ_μν (μν|P) [left (..|Q) right]_μν.
1093 : CALL gemm_context%gemm('N', 'N', b, b, b, 1.0_dp, left, b, Int_3c(:, :, q), b, &
1094 2856 : 0.0_dp, work, b)
1095 2856 : CALL gemm_context%gemm('N', 'N', b, b, b, 1.0_dp, work, b, right, b, 0.0_dp, weighted, b)
1096 207332 : DO p = 1, q
1097 204452 : B_PQ(p, q) = 0.0_dp
1098 1238860 : DO nu = 1, b
1099 7023228 : DO mu = 1, b
1100 6818776 : B_PQ(p, q) = B_PQ(p, q) + Int_3c(mu, nu, p)*weighted(mu, nu)
1101 : END DO
1102 : END DO
1103 207308 : B_PQ(q, p) = B_PQ(p, q)
1104 : END DO
1105 : END DO
1106 24 : CALL gemm_context%destroy()
1107 :
1108 24 : CALL timestop(handle)
1109 :
1110 48 : END SUBROUTINE contract_B_PQ
1111 :
1112 : ! **************************************************************************************************
1113 : !> \brief Adds separately normalized transition and exchange matrices.
1114 : !> \param X_Pi inverse metric factor, X_Pi X_Qi = M_PQ^+
1115 : !> \param B_PQ_transition transition matrix, replaced by the combined target
1116 : !> \param B_PQ_exchange unnormalized hard-window exchange covariance
1117 : !> \param target_norm number of nonzero normalized terms
1118 : ! **************************************************************************************************
1119 6 : SUBROUTINE combine_transition_and_exchange_targets( &
1120 6 : X_Pi, B_PQ_transition, B_PQ_exchange, target_norm)
1121 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN) :: X_Pi(:, :)
1122 : REAL(KIND=dp), INTENT(INOUT) :: B_PQ_transition(:, :)
1123 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN) :: B_PQ_exchange(:, :)
1124 : REAL(KIND=dp), INTENT(OUT) :: target_norm
1125 :
1126 : CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_transition_and_exchange_targets'
1127 :
1128 : INTEGER :: handle, k
1129 : REAL(KIND=dp) :: exchange_norm, transition_norm
1130 6 : REAL(KIND=dp), ALLOCATABLE :: work(:, :)
1131 :
1132 6 : CALL timeset(routineN, handle)
1133 :
1134 24 : ALLOCATE (work(SIZE(X_Pi, 1), SIZE(X_Pi, 2)))
1135 6 : work(:, :) = MATMUL(B_PQ_transition, X_Pi)
1136 6 : transition_norm = 0.0_dp
1137 482 : DO k = 1, SIZE(X_Pi, 2)
1138 45350 : transition_norm = transition_norm + DOT_PRODUCT(X_Pi(:, k), work(:, k))
1139 : END DO
1140 6 : work(:, :) = MATMUL(B_PQ_exchange, X_Pi)
1141 6 : exchange_norm = 0.0_dp
1142 482 : DO k = 1, SIZE(X_Pi, 2)
1143 45350 : exchange_norm = exchange_norm + DOT_PRODUCT(X_Pi(:, k), work(:, k))
1144 : END DO
1145 6 : target_norm = 0.0_dp
1146 : IF (.NOT. ieee_is_finite(transition_norm) .OR. &
1147 : .NOT. ieee_is_finite(exchange_norm) .OR. &
1148 6 : transition_norm < 0.0_dp .OR. exchange_norm < 0.0_dp) THEN
1149 0 : CPABORT("AUTO_RI target normalization is invalid")
1150 : END IF
1151 : ! B = B_ia / Tr(M^+ B_ia) + B_ni / Tr(M^+ B_ni), with zero-norm terms omitted.
1152 6 : IF (transition_norm > 0.0_dp) THEN
1153 45350 : B_PQ_transition(:, :) = B_PQ_transition/transition_norm
1154 6 : target_norm = target_norm + 1.0_dp
1155 : ELSE
1156 0 : B_PQ_transition(:, :) = 0.0_dp
1157 : END IF
1158 6 : IF (exchange_norm > 0.0_dp) THEN
1159 45350 : B_PQ_transition(:, :) = B_PQ_transition + B_PQ_exchange/exchange_norm
1160 6 : target_norm = target_norm + 1.0_dp
1161 : END IF
1162 45350 : IF (target_norm <= 0.0_dp .OR. .NOT. ALL(ieee_is_finite(B_PQ_transition))) THEN
1163 0 : CPABORT("AUTO_RI target normalization is invalid")
1164 : END IF
1165 :
1166 6 : CALL timestop(handle)
1167 :
1168 12 : END SUBROUTINE combine_transition_and_exchange_targets
1169 :
1170 : ! **************************************************************************************************
1171 : !> \brief Computes the atom-local matrix B_PQ^A in Eqs. (3)-(6).
1172 : !>
1173 : !> B_PQ^A = B_PQ^{A,C}/N_C^A + B_PQ^{A,X}/N_X^A.
1174 : !>
1175 : !> The molecular-orbital coefficient rows belonging to one atom are
1176 : !> collected on one MPI rank. No global AO-pair matrix is formed.
1177 : !> \param bs_env ...
1178 : !> \param iatom Atom whose AO rows are needed.
1179 : !> \param root_rank MPI rank performing the atomic contractions.
1180 : !> \param Int_3c atom-local three-center integrals (μν|P)
1181 : !> \param M_PQ atom-local metric of the large reference RI basis
1182 : !> \param B_PQ normalized atom-local matrix B_PQ^A
1183 : !> \param target_norm number of nonzero normalized contributions
1184 : ! **************************************************************************************************
1185 12 : SUBROUTINE compute_B_PQ_AA(bs_env, iatom, root_rank, Int_3c, M_PQ, &
1186 : B_PQ, target_norm)
1187 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1188 : INTEGER, INTENT(IN) :: iatom, root_rank
1189 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN) :: Int_3c(:, :, :), M_PQ(:, :)
1190 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: B_PQ(:, :)
1191 : REAL(KIND=dp), INTENT(OUT) :: target_norm
1192 :
1193 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_B_PQ_AA'
1194 :
1195 : INTEGER :: column, first, handle, i, j, nAO, nRI, &
1196 : nstates, rank, row, spin
1197 12 : INTEGER, POINTER :: cols(:), rows(:)
1198 : REAL(KIND=dp), ALLOCATABLE :: B_PQ_exchange(:, :), B_PQ_spin(:, :), &
1199 12 : C_mu_n(:, :), weights(:), X_Pi(:, :)
1200 :
1201 12 : CALL timeset(routineN, handle)
1202 12 : NULLIFY (rows, cols)
1203 12 : nAO = SIZE(Int_3c, 1)
1204 12 : nRI = SIZE(Int_3c, 3)
1205 12 : nstates = bs_env%n_mo_retained
1206 12 : CPASSERT(ALLOCATED(bs_env%eigenval_scf_Gamma))
1207 24 : first = 1 + SUM(bs_env%sizes_AO(:iatom - 1))
1208 48 : ALLOCATE (B_PQ(nRI, nRI), source=0.0_dp)
1209 96 : ALLOCATE (C_mu_n(nAO, nstates), weights(nstates), B_PQ_spin(nRI, nRI))
1210 12 : target_norm = 0.0_dp
1211 36 : ALLOCATE (B_PQ_exchange(nRI, nRI), source=0.0_dp)
1212 24 : DO spin = 1, bs_env%n_spin
1213 : CALL compute_orbital_weights(bs_env%eigenval_scf_Gamma(1:nstates, spin), &
1214 12 : bs_env%n_occ(spin), weights)
1215 12 : CALL cp_fm_get_info(bs_env%fm_mo_coeff_Gamma(spin), row_indices=rows, col_indices=cols)
1216 12 : C_mu_n(:, :) = 0.0_dp
1217 54 : DO i = 1, SIZE(rows)
1218 42 : row = rows(i) - first + 1
1219 42 : IF (row < 1 .OR. row > nAO) CYCLE
1220 124 : DO j = 1, SIZE(cols)
1221 98 : column = cols(j)
1222 98 : IF (column > nstates) CYCLE
1223 140 : C_mu_n(row, column) = bs_env%fm_mo_coeff_Gamma(spin)%local_data(i, j)
1224 : END DO
1225 : END DO
1226 12 : CALL bs_env%para_env%sum(C_mu_n, root_rank)
1227 12 : IF (bs_env%para_env%mepos /= root_rank) CYCLE
1228 6 : CALL compute_B_PQ_transition(Int_3c, C_mu_n, weights, bs_env%n_occ(spin), B_PQ_spin)
1229 45350 : B_PQ(:, :) = B_PQ + B_PQ_spin
1230 : CALL compute_B_PQ_exchange(Int_3c, C_mu_n, &
1231 : bs_env%eigenval_scf_Gamma(1:nstates, spin), &
1232 6 : bs_env%n_occ(spin), bs_env%auto_ri%occ_energy_window, B_PQ_spin)
1233 45368 : B_PQ_exchange(:, :) = B_PQ_exchange + B_PQ_spin
1234 : END DO
1235 12 : IF (bs_env%para_env%mepos /= root_rank) THEN
1236 6 : CALL timestop(handle)
1237 6 : RETURN
1238 : END IF
1239 6 : CALL compute_reference_M_PQ_inv_factor(M_PQ, X_Pi, rank)
1240 6 : CALL combine_transition_and_exchange_targets(X_Pi, B_PQ, B_PQ_exchange, target_norm)
1241 6 : CALL timestop(handle)
1242 :
1243 30 : END SUBROUTINE compute_B_PQ_AA
1244 :
1245 : ! **************************************************************************************************
1246 : !> \brief Normalizes each contraction coefficient vector U_Pp.
1247 : !>
1248 : !> U_Pp <- U_Pp/sqrt(Σ_P |U_Pp|^2).
1249 : !>
1250 : !> \param U_Pp contractions from the large RI basis P to optimized functions p
1251 : ! **************************************************************************************************
1252 18 : SUBROUTINE normalize_RI_coefficients(U_Pp)
1253 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: U_Pp
1254 :
1255 : CHARACTER(LEN=*), PARAMETER :: routineN = 'normalize_RI_coefficients'
1256 :
1257 : INTEGER :: handle, j
1258 : REAL(KIND=dp) :: column_norm
1259 :
1260 18 : CALL timeset(routineN, handle)
1261 :
1262 76 : DO j = 1, SIZE(U_Pp, 2)
1263 7742 : column_norm = NORM2(U_Pp(:, j))
1264 58 : IF (column_norm <= TINY(1.0_dp)) THEN
1265 0 : CPABORT("AUTO_RI produced a numerically zero localized column")
1266 : END IF
1267 7760 : U_Pp(:, j) = U_Pp(:, j)/column_norm
1268 : END DO
1269 :
1270 18 : CALL timestop(handle)
1271 :
1272 18 : END SUBROUTINE normalize_RI_coefficients
1273 :
1274 : ! **************************************************************************************************
1275 : !> \brief Creates a matrix in the existing molecular BLACS context.
1276 : !> \param qs_env ...
1277 : !> \param matrix ...
1278 : !> \param nr ...
1279 : !> \param nc ...
1280 : ! **************************************************************************************************
1281 24 : SUBROUTINE create_distributed_matrix(qs_env, matrix, nr, nc)
1282 : TYPE(qs_environment_type), POINTER :: qs_env
1283 : TYPE(cp_fm_type), INTENT(OUT) :: matrix
1284 : INTEGER, INTENT(IN) :: nr, nc
1285 :
1286 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_distributed_matrix'
1287 :
1288 : INTEGER :: handle
1289 : TYPE(cp_blacs_env_type), POINTER :: blacs
1290 : TYPE(cp_fm_struct_type), POINTER :: layout
1291 : TYPE(mp_para_env_type), POINTER :: para
1292 :
1293 24 : CALL timeset(routineN, handle)
1294 :
1295 24 : NULLIFY (blacs, para, layout)
1296 24 : CALL get_qs_env(qs_env, blacs_env=blacs, para_env=para)
1297 24 : CALL cp_fm_struct_create(layout, context=blacs, nrow_global=nr, ncol_global=nc, para_env=para)
1298 24 : CALL cp_fm_create(matrix, layout, set_zero=.TRUE.)
1299 24 : CALL cp_fm_struct_release(layout)
1300 :
1301 24 : CALL timestop(handle)
1302 :
1303 24 : END SUBROUTINE create_distributed_matrix
1304 :
1305 : ! **************************************************************************************************
1306 : !> \brief Computes and stores the two-center contractions U^AB in Eqs. (8)-(13).
1307 : !>
1308 : !> The selected AA space is projected out according to Eq. (8). The symmetric
1309 : !> pair density in Eq. (9) defines B^{AB,C} and B^{AB,X} in Eqs. (10) and (11).
1310 : !> Diagonalizing their sum gives U^AB in Eq. (12), which expands the contracted
1311 : !> pair functions in Eq. (13).
1312 : !> \param bs_env ...
1313 : !> \param U_Pp_A temporary atom-local contractions used to project the AA space
1314 : ! **************************************************************************************************
1315 4 : SUBROUTINE compute_AB_optimized_RI_basis(bs_env, U_Pp_A)
1316 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1317 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: U_Pp_A
1318 :
1319 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_AB_optimized_RI_basis'
1320 :
1321 : INTEGER :: AB_pair_count, center_index, current_atom, expected_RI_size, first_AB_pair, &
1322 : handle, handle_pair_batches, iatom, jatom, k, last_AB_pair, nAO_kind, nRI_A, nRI_AB_ref, &
1323 : nRI_B, RI_size_opt_AA, RI_size_opt_AB, RI_size_optimized, RI_size_ref, RI_size_requested, &
1324 : root_rank, row, size_opt_RI_AB, size_U_Pp_AB, size_U_Pp_AB_storage, spin
1325 4 : INTEGER, ALLOCATABLE :: AB_atom_A(:), AB_atom_B(:), AB_size_opt_RI_to_A(:), block_next(:), &
1326 4 : sizes_opt_RI_AA(:), sizes_opt_RI_AB(:), sizes_ref_RI(:), U_Pp_AB_offset(:)
1327 4 : REAL(KIND=dp), ALLOCATABLE :: B_PQ(:, :), C_mu_n_AB(:, :, :), C_mu_n_atom(:, :), &
1328 4 : eigenval(:), eigenval_AB(:, :), M_PQ(:, :), U_Pp_AA(:, :), U_Pp_AB(:, :), &
1329 4 : U_Pp_AB_packed(:)
1330 52 : TYPE(gw_3c_ctx_type) :: ctx
1331 4 : TYPE(gw_3c_ws_type) :: ws
1332 : TYPE(mp_para_env_type), POINTER :: para
1333 :
1334 4 : CALL timeset(routineN, handle)
1335 :
1336 4 : para => bs_env%para_env
1337 16 : ALLOCATE (sizes_ref_RI(bs_env%n_atom), sizes_opt_RI_AA(bs_env%n_atom))
1338 16 : DO iatom = 1, bs_env%n_atom
1339 16 : sizes_ref_RI(iatom) = get_ref_RI_size(bs_env, iatom)
1340 : END DO
1341 16 : sizes_opt_RI_AA(:) = bs_env%auto_ri%sizes_opt_RI
1342 16 : RI_size_ref = SUM(sizes_ref_RI)
1343 16 : RI_size_opt_AA = SUM(sizes_opt_RI_AA)
1344 :
1345 4 : CALL count_AB_neighbor_pairs(bs_env, AB_pair_count)
1346 : ALLOCATE (AB_atom_A(AB_pair_count), AB_atom_B(AB_pair_count), &
1347 : AB_size_opt_RI_to_A(AB_pair_count), sizes_opt_RI_AB(AB_pair_count), &
1348 18 : U_Pp_AB_offset(AB_pair_count))
1349 4 : AB_size_opt_RI_to_A = 0
1350 4 : U_Pp_AB_offset = 0
1351 26 : ALLOCATE (eigenval_AB(2*MAXVAL(sizes_ref_RI), AB_pair_count), source=0.0_dp)
1352 :
1353 : CALL gw_3c_ctx_create(ctx, bs_env, bs_env%ri_metric, &
1354 : basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, &
1355 4 : basis_i=bs_env%basis_set_RI)
1356 4 : CALL gw_3c_ws_create(ws, ctx)
1357 4 : k = 0
1358 12 : DO iatom = 1, bs_env%n_atom - 1
1359 24 : DO jatom = iatom + 1, bs_env%n_atom
1360 48 : IF (NORM2(bs_env%ri_rs%particle_set(iatom)%r - bs_env%ri_rs%particle_set(jatom)%r) > &
1361 : bs_env%auto_ri%neighbor_radius) CYCLE
1362 6 : k = k + 1
1363 6 : AB_atom_A(k) = iatom
1364 20 : AB_atom_B(k) = jatom
1365 : END DO
1366 : END DO
1367 : size_U_Pp_AB_storage = 0
1368 10 : DO k = 1, AB_pair_count
1369 6 : root_rank = MODULO(k - 1, para%num_pe)
1370 6 : IF (para%mepos /= root_rank) CYCLE
1371 3 : nRI_AB_ref = sizes_ref_RI(AB_atom_A(k)) + sizes_ref_RI(AB_atom_B(k))
1372 3 : U_Pp_AB_offset(k) = size_U_Pp_AB_storage + 1
1373 10 : size_U_Pp_AB_storage = size_U_Pp_AB_storage + nRI_AB_ref*nRI_AB_ref
1374 : END DO
1375 10 : ALLOCATE (U_Pp_AB_packed(size_U_Pp_AB_storage))
1376 4 : CALL timeset(routineN//'_pair_batches', handle_pair_batches)
1377 8 : DO first_AB_pair = 1, AB_pair_count, para%num_pe
1378 4 : last_AB_pair = MIN(AB_pair_count, first_AB_pair + para%num_pe - 1)
1379 : ! Gather only one pair's MO rows per rank, then solve all batch pairs concurrently.
1380 10 : DO k = first_AB_pair, last_AB_pair
1381 6 : root_rank = k - first_AB_pair
1382 6 : iatom = AB_atom_A(k)
1383 6 : jatom = AB_atom_B(k)
1384 6 : IF (para%mepos == root_rank) THEN
1385 0 : ALLOCATE (C_mu_n_AB( &
1386 : bs_env%sizes_AO(iatom) + bs_env%sizes_AO(jatom), &
1387 15 : bs_env%n_mo_retained, bs_env%n_spin))
1388 : END IF
1389 16 : DO spin = 1, bs_env%n_spin
1390 : row = 0
1391 24 : DO center_index = 1, 2
1392 12 : current_atom = iatom
1393 12 : IF (center_index == 2) current_atom = jatom
1394 12 : nAO_kind = bs_env%sizes_AO(current_atom)
1395 48 : ALLOCATE (C_mu_n_atom(nAO_kind, bs_env%n_mo_retained))
1396 : CALL cp_fm_get_submatrix(bs_env%fm_mo_coeff_Gamma(spin), C_mu_n_atom, &
1397 24 : start_row=1 + SUM(bs_env%sizes_AO(:current_atom - 1)))
1398 12 : IF (para%mepos == root_rank) THEN
1399 146 : C_mu_n_AB(row + 1:row + nAO_kind, :, spin) = C_mu_n_atom
1400 : END IF
1401 12 : DEALLOCATE (C_mu_n_atom)
1402 18 : row = row + nAO_kind
1403 : END DO
1404 : END DO
1405 : END DO
1406 4 : k = first_AB_pair + para%mepos
1407 8 : IF (k <= last_AB_pair) THEN
1408 3 : iatom = AB_atom_A(k)
1409 3 : jatom = AB_atom_B(k)
1410 3 : nRI_A = sizes_ref_RI(iatom)
1411 3 : nRI_B = sizes_ref_RI(jatom)
1412 3 : nRI_AB_ref = nRI_A + nRI_B
1413 : ! Eqs. (9)-(11): form the symmetric pair density and B^{AB,C}+B^{AB,X}.
1414 3 : CALL compute_B_PQ_AB(bs_env, ctx, ws, iatom, jatom, C_mu_n_AB, M_PQ, B_PQ)
1415 0 : ALLOCATE (U_Pp_AA(nRI_AB_ref, sizes_opt_RI_AA(iatom) + sizes_opt_RI_AA(jatom)), &
1416 12 : source=0.0_dp)
1417 : U_Pp_AA(1:nRI_A, 1:sizes_opt_RI_AA(iatom)) = &
1418 1865 : U_Pp_A(1:nRI_A, 1:sizes_opt_RI_AA(iatom), iatom)
1419 : U_Pp_AA(nRI_A + 1:, sizes_opt_RI_AA(iatom) + 1:) = &
1420 171 : U_Pp_A(1:nRI_B, 1:sizes_opt_RI_AA(jatom), jatom)
1421 : ! Eqs. (8) and (12): project out the AA space and diagonalize the remainder.
1422 3 : CALL solve_AB_eigenproblem(M_PQ, B_PQ, U_Pp_AA, U_Pp_AB, eigenval)
1423 17 : eigenval_AB(1:SIZE(eigenval), k) = eigenval
1424 3 : DEALLOCATE (eigenval)
1425 : U_Pp_AB_packed(U_Pp_AB_offset(k):U_Pp_AB_offset(k) + SIZE(U_Pp_AB) - 1) = &
1426 18 : RESHAPE(U_Pp_AB, [SIZE(U_Pp_AB)])
1427 3 : DEALLOCATE (U_Pp_AB)
1428 3 : DEALLOCATE (U_Pp_AA, B_PQ, M_PQ, C_mu_n_AB)
1429 : END IF
1430 : END DO
1431 4 : CALL timestop(handle_pair_batches)
1432 4 : k = AB_pair_count
1433 4 : CALL gw_3c_ws_release(ws)
1434 4 : CALL gw_3c_ctx_release(ctx)
1435 : RI_size_requested = NINT(MIN(bs_env%auto_ri%ri_ao_ratio, &
1436 : REAL(RI_size_ref, dp)/REAL(SUM(bs_env%sizes_AO), dp))* &
1437 16 : REAL(SUM(bs_env%sizes_AO), dp))
1438 4 : RI_size_requested = MAX(RI_size_requested, bs_env%n_atom)
1439 :
1440 4 : CALL para%sum(eigenval_AB)
1441 : CALL select_RI_sizes(eigenval_AB, RI_size_requested - SUM(sizes_opt_RI_AA), 0, &
1442 16 : sizes_opt_RI_AB)
1443 10 : RI_size_opt_AB = SUM(sizes_opt_RI_AB)
1444 16 : bs_env%auto_ri%sizes_opt_RI(:) = sizes_opt_RI_AA
1445 10 : DO k = 1, AB_pair_count
1446 6 : iatom = AB_atom_A(k)
1447 6 : jatom = AB_atom_B(k)
1448 6 : size_opt_RI_AB = sizes_opt_RI_AB(k)
1449 : AB_size_opt_RI_to_A(k) = balance_AB_columns(bs_env%auto_ri%sizes_opt_RI(iatom), &
1450 : bs_env%auto_ri%sizes_opt_RI(jatom), &
1451 6 : size_opt_RI_AB)
1452 : bs_env%auto_ri%sizes_opt_RI(iatom) = &
1453 6 : bs_env%auto_ri%sizes_opt_RI(iatom) + AB_size_opt_RI_to_A(k)
1454 : bs_env%auto_ri%sizes_opt_RI(jatom) = &
1455 10 : bs_env%auto_ri%sizes_opt_RI(jatom) + size_opt_RI_AB - AB_size_opt_RI_to_A(k)
1456 : END DO
1457 :
1458 16 : RI_size_optimized = SUM(bs_env%auto_ri%sizes_opt_RI)
1459 4 : expected_RI_size = RI_size_opt_AA + RI_size_opt_AB
1460 4 : CPASSERT(RI_size_optimized == expected_RI_size)
1461 16 : CPASSERT(ALL(bs_env%auto_ri%sizes_opt_RI > 0))
1462 4 : IF (RI_size_optimized > RI_size_ref) THEN
1463 0 : CPABORT("AUTO_RI request exceeds independent reference RI directions")
1464 : END IF
1465 4 : size_U_Pp_AB = 0
1466 16 : DO iatom = 1, bs_env%n_atom
1467 16 : size_U_Pp_AB = size_U_Pp_AB + sizes_ref_RI(iatom)*sizes_opt_RI_AA(iatom)
1468 : END DO
1469 10 : DO k = 1, AB_pair_count
1470 : size_U_Pp_AB = size_U_Pp_AB + &
1471 : (sizes_ref_RI(AB_atom_A(k)) + &
1472 10 : sizes_ref_RI(AB_atom_B(k)))*sizes_opt_RI_AB(k)
1473 : END DO
1474 30 : ALLOCATE (bs_env%auto_ri%AB_atom_A(bs_env%n_atom + AB_pair_count), source=0)
1475 30 : ALLOCATE (bs_env%auto_ri%AB_atom_B(bs_env%n_atom + AB_pair_count), source=0)
1476 30 : ALLOCATE (bs_env%auto_ri%AB_first_p_A(bs_env%n_atom + AB_pair_count), source=0)
1477 30 : ALLOCATE (bs_env%auto_ri%AB_first_p_B(bs_env%n_atom + AB_pair_count), source=0)
1478 30 : ALLOCATE (bs_env%auto_ri%AB_size_opt_RI_to_A(bs_env%n_atom + AB_pair_count), source=0)
1479 30 : ALLOCATE (bs_env%auto_ri%AB_size_ref_RI(bs_env%n_atom + AB_pair_count), source=0)
1480 30 : ALLOCATE (bs_env%auto_ri%AB_size_opt_RI(bs_env%n_atom + AB_pair_count), source=0)
1481 30 : ALLOCATE (bs_env%auto_ri%U_Pp_AB_offset(bs_env%n_atom + AB_pair_count), source=0)
1482 7696 : ALLOCATE (bs_env%auto_ri%U_Pp_AB(size_U_Pp_AB), source=0.0_dp)
1483 4 : bs_env%auto_ri%AB_block_count = 0
1484 :
1485 12 : ALLOCATE (block_next(bs_env%n_atom))
1486 4 : block_next(1) = 1
1487 12 : DO iatom = 2, bs_env%n_atom
1488 12 : block_next(iatom) = block_next(iatom - 1) + bs_env%auto_ri%sizes_opt_RI(iatom - 1)
1489 : END DO
1490 16 : DO iatom = 1, bs_env%n_atom
1491 12 : IF (sizes_opt_RI_AA(iatom) > 0) bs_env%auto_ri%AB_block_count = &
1492 12 : bs_env%auto_ri%AB_block_count + 1
1493 : CALL store_AA_columns(bs_env, iatom, &
1494 : U_Pp_A(1:sizes_ref_RI(iatom), 1:sizes_opt_RI_AA(iatom), iatom), &
1495 16 : block_next)
1496 : END DO
1497 :
1498 10 : DO k = 1, AB_pair_count
1499 6 : iatom = AB_atom_A(k)
1500 6 : jatom = AB_atom_B(k)
1501 6 : size_opt_RI_AB = sizes_opt_RI_AB(k)
1502 6 : IF (size_opt_RI_AB == 0) CYCLE
1503 6 : root_rank = MODULO(k - 1, para%num_pe)
1504 6 : nRI_AB_ref = sizes_ref_RI(iatom) + sizes_ref_RI(jatom)
1505 24 : ALLOCATE (U_Pp_AB(nRI_AB_ref, size_opt_RI_AB))
1506 6 : IF (para%mepos == root_rank) THEN
1507 : U_Pp_AB(:, :) = RESHAPE( &
1508 : U_Pp_AB_packed(U_Pp_AB_offset(k): &
1509 : U_Pp_AB_offset(k) + &
1510 : nRI_AB_ref*size_opt_RI_AB - 1), &
1511 9 : [nRI_AB_ref, size_opt_RI_AB])
1512 : END IF
1513 6 : CALL para%bcast(U_Pp_AB, root_rank)
1514 6 : bs_env%auto_ri%AB_block_count = bs_env%auto_ri%AB_block_count + 1
1515 : ! Eq. (13): store the selected neighbor-pair contractions in the global transform.
1516 6 : CALL store_AB_columns(bs_env, U_Pp_AB, iatom, jatom, AB_size_opt_RI_to_A(k), block_next)
1517 10 : DEALLOCATE (U_Pp_AB)
1518 : END DO
1519 :
1520 16 : DO iatom = 1, bs_env%n_atom
1521 36 : expected_RI_size = 1 + SUM(bs_env%auto_ri%sizes_opt_RI(:iatom))
1522 16 : CPASSERT(block_next(iatom) == expected_RI_size)
1523 : END DO
1524 0 : DEALLOCATE (block_next, AB_atom_A, AB_atom_B, AB_size_opt_RI_to_A, &
1525 4 : U_Pp_AB_offset, U_Pp_AB_packed)
1526 :
1527 4 : CALL timestop(handle)
1528 :
1529 16 : END SUBROUTINE compute_AB_optimized_RI_basis
1530 :
1531 : ! **************************************************************************************************
1532 : !> \brief Divides AB columns between atoms A and B for load balance.
1533 : !> \param load_a current number of columns owned by A
1534 : !> \param load_b current number of columns owned by B
1535 : !> \param RI_size_opt_AB number of new pair RI functions
1536 : !> \return ...
1537 : ! **************************************************************************************************
1538 6 : PURE INTEGER FUNCTION balance_AB_columns(load_a, load_b, RI_size_opt_AB) RESULT(RI_size_to_a)
1539 : INTEGER, INTENT(IN) :: load_a, load_b, RI_size_opt_AB
1540 :
1541 6 : RI_size_to_a = (load_b + RI_size_opt_AB - load_a + 1)/2
1542 6 : RI_size_to_a = MAX(0, MIN(RI_size_opt_AB, RI_size_to_a))
1543 :
1544 6 : END FUNCTION balance_AB_columns
1545 :
1546 : ! **************************************************************************************************
1547 : !> \brief Stores the AA contraction columns for one atom.
1548 : !> \param bs_env ...
1549 : !> \param atom ...
1550 : !> \param U_Pp_A atom-local contraction coefficients
1551 : !> \param block_next ...
1552 : ! **************************************************************************************************
1553 12 : SUBROUTINE store_AA_columns(bs_env, atom, U_Pp_A, block_next)
1554 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1555 : INTEGER, INTENT(IN) :: atom
1556 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: U_Pp_A
1557 : INTEGER, DIMENSION(:), INTENT(INOUT) :: block_next
1558 :
1559 : CHARACTER(LEN=*), PARAMETER :: routineN = 'store_AA_columns'
1560 :
1561 : INTEGER :: AB_block, first, handle, RI_size_opt_AA
1562 : REAL(KIND=dp), ALLOCATABLE :: U_Pp_A_normalized(:, :)
1563 :
1564 12 : CALL timeset(routineN, handle)
1565 :
1566 12 : first = block_next(atom)
1567 12 : AB_block = bs_env%auto_ri%AB_block_count
1568 12 : RI_size_opt_AA = SIZE(U_Pp_A, 2)
1569 48 : ALLOCATE (U_Pp_A_normalized(SIZE(U_Pp_A, 1), RI_size_opt_AA))
1570 12 : IF (RI_size_opt_AA > 0) THEN
1571 4588 : U_Pp_A_normalized(:, :) = U_Pp_A
1572 12 : CALL normalize_RI_coefficients(U_Pp_A_normalized)
1573 12 : bs_env%auto_ri%AB_atom_A(AB_block) = atom
1574 12 : bs_env%auto_ri%AB_atom_B(AB_block) = atom
1575 : bs_env%auto_ri%AB_first_p_A(AB_block) = &
1576 24 : first - SUM(bs_env%auto_ri%sizes_opt_RI(:atom - 1))
1577 12 : bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block) = RI_size_opt_AA
1578 12 : CALL store_U_Pp_AB(bs_env%auto_ri, AB_block, U_Pp_A_normalized)
1579 : END IF
1580 12 : DEALLOCATE (U_Pp_A_normalized)
1581 12 : block_next(atom) = first + RI_size_opt_AA
1582 :
1583 12 : CALL timestop(handle)
1584 :
1585 12 : END SUBROUTINE store_AA_columns
1586 :
1587 : ! **************************************************************************************************
1588 : !> \brief Stores an AB contraction in the blocks assigned to atoms A and B.
1589 : !> \param bs_env ...
1590 : !> \param U_Pp_AB neighbor-pair contractions
1591 : !> \param atom_A ...
1592 : !> \param atom_B ...
1593 : !> \param RI_size_to_a number of optimized pair RI functions assigned to atom A
1594 : !> \param block_next ...
1595 : ! **************************************************************************************************
1596 6 : SUBROUTINE store_AB_columns(bs_env, U_Pp_AB, atom_A, atom_B, RI_size_to_a, block_next)
1597 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1598 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1599 : INTENT(INOUT) :: U_Pp_AB
1600 : INTEGER, INTENT(IN) :: atom_A, atom_B, RI_size_to_a
1601 : INTEGER, DIMENSION(:), INTENT(INOUT) :: block_next
1602 :
1603 : CHARACTER(LEN=*), PARAMETER :: routineN = 'store_AB_columns'
1604 :
1605 : INTEGER :: AB_block, first_a, first_b, handle, &
1606 : RI_size_opt_AB
1607 :
1608 6 : CALL timeset(routineN, handle)
1609 :
1610 6 : RI_size_opt_AB = SIZE(U_Pp_AB, 2)
1611 6 : AB_block = bs_env%auto_ri%AB_block_count
1612 6 : first_a = block_next(atom_A)
1613 6 : first_b = block_next(atom_B)
1614 6 : CALL normalize_RI_coefficients(U_Pp_AB)
1615 :
1616 6 : bs_env%auto_ri%AB_atom_A(AB_block) = atom_A
1617 6 : bs_env%auto_ri%AB_atom_B(AB_block) = atom_B
1618 8 : bs_env%auto_ri%AB_first_p_A(AB_block) = first_a - SUM(bs_env%auto_ri%sizes_opt_RI(:atom_A - 1))
1619 16 : bs_env%auto_ri%AB_first_p_B(AB_block) = first_b - SUM(bs_env%auto_ri%sizes_opt_RI(:atom_B - 1))
1620 6 : bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block) = RI_size_to_a
1621 6 : CALL store_U_Pp_AB(bs_env%auto_ri, AB_block, U_Pp_AB)
1622 6 : block_next(atom_A) = first_a + RI_size_to_a
1623 6 : block_next(atom_B) = first_b + RI_size_opt_AB - RI_size_to_a
1624 :
1625 6 : CALL timestop(handle)
1626 :
1627 6 : END SUBROUTINE store_AB_columns
1628 :
1629 : ! **************************************************************************************************
1630 : !> \brief Packs one AB U_Pp matrix into the persistent AUTO_RI data; B=A denotes an AA block.
1631 : !> \param auto_ri ...
1632 : !> \param AB_block ...
1633 : !> \param U_Pp ...
1634 : ! **************************************************************************************************
1635 18 : SUBROUTINE store_U_Pp_AB(auto_ri, AB_block, U_Pp)
1636 : TYPE(auto_ri_type), INTENT(INOUT) :: auto_ri
1637 : INTEGER, INTENT(IN) :: AB_block
1638 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: U_Pp
1639 :
1640 : INTEGER :: first, last
1641 :
1642 18 : first = 1
1643 18 : IF (AB_block > 1) THEN
1644 : first = auto_ri%U_Pp_AB_offset(AB_block - 1) + &
1645 : auto_ri%AB_size_ref_RI(AB_block - 1)* &
1646 14 : auto_ri%AB_size_opt_RI(AB_block - 1)
1647 : END IF
1648 54 : last = first + SIZE(U_Pp) - 1
1649 18 : CPASSERT(last <= SIZE(auto_ri%U_Pp_AB))
1650 18 : auto_ri%U_Pp_AB_offset(AB_block) = first
1651 18 : auto_ri%AB_size_ref_RI(AB_block) = SIZE(U_Pp, 1)
1652 18 : auto_ri%AB_size_opt_RI(AB_block) = SIZE(U_Pp, 2)
1653 7756 : auto_ri%U_Pp_AB(first:last) = RESHAPE(U_Pp, [SIZE(U_Pp)])
1654 :
1655 18 : END SUBROUTINE store_U_Pp_AB
1656 :
1657 : ! **************************************************************************************************
1658 : !> \brief Computes the two-center matrix B_PQ^AB in Eqs. (9)-(12).
1659 : !>
1660 : !> Its numerator uses the symmetric off-diagonal density
1661 : !> ρ̅_mn^AB=ρ_mn^AB+ρ_mn^BA after projection from the retained AA
1662 : !> space. Its normalization uses the complete two-atom density.
1663 : !> \param bs_env ...
1664 : !> \param ctx ...
1665 : !> \param ws ...
1666 : !> \param iatom ...
1667 : !> \param jatom ...
1668 : !> \param C_mu_n_AB molecular-orbital coefficients on the AO rows of atoms A and B
1669 : !> \param M_PQ combined metric of the large reference RI bases on A and B
1670 : !> \param B_PQ normalized two-center matrix B_PQ^AB
1671 : ! **************************************************************************************************
1672 3 : SUBROUTINE compute_B_PQ_AB(bs_env, ctx, ws, iatom, jatom, C_mu_n_AB, M_PQ, B_PQ)
1673 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1674 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
1675 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
1676 : INTEGER, INTENT(IN) :: iatom, jatom
1677 : REAL(KIND=dp), INTENT(IN) :: C_mu_n_AB(:, :, :)
1678 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: M_PQ(:, :), B_PQ(:, :)
1679 :
1680 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_B_PQ_AB'
1681 :
1682 : INTEGER :: atoms(2), handle, homo, i, j, k, kind_a, kind_b, nAO_AB, nRI_AB_ref, &
1683 : offset_AO(2), offset_ref_RI(2), rank, sizes_AO(2), sizes_ref_RI(2), spin
1684 : LOGICAL :: screened
1685 : REAL(KIND=dp) :: norm_exchange, norm_transition
1686 3 : REAL(KIND=dp), ALLOCATABLE :: B_Pi(:, :), B_PQ_exchange_AB(:, :), B_PQ_exchange_full(:, :), &
1687 3 : B_PQ_spin(:, :), B_PQ_transition_AB(:, :), B_PQ_transition_full(:, :), &
1688 3 : Int_3c_AB(:, :, :), Int_3c_all(:, :, :), M_PQ_block(:, :), weights(:), X_Pi(:, :)
1689 :
1690 3 : CALL timeset(routineN, handle)
1691 :
1692 3 : atoms(1) = iatom
1693 3 : atoms(2) = jatom
1694 9 : sizes_AO = bs_env%sizes_AO(atoms)
1695 9 : DO i = 1, 2
1696 9 : sizes_ref_RI(i) = get_ref_RI_size(bs_env, atoms(i))
1697 : END DO
1698 3 : offset_AO(1) = 0
1699 3 : offset_AO(2) = sizes_AO(1)
1700 3 : offset_ref_RI(1) = 0
1701 3 : offset_ref_RI(2) = sizes_ref_RI(1)
1702 9 : nAO_AB = SUM(sizes_AO)
1703 9 : nRI_AB_ref = SUM(sizes_ref_RI)
1704 0 : ALLOCATE (B_PQ(nRI_AB_ref, nRI_AB_ref), &
1705 18 : M_PQ(nRI_AB_ref, nRI_AB_ref), source=0.0_dp)
1706 15 : ALLOCATE (Int_3c_all(nAO_AB, nAO_AB, nRI_AB_ref), source=0.0_dp)
1707 9 : DO i = 1, 2
1708 21 : DO j = 1, 2
1709 : CALL compute_M_PQ_AB( &
1710 12 : bs_env, atoms(i), atoms(j), M_PQ_block, bs_env%ri_metric)
1711 : M_PQ(offset_ref_RI(i) + 1:offset_ref_RI(i) + sizes_ref_RI(i), &
1712 80042 : offset_ref_RI(j) + 1:offset_ref_RI(j) + sizes_ref_RI(j)) = M_PQ_block
1713 12 : DEALLOCATE (M_PQ_block)
1714 42 : DO k = 1, 2
1715 : CALL build_3c_integral_block_ctx( &
1716 : Int_3c_all, ctx, ws, atom_j=atoms(j), atom_k=atoms(k), &
1717 : atom_i=atoms(i), j_offset=offset_AO(j), k_offset=offset_AO(k), &
1718 36 : i_offset=offset_ref_RI(i), screened=screened)
1719 : END DO
1720 : END DO
1721 : END DO
1722 3 : kind_a = bs_env%ri_rs%particle_set(iatom)%atomic_kind%kind_number
1723 3 : kind_b = bs_env%ri_rs%particle_set(jatom)%atomic_kind%kind_number
1724 : CALL add_RI_metric_regularization(bs_env%basis_set_RI(kind_a)%gto_basis_set, &
1725 3 : bs_env%regularization_RI, M_PQ)
1726 : CALL add_RI_metric_regularization(bs_env%basis_set_RI(kind_b)%gto_basis_set, &
1727 : bs_env%regularization_RI, M_PQ, &
1728 3 : offset=sizes_ref_RI(1))
1729 16523 : ALLOCATE (Int_3c_AB, source=Int_3c_all)
1730 11679 : Int_3c_AB(1:sizes_AO(1), 1:sizes_AO(1), :) = 0.0_dp
1731 1431 : Int_3c_AB(sizes_AO(1) + 1:, sizes_AO(1) + 1:, :) = 0.0_dp
1732 15 : ALLOCATE (weights(bs_env%n_mo_retained), B_PQ_spin(nRI_AB_ref, nRI_AB_ref))
1733 : ALLOCATE (B_PQ_transition_AB(nRI_AB_ref, nRI_AB_ref), &
1734 : B_PQ_exchange_AB(nRI_AB_ref, nRI_AB_ref), &
1735 : B_PQ_transition_full(nRI_AB_ref, nRI_AB_ref), &
1736 : B_PQ_exchange_full(nRI_AB_ref, nRI_AB_ref), &
1737 27 : source=0.0_dp)
1738 6 : DO spin = 1, bs_env%n_spin
1739 3 : homo = bs_env%n_occ(spin)
1740 : CALL compute_orbital_weights( &
1741 : eigenval_MO=bs_env%eigenval_scf_Gamma(1:bs_env%n_mo_retained, spin), &
1742 3 : homo=homo, weights=weights)
1743 122 : CALL compute_B_PQ_transition(Int_3c_all, C_mu_n_AB(:, :, spin), weights, homo, B_PQ_spin)
1744 79557 : B_PQ_transition_full(:, :) = B_PQ_transition_full + B_PQ_spin
1745 122 : CALL compute_B_PQ_transition(Int_3c_AB, C_mu_n_AB(:, :, spin), weights, homo, B_PQ_spin)
1746 79557 : B_PQ_transition_AB(:, :) = B_PQ_transition_AB + B_PQ_spin
1747 : CALL compute_B_PQ_exchange(Int_3c_all, C_mu_n_AB(:, :, spin), &
1748 : bs_env%eigenval_scf_Gamma(1:bs_env%n_mo_retained, spin), &
1749 122 : homo, bs_env%auto_ri%occ_energy_window, B_PQ_spin)
1750 79557 : B_PQ_exchange_full(:, :) = B_PQ_exchange_full + B_PQ_spin
1751 : CALL compute_B_PQ_exchange(Int_3c_AB, C_mu_n_AB(:, :, spin), &
1752 : bs_env%eigenval_scf_Gamma(1:bs_env%n_mo_retained, spin), &
1753 122 : homo, bs_env%auto_ri%occ_energy_window, B_PQ_spin)
1754 79560 : B_PQ_exchange_AB(:, :) = B_PQ_exchange_AB + B_PQ_spin
1755 : END DO
1756 3 : CALL compute_reference_M_PQ_inv_factor(M_PQ, X_Pi, rank)
1757 3 : DEALLOCATE (B_PQ_spin)
1758 12 : ALLOCATE (B_Pi(nRI_AB_ref, rank))
1759 3 : B_Pi(:, :) = MATMUL(B_PQ_transition_full, X_Pi)
1760 79557 : norm_transition = SUM(X_Pi*B_Pi)
1761 3 : B_Pi(:, :) = MATMUL(B_PQ_exchange_full, X_Pi)
1762 79557 : norm_exchange = SUM(X_Pi*B_Pi)
1763 : ! Eqs. (10)-(12): B_PQ^AB = (B_PQ^{AB,C}/N_C^{AB}+B_PQ^{AB,X}/N_X^{AB})/2.
1764 3 : IF (norm_transition > 0.0_dp) THEN
1765 79557 : B_PQ(:, :) = B_PQ + 0.5_dp*B_PQ_transition_AB/norm_transition
1766 : END IF
1767 3 : IF (norm_exchange > 0.0_dp) THEN
1768 79557 : B_PQ(:, :) = B_PQ + 0.5_dp*B_PQ_exchange_AB/norm_exchange
1769 : END IF
1770 :
1771 3 : CALL timestop(handle)
1772 :
1773 6 : END SUBROUTINE compute_B_PQ_AB
1774 :
1775 : ! **************************************************************************************************
1776 : !> \brief Builds lists of AB blocks whose columns belong to each atom; B=A denotes an AA block.
1777 : !> \param bs_env ...
1778 : !> \param lists ...
1779 : ! **************************************************************************************************
1780 4 : SUBROUTINE build_AB_block_lists(bs_env, lists)
1781 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1782 : TYPE(AB_block_list_type), ALLOCATABLE, INTENT(OUT) :: lists(:)
1783 :
1784 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_AB_block_lists'
1785 :
1786 : INTEGER :: AB_block, atom_A, atom_B, handle, &
1787 : position, RI_size_to_a, size_opt_RI_AB
1788 4 : INTEGER, ALLOCATABLE :: n_AB_blocks_by_atom(:)
1789 :
1790 4 : CALL timeset(routineN, handle)
1791 :
1792 24 : ALLOCATE (lists(bs_env%n_atom))
1793 12 : ALLOCATE (n_AB_blocks_by_atom(bs_env%n_atom), source=0)
1794 22 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
1795 18 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
1796 18 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
1797 18 : size_opt_RI_AB = bs_env%auto_ri%AB_size_opt_RI(AB_block)
1798 18 : RI_size_to_a = bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
1799 18 : IF (RI_size_to_a > 0) n_AB_blocks_by_atom(atom_A) = n_AB_blocks_by_atom(atom_A) + 1
1800 22 : IF (atom_B /= atom_A .AND. RI_size_to_a < size_opt_RI_AB) THEN
1801 4 : n_AB_blocks_by_atom(atom_B) = n_AB_blocks_by_atom(atom_B) + 1
1802 : END IF
1803 : END DO
1804 :
1805 16 : DO atom_A = 1, bs_env%n_atom
1806 36 : ALLOCATE (lists(atom_A)%block(n_AB_blocks_by_atom(atom_A)))
1807 24 : ALLOCATE (lists(atom_A)%first_p_AB(n_AB_blocks_by_atom(atom_A)))
1808 24 : ALLOCATE (lists(atom_A)%first_p_atom(n_AB_blocks_by_atom(atom_A)))
1809 28 : ALLOCATE (lists(atom_A)%size_opt_RI(n_AB_blocks_by_atom(atom_A)))
1810 : END DO
1811 4 : n_AB_blocks_by_atom = 0
1812 22 : DO AB_block = 1, bs_env%auto_ri%AB_block_count
1813 18 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
1814 18 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
1815 18 : size_opt_RI_AB = bs_env%auto_ri%AB_size_opt_RI(AB_block)
1816 18 : RI_size_to_a = bs_env%auto_ri%AB_size_opt_RI_to_A(AB_block)
1817 18 : IF (RI_size_to_a > 0) THEN
1818 14 : position = n_AB_blocks_by_atom(atom_A) + 1
1819 14 : n_AB_blocks_by_atom(atom_A) = position
1820 14 : lists(atom_A)%block(position) = AB_block
1821 14 : lists(atom_A)%first_p_AB(position) = 1
1822 14 : lists(atom_A)%first_p_atom(position) = bs_env%auto_ri%AB_first_p_A(AB_block)
1823 14 : lists(atom_A)%size_opt_RI(position) = RI_size_to_a
1824 : END IF
1825 22 : IF (atom_B /= atom_A .AND. RI_size_to_a < size_opt_RI_AB) THEN
1826 4 : position = n_AB_blocks_by_atom(atom_B) + 1
1827 4 : n_AB_blocks_by_atom(atom_B) = position
1828 4 : lists(atom_B)%block(position) = AB_block
1829 4 : lists(atom_B)%first_p_AB(position) = RI_size_to_a + 1
1830 : lists(atom_B)%first_p_atom(position) = &
1831 4 : bs_env%auto_ri%AB_first_p_B(AB_block)
1832 4 : lists(atom_B)%size_opt_RI(position) = size_opt_RI_AB - RI_size_to_a
1833 : END IF
1834 : END DO
1835 16 : DO atom_A = 1, bs_env%n_atom
1836 34 : IF (SUM(lists(atom_A)%size_opt_RI) /= bs_env%auto_ri%sizes_opt_RI(atom_A)) THEN
1837 0 : CPABORT("AUTO_RI AB-block columns do not fill the optimized atom block")
1838 : END IF
1839 : END DO
1840 4 : DEALLOCATE (n_AB_blocks_by_atom)
1841 :
1842 4 : CALL timestop(handle)
1843 :
1844 4 : END SUBROUTINE build_AB_block_lists
1845 :
1846 : ! **************************************************************************************************
1847 : !> \brief Returns the location of atom in atoms(:count), or zero if it is absent.
1848 : !> \param atom ...
1849 : !> \param atoms ...
1850 : !> \param count ...
1851 : !> \return ...
1852 : ! **************************************************************************************************
1853 128 : PURE INTEGER FUNCTION find_atom_position(atom, atoms, count) RESULT(position)
1854 : INTEGER, INTENT(IN) :: atom
1855 : INTEGER, DIMENSION(:), INTENT(IN) :: atoms
1856 : INTEGER, INTENT(IN) :: count
1857 :
1858 : INTEGER :: i
1859 :
1860 128 : position = 0
1861 172 : DO i = 1, count
1862 172 : IF (atoms(i) == atom) THEN
1863 128 : position = i
1864 : RETURN
1865 : END IF
1866 : END DO
1867 :
1868 : END FUNCTION find_atom_position
1869 :
1870 : ! **************************************************************************************************
1871 : !> \brief Collects U_Pp from the large reference basis for one optimized atom block.
1872 : !> \param bs_env ...
1873 : !> \param atom ...
1874 : !> \param list ...
1875 : !> \param atom_U ...
1876 : ! **************************************************************************************************
1877 32 : SUBROUTINE build_optimized_atom_U(bs_env, atom, list, atom_U)
1878 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1879 : INTEGER, INTENT(IN) :: atom
1880 : TYPE(AB_block_list_type), INTENT(IN) :: list
1881 : TYPE(optimized_atom_U_type), INTENT(OUT) :: atom_U
1882 :
1883 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_optimized_atom_U'
1884 :
1885 : INTEGER :: AB_block, atom_A, atom_B, first, first_p_AB, first_p_atom, first_row, handle, i, &
1886 : last, nRI_ref_A, nRI_ref_block, nsupport, position, size_opt_RI_AB
1887 : INTEGER, ALLOCATABLE :: support_atoms(:)
1888 32 : REAL(KIND=dp), ALLOCATABLE :: U_Pp_AB(:, :)
1889 :
1890 32 : CALL timeset(routineN, handle)
1891 :
1892 96 : ALLOCATE (support_atoms(2*SIZE(list%block)))
1893 32 : nsupport = 0
1894 80 : DO i = 1, SIZE(list%block)
1895 48 : AB_block = list%block(i)
1896 48 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
1897 48 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
1898 48 : position = find_atom_position(atom_A, support_atoms, nsupport)
1899 48 : IF (position == 0) THEN
1900 42 : nsupport = nsupport + 1
1901 42 : support_atoms(nsupport) = atom_A
1902 : END IF
1903 80 : IF (atom_B /= atom_A) THEN
1904 16 : position = find_atom_position(atom_B, support_atoms, nsupport)
1905 16 : IF (position == 0) THEN
1906 6 : nsupport = nsupport + 1
1907 6 : support_atoms(nsupport) = atom_B
1908 : END IF
1909 : END IF
1910 : END DO
1911 32 : CPASSERT(nsupport > 0)
1912 128 : ALLOCATE (atom_U%ref_atom(nsupport), atom_U%first_P(nsupport))
1913 80 : atom_U%ref_atom(:) = support_atoms(:nsupport)
1914 : first_row = 1
1915 80 : DO i = 1, nsupport
1916 48 : atom_U%first_P(i) = first_row
1917 80 : first_row = first_row + get_ref_RI_size(bs_env, atom_U%ref_atom(i))
1918 : END DO
1919 0 : ALLOCATE (atom_U%U_Pp(first_row - 1, bs_env%auto_ri%sizes_opt_RI(atom)), &
1920 25356 : source=0.0_dp)
1921 :
1922 80 : DO i = 1, SIZE(list%block)
1923 48 : AB_block = list%block(i)
1924 48 : atom_A = bs_env%auto_ri%AB_atom_A(AB_block)
1925 48 : atom_B = bs_env%auto_ri%AB_atom_B(AB_block)
1926 48 : nRI_ref_block = bs_env%auto_ri%AB_size_ref_RI(AB_block)
1927 48 : size_opt_RI_AB = bs_env%auto_ri%AB_size_opt_RI(AB_block)
1928 48 : first = bs_env%auto_ri%U_Pp_AB_offset(AB_block)
1929 48 : last = first + nRI_ref_block*size_opt_RI_AB - 1
1930 192 : ALLOCATE (U_Pp_AB(nRI_ref_block, size_opt_RI_AB))
1931 : U_Pp_AB(:, :) = RESHAPE(bs_env%auto_ri%U_Pp_AB(first:last), &
1932 144 : [nRI_ref_block, size_opt_RI_AB])
1933 48 : first_p_AB = list%first_p_AB(i)
1934 48 : first_p_atom = list%first_p_atom(i)
1935 48 : size_opt_RI_AB = list%size_opt_RI(i)
1936 48 : IF (first_p_AB < 1 .OR. &
1937 : first_p_AB + size_opt_RI_AB - 1 > SIZE(U_Pp_AB, 2)) THEN
1938 0 : CPABORT("AUTO_RI AB-block column range is invalid")
1939 : END IF
1940 48 : IF (first_p_atom < 1 .OR. &
1941 : first_p_atom + size_opt_RI_AB - 1 > bs_env%auto_ri%sizes_opt_RI(atom)) THEN
1942 0 : CPABORT("AUTO_RI optimized atom-block column range is invalid")
1943 : END IF
1944 :
1945 48 : nRI_ref_A = get_ref_RI_size(bs_env, atom_A)
1946 48 : position = find_atom_position(atom_A, atom_U%ref_atom, nsupport)
1947 48 : CPASSERT(position > 0)
1948 48 : first_row = atom_U%first_P(position)
1949 : atom_U%U_Pp(first_row:first_row + nRI_ref_A - 1, &
1950 : first_p_atom:first_p_atom + size_opt_RI_AB - 1) = &
1951 19048 : U_Pp_AB(1:nRI_ref_A, first_p_AB:first_p_AB + size_opt_RI_AB - 1)
1952 48 : IF (atom_B /= atom_A) THEN
1953 16 : position = find_atom_position(atom_B, atom_U%ref_atom, nsupport)
1954 16 : CPASSERT(position > 0)
1955 16 : first_row = atom_U%first_P(position)
1956 : atom_U%U_Pp(first_row:first_row + get_ref_RI_size(bs_env, atom_B) - 1, &
1957 : first_p_atom:first_p_atom + size_opt_RI_AB - 1) = &
1958 2592 : U_Pp_AB(nRI_ref_A + 1:, first_p_AB:first_p_AB + size_opt_RI_AB - 1)
1959 : END IF
1960 80 : DEALLOCATE (U_Pp_AB)
1961 : END DO
1962 32 : DEALLOCATE (support_atoms)
1963 :
1964 32 : CALL timestop(handle)
1965 :
1966 64 : END SUBROUTINE build_optimized_atom_U
1967 :
1968 : ! **************************************************************************************************
1969 : !> \brief Adds REGULARIZATION_RI to one atom-local reference metric block.
1970 : !> \param bs_env ...
1971 : !> \param atom ...
1972 : !> \param M_PQ atom-local metric block in the large reference RI basis
1973 : ! **************************************************************************************************
1974 0 : SUBROUTINE regularize_M_PQ_block(bs_env, atom, M_PQ)
1975 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1976 : INTEGER, INTENT(IN) :: atom
1977 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: M_PQ
1978 :
1979 : CHARACTER(LEN=*), PARAMETER :: routineN = 'regularize_M_PQ_block'
1980 :
1981 : INTEGER :: handle, i, iset, kind, nset
1982 0 : INTEGER, DIMENSION(:), POINTER :: npgf, nsgf_set
1983 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1984 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
1985 :
1986 0 : CALL timeset(routineN, handle)
1987 :
1988 0 : CPASSERT(SIZE(M_PQ, 1) == SIZE(M_PQ, 2))
1989 0 : kind = bs_env%ri_rs%particle_set(atom)%atomic_kind%kind_number
1990 : CALL get_gto_basis_set(bs_env%basis_set_RI(kind)%gto_basis_set, &
1991 : first_sgf=first_sgf, npgf=npgf, nset=nset, &
1992 0 : nsgf_set=nsgf_set, zet=zet)
1993 0 : DO iset = 1, nset
1994 0 : DO i = first_sgf(1, iset), first_sgf(1, iset) + nsgf_set(iset) - 1
1995 : M_PQ(i, i) = M_PQ(i, i) + &
1996 : bs_env%regularization_RI* &
1997 0 : MAX(1.0_dp, 1.0_dp/MINVAL(zet(1:npgf(iset), iset)))
1998 : END DO
1999 : END DO
2000 :
2001 0 : CALL timestop(handle)
2002 :
2003 0 : END SUBROUTINE regularize_M_PQ_block
2004 :
2005 : ! **************************************************************************************************
2006 : !> \brief Directly contracts one optimized atom-pair metric block.
2007 : !>
2008 : !> M_pq^{AB} = Σ_IJ (U_Ip^A)^T M_IJ U_Jq^B. Each reference
2009 : !> atom-pair block M_IJ is discarded immediately after this contraction.
2010 : !> \param bs_env ...
2011 : !> \param atom_U_a contractions U_Ip for optimized block A
2012 : !> \param atom_U_b contractions U_Jq for optimized block B
2013 : !> \param coulomb true for Coulomb; false for the GW fitting operator
2014 : !> \param symmetric true when both transforms describe the same optimized block
2015 : !> \param matrix_optimized contracted M_pq or V_pq block in the optimized RI basis
2016 : ! **************************************************************************************************
2017 24 : SUBROUTINE compute_optimized_RI_matrix_block( &
2018 : bs_env, atom_U_a, atom_U_b, coulomb, symmetric, matrix_optimized)
2019 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2020 : TYPE(optimized_atom_U_type), INTENT(IN) :: atom_U_a, atom_U_b
2021 : LOGICAL, INTENT(IN) :: coulomb, symmetric
2022 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
2023 : INTENT(OUT) :: matrix_optimized
2024 :
2025 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_optimized_RI_matrix_block'
2026 :
2027 : INTEGER :: atom_A, atom_B, first_a, first_b, &
2028 : handle, i, j, last_a, last_b, &
2029 : nRI_small_A, nRI_small_B
2030 24 : REAL(KIND=dp), ALLOCATABLE :: matrix_pq_contribution(:, :), matrix_reference(:, :), &
2031 24 : matrix_reference_times_U_b(:, :)
2032 :
2033 24 : CALL timeset(routineN, handle)
2034 :
2035 24 : nRI_small_A = SIZE(atom_U_a%U_Pp, 2)
2036 24 : nRI_small_B = SIZE(atom_U_b%U_Pp, 2)
2037 24 : IF (symmetric) THEN
2038 12 : CPASSERT(nRI_small_A == nRI_small_B)
2039 : END IF
2040 96 : ALLOCATE (matrix_optimized(nRI_small_A, nRI_small_B), source=0.0_dp)
2041 62 : DO i = 1, SIZE(atom_U_a%ref_atom)
2042 38 : atom_A = atom_U_a%ref_atom(i)
2043 38 : first_a = atom_U_a%first_P(i)
2044 38 : last_a = first_a + get_ref_RI_size(bs_env, atom_A) - 1
2045 124 : DO j = 1, SIZE(atom_U_b%ref_atom)
2046 62 : IF (symmetric .AND. j > i) CYCLE
2047 54 : atom_B = atom_U_b%ref_atom(j)
2048 54 : first_b = atom_U_b%first_P(j)
2049 54 : last_b = first_b + get_ref_RI_size(bs_env, atom_B) - 1
2050 54 : IF (coulomb) THEN
2051 27 : CALL compute_M_PQ_AB(bs_env, atom_A, atom_B, matrix_reference)
2052 : ELSE
2053 : CALL compute_M_PQ_AB( &
2054 27 : bs_env, atom_A, atom_B, matrix_reference, bs_env%ri_metric)
2055 27 : IF (atom_A == atom_B .AND. bs_env%regularization_RI > 0.0_dp) THEN
2056 0 : CALL regularize_M_PQ_block(bs_env, atom_A, matrix_reference)
2057 : END IF
2058 : END IF
2059 216 : ALLOCATE (matrix_reference_times_U_b(SIZE(matrix_reference, 1), nRI_small_B))
2060 : matrix_reference_times_U_b(:, :) = &
2061 694010 : MATMUL(matrix_reference, atom_U_b%U_Pp(first_b:last_b, :))
2062 162 : ALLOCATE (matrix_pq_contribution(nRI_small_A, nRI_small_B))
2063 : matrix_pq_contribution(:, :) = &
2064 54 : MATMUL(TRANSPOSE(atom_U_a%U_Pp(first_a:last_a, :)), &
2065 153104 : matrix_reference_times_U_b)
2066 2008 : matrix_optimized(:, :) = matrix_optimized + matrix_pq_contribution
2067 54 : IF (symmetric .AND. i /= j) THEN
2068 320 : matrix_optimized(:, :) = matrix_optimized + TRANSPOSE(matrix_pq_contribution)
2069 : END IF
2070 100 : DEALLOCATE (matrix_pq_contribution, matrix_reference, matrix_reference_times_U_b)
2071 : END DO
2072 : END DO
2073 798 : IF (.NOT. ALL(ieee_is_finite(matrix_optimized))) THEN
2074 0 : CPABORT("AUTO_RI optimized metric block contains invalid values")
2075 : END IF
2076 :
2077 24 : CALL timestop(handle)
2078 :
2079 48 : END SUBROUTINE compute_optimized_RI_matrix_block
2080 :
2081 : ! **************************************************************************************************
2082 : !> \brief Builds M_opt or V_opt directly from contracted atom-pair blocks.
2083 : !> \param qs_env ...
2084 : !> \param bs_env ...
2085 : !> \param lists ...
2086 : !> \param matrix_pq contracted M_pq or V_pq in the optimized RI basis
2087 : !> \param coulomb true for Coulomb; false for the GW fitting operator
2088 : ! **************************************************************************************************
2089 8 : SUBROUTINE compute_contracted_RI_matrix(qs_env, bs_env, lists, matrix_pq, coulomb)
2090 : TYPE(qs_environment_type), POINTER :: qs_env
2091 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2092 : TYPE(AB_block_list_type), DIMENSION(:), INTENT(IN) :: lists
2093 : TYPE(cp_fm_type), INTENT(OUT) :: matrix_pq
2094 : LOGICAL, INTENT(IN) :: coulomb
2095 :
2096 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_contracted_RI_matrix'
2097 :
2098 : INTEGER :: handle, iatom, jatom, npcol, nprow, &
2099 : pair_index, RI_size_opt
2100 8 : INTEGER, POINTER :: col_dist(:), row_dist(:), sizes_opt_RI(:)
2101 : LOGICAL :: atom_U_a_ready
2102 8 : REAL(KIND=dp), ALLOCATABLE :: matrix_pq_block(:, :)
2103 : TYPE(dbcsr_distribution_type) :: dist
2104 : TYPE(dbcsr_distribution_type), POINTER :: template
2105 : TYPE(dbcsr_type) :: blocks
2106 8 : TYPE(optimized_atom_U_type) :: atom_U_a, atom_U_b
2107 :
2108 8 : CALL timeset(routineN, handle)
2109 :
2110 40 : ALLOCATE (sizes_opt_RI(bs_env%n_atom), row_dist(bs_env%n_atom), col_dist(bs_env%n_atom))
2111 64 : sizes_opt_RI(:) = bs_env%auto_ri%sizes_opt_RI
2112 8 : CALL get_qs_env(qs_env, dbcsr_dist=template)
2113 8 : CALL dbcsr_distribution_get(template, nprows=nprow, npcols=npcol)
2114 32 : DO iatom = 1, SIZE(sizes_opt_RI)
2115 24 : row_dist(iatom) = MOD(iatom - 1, nprow)
2116 32 : col_dist(iatom) = MOD(iatom - 1, npcol)
2117 : END DO
2118 32 : RI_size_opt = SUM(sizes_opt_RI)
2119 8 : CALL create_distributed_matrix(qs_env, matrix_pq, RI_size_opt, RI_size_opt)
2120 8 : CALL dbcsr_distribution_new(dist, template=template, row_dist=row_dist, col_dist=col_dist)
2121 : CALL dbcsr_create(blocks, name='AUTO_RI contracted metric', dist=dist, &
2122 : matrix_type=dbcsr_type_symmetric, &
2123 8 : row_blk_size=sizes_opt_RI, col_blk_size=sizes_opt_RI)
2124 :
2125 8 : pair_index = 0
2126 32 : DO iatom = 1, SIZE(sizes_opt_RI)
2127 24 : atom_U_a_ready = .FALSE.
2128 80 : DO jatom = 1, iatom
2129 48 : pair_index = pair_index + 1
2130 48 : IF (MOD(pair_index - 1, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
2131 24 : IF (.NOT. atom_U_a_ready) THEN
2132 20 : CALL build_optimized_atom_U(bs_env, iatom, lists(iatom), atom_U_a)
2133 20 : atom_U_a_ready = .TRUE.
2134 : END IF
2135 24 : IF (iatom == jatom) THEN
2136 : CALL compute_optimized_RI_matrix_block( &
2137 : bs_env, atom_U_a, atom_U_a, coulomb, .TRUE., &
2138 12 : matrix_pq_block)
2139 912 : matrix_pq_block(:, :) = 0.5_dp*(matrix_pq_block + TRANSPOSE(matrix_pq_block))
2140 : ELSE
2141 12 : CALL build_optimized_atom_U(bs_env, jatom, lists(jatom), atom_U_b)
2142 : CALL compute_optimized_RI_matrix_block( &
2143 : bs_env, atom_U_a, atom_U_b, coulomb, .FALSE., &
2144 12 : matrix_pq_block)
2145 : END IF
2146 24 : CALL dbcsr_put_block(blocks, iatom, jatom, matrix_pq_block)
2147 72 : DEALLOCATE (matrix_pq_block)
2148 : END DO
2149 : END DO
2150 8 : CALL dbcsr_finalize(blocks)
2151 8 : CALL copy_dbcsr_to_fm(blocks, matrix_pq)
2152 8 : CALL dbcsr_release(blocks)
2153 8 : CALL dbcsr_distribution_release(dist)
2154 8 : DEALLOCATE (sizes_opt_RI, row_dist, col_dist)
2155 :
2156 8 : CALL timestop(handle)
2157 :
2158 32 : END SUBROUTINE compute_contracted_RI_matrix
2159 :
2160 : ! **************************************************************************************************
2161 : !> \brief Computes the pseudoinverse of the optimized fitting metric M_opt.
2162 : !>
2163 : !> The generalized problem M_opt X=M_AA X λ separates linearly
2164 : !> dependent directions. The result is M_opt^+=X diag(1/λ) X^T.
2165 : !> \param qs_env ...
2166 : !> \param bs_env ...
2167 : !> \param M_pq optimized RI fitting metric
2168 : ! **************************************************************************************************
2169 4 : SUBROUTINE compute_M_pq_inv(qs_env, bs_env, M_pq)
2170 : TYPE(qs_environment_type), POINTER :: qs_env
2171 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2172 : TYPE(cp_fm_type), INTENT(INOUT) :: M_pq
2173 :
2174 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_M_pq_inv'
2175 :
2176 : INTEGER :: handle, ncol_local, nrow_local, &
2177 : RI_size_opt
2178 4 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2179 4 : REAL(KIND=dp), ALLOCATABLE :: eigenval(:), factors(:)
2180 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2181 4 : POINTER :: M_pq_AA_local, M_pq_local
2182 : TYPE(cp_fm_type) :: M_pq_AA, work, X_pq
2183 :
2184 4 : CALL timeset(routineN, handle)
2185 :
2186 16 : RI_size_opt = SUM(bs_env%auto_ri%sizes_opt_RI)
2187 : CALL create_distributed_matrix(qs_env, bs_env%auto_ri%M_pq_inv, &
2188 4 : RI_size_opt, RI_size_opt)
2189 4 : CALL create_distributed_matrix(qs_env, X_pq, RI_size_opt, RI_size_opt)
2190 4 : CALL create_distributed_matrix(qs_env, M_pq_AA, RI_size_opt, RI_size_opt)
2191 4 : CALL create_distributed_matrix(qs_env, work, RI_size_opt, RI_size_opt)
2192 : CALL cp_fm_get_info(M_pq, nrow_local=nrow_local, ncol_local=ncol_local, &
2193 : row_indices=row_indices, col_indices=col_indices, &
2194 4 : local_data=M_pq_local)
2195 4 : CALL cp_fm_get_info(M_pq_AA, local_data=M_pq_AA_local)
2196 4 : M_pq_AA_local = 0.0_dp
2197 : CALL copy_M_pq_AA_blocks(M_pq_local(:nrow_local, :ncol_local), &
2198 : M_pq_AA_local(:nrow_local, :ncol_local), &
2199 : row_indices(:nrow_local), col_indices(:ncol_local), &
2200 4 : bs_env%auto_ri%sizes_opt_RI)
2201 :
2202 12 : ALLOCATE (eigenval(RI_size_opt))
2203 4 : CALL cp_fm_geeig(M_pq, M_pq_AA, X_pq, eigenval, work)
2204 4 : CALL compute_optimized_M_pq_inv_factors(eigenval, factors)
2205 4 : CALL cp_fm_column_scale(X_pq, factors)
2206 : CALL cp_fm_syrk("U", "N", RI_size_opt, 1.0_dp, X_pq, 1, 1, 0.0_dp, &
2207 4 : bs_env%auto_ri%M_pq_inv)
2208 4 : CALL cp_fm_uplo_to_full(bs_env%auto_ri%M_pq_inv, work)
2209 :
2210 4 : DEALLOCATE (eigenval, factors)
2211 4 : CALL cp_fm_release(X_pq)
2212 4 : CALL cp_fm_release(M_pq_AA)
2213 4 : CALL cp_fm_release(work)
2214 :
2215 4 : CALL timestop(handle)
2216 :
2217 8 : END SUBROUTINE compute_M_pq_inv
2218 :
2219 : ! **************************************************************************************************
2220 : !> \brief Computes M_opt and V_opt without constructing M_ref or V_ref.
2221 : !> \param qs_env ...
2222 : !> \param bs_env ...
2223 : ! **************************************************************************************************
2224 4 : SUBROUTINE compute_M_pq_inv_and_V_pq(qs_env, bs_env)
2225 : TYPE(qs_environment_type), POINTER :: qs_env
2226 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2227 :
2228 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_M_pq_inv_and_V_pq'
2229 :
2230 : INTEGER :: handle
2231 4 : TYPE(AB_block_list_type), ALLOCATABLE :: lists(:)
2232 : TYPE(cp_fm_type) :: M_pq
2233 :
2234 4 : CALL timeset(routineN, handle)
2235 :
2236 4 : CALL build_AB_block_lists(bs_env, lists)
2237 4 : CALL initialize_auto_ri_integrals(bs_env)
2238 4 : CALL compute_contracted_RI_matrix(qs_env, bs_env, lists, M_pq, .FALSE.)
2239 4 : CALL compute_M_pq_inv(qs_env, bs_env, M_pq)
2240 4 : CALL cp_fm_release(M_pq)
2241 4 : CALL compute_contracted_RI_matrix(qs_env, bs_env, lists, bs_env%auto_ri%V_pq, .TRUE.)
2242 :
2243 : ! Activate the optimized per-atom RI block sizes for RI-RS.
2244 4 : CALL set_optimized_RI_basis_sizes(bs_env)
2245 :
2246 4 : CALL timestop(handle)
2247 :
2248 20 : END SUBROUTINE compute_M_pq_inv_and_V_pq
2249 :
2250 : ! **************************************************************************************************
2251 : !> \brief Computes the atom-local two-center RI metric (P|Q) and the
2252 : !> three-center integrals (μν|P) used in Eqs. (3) and (4).
2253 : !> \param ao_basis ...
2254 : !> \param ri_basis ...
2255 : !> \param potential fitting operator used for both integral types
2256 : !> \param M_PQ two-center metric of the large reference RI basis
2257 : !> \param Int_3c three-center integrals (μν|P)
2258 : ! **************************************************************************************************
2259 8 : SUBROUTINE compute_int_3c(ao_basis, ri_basis, potential, M_PQ, Int_3c)
2260 : TYPE(gto_basis_set_type), INTENT(IN) :: ao_basis, ri_basis
2261 : TYPE(coulomb_operator_type), INTENT(IN) :: potential
2262 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
2263 : INTENT(OUT) :: M_PQ
2264 : REAL(KIND=dp), INTENT(OUT) :: Int_3c(:, :, :)
2265 :
2266 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_int_3c'
2267 :
2268 : INTEGER :: handle, iset, jset, kset, max_ao_l, &
2269 : max_l, max_ri_l, nAO, ncoa, ncob, &
2270 : ncoc, nRI, nsgfa, nsgfb, nsgfc, sgfa, &
2271 : sgfb, sgfc
2272 8 : INTEGER, DIMENSION(:), POINTER :: lmax_a, lmax_c, lmin_a, lmin_c, npgf_a, &
2273 8 : npgf_c, nsgf_set_a, nsgf_set_c
2274 8 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_a, first_sgf_c
2275 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Int_2c_cart, Int_2c_spherical
2276 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Int_3c_cart, Int_3c_spherical_all, &
2277 8 : Int_3c_spherical_block
2278 : REAL(KIND=dp), DIMENSION(3) :: center
2279 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: radius_a, radius_b, radius_c
2280 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius_a, pgf_radius_c, sphi_a, &
2281 8 : sphi_c, zeta_a, zeta_c
2282 : TYPE(coulomb_integral_context_type) :: context
2283 :
2284 8 : CALL timeset(routineN, handle)
2285 :
2286 : CALL get_gto_basis_set(ao_basis, nsgf=nAO, lmax=lmax_a, lmin=lmin_a, &
2287 : npgf=npgf_a, nsgf_set=nsgf_set_a, first_sgf=first_sgf_a, &
2288 8 : pgf_radius=pgf_radius_a, sphi=sphi_a, zet=zeta_a)
2289 : CALL get_gto_basis_set(ri_basis, nsgf=nRI, lmax=lmax_c, lmin=lmin_c, &
2290 : npgf=npgf_c, nsgf_set=nsgf_set_c, first_sgf=first_sgf_c, &
2291 8 : pgf_radius=pgf_radius_c, sphi=sphi_c, zet=zeta_c)
2292 8 : CPASSERT(nAO > 0 .AND. nRI > 0)
2293 20 : max_ao_l = MAXVAL(lmax_a)
2294 180 : max_ri_l = MAXVAL(lmax_c)
2295 8 : max_l = MAX(max_ao_l, max_ri_l)
2296 8 : CALL init_orbital_pointers(MAX(max_l, 2*max_ao_l))
2297 8 : CALL init_spherical_harmonics(max_l, -1)
2298 8 : IF (potential%potential_type == do_potential_truncated) THEN
2299 8 : IF (active_integral_library%coulomb2_library /= library_libint .OR. &
2300 : active_integral_library%coulomb3_library /= library_libint) THEN
2301 : CALL cp_abort(__LOCATION__, &
2302 : "AUTO_RI TRUNCATED fitting requires LIBINT for "// &
2303 0 : "COULOMB_2C and COULOMB_3C")
2304 : END IF
2305 : END IF
2306 : ! Libint initialization takes the maximum angular momentum of one shell,
2307 : ! not the sum of angular momenta entering an integral.
2308 8 : CALL coulomb_integral_init(context, max_am_2c=max_ri_l, max_am_3c=max_l)
2309 8 : center = 0.0_dp
2310 :
2311 64 : ALLOCATE (M_PQ(nRI, nRI), Int_3c_spherical_all(nAO, nAO, nRI))
2312 8 : M_PQ = 0.0_dp
2313 8 : Int_3c_spherical_all = 0.0_dp
2314 :
2315 180 : DO iset = 1, ri_basis%nset
2316 172 : ncoa = npgf_c(iset)*ncoset(lmax_c(iset))
2317 172 : sgfa = first_sgf_c(1, iset)
2318 172 : nsgfa = nsgf_set_c(iset)
2319 172 : radius_a => pgf_radius_c(1:npgf_c(iset), iset)
2320 4216 : DO jset = 1, ri_basis%nset
2321 4036 : ncob = npgf_c(jset)*ncoset(lmax_c(jset))
2322 4036 : sgfb = first_sgf_c(1, jset)
2323 4036 : nsgfb = nsgf_set_c(jset)
2324 4036 : radius_b => pgf_radius_c(1:npgf_c(jset), jset)
2325 28252 : ALLOCATE (Int_2c_cart(ncoa, ncob), Int_2c_spherical(nsgfa, nsgfb))
2326 4036 : Int_2c_cart = 0.0_dp
2327 : CALL compute_coulomb_2c(context, lmin_c(iset), lmax_c(iset), &
2328 : lmin_c(jset), lmax_c(jset), &
2329 : npgf_c(iset), npgf_c(jset), zeta_c(1:npgf_c(iset), iset), &
2330 : zeta_c(1:npgf_c(jset), jset), &
2331 4036 : radius_a, radius_b, center, center, Int_2c_cart, potential)
2332 : CALL ab_contract(Int_2c_spherical, Int_2c_cart, sphi_c(:, sgfa:), sphi_c(:, sgfb:), &
2333 4036 : ncoa, ncob, nsgfa, nsgfb)
2334 99308 : M_PQ(sgfa:sgfa + nsgfa - 1, sgfb:sgfb + nsgfb - 1) = Int_2c_spherical
2335 4208 : DEALLOCATE (Int_2c_spherical, Int_2c_cart)
2336 : END DO
2337 : END DO
2338 :
2339 20 : DO iset = 1, ao_basis%nset
2340 12 : ncoa = npgf_a(iset)*ncoset(lmax_a(iset))
2341 12 : sgfa = first_sgf_a(1, iset)
2342 12 : nsgfa = nsgf_set_a(iset)
2343 12 : radius_a => pgf_radius_a(1:npgf_a(iset), iset)
2344 40 : DO jset = 1, ao_basis%nset
2345 20 : ncob = npgf_a(jset)*ncoset(lmax_a(jset))
2346 20 : sgfb = first_sgf_a(1, jset)
2347 20 : nsgfb = nsgf_set_a(jset)
2348 20 : radius_b => pgf_radius_a(1:npgf_a(jset), jset)
2349 540 : DO kset = 1, ri_basis%nset
2350 508 : ncoc = npgf_c(kset)*ncoset(lmax_c(kset))
2351 508 : sgfc = first_sgf_c(1, kset)
2352 508 : nsgfc = nsgf_set_c(kset)
2353 508 : radius_c => pgf_radius_c(1:npgf_c(kset), kset)
2354 0 : ALLOCATE (Int_3c_cart(ncoa, ncob, ncoc), &
2355 4572 : Int_3c_spherical_block(nsgfa, nsgfb, nsgfc))
2356 508 : Int_3c_cart = 0.0_dp
2357 : CALL compute_coulomb_3c(context, lmin_a(iset), lmax_a(iset), &
2358 : lmin_a(jset), lmax_a(jset), &
2359 : lmin_c(kset), lmax_c(kset), npgf_a(iset), npgf_a(jset), &
2360 : npgf_c(kset), zeta_a(1:npgf_a(iset), iset), &
2361 : zeta_a(1:npgf_a(jset), jset), zeta_c(1:npgf_c(kset), kset), &
2362 : radius_a, radius_b, radius_c, center, center, center, &
2363 508 : Int_3c_cart, potential)
2364 : CALL abc_contract(Int_3c_spherical_block, Int_3c_cart, &
2365 : sphi_a(:, sgfa:), sphi_a(:, sgfb:), &
2366 508 : sphi_c(:, sgfc:), ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
2367 : Int_3c_spherical_all(sgfa:sgfa + nsgfa - 1, sgfb:sgfb + nsgfb - 1, &
2368 21136 : sgfc:sgfc + nsgfc - 1) = Int_3c_spherical_block
2369 528 : DEALLOCATE (Int_3c_spherical_block, Int_3c_cart)
2370 : END DO
2371 : END DO
2372 : END DO
2373 :
2374 16540 : Int_3c(:, :, :) = Int_3c_spherical_all
2375 8 : CALL coulomb_integral_cleanup(context)
2376 78376 : IF (.NOT. ALL(ieee_is_finite(M_PQ))) THEN
2377 0 : CPABORT("AUTO_RI metric contains invalid values")
2378 : END IF
2379 16540 : IF (.NOT. ALL(ieee_is_finite(Int_3c))) THEN
2380 0 : CPABORT("AUTO_RI three-center integrals contain invalid values")
2381 : END IF
2382 :
2383 8 : CALL timestop(handle)
2384 :
2385 16 : END SUBROUTINE compute_int_3c
2386 :
2387 : ! **************************************************************************************************
2388 : !> \brief Initializes the integral kernels used by the AUTO_RI optimization.
2389 : !> \param bs_env ...
2390 : ! **************************************************************************************************
2391 8 : SUBROUTINE initialize_auto_ri_integrals(bs_env)
2392 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2393 :
2394 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_auto_ri_integrals'
2395 :
2396 : CHARACTER(LEN=default_path_length) :: table_file
2397 : INTEGER :: handle, ikind, max_ao_l, max_l, &
2398 : max_ri_l, table_lmax, unit_id
2399 :
2400 8 : CALL timeset(routineN, handle)
2401 :
2402 8 : CPASSERT(ALLOCATED(bs_env%basis_set_AO))
2403 8 : CPASSERT(ALLOCATED(bs_env%basis_set_RI))
2404 8 : CPASSERT(SIZE(bs_env%basis_set_AO) == SIZE(bs_env%basis_set_RI))
2405 8 : CPASSERT(SIZE(bs_env%basis_set_RI) > 0)
2406 8 : max_ao_l = 0
2407 8 : max_ri_l = 0
2408 24 : DO ikind = 1, SIZE(bs_env%basis_set_RI)
2409 16 : CPASSERT(ASSOCIATED(bs_env%basis_set_AO(ikind)%gto_basis_set))
2410 16 : CPASSERT(ASSOCIATED(bs_env%basis_set_RI(ikind)%gto_basis_set))
2411 40 : max_ao_l = MAX(max_ao_l, MAXVAL(bs_env%basis_set_AO(ikind)%gto_basis_set%lmax))
2412 368 : max_ri_l = MAX(max_ri_l, MAXVAL(bs_env%basis_set_RI(ikind)%gto_basis_set%lmax))
2413 : END DO
2414 8 : max_l = MAX(max_ao_l, max_ri_l)
2415 : ! The OS two-center kernel accesses l+1 even without force evaluation.
2416 8 : CALL init_orbital_pointers(max_l + 1)
2417 8 : CALL init_spherical_harmonics(max_l, -1)
2418 :
2419 8 : IF (bs_env%ri_metric%potential_type == do_potential_truncated) THEN
2420 : ! Cover both the RI-RI two-center and AO-AO-RI three-center kernels.
2421 8 : table_lmax = MAX(2*max_ri_l + 1, 2*max_ao_l + max_ri_l)
2422 8 : IF (get_lmax_init() < table_lmax) THEN
2423 2 : unit_id = -1
2424 2 : IF (bs_env%para_env%is_source()) THEN
2425 1 : table_file = discover_file(bs_env%ri_metric%filename)
2426 1 : CALL open_file(unit_number=unit_id, file_name=TRIM(table_file))
2427 : END IF
2428 2 : CALL init_t_c_g0(table_lmax, unit_id, bs_env%para_env%mepos, bs_env%para_env)
2429 2 : IF (bs_env%para_env%is_source()) CALL close_file(unit_id)
2430 : END IF
2431 : END IF
2432 :
2433 8 : CALL timestop(handle)
2434 :
2435 8 : END SUBROUTINE initialize_auto_ri_integrals
2436 :
2437 : ! **************************************************************************************************
2438 : !> \brief Computes the atom-pair two-center RI metric matrix (P_A|Q_B).
2439 : !> \param bs_env ...
2440 : !> \param atom_A ...
2441 : !> \param atom_B ...
2442 : !> \param M_PQ metric block in the large reference RI basis
2443 : !> \param ri_potential ...
2444 : ! **************************************************************************************************
2445 66 : SUBROUTINE compute_M_PQ_AB(bs_env, atom_A, atom_B, M_PQ, ri_potential)
2446 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2447 : INTEGER, INTENT(IN) :: atom_A, atom_B
2448 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
2449 : INTENT(OUT) :: M_PQ
2450 : TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: ri_potential
2451 :
2452 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_M_PQ_AB'
2453 :
2454 : INTEGER :: atom_1, atom_2, handle, kind_1, kind_2
2455 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: M_PQ_canonical
2456 : TYPE(libint_potential_type) :: potential
2457 :
2458 66 : CALL timeset(routineN, handle)
2459 :
2460 66 : CPASSERT(ASSOCIATED(bs_env%ri_rs%particle_set))
2461 66 : CPASSERT(atom_A >= 1 .AND. atom_A <= bs_env%n_atom)
2462 66 : CPASSERT(atom_B >= 1 .AND. atom_B <= bs_env%n_atom)
2463 66 : IF (atom_A >= atom_B) THEN
2464 : atom_1 = atom_A
2465 : atom_2 = atom_B
2466 : ELSE
2467 11 : atom_1 = atom_B
2468 11 : atom_2 = atom_A
2469 : END IF
2470 66 : kind_1 = bs_env%ri_rs%particle_set(atom_1)%atomic_kind%kind_number
2471 66 : kind_2 = bs_env%ri_rs%particle_set(atom_2)%atomic_kind%kind_number
2472 :
2473 66 : IF (PRESENT(ri_potential)) THEN
2474 39 : CPASSERT(ri_potential%potential_type == do_potential_truncated)
2475 39 : potential = ri_potential
2476 : ELSE
2477 : potential%potential_type = do_potential_coulomb
2478 : END IF
2479 : CALL compute_M_PQ( &
2480 : basis_a=bs_env%basis_set_RI(kind_1)%gto_basis_set, &
2481 : basis_b=bs_env%basis_set_RI(kind_2)%gto_basis_set, &
2482 : center_a=bs_env%ri_rs%particle_set(atom_1)%r, &
2483 : center_b=bs_env%ri_rs%particle_set(atom_2)%r, &
2484 66 : potential=potential, M_PQ=M_PQ_canonical)
2485 :
2486 66 : IF (atom_A >= atom_B) THEN
2487 55 : CALL MOVE_ALLOC(M_PQ_canonical, M_PQ)
2488 : ELSE
2489 44 : ALLOCATE (M_PQ(SIZE(M_PQ_canonical, 2), SIZE(M_PQ_canonical, 1)))
2490 74041 : M_PQ(:, :) = TRANSPOSE(M_PQ_canonical)
2491 11 : DEALLOCATE (M_PQ_canonical)
2492 : END IF
2493 :
2494 66 : CALL timestop(handle)
2495 :
2496 132 : END SUBROUTINE compute_M_PQ_AB
2497 :
2498 : ! **************************************************************************************************
2499 : !> \brief Computes one contracted two-center RI metric block (P_A|Q_B).
2500 : !> \param basis_a ...
2501 : !> \param basis_b ...
2502 : !> \param center_a ...
2503 : !> \param center_b ...
2504 : !> \param potential ...
2505 : !> \param M_PQ metric block in the large reference RI basis
2506 : ! **************************************************************************************************
2507 66 : SUBROUTINE compute_M_PQ(basis_a, basis_b, center_a, center_b, potential, M_PQ)
2508 : TYPE(gto_basis_set_type), POINTER :: basis_a, basis_b
2509 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: center_a, center_b
2510 : TYPE(libint_potential_type), INTENT(IN) :: potential
2511 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
2512 : INTENT(OUT) :: M_PQ
2513 :
2514 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_M_PQ'
2515 :
2516 : INTEGER :: handle
2517 :
2518 66 : CALL timeset(routineN, handle)
2519 :
2520 264 : ALLOCATE (M_PQ(basis_a%nsgf, basis_b%nsgf), source=0.0_dp)
2521 66 : SELECT CASE (potential%potential_type)
2522 : CASE (do_potential_truncated)
2523 : CALL int_operators_r12_ab_os(operator_truncated, M_PQ, rab=center_b - center_a, &
2524 : fba=basis_a, fbb=basis_b, &
2525 156 : r_cutoff=potential%cutoff_radius, calculate_forces=.FALSE.)
2526 : CASE (do_potential_coulomb)
2527 : CALL int_operators_r12_ab_os(operator_coulomb, M_PQ, rab=center_b - center_a, &
2528 108 : fba=basis_a, fbb=basis_b, calculate_forces=.FALSE.)
2529 : CASE DEFAULT
2530 66 : CPABORT("AUTO_RI two-center metric requires COULOMB or TRUNCATED")
2531 : END SELECT
2532 :
2533 66 : CALL timestop(handle)
2534 :
2535 66 : END SUBROUTINE compute_M_PQ
2536 :
2537 : ! **************************************************************************************************
2538 : !> \brief Returns the reference RI basis size of one atom.
2539 : !> \param bs_env ...
2540 : !> \param iatom ...
2541 : !> \return ...
2542 : ! **************************************************************************************************
2543 206 : INTEGER FUNCTION get_ref_RI_size(bs_env, iatom) RESULT(nRI_ref)
2544 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
2545 : INTEGER, INTENT(IN) :: iatom
2546 :
2547 : INTEGER :: ikind
2548 :
2549 206 : ikind = bs_env%ri_rs%particle_set(iatom)%atomic_kind%kind_number
2550 206 : nRI_ref = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
2551 :
2552 206 : END FUNCTION get_ref_RI_size
2553 :
2554 78 : END MODULE gw_auto_ri_generate_RI_basis
|
