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 Calculation of contracted, spherical Gaussian integrals using the solid harmonic
10 : !> Gaussian (SHG) integral scheme. Routines for the following two-center integrals:
11 : !> i) (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc.
12 : !> ii) (aba) and (abb) s-overlaps
13 : !> \par Literature
14 : !> T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008)
15 : !> T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure
16 : !> Theory, Wiley
17 : !> \par History
18 : !> created [05.2016]
19 : !> \author Dorothea Golze
20 : ! **************************************************************************************************
21 : MODULE generic_shg_integrals
22 : USE basis_set_types, ONLY: gto_basis_set_type
23 : USE constants_operator, ONLY: operator_coulomb,&
24 : operator_gauss,&
25 : operator_verf,&
26 : operator_verfc,&
27 : operator_vgauss
28 : USE construct_shg, ONLY: &
29 : construct_dev_shg_ab, construct_int_shg_ab, construct_overlap_shg_aba, &
30 : construct_overlap_shg_abb, dev_overlap_shg_aba, dev_overlap_shg_abb, get_W_matrix, &
31 : get_dW_matrix, get_real_scaled_solid_harmonic
32 : USE kinds, ONLY: dp
33 : USE orbital_pointers, ONLY: nsoset
34 : USE s_contract_shg, ONLY: &
35 : contract_s_overlap_aba, contract_s_overlap_abb, contract_s_ra2m_ab, &
36 : contract_sint_ab_chigh, contract_sint_ab_clow, s_coulomb_ab, s_gauss_ab, s_overlap_ab, &
37 : s_overlap_abb, s_ra2m_ab, s_verf_ab, s_verfc_ab, s_vgauss_ab
38 : #include "../base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals'
45 :
46 : PUBLIC :: int_operators_r12_ab_shg, int_overlap_ab_shg, &
47 : int_ra2m_ab_shg, int_overlap_aba_shg, int_overlap_abb_shg, &
48 : get_abb_same_kind, lri_precalc_angular_shg_part, &
49 : int_overlap_ab_shg_low, int_ra2m_ab_shg_low, int_overlap_aba_shg_low, &
50 : int_overlap_abb_shg_low
51 :
52 : ABSTRACT INTERFACE
53 : ! **************************************************************************************************
54 : !> \brief Interface for the calculation of integrals over s-functions and their scalar derivatives
55 : !> with respect to rab2
56 : !> \param la_max ...
57 : !> \param npgfa ...
58 : !> \param zeta ...
59 : !> \param lb_max ...
60 : !> \param npgfb ...
61 : !> \param zetb ...
62 : !> \param omega ...
63 : !> \param rab ...
64 : !> \param v matrix storing the integrals and scalar derivatives
65 : !> \param calculate_forces ...
66 : ! **************************************************************************************************
67 : SUBROUTINE ab_sint_shg(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
68 : USE kinds, ONLY: dp
69 : INTEGER, INTENT(IN) :: la_max, npgfa
70 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
71 : INTEGER, INTENT(IN) :: lb_max, npgfb
72 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
73 : REAL(KIND=dp), INTENT(IN) :: omega
74 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
75 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
76 : LOGICAL, INTENT(IN) :: calculate_forces
77 :
78 : END SUBROUTINE ab_sint_shg
79 : END INTERFACE
80 :
81 : ! **************************************************************************************************
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the SHG scheme
87 : !> \param r12_operator the integral operator, which depends on r12=|r1-r2|
88 : !> \param vab integral matrix of spherical contracted Gaussian functions
89 : !> \param dvab derivative of the integrals
90 : !> \param rab distance vector between center A and B
91 : !> \param fba basis at center A
92 : !> \param fbb basis at center B
93 : !> \param scona_shg SHG contraction matrix for A
94 : !> \param sconb_shg SHG contraction matrix for B
95 : !> \param omega parameter in the operator
96 : !> \param calculate_forces ...
97 : ! **************************************************************************************************
98 3461452 : SUBROUTINE int_operators_r12_ab_shg(r12_operator, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
99 : omega, calculate_forces)
100 :
101 : INTEGER, INTENT(IN) :: r12_operator
102 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
103 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
104 : OPTIONAL :: dvab
105 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
106 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
107 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg
108 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: omega
109 : LOGICAL, INTENT(IN) :: calculate_forces
110 :
111 : INTEGER :: la_max, lb_max
112 : REAL(KIND=dp) :: my_omega
113 3461452 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
114 3461452 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
115 :
116 : PROCEDURE(ab_sint_shg), POINTER :: s_operator_ab
117 :
118 3461452 : NULLIFY (s_operator_ab)
119 :
120 9566628 : la_max = MAXVAL(fba%lmax)
121 9567428 : lb_max = MAXVAL(fbb%lmax)
122 :
123 : CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
124 3461452 : my_omega = 1.0_dp
125 :
126 6922776 : SELECT CASE (r12_operator)
127 : CASE (operator_coulomb)
128 3461324 : s_operator_ab => s_coulomb_ab
129 : CASE (operator_verf)
130 32 : s_operator_ab => s_verf_ab
131 32 : IF (PRESENT(omega)) my_omega = omega
132 : CASE (operator_verfc)
133 32 : s_operator_ab => s_verfc_ab
134 32 : IF (PRESENT(omega)) my_omega = omega
135 : CASE (operator_vgauss)
136 32 : s_operator_ab => s_vgauss_ab
137 32 : IF (PRESENT(omega)) my_omega = omega
138 : CASE (operator_gauss)
139 32 : s_operator_ab => s_gauss_ab
140 32 : IF (PRESENT(omega)) my_omega = omega
141 : CASE DEFAULT
142 3461452 : CPABORT("Operator not available")
143 : END SELECT
144 :
145 : CALL int_operator_ab_shg_low(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
146 6922744 : my_omega, Waux_mat, dWaux_mat, calculate_forces)
147 :
148 3461452 : DEALLOCATE (Waux_mat, dWaux_mat)
149 :
150 3461452 : END SUBROUTINE int_operators_r12_ab_shg
151 :
152 : ! **************************************************************************************************
153 : !> \brief calculate overlap integrals (a,b)
154 : !> \param vab integral (a,b)
155 : !> \param dvab derivative of sab
156 : !> \param rab distance vector
157 : !> \param fba basis at center A
158 : !> \param fbb basis at center B
159 : !> \param scona_shg contraction matrix A
160 : !> \param sconb_shg contraxtion matrix B
161 : !> \param calculate_forces ...
162 : ! **************************************************************************************************
163 32 : SUBROUTINE int_overlap_ab_shg(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
164 : calculate_forces)
165 :
166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
167 : INTENT(INOUT) :: vab
168 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
169 : INTENT(INOUT) :: dvab
170 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
171 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
172 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg
173 : LOGICAL, INTENT(IN) :: calculate_forces
174 :
175 : INTEGER :: la_max, lb_max
176 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
177 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
178 :
179 352 : la_max = MAXVAL(fba%lmax)
180 512 : lb_max = MAXVAL(fbb%lmax)
181 :
182 : CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
183 :
184 : CALL int_overlap_ab_shg_low(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
185 32 : Waux_mat, dWaux_mat, .TRUE., calculate_forces, contraction_high=.TRUE.)
186 :
187 32 : DEALLOCATE (Waux_mat, dWaux_mat)
188 :
189 32 : END SUBROUTINE int_overlap_ab_shg
190 :
191 : ! **************************************************************************************************
192 : !> \brief Calcululates the two-center integrals of the type (a|(r-Ra)^(2m)|b) using the SHG scheme
193 : !> \param vab integral matrix of spherical contracted Gaussian functions
194 : !> \param dvab derivative of the integrals
195 : !> \param rab distance vector between center A and B
196 : !> \param fba basis at center A
197 : !> \param fbb basis at center B
198 : !> \param scon_ra2m contraction matrix for A including the combinatorial factors
199 : !> \param sconb_shg SHG contraction matrix for B
200 : !> \param m exponent in (r-Ra)^(2m) operator
201 : !> \param calculate_forces ...
202 : ! **************************************************************************************************
203 32 : SUBROUTINE int_ra2m_ab_shg(vab, dvab, rab, fba, fbb, scon_ra2m, sconb_shg, &
204 : m, calculate_forces)
205 :
206 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
207 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dvab
208 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
209 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
210 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scon_ra2m
211 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: sconb_shg
212 : INTEGER, INTENT(IN) :: m
213 : LOGICAL, INTENT(IN) :: calculate_forces
214 :
215 : INTEGER :: la_max, lb_max
216 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
217 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
218 :
219 352 : la_max = MAXVAL(fba%lmax)
220 512 : lb_max = MAXVAL(fbb%lmax)
221 :
222 : CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
223 : CALL int_ra2m_ab_shg_low(vab, dvab, rab, fba, fbb, sconb_shg, scon_ra2m, &
224 32 : m, Waux_mat, dWaux_mat, calculate_forces)
225 :
226 32 : DEALLOCATE (Waux_mat, dWaux_mat)
227 :
228 32 : END SUBROUTINE int_ra2m_ab_shg
229 :
230 : ! **************************************************************************************************
231 : !> \brief calculate integrals (a,b,fa)
232 : !> \param saba integral [aba]
233 : !> \param dsaba derivative of [aba]
234 : !> \param rab distance vector between A and B
235 : !> \param oba orbital basis at center A
236 : !> \param obb orbital basis at center B
237 : !> \param fba auxiliary basis set at center A
238 : !> \param scon_obb contraction matrix for orb bas on B
239 : !> \param scona_mix mixed contraction matrix orb + ri basis on A
240 : !> \param oba_index orbital basis index for scona_mix
241 : !> \param fba_index ri basis index for scona_mix
242 : !> \param cg_coeff Clebsch-Gordon coefficients
243 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
244 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
245 : !> \param calculate_forces ...
246 : ! **************************************************************************************************
247 178 : SUBROUTINE int_overlap_aba_shg(saba, dsaba, rab, oba, obb, fba, scon_obb, &
248 178 : scona_mix, oba_index, fba_index, &
249 178 : cg_coeff, cg_none0_list, ncg_none0, &
250 : calculate_forces)
251 :
252 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
253 : INTENT(INOUT) :: saba
254 : REAL(KIND=dp), ALLOCATABLE, &
255 : DIMENSION(:, :, :, :), INTENT(INOUT) :: dsaba
256 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
257 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fba
258 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_obb
259 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix
260 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: oba_index, fba_index
261 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
262 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
263 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
264 : LOGICAL, INTENT(IN) :: calculate_forces
265 :
266 : INTEGER :: laa_max, lb_max
267 178 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
268 178 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
269 :
270 2610 : laa_max = MAXVAL(oba%lmax) + MAXVAL(fba%lmax)
271 400 : lb_max = MAXVAL(obb%lmax)
272 :
273 2307190 : saba = 0.0_dp
274 952018 : IF (calculate_forces) dsaba = 0.0_dp
275 : CALL precalc_angular_shg_part(laa_max, lb_max, rab, Waux_mat, dWaux_mat, &
276 : calculate_forces)
277 : CALL int_overlap_aba_shg_low(saba, dsaba, rab, oba, obb, fba, &
278 : scon_obb, scona_mix, oba_index, fba_index, &
279 : cg_coeff, cg_none0_list, ncg_none0, &
280 178 : Waux_mat, dWaux_mat, .TRUE., calculate_forces)
281 :
282 178 : DEALLOCATE (Waux_mat, dWaux_mat)
283 :
284 178 : END SUBROUTINE int_overlap_aba_shg
285 :
286 : ! **************************************************************************************************
287 : !> \brief calculate integrals (a,b,fb)
288 : !> \param sabb integral [abb]
289 : !> \param dsabb derivative of [abb]
290 : !> \param rab distance vector between A and B
291 : !> \param oba orbital basis at center A
292 : !> \param obb orbital basis at center B
293 : !> \param fbb auxiliary basis set at center B
294 : !> \param scon_oba contraction matrix for orb bas on A
295 : !> \param sconb_mix mixed contraction matrix orb + ri basis on B
296 : !> \param obb_index orbital basis index for sconb_mix
297 : !> \param fbb_index ri basis index for sconb_mix
298 : !> \param cg_coeff Clebsch-Gordon coefficients
299 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
300 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
301 : !> \param calculate_forces ...
302 : ! **************************************************************************************************
303 16 : SUBROUTINE int_overlap_abb_shg(sabb, dsabb, rab, oba, obb, fbb, scon_oba, &
304 16 : sconb_mix, obb_index, fbb_index, &
305 16 : cg_coeff, cg_none0_list, ncg_none0, &
306 : calculate_forces)
307 :
308 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
309 : INTENT(INOUT) :: sabb
310 : REAL(KIND=dp), ALLOCATABLE, &
311 : DIMENSION(:, :, :, :), INTENT(INOUT) :: dsabb
312 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
313 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fbb
314 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_oba
315 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix
316 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: obb_index, fbb_index
317 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
318 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
319 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
320 : LOGICAL, INTENT(IN) :: calculate_forces
321 :
322 : INTEGER :: la_max, lbb_max
323 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
324 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
325 :
326 32 : la_max = MAXVAL(oba%lmax)
327 272 : lbb_max = MAXVAL(obb%lmax) + MAXVAL(fbb%lmax)
328 :
329 317280 : sabb = 0.0_dp
330 951856 : IF (calculate_forces) dsabb = 0.0_dp
331 : CALL precalc_angular_shg_part(lbb_max, la_max, rab, Waux_mat, dWaux_mat, &
332 : calculate_forces)
333 : CALL int_overlap_abb_shg_low(sabb, dsabb, rab, oba, obb, fbb, &
334 : scon_oba, sconb_mix, obb_index, fbb_index, &
335 : cg_coeff, cg_none0_list, ncg_none0, &
336 16 : Waux_mat, dWaux_mat, .TRUE., calculate_forces)
337 :
338 16 : DEALLOCATE (Waux_mat, dWaux_mat)
339 :
340 16 : END SUBROUTINE int_overlap_abb_shg
341 :
342 : ! **************************************************************************************************
343 : !> \brief precalculates the angular part of the SHG integrals
344 : !> \param la_max ...
345 : !> \param lb_max ...
346 : !> \param rab distance vector between a and b
347 : !> \param Waux_mat W matrix that contains the angular-dependent part
348 : !> \param dWaux_mat derivative of the W matrix
349 : !> \param calculate_forces ...
350 : ! **************************************************************************************************
351 3461710 : SUBROUTINE precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
352 :
353 : INTEGER, INTENT(IN) :: la_max, lb_max
354 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
355 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
356 : INTENT(OUT) :: Waux_mat
357 : REAL(KIND=dp), ALLOCATABLE, &
358 : DIMENSION(:, :, :, :), INTENT(OUT) :: dWaux_mat
359 : LOGICAL, INTENT(IN) :: calculate_forces
360 :
361 : INTEGER :: lmax, mdim(3)
362 : INTEGER, DIMENSION(:), POINTER :: la_max_all
363 : REAL(KIND=dp) :: rab2
364 3461710 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Rc, Rs
365 :
366 3461710 : NULLIFY (la_max_all)
367 3461710 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
368 3461710 : lmax = MAX(la_max, lb_max)
369 :
370 10385130 : ALLOCATE (la_max_all(0:lb_max))
371 20770260 : ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
372 47115174 : Rc = 0._dp
373 47115174 : Rs = 0._dp
374 3461710 : mdim(1) = MIN(la_max, lb_max) + 1
375 3461710 : mdim(2) = nsoset(la_max) + 1
376 3461710 : mdim(3) = nsoset(lb_max) + 1
377 17308550 : ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
378 17308550 : ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3)))
379 :
380 8570598 : la_max_all(0:lb_max) = la_max
381 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
382 13846840 : CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
383 3461710 : CALL get_W_matrix(la_max_all, lb_max, lmax, Rc, Rs, Waux_mat)
384 3461710 : IF (calculate_forces) THEN
385 256 : CALL get_dW_matrix(la_max_all, lb_max, Waux_mat, dWaux_mat)
386 : END IF
387 :
388 3461710 : DEALLOCATE (Rc, Rs, la_max_all)
389 :
390 3461710 : END SUBROUTINE precalc_angular_shg_part
391 :
392 : ! **************************************************************************************************
393 : !> \brief calculate integrals (a|O(r12)|b)
394 : !> \param s_operator_ab procedure pointer for the respective operator. The integral evaluation
395 : !> differs only in the calculation of the [s|O(r12)|s] integrals and their scalar
396 : !> derivatives.
397 : !> \param vab integral matrix of spherical contracted Gaussian functions
398 : !> \param dvab derivative of the integrals
399 : !> \param rab distance vector between center A and B
400 : !> \param fba basis at center A
401 : !> \param fbb basis at center B
402 : !> \param scona_shg SHG contraction matrix for A
403 : !> \param sconb_shg SHG contraction matrix for B
404 : !> \param omega parameter in the operator
405 : !> \param Waux_mat W matrix that contains the angular-dependent part
406 : !> \param dWaux_mat derivative of the W matrix
407 : !> \param calculate_forces ...
408 : ! **************************************************************************************************
409 3461452 : SUBROUTINE int_operator_ab_shg_low(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
410 3461452 : omega, Waux_mat, dWaux_mat, calculate_forces)
411 :
412 : PROCEDURE(ab_sint_shg), POINTER :: s_operator_ab
413 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
414 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, INTENT(INOUT) :: dvab
415 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
416 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
417 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: scona_shg, sconb_shg
418 : REAL(KIND=dp), INTENT(IN) :: omega
419 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
420 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
421 : LOGICAL, INTENT(IN) :: calculate_forces
422 :
423 : INTEGER :: iset, jset, la_max_set, lb_max_set, ndev, nds, nds_max, npgfa_set, &
424 : npgfb_set, nseta, nsetb, nsgfa_set, nsgfb_set, nshella_set, nshellb_set
425 3461452 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
426 3461452 : nsgfb, nshella, nshellb
427 3461452 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb
428 : REAL(KIND=dp) :: dab
429 3461452 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
430 3461452 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
431 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: swork, swork_cont
432 :
433 3461452 : NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
434 3461452 : set_radius_b, zeta, zetb)
435 :
436 : ! basis ikind
437 3461452 : first_sgfa => fba%first_sgf
438 3461452 : la_max => fba%lmax
439 3461452 : la => fba%l
440 3461452 : npgfa => fba%npgf
441 3461452 : nsgfa => fba%nsgf_set
442 3461452 : nseta = fba%nset
443 3461452 : set_radius_a => fba%set_radius
444 3461452 : zeta => fba%zet
445 3461452 : nshella => fba%nshell
446 : ! basis jkind
447 3461452 : first_sgfb => fbb%first_sgf
448 3461452 : lb_max => fbb%lmax
449 3461452 : lb => fbb%l
450 3461452 : npgfb => fbb%npgf
451 3461452 : nsgfb => fbb%nsgf_set
452 3461452 : nsetb = fbb%nset
453 3461452 : set_radius_b => fbb%set_radius
454 3461452 : zetb => fbb%zet
455 3461452 : nshellb => fbb%nshell
456 :
457 13845808 : dab = SQRT(SUM(rab**2))
458 :
459 9566628 : la_max_set = MAXVAL(la_max)
460 9567428 : lb_max_set = MAXVAL(lb_max)
461 :
462 : ! allocate some work matrices
463 9566628 : npgfa_set = MAXVAL(npgfa)
464 9567428 : npgfb_set = MAXVAL(npgfb)
465 9566628 : nshella_set = MAXVAL(nshella)
466 9567428 : nshellb_set = MAXVAL(nshellb)
467 9566628 : nsgfa_set = MAXVAL(nsgfa)
468 9567428 : nsgfb_set = MAXVAL(nsgfb)
469 3461452 : ndev = 0
470 3461452 : IF (calculate_forces) ndev = 1
471 3461452 : nds_max = la_max_set + lb_max_set + ndev + 1
472 17307260 : ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
473 17307260 : ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
474 :
475 164848664 : vab = 0.0_dp
476 17919692 : IF (calculate_forces) dvab = 0.0_dp
477 :
478 9566628 : DO iset = 1, nseta
479 :
480 31567288 : DO jset = 1, nsetb
481 :
482 22000660 : nds = la_max(iset) + lb_max(jset) + ndev + 1
483 187390072 : swork(1:npgfa(iset), 1:npgfb(jset), 1:nds) = 0.0_dp
484 : CALL s_operator_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
485 : lb_max(jset), npgfb(jset), zetb(:, jset), &
486 22000660 : omega, rab, swork, calculate_forces)
487 : CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), &
488 : scona_shg(1:npgfa(iset), 1:nshella(iset), iset), &
489 : npgfb(jset), nshellb(jset), &
490 : sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
491 : nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), &
492 22000660 : swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)))
493 : CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
494 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
495 22000660 : swork_cont, Waux_mat, vab)
496 28105836 : IF (calculate_forces) THEN
497 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
498 : CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
499 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
500 96000 : -rab, swork_cont, Waux_mat, dWaux_mat, dvab)
501 : END IF
502 : END DO
503 : END DO
504 :
505 3461452 : DEALLOCATE (swork, swork_cont)
506 :
507 3461452 : END SUBROUTINE int_operator_ab_shg_low
508 :
509 : ! **************************************************************************************************
510 : !> \brief calculate overlap integrals (a,b); requires angular-dependent part as input
511 : !> \param sab integral (a,b)
512 : !> \param dsab derivative of sab
513 : !> \param rab distance vector
514 : !> \param fba basis at center A
515 : !> \param fbb basis at center B
516 : !> \param scona_shg contraction matrix A
517 : !> \param sconb_shg contraxtion matrix B
518 : !> \param Waux_mat W matrix that contains the angular-dependent part
519 : !> \param dWaux_mat derivative of the W matrix
520 : !> \param calculate_ints ...
521 : !> \param calculate_forces ...
522 : !> \param contraction_high ...
523 : ! **************************************************************************************************
524 30578 : SUBROUTINE int_overlap_ab_shg_low(sab, dsab, rab, fba, fbb, scona_shg, sconb_shg, Waux_mat, dWaux_mat, &
525 : calculate_ints, calculate_forces, contraction_high)
526 :
527 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
528 : INTENT(INOUT) :: sab
529 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
530 : INTENT(INOUT) :: dsab
531 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
532 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
533 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg, Waux_mat
534 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
535 : LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
536 : LOGICAL, INTENT(IN), OPTIONAL :: contraction_high
537 :
538 : INTEGER :: iset, jset, la_max_set, lb_max_set, &
539 : ndev, nds, nds_max, npgfa_set, &
540 : npgfb_set, nseta, nsetb, nshella_set, &
541 : nshellb_set
542 30578 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nshella, &
543 30578 : nshellb
544 30578 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb
545 : LOGICAL :: my_contraction_high
546 : REAL(KIND=dp) :: dab
547 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: swork, swork_cont
548 30578 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
549 30578 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
550 :
551 30578 : NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
552 30578 : set_radius_b, zeta, zetb)
553 :
554 : ! basis ikind
555 30578 : first_sgfa => fba%first_sgf
556 30578 : la_max => fba%lmax
557 30578 : la => fba%l
558 30578 : npgfa => fba%npgf
559 30578 : nseta = fba%nset
560 30578 : set_radius_a => fba%set_radius
561 30578 : zeta => fba%zet
562 30578 : nshella => fba%nshell
563 : ! basis jkind
564 30578 : first_sgfb => fbb%first_sgf
565 30578 : lb_max => fbb%lmax
566 30578 : lb => fbb%l
567 30578 : npgfb => fbb%npgf
568 30578 : nsetb = fbb%nset
569 30578 : set_radius_b => fbb%set_radius
570 30578 : zetb => fbb%zet
571 30578 : nshellb => fbb%nshell
572 :
573 122312 : dab = SQRT(SUM(rab**2))
574 :
575 82433 : la_max_set = MAXVAL(la_max)
576 82593 : lb_max_set = MAXVAL(lb_max)
577 :
578 : ! allocate some work matrices
579 82433 : npgfa_set = MAXVAL(npgfa)
580 82593 : npgfb_set = MAXVAL(npgfb)
581 82433 : nshella_set = MAXVAL(nshella)
582 82593 : nshellb_set = MAXVAL(nshellb)
583 30578 : ndev = 0
584 30578 : IF (calculate_forces) ndev = 1
585 30578 : nds_max = la_max_set + lb_max_set + ndev + 1
586 152890 : ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
587 152890 : ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
588 :
589 15202264 : IF (calculate_ints) sab = 0.0_dp
590 27663425 : IF (calculate_forces) dsab = 0.0_dp
591 :
592 30578 : my_contraction_high = .TRUE.
593 30578 : IF (PRESENT(contraction_high)) my_contraction_high = contraction_high
594 :
595 82433 : DO iset = 1, nseta
596 :
597 333888 : DO jset = 1, nsetb
598 :
599 251455 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
600 :
601 210110 : nds = la_max(iset) + lb_max(jset) + ndev + 1
602 7948160 : swork(1:npgfa(iset), 1:npgfb(jset), 1:nds) = 0.0_dp
603 : CALL s_overlap_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
604 : lb_max(jset), npgfb(jset), zetb(:, jset), &
605 210110 : rab, swork, calculate_forces)
606 210110 : IF (my_contraction_high) THEN
607 : CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), &
608 : scona_shg(1:npgfa(iset), 1:nshella(iset), iset), &
609 : npgfb(jset), nshellb(jset), &
610 : sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
611 : nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), &
612 32555 : swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)))
613 : ELSE
614 : CALL contract_sint_ab_clow(la(:, iset), npgfa(iset), nshella(iset), &
615 : scona_shg(:, :, iset), &
616 : lb(:, jset), npgfb(jset), nshellb(jset), &
617 : sconb_shg(:, :, jset), &
618 177555 : swork, swork_cont, calculate_forces)
619 : END IF
620 210110 : IF (calculate_ints) THEN
621 : CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
622 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
623 126764 : swork_cont, Waux_mat, sab)
624 : END IF
625 261965 : IF (calculate_forces) THEN
626 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
627 : CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
628 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
629 352584 : -rab, swork_cont, Waux_mat, dWaux_mat, dsab)
630 : END IF
631 : END DO
632 : END DO
633 :
634 30578 : DEALLOCATE (swork, swork_cont)
635 :
636 30578 : END SUBROUTINE int_overlap_ab_shg_low
637 :
638 : ! **************************************************************************************************
639 : !> \brief calculate integrals (a|ra^2m)|b); requires angular-dependent part as input
640 : !> \param vab integral matrix of spherical contracted Gaussian functions
641 : !> \param dvab derivative of the integrals
642 : !> \param rab distance vector between center A and B
643 : !> \param fba basis at center A
644 : !> \param fbb basis at center B
645 : !> \param sconb_shg SHG contraction matrix for B
646 : !> \param scon_ra2m contraction matrix for A including the combinatorial factors
647 : !> \param m exponent in (r-Ra)^(2m) operator
648 : !> \param Waux_mat W matrix that contains the angular-dependent part
649 : !> \param dWaux_mat derivative of the W matrix
650 : !> \param calculate_forces ...
651 : ! **************************************************************************************************
652 32 : SUBROUTINE int_ra2m_ab_shg_low(vab, dvab, rab, fba, fbb, sconb_shg, scon_ra2m, m, Waux_mat, dWaux_mat, &
653 : calculate_forces)
654 :
655 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
656 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dvab
657 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
658 : TYPE(gto_basis_set_type), POINTER :: fba, fbb
659 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: sconb_shg
660 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scon_ra2m
661 : INTEGER, INTENT(IN) :: m
662 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
663 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
664 : LOGICAL, INTENT(IN) :: calculate_forces
665 :
666 : INTEGER :: iset, jset, la_max_set, lb_max_set, &
667 : ndev, nds, nds_max, npgfa_set, &
668 : npgfb_set, nseta, nsetb, nshella_set, &
669 : nshellb_set
670 32 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
671 32 : nsgfb, nshella, nshellb
672 32 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb
673 : REAL(KIND=dp) :: dab
674 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: swork_cont
675 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: swork
676 32 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
677 :
678 32 : NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, zeta, zetb)
679 :
680 : ! basis ikind
681 32 : first_sgfa => fba%first_sgf
682 32 : la_max => fba%lmax
683 32 : la => fba%l
684 32 : npgfa => fba%npgf
685 32 : nsgfa => fba%nsgf_set
686 32 : nseta = fba%nset
687 32 : zeta => fba%zet
688 32 : nshella => fba%nshell
689 : ! basis jkind
690 32 : first_sgfb => fbb%first_sgf
691 32 : lb_max => fbb%lmax
692 32 : lb => fbb%l
693 32 : npgfb => fbb%npgf
694 32 : nsgfb => fbb%nsgf_set
695 32 : nsetb = fbb%nset
696 32 : zetb => fbb%zet
697 32 : nshellb => fbb%nshell
698 :
699 128 : dab = SQRT(SUM(rab**2))
700 :
701 352 : la_max_set = MAXVAL(la_max)
702 512 : lb_max_set = MAXVAL(lb_max)
703 :
704 : ! allocate some work matrices
705 352 : npgfa_set = MAXVAL(npgfa)
706 512 : npgfb_set = MAXVAL(npgfb)
707 352 : nshella_set = MAXVAL(nshella)
708 512 : nshellb_set = MAXVAL(nshellb)
709 32 : ndev = 0
710 32 : IF (calculate_forces) ndev = 1
711 32 : nds_max = la_max_set + lb_max_set + ndev + 1
712 192 : ALLOCATE (swork(npgfa_set, npgfb_set, 1:m + 1, nds_max))
713 160 : ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
714 :
715 963872 : vab = 0.0_dp
716 2891680 : IF (calculate_forces) dvab = 0.0_dp
717 :
718 352 : DO iset = 1, nseta
719 :
720 5152 : DO jset = 1, nsetb
721 :
722 4800 : nds = la_max(iset) + lb_max(jset) + ndev + 1
723 447840 : swork(1:npgfa(iset), 1:npgfb(jset), 1:m + 1, 1:nds) = 0.0_dp
724 : CALL s_ra2m_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
725 : lb_max(jset), npgfb(jset), zetb(:, jset), &
726 4800 : m, rab, swork, calculate_forces)
727 : CALL contract_s_ra2m_ab(npgfa(iset), nshella(iset), &
728 : scon_ra2m(1:npgfa(iset), 1:m + 1, 1:nshella(iset), iset), &
729 : npgfb(jset), nshellb(jset), &
730 : sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
731 : swork(1:npgfa(iset), 1:npgfb(jset), 1:m + 1, 1:nds), &
732 : swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)), &
733 4800 : m, nds)
734 : CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
735 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
736 4800 : swork_cont, Waux_mat, vab)
737 5120 : IF (calculate_forces) THEN
738 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
739 : CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
740 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
741 19200 : -rab, swork_cont, Waux_mat, dWaux_mat, dvab)
742 : END IF
743 : END DO
744 : END DO
745 :
746 32 : DEALLOCATE (swork, swork_cont)
747 :
748 32 : END SUBROUTINE int_ra2m_ab_shg_low
749 : ! **************************************************************************************************
750 : !> \brief calculate integrals (a,b,fb); requires angular-dependent part as input
751 : !> \param abbint integral (a,b,fb)
752 : !> \param dabbint derivative of abbint
753 : !> \param rab distance vector between A and B
754 : !> \param oba orbital basis at center A
755 : !> \param obb orbital basis at center B
756 : !> \param fbb auxiliary basis set at center B
757 : !> \param scon_oba contraction matrix for orb bas on A
758 : !> \param sconb_mix mixed contraction matrix orb + ri basis on B
759 : !> \param obb_index orbital basis index for sconb_mix
760 : !> \param fbb_index ri basis index for sconb_mix
761 : !> \param cg_coeff Clebsch-Gordon coefficients
762 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
763 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
764 : !> \param Waux_mat W matrix that contains the angular-dependent part
765 : !> \param dWaux_mat derivative of the W matrix
766 : !> \param calculate_ints ...
767 : !> \param calculate_forces ...
768 : ! **************************************************************************************************
769 760 : SUBROUTINE int_overlap_abb_shg_low(abbint, dabbint, rab, oba, obb, fbb, scon_oba, sconb_mix, &
770 760 : obb_index, fbb_index, cg_coeff, cg_none0_list, &
771 760 : ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
772 :
773 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
774 : INTENT(INOUT) :: abbint
775 : REAL(KIND=dp), ALLOCATABLE, &
776 : DIMENSION(:, :, :, :), INTENT(INOUT) :: dabbint
777 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
778 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fbb
779 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_oba
780 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix
781 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: obb_index, fbb_index
782 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
783 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
784 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
785 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
786 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
787 : LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
788 :
789 : INTEGER :: iset, jset, kset, la_max_set, lb_max_set, lbb_max, lbb_max_set, lcb_max_set, na, &
790 : nb, ncb, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfcb_set, nseta, nsetb, &
791 : nsetcb, nshella_set, nshellb_set, nshellcb_set, sgfa, sgfb, sgfcb
792 760 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lcb_max, npgfa, npgfb, &
793 760 : npgfcb, nsgfa, nsgfb, nsgfcb, nshella, &
794 760 : nshellb, nshellcb
795 760 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfcb, la, &
796 760 : lb, lcb
797 : REAL(KIND=dp) :: dab, rab2
798 : REAL(KIND=dp), ALLOCATABLE, &
799 : DIMENSION(:, :, :, :, :) :: swork, swork_cont
800 760 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_cb
801 760 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetcb
802 :
803 760 : NULLIFY (la_max, lb_max, lcb_max, npgfa, npgfb, npgfcb)
804 760 : NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
805 760 : set_radius_cb, zeta, zetb, zetcb)
806 :
807 : ! basis ikind
808 760 : first_sgfa => oba%first_sgf
809 760 : la_max => oba%lmax
810 760 : la => oba%l
811 760 : nsgfa => oba%nsgf_set
812 760 : npgfa => oba%npgf
813 760 : nshella => oba%nshell
814 760 : nseta = oba%nset
815 760 : set_radius_a => oba%set_radius
816 760 : zeta => oba%zet
817 : ! basis jkind
818 760 : first_sgfb => obb%first_sgf
819 760 : lb_max => obb%lmax
820 760 : lb => obb%l
821 760 : nsgfb => obb%nsgf_set
822 760 : npgfb => obb%npgf
823 760 : nshellb => obb%nshell
824 760 : nsetb = obb%nset
825 760 : set_radius_b => obb%set_radius
826 760 : zetb => obb%zet
827 :
828 : ! basis RI on B
829 760 : first_sgfcb => fbb%first_sgf
830 760 : lcb_max => fbb%lmax
831 760 : lcb => fbb%l
832 760 : nsgfcb => fbb%nsgf_set
833 760 : npgfcb => fbb%npgf
834 760 : nshellcb => fbb%nshell
835 760 : nsetcb = fbb%nset
836 760 : set_radius_cb => fbb%set_radius
837 760 : zetcb => fbb%zet
838 :
839 3040 : dab = SQRT(SUM(rab**2))
840 760 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
841 :
842 1748 : la_max_set = MAXVAL(la_max)
843 1748 : lb_max_set = MAXVAL(lb_max)
844 7546 : lcb_max_set = MAXVAL(lcb_max)
845 1748 : npgfa_set = MAXVAL(npgfa)
846 1748 : npgfb_set = MAXVAL(npgfb)
847 7546 : npgfcb_set = MAXVAL(npgfcb)
848 1748 : nshella_set = MAXVAL(nshella)
849 1748 : nshellb_set = MAXVAL(nshellb)
850 7546 : nshellcb_set = MAXVAL(nshellcb)
851 : !*** for forces: derivative+1 in auxiliary vector required
852 760 : ndev = 0
853 760 : IF (calculate_forces) ndev = 1
854 :
855 760 : lbb_max_set = lb_max_set + lcb_max_set
856 :
857 : ! allocate some work storage....
858 760 : nds_max = la_max_set + lbb_max_set + ndev + 1
859 760 : nl_set = INT((lbb_max_set)/2)
860 5320 : ALLOCATE (swork(npgfa_set, npgfb_set, npgfcb_set, nl_set + 1, nds_max))
861 5320 : ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellcb_set))
862 :
863 1748 : DO iset = 1, nseta
864 :
865 3192 : DO jset = 1, nsetb
866 :
867 1444 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
868 :
869 8984 : DO kset = 1, nsetcb
870 :
871 7104 : IF (set_radius_a(iset) + set_radius_cb(kset) < dab) CYCLE
872 :
873 4292 : lbb_max = lb_max(jset) + lcb_max(kset)
874 4292 : nds = la_max(iset) + lbb_max + ndev + 1
875 4292 : nl = INT((lbb_max)/2) + 1
876 4972648 : swork(1:npgfa(iset), 1:npgfb(jset), 1:npgfcb(kset), 1:nl, 1:nds) = 0.0_dp
877 : CALL s_overlap_abb(la_max(iset), npgfa(iset), zeta(:, iset), &
878 : lb_max(jset), npgfb(jset), zetb(:, jset), &
879 : lcb_max(kset), npgfcb(kset), zetcb(:, kset), &
880 4292 : rab, swork, calculate_forces)
881 :
882 : CALL contract_s_overlap_abb(la(:, iset), npgfa(iset), nshella(iset), &
883 : scon_oba(1:npgfa(iset), 1:nshella(iset), iset), &
884 : lb(:, jset), npgfb(jset), nshellb(jset), &
885 : lcb(:, kset), npgfcb(kset), nshellcb(kset), &
886 : obb_index(:, :, jset), fbb_index(:, :, kset), sconb_mix, nl, nds, &
887 : swork(1:npgfa(iset), 1:npgfb(jset), 1:npgfcb(kset), 1:nl, 1:nds), &
888 4292 : swork_cont, calculate_forces)
889 :
890 4292 : IF (calculate_ints) THEN
891 : CALL construct_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
892 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
893 : lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
894 : cg_coeff, cg_none0_list, &
895 3421 : ncg_none0, swork_cont, Waux_mat, abbint)
896 : END IF
897 4292 : IF (calculate_forces) THEN
898 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
899 : CALL dev_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
900 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
901 : lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
902 : cg_coeff, cg_none0_list, ncg_none0, -rab, swork_cont, &
903 4444 : Waux_mat, dWaux_mat, dabbint)
904 : END IF
905 : ! max value of integrals in this set triple
906 4292 : sgfa = first_sgfa(1, iset)
907 4292 : na = sgfa + nsgfa(iset) - 1
908 4292 : sgfb = first_sgfb(1, jset)
909 4292 : nb = sgfb + nsgfb(jset) - 1
910 4292 : sgfcb = first_sgfcb(1, kset)
911 8548 : ncb = sgfcb + nsgfcb(kset) - 1
912 : END DO
913 : END DO
914 : END DO
915 :
916 760 : DEALLOCATE (swork_cont)
917 760 : DEALLOCATE (swork)
918 :
919 760 : END SUBROUTINE int_overlap_abb_shg_low
920 : ! **************************************************************************************************
921 : !> \brief obtain integrals (a,b,fb) by symmetry relations from (a,b,fa) if basis sets at a and
922 : !> b are of the same kind, i.e. a and b are same atom type
923 : !> \param abbint integral (a,b,fb)
924 : !> \param dabbint derivative of abbint
925 : !> \param abaint integral (a,b,fa)
926 : !> \param dabdaint derivative of abaint
927 : !> \param rab distance vector between A and B
928 : !> \param oba orbital basis at center A
929 : !> \param fba auxiliary basis set at center A
930 : !> \param calculate_ints ...
931 : !> \param calculate_forces ...
932 : ! **************************************************************************************************
933 26725 : SUBROUTINE get_abb_same_kind(abbint, dabbint, abaint, dabdaint, rab, oba, fba, &
934 : calculate_ints, calculate_forces)
935 :
936 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
937 : INTENT(INOUT) :: abbint
938 : REAL(KIND=dp), ALLOCATABLE, &
939 : DIMENSION(:, :, :, :), INTENT(INOUT) :: dabbint
940 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
941 : INTENT(INOUT) :: abaint
942 : REAL(KIND=dp), ALLOCATABLE, &
943 : DIMENSION(:, :, :, :), INTENT(INOUT) :: dabdaint
944 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
945 : TYPE(gto_basis_set_type), POINTER :: oba, fba
946 : LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
947 :
948 : INTEGER :: i, iend, iset, isgfa, ishella, istart, jend, jset, jsgfa, jshella, jstart, kend, &
949 : kset, ksgfa, kshella, kstart, lai, laj, lak, nseta, nsetca, nsgfa, nsgfca, sgfa_end, &
950 : sgfa_start
951 26725 : INTEGER, ALLOCATABLE, DIMENSION(:) :: lai_set, laj_set, lak_set
952 26725 : INTEGER, DIMENSION(:), POINTER :: nsgfa_set, nsgfca_set, nshella, nshellca
953 26725 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfca, la, lca
954 : REAL(KIND=dp) :: dab
955 26725 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_ca
956 :
957 26725 : NULLIFY (nshellca, first_sgfa, first_sgfca, lca, set_radius_a, &
958 26725 : set_radius_ca)
959 :
960 : ! basis ikind
961 26725 : first_sgfa => oba%first_sgf
962 26725 : set_radius_a => oba%set_radius
963 26725 : nseta = oba%nset
964 26725 : nsgfa = oba%nsgf
965 26725 : nsgfa_set => oba%nsgf_set
966 26725 : nshella => oba%nshell
967 26725 : la => oba%l
968 :
969 : ! basis RI
970 26725 : first_sgfca => fba%first_sgf
971 26725 : set_radius_ca => fba%set_radius
972 26725 : nsetca = fba%nset
973 26725 : nshellca => fba%nshell
974 26725 : nsgfca = fba%nsgf
975 26725 : nsgfca_set => fba%nsgf_set
976 26725 : lca => fba%l
977 :
978 80175 : ALLOCATE (lai_set(nsgfa))
979 53450 : ALLOCATE (laj_set(nsgfa))
980 80175 : ALLOCATE (lak_set(nsgfca))
981 :
982 106900 : dab = SQRT(SUM(rab**2))
983 53966 : DO iset = 1, nseta
984 :
985 158748 : DO ishella = 1, nshella(iset)
986 131507 : sgfa_start = first_sgfa(ishella, iset)
987 131507 : sgfa_end = sgfa_start + 2*la(ishella, iset)
988 500795 : lai_set(sgfa_start:sgfa_end) = la(ishella, iset)
989 : END DO
990 27241 : istart = first_sgfa(1, iset)
991 27241 : iend = istart + nsgfa_set(iset) - 1
992 :
993 82239 : DO jset = 1, nseta
994 :
995 28273 : IF (set_radius_a(iset) + set_radius_a(jset) < dab) CYCLE
996 157624 : DO jshella = 1, nshella(jset)
997 130745 : sgfa_start = first_sgfa(jshella, jset)
998 130745 : sgfa_end = sgfa_start + 2*la(jshella, jset)
999 494245 : laj_set(sgfa_start:sgfa_end) = la(jshella, jset)
1000 : END DO
1001 26879 : jstart = first_sgfa(1, jset)
1002 26879 : jend = jstart + nsgfa_set(jset) - 1
1003 :
1004 292644 : DO kset = 1, nsetca
1005 :
1006 238524 : IF (set_radius_a(iset) + set_radius_ca(kset) < dab) CYCLE
1007 :
1008 403696 : DO kshella = 1, nshellca(kset)
1009 288548 : sgfa_start = first_sgfca(kshella, kset)
1010 288548 : sgfa_end = sgfa_start + 2*lca(kshella, kset)
1011 1396394 : lak_set(sgfa_start:sgfa_end) = lca(kshella, kset)
1012 : END DO
1013 115148 : kstart = first_sgfca(1, kset)
1014 115148 : kend = kstart + nsgfca_set(kset) - 1
1015 1136119 : DO ksgfa = kstart, kend
1016 992698 : lak = lak_set(ksgfa)
1017 13744064 : DO jsgfa = jstart, jend
1018 12512842 : laj = laj_set(jsgfa)
1019 174158670 : DO isgfa = istart, iend
1020 160653130 : lai = lai_set(isgfa)
1021 173165972 : IF (MODULO((lai + laj + lak), 2) /= 0) THEN
1022 80305489 : IF (calculate_ints) THEN
1023 : abbint(isgfa, jsgfa, ksgfa) = &
1024 43130519 : -abaint(jsgfa, isgfa, ksgfa)
1025 : END IF
1026 80305489 : IF (calculate_forces) THEN
1027 148699880 : DO i = 1, 3
1028 : dabbint(isgfa, jsgfa, ksgfa, i) = &
1029 148699880 : -dabdaint(jsgfa, isgfa, ksgfa, i)
1030 : END DO
1031 : END IF
1032 : ELSE
1033 80347641 : IF (calculate_ints) THEN
1034 : abbint(isgfa, jsgfa, ksgfa) = &
1035 43149633 : abaint(jsgfa, isgfa, ksgfa)
1036 : END IF
1037 80347641 : IF (calculate_forces) THEN
1038 148792032 : DO i = 1, 3
1039 : dabbint(isgfa, jsgfa, ksgfa, i) = &
1040 148792032 : dabdaint(jsgfa, isgfa, ksgfa, i)
1041 : END DO
1042 : END IF
1043 : END IF
1044 : END DO
1045 : END DO
1046 : END DO
1047 : END DO
1048 : END DO
1049 : END DO
1050 :
1051 26725 : DEALLOCATE (lai_set, laj_set, lak_set)
1052 :
1053 26725 : END SUBROUTINE get_abb_same_kind
1054 :
1055 : ! **************************************************************************************************
1056 : !> \brief calculate integrals (a,b,fa); requires angular-dependent part as input
1057 : !> \param abaint integral (a,b,fa)
1058 : !> \param dabdaint ...
1059 : !> \param rab distance vector between A and B
1060 : !> \param oba orbital basis at center A
1061 : !> \param obb orbital basis at center B
1062 : !> \param fba auxiliary basis set at center A
1063 : !> \param scon_obb contraction matrix for orb bas on B
1064 : !> \param scona_mix mixed contraction matrix orb + ri basis on A
1065 : !> \param oba_index orbital basis index for scona_mix
1066 : !> \param fba_index ri basis index for scona_mix
1067 : !> \param cg_coeff Clebsch-Gordon coefficients
1068 : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
1069 : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
1070 : !> \param Waux_mat W matrix that contains the angular-dependent part
1071 : !> \param dWaux_mat derivative of the W matrix
1072 : !> \param calculate_ints ...
1073 : !> \param calculate_forces ...
1074 : ! **************************************************************************************************
1075 27647 : SUBROUTINE int_overlap_aba_shg_low(abaint, dabdaint, rab, oba, obb, fba, scon_obb, scona_mix, &
1076 27647 : oba_index, fba_index, cg_coeff, cg_none0_list, &
1077 27647 : ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
1078 :
1079 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1080 : INTENT(INOUT) :: abaint
1081 : REAL(KIND=dp), ALLOCATABLE, &
1082 : DIMENSION(:, :, :, :), INTENT(INOUT) :: dabdaint
1083 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
1084 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fba
1085 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_obb
1086 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix
1087 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: oba_index, fba_index
1088 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
1089 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
1090 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
1091 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Waux_mat
1092 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dWaux_mat
1093 : LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
1094 :
1095 : INTEGER :: iset, jset, kset, la_max_set, laa_max, laa_max_set, lb_max_set, lca_max_set, na, &
1096 : nb, nca, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfca_set, nseta, nsetb, &
1097 : nsetca, nshella_set, nshellb_set, nshellca_set, sgfa, sgfb, sgfca
1098 27647 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lca_max, npgfa, npgfb, &
1099 27647 : npgfca, nsgfa, nsgfb, nsgfca, nshella, &
1100 27647 : nshellb, nshellca
1101 27647 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfca, la, &
1102 27647 : lb, lca
1103 : REAL(KIND=dp) :: dab, rab2
1104 : REAL(KIND=dp), ALLOCATABLE, &
1105 : DIMENSION(:, :, :, :, :) :: swork, swork_cont
1106 27647 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_ca
1107 27647 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetca
1108 :
1109 27647 : NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca)
1110 27647 : NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
1111 27647 : set_radius_ca, zeta, zetb, zetca)
1112 :
1113 : ! basis ikind
1114 27647 : first_sgfa => oba%first_sgf
1115 27647 : la_max => oba%lmax
1116 27647 : la => oba%l
1117 27647 : nsgfa => oba%nsgf_set
1118 27647 : npgfa => oba%npgf
1119 27647 : nshella => oba%nshell
1120 27647 : nseta = oba%nset
1121 27647 : set_radius_a => oba%set_radius
1122 27647 : zeta => oba%zet
1123 : ! basis jkind
1124 27647 : first_sgfb => obb%first_sgf
1125 27647 : lb_max => obb%lmax
1126 27647 : lb => obb%l
1127 27647 : nsgfb => obb%nsgf_set
1128 27647 : npgfb => obb%npgf
1129 27647 : nshellb => obb%nshell
1130 27647 : nsetb = obb%nset
1131 27647 : set_radius_b => obb%set_radius
1132 27647 : zetb => obb%zet
1133 :
1134 : ! basis RI A
1135 27647 : first_sgfca => fba%first_sgf
1136 27647 : lca_max => fba%lmax
1137 27647 : lca => fba%l
1138 27647 : nsgfca => fba%nsgf_set
1139 27647 : npgfca => fba%npgf
1140 27647 : nshellca => fba%nshell
1141 27647 : nsetca = fba%nset
1142 27647 : set_radius_ca => fba%set_radius
1143 27647 : zetca => fba%zet
1144 :
1145 110588 : dab = SQRT(SUM(rab**2))
1146 27647 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1147 :
1148 56082 : la_max_set = MAXVAL(la_max)
1149 56082 : lb_max_set = MAXVAL(lb_max)
1150 274583 : lca_max_set = MAXVAL(lca_max)
1151 56082 : npgfa_set = MAXVAL(npgfa)
1152 56082 : npgfb_set = MAXVAL(npgfb)
1153 274583 : npgfca_set = MAXVAL(npgfca)
1154 56082 : nshella_set = MAXVAL(nshella)
1155 56082 : nshellb_set = MAXVAL(nshellb)
1156 274583 : nshellca_set = MAXVAL(nshellca)
1157 : !*** for forces: derivative+1 in auxiliary vector required
1158 27647 : ndev = 0
1159 27647 : IF (calculate_forces) ndev = 1
1160 :
1161 27647 : laa_max_set = la_max_set + lca_max_set
1162 :
1163 : ! allocate some work storage....
1164 27647 : nds_max = laa_max_set + lb_max_set + ndev + 1
1165 27647 : nl_set = INT((laa_max_set)/2)
1166 193529 : ALLOCATE (swork(npgfb_set, npgfa_set, npgfca_set, nl_set + 1, nds_max))
1167 193529 : ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellca_set))
1168 :
1169 56082 : DO iset = 1, nseta
1170 :
1171 86189 : DO jset = 1, nsetb
1172 :
1173 30107 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
1174 :
1175 306534 : DO kset = 1, nsetca
1176 :
1177 249938 : IF (set_radius_b(jset) + set_radius_ca(kset) < dab) CYCLE
1178 :
1179 : !*** calculate s_baa here
1180 123750 : laa_max = la_max(iset) + lca_max(kset)
1181 123750 : nds = laa_max + lb_max(jset) + ndev + 1
1182 123750 : nl = INT(laa_max/2) + 1
1183 97515203 : swork(1:npgfb(jset), 1:npgfa(iset), 1:npgfca(kset), 1:nl, 1:nds) = 0.0_dp
1184 : CALL s_overlap_abb(lb_max(jset), npgfb(jset), zetb(:, jset), &
1185 : la_max(iset), npgfa(iset), zeta(:, iset), &
1186 : lca_max(kset), npgfca(kset), zetca(:, kset), &
1187 123750 : rab, swork, calculate_forces)
1188 :
1189 : CALL contract_s_overlap_aba(la(:, iset), npgfa(iset), nshella(iset), &
1190 : lb(:, jset), npgfb(jset), nshellb(jset), &
1191 : scon_obb(1:npgfb(jset), 1:nshellb(jset), jset), &
1192 : lca(:, kset), npgfca(kset), nshellca(kset), &
1193 : oba_index(:, :, iset), fba_index(:, :, kset), scona_mix, nl, nds, &
1194 : swork(1:npgfb(jset), 1:npgfa(iset), 1:npgfca(kset), 1:nl, 1:nds), &
1195 123750 : swork_cont, calculate_forces)
1196 123750 : IF (calculate_ints) THEN
1197 : CALL construct_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
1198 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
1199 : lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
1200 : cg_coeff, cg_none0_list, ncg_none0, &
1201 70378 : swork_cont, Waux_mat, abaint)
1202 : END IF
1203 123750 : IF (calculate_forces) THEN
1204 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
1205 : CALL dev_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
1206 : lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
1207 : lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
1208 : cg_coeff, cg_none0_list, ncg_none0, &
1209 214448 : -rab, swork_cont, Waux_mat, dWaux_mat, dabdaint)
1210 : END IF
1211 : ! max value of integrals in this set triple
1212 123750 : sgfa = first_sgfa(1, iset)
1213 123750 : na = sgfa + nsgfa(iset) - 1
1214 123750 : sgfb = first_sgfb(1, jset)
1215 123750 : nb = sgfb + nsgfb(jset) - 1
1216 123750 : sgfca = first_sgfca(1, kset)
1217 280045 : nca = sgfca + nsgfca(kset) - 1
1218 :
1219 : END DO
1220 : END DO
1221 : END DO
1222 :
1223 27647 : DEALLOCATE (swork_cont)
1224 27647 : DEALLOCATE (swork)
1225 :
1226 27647 : END SUBROUTINE int_overlap_aba_shg_low
1227 :
1228 : ! **************************************************************************************************
1229 : !> \brief precalculates the angular part of the SHG integrals for the matrices
1230 : !> (fa,fb), (a,b), (a,b,fa) and (b,fb,a); the same Waux_mat can then be used for all
1231 : !> for integrals; specific for LRIGPW
1232 : !> \param oba orbital basis on a
1233 : !> \param obb orbital basis on b
1234 : !> \param fba aux basis on a
1235 : !> \param fbb aux basis on b
1236 : !> \param rab distance vector between a and b
1237 : !> \param Waux_mat W matrix that contains the angular-dependent part
1238 : !> \param dWaux_mat derivative of the W matrix
1239 : !> \param calculate_forces ...
1240 : ! **************************************************************************************************
1241 27469 : SUBROUTINE lri_precalc_angular_shg_part(oba, obb, fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces)
1242 :
1243 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fba, fbb
1244 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
1245 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1246 : INTENT(OUT) :: Waux_mat
1247 : REAL(KIND=dp), ALLOCATABLE, &
1248 : DIMENSION(:, :, :, :), INTENT(OUT) :: dWaux_mat
1249 : LOGICAL, INTENT(IN) :: calculate_forces
1250 :
1251 : INTEGER :: i, isize, j, k, la_max, laa_max, lb_max, &
1252 : lbb_max, lca_max, lcb_max, li_max, &
1253 : lj_max, lmax, mdim(3), size_int(4, 2), &
1254 : temp
1255 27469 : INTEGER, DIMENSION(:), POINTER :: li_max_all
1256 : REAL(KIND=dp) :: rab2
1257 27469 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Rc, Rs
1258 :
1259 27469 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1260 :
1261 : !*** 1 Waux_mat of size (li_max,lj_max) for elements
1262 : ! i j
1263 : ! [aab] --> (laa_max, lb_max)
1264 : ! [bba] --> (lbb_max, la_max) --> use for [abb]
1265 : ! [ab] ri --> (lca_max, lcb_max)
1266 : ! [ab] orb --> (la_max , lb_max)
1267 :
1268 55682 : la_max = MAXVAL(oba%lmax)
1269 55682 : lb_max = MAXVAL(obb%lmax)
1270 272195 : lca_max = MAXVAL(fba%lmax)
1271 272195 : lcb_max = MAXVAL(fbb%lmax)
1272 :
1273 27469 : laa_max = la_max + lca_max
1274 27469 : lbb_max = lb_max + lcb_max
1275 27469 : li_max = MAX(laa_max, lbb_max)
1276 27469 : lj_max = MAX(la_max, lb_max, lcb_max)
1277 27469 : lmax = li_max
1278 :
1279 82407 : ALLOCATE (li_max_all(0:lj_max))
1280 164814 : ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
1281 4219853 : Rc = 0._dp
1282 4219853 : Rs = 0._dp
1283 27469 : mdim(1) = li_max + lj_max + 1
1284 27469 : mdim(2) = nsoset(li_max) + 1
1285 27469 : mdim(3) = nsoset(lj_max) + 1
1286 137345 : ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
1287 137345 : ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3)))
1288 : !Waux_mat = 0._dp !.. takes time
1289 : !dWaux_mat =0._dp !.. takes time
1290 :
1291 : !*** Waux_mat (li_max,lj_max) contains elements not needed,
1292 : !*** make indixing so that only required ones are computed
1293 : !*** li_max_all(j) --> li_max dependent on j
1294 82407 : size_int(1, :) = [laa_max, lb_max]
1295 82407 : size_int(2, :) = [lbb_max, la_max]
1296 82407 : size_int(3, :) = [lca_max, lcb_max]
1297 82407 : size_int(4, :) = [la_max, lb_max]
1298 :
1299 139127 : li_max_all(:) = 0
1300 137345 : DO isize = 1, 4
1301 109876 : i = size_int(isize, 1)
1302 109876 : j = size_int(isize, 2)
1303 109876 : k = li_max_all(j)
1304 137345 : IF (k < i) li_max_all(j) = i
1305 : END DO
1306 27469 : temp = li_max_all(lj_max)
1307 139127 : DO j = lj_max, 0, -1
1308 139127 : IF (li_max_all(j) < temp) THEN
1309 56040 : li_max_all(j) = temp
1310 : ELSE
1311 : temp = li_max_all(j)
1312 : END IF
1313 : END DO
1314 :
1315 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
1316 109876 : CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
1317 27469 : CALL get_W_matrix(li_max_all, lj_max, lmax, Rc, Rs, Waux_mat)
1318 27469 : IF (calculate_forces) THEN
1319 11948 : CALL get_dW_matrix(li_max_all, lj_max, Waux_mat, dWaux_mat)
1320 : END IF
1321 :
1322 27469 : DEALLOCATE (Rc, Rs, li_max_all)
1323 :
1324 27469 : END SUBROUTINE lri_precalc_angular_shg_part
1325 :
1326 : END MODULE generic_shg_integrals
|
