Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 905220 : 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 905220 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
114 905220 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
115 :
116 : PROCEDURE(ab_sint_shg), POINTER :: s_operator_ab
117 :
118 905220 : NULLIFY (s_operator_ab)
119 :
120 3684604 : la_max = MAXVAL(fba%lmax)
121 3685404 : lb_max = MAXVAL(fbb%lmax)
122 :
123 : CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
124 905220 : my_omega = 1.0_dp
125 :
126 1810312 : SELECT CASE (r12_operator)
127 : CASE (operator_coulomb)
128 905092 : 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 905220 : 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 1810280 : my_omega, Waux_mat, dWaux_mat, calculate_forces)
147 :
148 905220 : DEALLOCATE (Waux_mat, dWaux_mat)
149 :
150 905220 : 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 144 : SUBROUTINE int_overlap_aba_shg(saba, dsaba, rab, oba, obb, fba, scon_obb, &
248 144 : scona_mix, oba_index, fba_index, &
249 144 : 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 144 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
268 144 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
269 :
270 1948 : laa_max = MAXVAL(oba%lmax) + MAXVAL(fba%lmax)
271 332 : lb_max = MAXVAL(obb%lmax)
272 :
273 1884366 : saba = 0.0_dp
274 951984 : 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 144 : Waux_mat, dWaux_mat, .TRUE., calculate_forces)
281 :
282 144 : DEALLOCATE (Waux_mat, dWaux_mat)
283 :
284 144 : 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 905444 : 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 905444 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Rc, Rs
365 :
366 905444 : NULLIFY (la_max_all)
367 905444 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
368 905444 : lmax = MAX(la_max, lb_max)
369 :
370 2716332 : ALLOCATE (la_max_all(0:lb_max))
371 6338108 : ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
372 20272258 : Rc = 0._dp
373 20272258 : Rs = 0._dp
374 905444 : mdim(1) = MIN(la_max, lb_max) + 1
375 905444 : mdim(2) = nsoset(la_max) + 1
376 905444 : mdim(3) = nsoset(lb_max) + 1
377 4527220 : ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
378 4527220 : ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3)))
379 :
380 2626414 : la_max_all(0:lb_max) = la_max
381 : !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
382 3621776 : CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
383 905444 : CALL get_W_matrix(la_max_all, lb_max, lmax, Rc, Rs, Waux_mat)
384 905444 : IF (calculate_forces) THEN
385 256 : CALL get_dW_matrix(la_max_all, lb_max, Waux_mat, dWaux_mat)
386 : END IF
387 :
388 905444 : DEALLOCATE (Rc, Rs, la_max_all)
389 :
390 905444 : 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 905220 : SUBROUTINE int_operator_ab_shg_low(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
410 905220 : 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 905220 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
426 905220 : nsgfb, nshella, nshellb
427 905220 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb
428 : REAL(KIND=dp) :: dab
429 905220 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
430 905220 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
431 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: swork, swork_cont
432 :
433 905220 : NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
434 905220 : set_radius_b, zeta, zetb)
435 :
436 : ! basis ikind
437 905220 : first_sgfa => fba%first_sgf
438 905220 : la_max => fba%lmax
439 905220 : la => fba%l
440 905220 : npgfa => fba%npgf
441 905220 : nsgfa => fba%nsgf_set
442 905220 : nseta = fba%nset
443 905220 : set_radius_a => fba%set_radius
444 905220 : zeta => fba%zet
445 905220 : nshella => fba%nshell
446 : ! basis jkind
447 905220 : first_sgfb => fbb%first_sgf
448 905220 : lb_max => fbb%lmax
449 905220 : lb => fbb%l
450 905220 : npgfb => fbb%npgf
451 905220 : nsgfb => fbb%nsgf_set
452 905220 : nsetb = fbb%nset
453 905220 : set_radius_b => fbb%set_radius
454 905220 : zetb => fbb%zet
455 905220 : nshellb => fbb%nshell
456 :
457 3620880 : dab = SQRT(SUM(rab**2))
458 :
459 3684604 : la_max_set = MAXVAL(la_max)
460 3685404 : lb_max_set = MAXVAL(lb_max)
461 :
462 : ! allocate some work matrices
463 3684604 : npgfa_set = MAXVAL(npgfa)
464 3685404 : npgfb_set = MAXVAL(npgfb)
465 3684604 : nshella_set = MAXVAL(nshella)
466 3685404 : nshellb_set = MAXVAL(nshellb)
467 3684604 : nsgfa_set = MAXVAL(nsgfa)
468 3685404 : nsgfb_set = MAXVAL(nsgfb)
469 905220 : ndev = 0
470 905220 : IF (calculate_forces) ndev = 1
471 905220 : nds_max = la_max_set + lb_max_set + ndev + 1
472 4526100 : ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
473 4526100 : ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
474 :
475 148055868 : vab = 0.0_dp
476 15363460 : IF (calculate_forces) dvab = 0.0_dp
477 :
478 3684604 : DO iset = 1, nseta
479 :
480 21440796 : DO jset = 1, nsetb
481 :
482 17756192 : nds = la_max(iset) + lb_max(jset) + ndev + 1
483 164330720 : 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 17756192 : 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 17756192 : 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 17756192 : swork_cont, Waux_mat, vab)
496 20535576 : 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 905220 : DEALLOCATE (swork, swork_cont)
506 :
507 905220 : 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 16834 : 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 16834 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nshella, &
543 16834 : nshellb
544 16834 : 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 16834 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
549 16834 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
550 :
551 16834 : NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
552 16834 : set_radius_b, zeta, zetb)
553 :
554 : ! basis ikind
555 16834 : first_sgfa => fba%first_sgf
556 16834 : la_max => fba%lmax
557 16834 : la => fba%l
558 16834 : npgfa => fba%npgf
559 16834 : nseta = fba%nset
560 16834 : set_radius_a => fba%set_radius
561 16834 : zeta => fba%zet
562 16834 : nshella => fba%nshell
563 : ! basis jkind
564 16834 : first_sgfb => fbb%first_sgf
565 16834 : lb_max => fbb%lmax
566 16834 : lb => fbb%l
567 16834 : npgfb => fbb%npgf
568 16834 : nsetb = fbb%nset
569 16834 : set_radius_b => fbb%set_radius
570 16834 : zetb => fbb%zet
571 16834 : nshellb => fbb%nshell
572 :
573 67336 : dab = SQRT(SUM(rab**2))
574 :
575 46529 : la_max_set = MAXVAL(la_max)
576 46689 : lb_max_set = MAXVAL(lb_max)
577 :
578 : ! allocate some work matrices
579 46529 : npgfa_set = MAXVAL(npgfa)
580 46689 : npgfb_set = MAXVAL(npgfb)
581 46529 : nshella_set = MAXVAL(nshella)
582 46689 : nshellb_set = MAXVAL(nshellb)
583 16834 : ndev = 0
584 16834 : IF (calculate_forces) ndev = 1
585 16834 : nds_max = la_max_set + lb_max_set + ndev + 1
586 84170 : ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
587 84170 : ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
588 :
589 10392820 : IF (calculate_ints) sab = 0.0_dp
590 16449973 : IF (calculate_forces) dsab = 0.0_dp
591 :
592 16834 : my_contraction_high = .TRUE.
593 16834 : IF (PRESENT(contraction_high)) my_contraction_high = contraction_high
594 :
595 46529 : DO iset = 1, nseta
596 :
597 193520 : DO jset = 1, nsetb
598 :
599 146991 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
600 :
601 120606 : nds = la_max(iset) + lb_max(jset) + ndev + 1
602 4407240 : 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 120606 : rab, swork, calculate_forces)
606 120606 : 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 19819 : 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 100787 : swork, swork_cont, calculate_forces)
619 : END IF
620 120606 : 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 76964 : swork_cont, Waux_mat, sab)
624 : END IF
625 150301 : 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 193768 : -rab, swork_cont, Waux_mat, dWaux_mat, dsab)
630 : END IF
631 : END DO
632 : END DO
633 :
634 16834 : DEALLOCATE (swork, swork_cont)
635 :
636 16834 : 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 744 : SUBROUTINE int_overlap_abb_shg_low(abbint, dabbint, rab, oba, obb, fbb, scon_oba, sconb_mix, &
770 744 : obb_index, fbb_index, cg_coeff, cg_none0_list, &
771 744 : 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 744 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lcb_max, npgfa, npgfb, &
793 744 : npgfcb, nsgfa, nsgfb, nsgfcb, nshella, &
794 744 : nshellb, nshellcb
795 744 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfcb, la, &
796 744 : lb, lcb
797 : REAL(KIND=dp) :: dab, rab2
798 : REAL(KIND=dp), ALLOCATABLE, &
799 : DIMENSION(:, :, :, :, :) :: swork, swork_cont
800 744 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_cb
801 744 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetcb
802 :
803 744 : NULLIFY (la_max, lb_max, lcb_max, npgfa, npgfb, npgfcb)
804 744 : NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
805 744 : set_radius_cb, zeta, zetb, zetcb)
806 :
807 : ! basis ikind
808 744 : first_sgfa => oba%first_sgf
809 744 : la_max => oba%lmax
810 744 : la => oba%l
811 744 : nsgfa => oba%nsgf_set
812 744 : npgfa => oba%npgf
813 744 : nshella => oba%nshell
814 744 : nseta = oba%nset
815 744 : set_radius_a => oba%set_radius
816 744 : zeta => oba%zet
817 : ! basis jkind
818 744 : first_sgfb => obb%first_sgf
819 744 : lb_max => obb%lmax
820 744 : lb => obb%l
821 744 : nsgfb => obb%nsgf_set
822 744 : npgfb => obb%npgf
823 744 : nshellb => obb%nshell
824 744 : nsetb = obb%nset
825 744 : set_radius_b => obb%set_radius
826 744 : zetb => obb%zet
827 :
828 : ! basis RI on B
829 744 : first_sgfcb => fbb%first_sgf
830 744 : lcb_max => fbb%lmax
831 744 : lcb => fbb%l
832 744 : nsgfcb => fbb%nsgf_set
833 744 : npgfcb => fbb%npgf
834 744 : nshellcb => fbb%nshell
835 744 : nsetcb = fbb%nset
836 744 : set_radius_cb => fbb%set_radius
837 744 : zetcb => fbb%zet
838 :
839 2976 : dab = SQRT(SUM(rab**2))
840 744 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
841 :
842 1716 : la_max_set = MAXVAL(la_max)
843 1716 : lb_max_set = MAXVAL(lb_max)
844 7242 : lcb_max_set = MAXVAL(lcb_max)
845 1716 : npgfa_set = MAXVAL(npgfa)
846 1716 : npgfb_set = MAXVAL(npgfb)
847 7242 : npgfcb_set = MAXVAL(npgfcb)
848 1716 : nshella_set = MAXVAL(nshella)
849 1716 : nshellb_set = MAXVAL(nshellb)
850 7242 : nshellcb_set = MAXVAL(nshellcb)
851 : !*** for forces: derivative+1 in auxiliary vector required
852 744 : ndev = 0
853 744 : IF (calculate_forces) ndev = 1
854 :
855 744 : lbb_max_set = lb_max_set + lcb_max_set
856 :
857 : ! allocate some work storage....
858 744 : nds_max = la_max_set + lbb_max_set + ndev + 1
859 744 : nl_set = INT((lbb_max_set)/2)
860 5208 : ALLOCATE (swork(npgfa_set, npgfb_set, npgfcb_set, nl_set + 1, nds_max))
861 5208 : ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellcb_set))
862 :
863 1716 : DO iset = 1, nseta
864 :
865 3144 : DO jset = 1, nsetb
866 :
867 1428 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
868 :
869 8664 : DO kset = 1, nsetcb
870 :
871 6816 : IF (set_radius_a(iset) + set_radius_cb(kset) < dab) CYCLE
872 :
873 4004 : lbb_max = lb_max(jset) + lcb_max(kset)
874 4004 : nds = la_max(iset) + lbb_max + ndev + 1
875 4004 : nl = INT((lbb_max)/2) + 1
876 4749072 : 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 4004 : 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 4004 : swork_cont, calculate_forces)
889 :
890 4004 : 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 3277 : ncg_none0, swork_cont, Waux_mat, abbint)
896 : END IF
897 4004 : 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 3868 : Waux_mat, dWaux_mat, dabbint)
904 : END IF
905 : ! max value of integrals in this set triple
906 4004 : sgfa = first_sgfa(1, iset)
907 4004 : na = sgfa + nsgfa(iset) - 1
908 4004 : sgfb = first_sgfb(1, jset)
909 4004 : nb = sgfb + nsgfb(jset) - 1
910 4004 : sgfcb = first_sgfcb(1, kset)
911 8244 : ncb = sgfcb + nsgfcb(kset) - 1
912 : END DO
913 : END DO
914 : END DO
915 :
916 744 : DEALLOCATE (swork_cont)
917 744 : DEALLOCATE (swork)
918 :
919 744 : 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 14005 : 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 14005 : INTEGER, ALLOCATABLE, DIMENSION(:) :: lai_set, laj_set, lak_set
952 14005 : INTEGER, DIMENSION(:), POINTER :: nsgfa_set, nsgfca_set, nshella, nshellca
953 14005 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfca, la, lca
954 : REAL(KIND=dp) :: dab
955 14005 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_ca
956 :
957 14005 : NULLIFY (nshellca, first_sgfa, first_sgfca, lca, set_radius_a, &
958 14005 : set_radius_ca)
959 :
960 : ! basis ikind
961 14005 : first_sgfa => oba%first_sgf
962 14005 : set_radius_a => oba%set_radius
963 14005 : nseta = oba%nset
964 14005 : nsgfa = oba%nsgf
965 14005 : nsgfa_set => oba%nsgf_set
966 14005 : nshella => oba%nshell
967 14005 : la => oba%l
968 :
969 : ! basis RI
970 14005 : first_sgfca => fba%first_sgf
971 14005 : set_radius_ca => fba%set_radius
972 14005 : nsetca = fba%nset
973 14005 : nshellca => fba%nshell
974 14005 : nsgfca = fba%nsgf
975 14005 : nsgfca_set => fba%nsgf_set
976 14005 : lca => fba%l
977 :
978 42015 : ALLOCATE (lai_set(nsgfa))
979 42015 : ALLOCATE (laj_set(nsgfa))
980 42015 : ALLOCATE (lak_set(nsgfca))
981 :
982 56020 : dab = SQRT(SUM(rab**2))
983 28526 : DO iset = 1, nseta
984 :
985 82460 : DO ishella = 1, nshella(iset)
986 67939 : sgfa_start = first_sgfa(ishella, iset)
987 67939 : sgfa_end = sgfa_start + 2*la(ishella, iset)
988 259275 : lai_set(sgfa_start:sgfa_end) = la(ishella, iset)
989 : END DO
990 14521 : istart = first_sgfa(1, iset)
991 14521 : iend = istart + nsgfa_set(iset) - 1
992 :
993 44079 : DO jset = 1, nseta
994 :
995 15553 : IF (set_radius_a(iset) + set_radius_a(jset) < dab) CYCLE
996 81336 : DO jshella = 1, nshella(jset)
997 67177 : sgfa_start = first_sgfa(jshella, jset)
998 67177 : sgfa_end = sgfa_start + 2*la(jshella, jset)
999 252725 : laj_set(sgfa_start:sgfa_end) = la(jshella, jset)
1000 : END DO
1001 14159 : jstart = first_sgfa(1, jset)
1002 14159 : jend = jstart + nsgfa_set(jset) - 1
1003 :
1004 152516 : DO kset = 1, nsetca
1005 :
1006 123836 : IF (set_radius_a(iset) + set_radius_ca(kset) < dab) CYCLE
1007 :
1008 210208 : DO kshella = 1, nshellca(kset)
1009 150740 : sgfa_start = first_sgfca(kshella, kset)
1010 150740 : sgfa_end = sgfa_start + 2*lca(kshella, kset)
1011 732074 : lak_set(sgfa_start:sgfa_end) = lca(kshella, kset)
1012 : END DO
1013 59468 : kstart = first_sgfca(1, kset)
1014 59468 : kend = kstart + nsgfca_set(kset) - 1
1015 596887 : DO ksgfa = kstart, kend
1016 521866 : lak = lak_set(ksgfa)
1017 7050400 : DO jsgfa = jstart, jend
1018 6404698 : laj = laj_set(jsgfa)
1019 88237182 : DO isgfa = istart, iend
1020 81310618 : lai = lai_set(isgfa)
1021 87715316 : IF (MODULO((lai + laj + lak), 2) /= 0) THEN
1022 40645729 : IF (calculate_ints) THEN
1023 : abbint(isgfa, jsgfa, ksgfa) = &
1024 22047455 : -abaint(jsgfa, isgfa, ksgfa)
1025 : END IF
1026 40645729 : IF (calculate_forces) THEN
1027 74393096 : DO i = 1, 3
1028 : dabbint(isgfa, jsgfa, ksgfa, i) = &
1029 74393096 : -dabdaint(jsgfa, isgfa, ksgfa, i)
1030 : END DO
1031 : END IF
1032 : ELSE
1033 40664889 : IF (calculate_ints) THEN
1034 : abbint(isgfa, jsgfa, ksgfa) = &
1035 22054833 : abaint(jsgfa, isgfa, ksgfa)
1036 : END IF
1037 40664889 : IF (calculate_forces) THEN
1038 74440224 : DO i = 1, 3
1039 : dabbint(isgfa, jsgfa, ksgfa, i) = &
1040 74440224 : 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 14005 : DEALLOCATE (lai_set, laj_set, lak_set)
1052 :
1053 14005 : 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 14877 : SUBROUTINE int_overlap_aba_shg_low(abaint, dabdaint, rab, oba, obb, fba, scon_obb, scona_mix, &
1076 14877 : oba_index, fba_index, cg_coeff, cg_none0_list, &
1077 14877 : 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 14877 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lca_max, npgfa, npgfb, &
1099 14877 : npgfca, nsgfa, nsgfb, nsgfca, nshella, &
1100 14877 : nshellb, nshellca
1101 14877 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfca, la, &
1102 14877 : lb, lca
1103 : REAL(KIND=dp) :: dab, rab2
1104 : REAL(KIND=dp), ALLOCATABLE, &
1105 : DIMENSION(:, :, :, :, :) :: swork, swork_cont
1106 14877 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_ca
1107 14877 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetca
1108 :
1109 14877 : NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca)
1110 14877 : NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
1111 14877 : set_radius_ca, zeta, zetb, zetca)
1112 :
1113 : ! basis ikind
1114 14877 : first_sgfa => oba%first_sgf
1115 14877 : la_max => oba%lmax
1116 14877 : la => oba%l
1117 14877 : nsgfa => oba%nsgf_set
1118 14877 : npgfa => oba%npgf
1119 14877 : nshella => oba%nshell
1120 14877 : nseta = oba%nset
1121 14877 : set_radius_a => oba%set_radius
1122 14877 : zeta => oba%zet
1123 : ! basis jkind
1124 14877 : first_sgfb => obb%first_sgf
1125 14877 : lb_max => obb%lmax
1126 14877 : lb => obb%l
1127 14877 : nsgfb => obb%nsgf_set
1128 14877 : npgfb => obb%npgf
1129 14877 : nshellb => obb%nshell
1130 14877 : nsetb = obb%nset
1131 14877 : set_radius_b => obb%set_radius
1132 14877 : zetb => obb%zet
1133 :
1134 : ! basis RI A
1135 14877 : first_sgfca => fba%first_sgf
1136 14877 : lca_max => fba%lmax
1137 14877 : lca => fba%l
1138 14877 : nsgfca => fba%nsgf_set
1139 14877 : npgfca => fba%npgf
1140 14877 : nshellca => fba%nshell
1141 14877 : nsetca = fba%nset
1142 14877 : set_radius_ca => fba%set_radius
1143 14877 : zetca => fba%zet
1144 :
1145 59508 : dab = SQRT(SUM(rab**2))
1146 14877 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1147 :
1148 30542 : la_max_set = MAXVAL(la_max)
1149 30542 : lb_max_set = MAXVAL(lb_max)
1150 146243 : lca_max_set = MAXVAL(lca_max)
1151 30542 : npgfa_set = MAXVAL(npgfa)
1152 30542 : npgfb_set = MAXVAL(npgfb)
1153 146243 : npgfca_set = MAXVAL(npgfca)
1154 30542 : nshella_set = MAXVAL(nshella)
1155 30542 : nshellb_set = MAXVAL(nshellb)
1156 146243 : nshellca_set = MAXVAL(nshellca)
1157 : !*** for forces: derivative+1 in auxiliary vector required
1158 14877 : ndev = 0
1159 14877 : IF (calculate_forces) ndev = 1
1160 :
1161 14877 : laa_max_set = la_max_set + lca_max_set
1162 :
1163 : ! allocate some work storage....
1164 14877 : nds_max = laa_max_set + lb_max_set + ndev + 1
1165 14877 : nl_set = INT((laa_max_set)/2)
1166 104139 : ALLOCATE (swork(npgfb_set, npgfa_set, npgfca_set, nl_set + 1, nds_max))
1167 104139 : ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellca_set))
1168 :
1169 30542 : DO iset = 1, nseta
1170 :
1171 47879 : DO jset = 1, nsetb
1172 :
1173 17337 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
1174 :
1175 165424 : DO kset = 1, nsetca
1176 :
1177 134368 : IF (set_radius_b(jset) + set_radius_ca(kset) < dab) CYCLE
1178 :
1179 : !*** calculate s_baa here
1180 67188 : laa_max = la_max(iset) + lca_max(kset)
1181 67188 : nds = laa_max + lb_max(jset) + ndev + 1
1182 67188 : nl = INT(laa_max/2) + 1
1183 55180887 : 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 67188 : 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 67188 : swork_cont, calculate_forces)
1196 67188 : 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 40040 : swork_cont, Waux_mat, abaint)
1202 : END IF
1203 67188 : 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 109552 : -rab, swork_cont, Waux_mat, dWaux_mat, dabdaint)
1210 : END IF
1211 : ! max value of integrals in this set triple
1212 67188 : sgfa = first_sgfa(1, iset)
1213 67188 : na = sgfa + nsgfa(iset) - 1
1214 67188 : sgfb = first_sgfb(1, jset)
1215 67188 : nb = sgfb + nsgfb(jset) - 1
1216 67188 : sgfca = first_sgfca(1, kset)
1217 151705 : nca = sgfca + nsgfca(kset) - 1
1218 :
1219 : END DO
1220 : END DO
1221 : END DO
1222 :
1223 14877 : DEALLOCATE (swork_cont)
1224 14877 : DEALLOCATE (swork)
1225 :
1226 14877 : 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 14733 : 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 14733 : INTEGER, DIMENSION(:), POINTER :: li_max_all
1256 : REAL(KIND=dp) :: rab2
1257 14733 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Rc, Rs
1258 :
1259 14733 : 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 30210 : la_max = MAXVAL(oba%lmax)
1269 30210 : lb_max = MAXVAL(obb%lmax)
1270 144483 : lca_max = MAXVAL(fba%lmax)
1271 144483 : lcb_max = MAXVAL(fbb%lmax)
1272 :
1273 14733 : laa_max = la_max + lca_max
1274 14733 : lbb_max = lb_max + lcb_max
1275 14733 : li_max = MAX(laa_max, lbb_max)
1276 14733 : lj_max = MAX(la_max, lb_max, lcb_max)
1277 14733 : lmax = li_max
1278 :
1279 44199 : ALLOCATE (li_max_all(0:lj_max))
1280 103131 : ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
1281 2336237 : Rc = 0._dp
1282 2336237 : Rs = 0._dp
1283 14733 : mdim(1) = li_max + lj_max + 1
1284 14733 : mdim(2) = nsoset(li_max) + 1
1285 14733 : mdim(3) = nsoset(lj_max) + 1
1286 73665 : ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
1287 73665 : 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 44199 : size_int(1, :) = (/laa_max, lb_max/)
1295 44199 : size_int(2, :) = (/lbb_max, la_max/)
1296 44199 : size_int(3, :) = (/lca_max, lcb_max/)
1297 44199 : size_int(4, :) = (/la_max, lb_max/)
1298 :
1299 75471 : li_max_all(:) = 0
1300 73665 : DO isize = 1, 4
1301 58932 : i = size_int(isize, 1)
1302 58932 : j = size_int(isize, 2)
1303 58932 : k = li_max_all(j)
1304 73665 : IF (k < i) li_max_all(j) = i
1305 : END DO
1306 14733 : temp = li_max_all(lj_max)
1307 75471 : DO j = lj_max, 0, -1
1308 75471 : IF (li_max_all(j) < temp) THEN
1309 30600 : 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 58932 : CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
1317 14733 : CALL get_W_matrix(li_max_all, lj_max, lmax, Rc, Rs, Waux_mat)
1318 14733 : IF (calculate_forces) THEN
1319 6140 : CALL get_dW_matrix(li_max_all, lj_max, Waux_mat, dWaux_mat)
1320 : END IF
1321 :
1322 14733 : DEALLOCATE (Rc, Rs, li_max_all)
1323 :
1324 14733 : END SUBROUTINE lri_precalc_angular_shg_part
1325 :
1326 : END MODULE generic_shg_integrals
|
