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 Routines for calculating the s-integrals and their scalar derivatives with respect to rab2
10 : !> over solid harmonic Gaussian (SHG) functions + contraction routines for these integrals
11 : !> i) (s|O(r12)|s) where O(r12) is the overlap, coulomb operator etc.
12 : !> ii) (aba) and (abb) s-overlaps
13 : !> \par Literature (partly)
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 : !> R. Ahlrichs, PCCP, 8, 3072 (2006)
18 : !> \par History
19 : !> created [05.2016]
20 : !> \author Dorothea Golze
21 : ! **************************************************************************************************
22 : MODULE s_contract_shg
23 : USE gamma, ONLY: fgamma => fgamma_0
24 : USE kinds, ONLY: dp
25 : USE mathconstants, ONLY: dfac,&
26 : fac,&
27 : pi
28 : #include "../base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 :
34 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 's_contract_shg'
35 :
36 : ! *** Public subroutines ***
37 : PUBLIC :: s_overlap_ab, s_overlap_abb, s_coulomb_ab, s_verf_ab, s_verfc_ab, &
38 : s_vgauss_ab, s_gauss_ab, s_ra2m_ab
39 :
40 : PUBLIC :: contract_sint_ab_clow, contract_sint_ab_chigh, contract_s_overlap_aba, &
41 : contract_s_overlap_abb, contract_s_ra2m_ab
42 :
43 : CONTAINS
44 :
45 : ! **************************************************************************************************
46 : !> \brief calculates the uncontracted, not normalized [s|s] overlap
47 : !> \param la_max maximal l quantum number on a
48 : !> \param npgfa number of primitive Gaussian on a
49 : !> \param zeta set of exponents on a
50 : !> \param lb_max maximal l quantum number on b
51 : !> \param npgfb number of primitive Gaussian on a
52 : !> \param zetb set of exponents on a
53 : !> \param rab distance vector between a and b
54 : !> \param s uncontracted overlap of s functions
55 : !> \param calculate_forces ...
56 : ! **************************************************************************************************
57 120606 : SUBROUTINE s_overlap_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, rab, s, calculate_forces)
58 :
59 : INTEGER, INTENT(IN) :: la_max, npgfa
60 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
61 : INTEGER, INTENT(IN) :: lb_max, npgfb
62 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
63 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
64 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
65 : LOGICAL, INTENT(IN) :: calculate_forces
66 :
67 : INTEGER :: ids, ipgfa, jpgfb, ndev
68 : REAL(KIND=dp) :: a, b, rab2, xhi, zet
69 :
70 : ! Distance of the centers a and b
71 120606 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
72 120606 : ndev = 0
73 120606 : IF (calculate_forces) ndev = 1
74 : ! Loops over all pairs of primitive Gaussian-type functions
75 305678 : DO ipgfa = 1, npgfa
76 833702 : DO jpgfb = 1, npgfb
77 :
78 : ! Distance Screening !maybe later
79 :
80 : ! Calculate some prefactors
81 528024 : a = zeta(ipgfa)
82 528024 : b = zetb(jpgfb)
83 528024 : zet = a + b
84 528024 : xhi = a*b/zet
85 :
86 : ! [s|s] integral
87 528024 : s(ipgfa, jpgfb, 1) = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
88 :
89 2810700 : DO ids = 2, la_max + lb_max + ndev + 1
90 2625628 : s(ipgfa, jpgfb, ids) = -xhi*s(ipgfa, jpgfb, ids - 1)
91 : END DO
92 :
93 : END DO
94 : END DO
95 :
96 120606 : END SUBROUTINE s_overlap_ab
97 :
98 : ! **************************************************************************************************
99 : !> \brief calculates [s|ra^n|s] integrals for [aba] and the [s|rb^n|s]
100 : !> integrals for [abb]
101 : !> \param la_max maximal l quantum number on a, orbital basis
102 : !> \param npgfa number of primitive Gaussian on a, orbital basis
103 : !> \param zeta set of exponents on a, orbital basis
104 : !> \param lb_max maximal l quantum number on b, orbital basis
105 : !> \param npgfb number of primitive Gaussian on a, orbital basis
106 : !> \param zetb set of exponents on b, orbital basis
107 : !> \param lcb_max maximal l quantum number of aux basis on b
108 : !> \param npgfcb number of primitive Gaussian on b. aux basis
109 : !> \param zetcb set of exponents on b, aux basis
110 : !> \param rab distance vector between a and b
111 : !> \param s uncontracted [s|r^n|s] integrals
112 : !> \param calculate_forces ...
113 : ! **************************************************************************************************
114 71192 : SUBROUTINE s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, &
115 71192 : rab, s, calculate_forces)
116 :
117 : INTEGER, INTENT(IN) :: la_max, npgfa
118 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
119 : INTEGER, INTENT(IN) :: lb_max, npgfb
120 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
121 : INTEGER, INTENT(IN) :: lcb_max, npgfcb
122 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetcb
123 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
124 : REAL(KIND=dp), DIMENSION(:, :, :, :, :), &
125 : INTENT(INOUT) :: s
126 : LOGICAL, INTENT(IN) :: calculate_forces
127 :
128 : INTEGER :: i, ids, il, ipgfa, j, jpgfb, kpgfb, &
129 : lbb_max, lmax, n, ndev, nds, nfac, nl
130 : REAL(KIND=dp) :: a, b, dfsr_int, exp_rab2, k, pfac, &
131 : prefac, rab2, sqrt_pi3, sqrt_zet, &
132 : sr_int, temp, xhi, zet
133 71192 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dsr_int, dtemp
134 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coeff_srs
135 :
136 : ! Distance of the centers a and b
137 71192 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
138 71192 : ndev = 0
139 71192 : IF (calculate_forces) ndev = 1
140 :
141 71192 : lbb_max = lb_max + lcb_max
142 71192 : nl = INT(lbb_max/2)
143 71192 : IF (lb_max == 0 .OR. lcb_max == 0) nl = 0
144 71192 : lmax = la_max + lbb_max
145 :
146 355960 : ALLOCATE (dtemp(nl + 1), dsr_int(nl + 1))
147 355960 : ALLOCATE (coeff_srs(nl + 1, nl + 1, nl + 1))
148 71192 : IF (nl > 5) CALL get_prefac_sabb(nl, coeff_srs)
149 :
150 71192 : sqrt_pi3 = SQRT(pi**3)
151 :
152 : ! Loops over all pairs of primitive Gaussian-type functions
153 156950 : DO kpgfb = 1, npgfcb
154 622648 : DO jpgfb = 1, npgfb
155 3161930 : DO ipgfa = 1, npgfa
156 :
157 : !Calculate some prefactors
158 2610474 : a = zeta(ipgfa)
159 2610474 : b = zetb(jpgfb) + zetcb(kpgfb)
160 :
161 2610474 : zet = a + b
162 2610474 : xhi = a*b/zet
163 2610474 : exp_rab2 = EXP(-xhi*rab2)
164 :
165 2610474 : pfac = a**2/zet
166 2610474 : sqrt_zet = SQRT(zet)
167 :
168 9686688 : DO il = 0, nl
169 6610516 : nds = lmax - 2*il + ndev + 1
170 2610474 : SELECT CASE (il)
171 : CASE (0)
172 : ! [s|s] integral
173 2610474 : s(ipgfa, jpgfb, kpgfb, il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2
174 16569083 : DO ids = 2, nds
175 13958609 : n = ids - 1
176 16569083 : s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, kpgfb, il + 1, 1)
177 : END DO
178 : CASE (1)
179 : ![s|r^2|s] integral
180 2056110 : sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp
181 2056110 : s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
182 2056110 : k = sqrt_pi3*a**2/sqrt_zet**7
183 10717201 : DO ids = 2, nds
184 8661091 : n = ids - 1
185 : s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
186 10717201 : + n*(-xhi)**(n - 1)*k*exp_rab2
187 : END DO
188 : CASE (2)
189 : ![s|r^4|s] integral
190 1793960 : prefac = sqrt_pi3/4.0_dp/sqrt_zet**7
191 1793960 : temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2
192 1793960 : sr_int = prefac*temp
193 1793960 : s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
194 : !** derivatives
195 1793960 : k = sqrt_pi3*a**4/sqrt_zet**11
196 1793960 : dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
197 6414440 : DO ids = 2, nds
198 4620480 : n = ids - 1
199 4620480 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
200 4620480 : dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
201 4620480 : dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2
202 6414440 : s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3)
203 : END DO
204 : CASE (3)
205 : ![s|r^6|s] integral
206 148048 : prefac = sqrt_pi3/8.0_dp/sqrt_zet**9
207 148048 : temp = 105.0_dp + 210.0_dp*pfac*rab2
208 148048 : temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3
209 148048 : sr_int = prefac*temp
210 148048 : s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
211 : !** derivatives
212 148048 : k = sqrt_pi3*a**6/sqrt_zet**15
213 : dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 &
214 148048 : + 24.0_dp*pfac**3*rab2**2)
215 148048 : dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2)
216 322848 : DO ids = 2, nds
217 174800 : n = ids - 1
218 174800 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
219 174800 : dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
220 : dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
221 174800 : *exp_rab2*dsr_int(2)
222 174800 : dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)*(-xhi)**(n - 3)*k*exp_rab2
223 : s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) &
224 322848 : + dtemp(3) + dtemp(4)
225 : END DO
226 : CASE (4)
227 : ![s|r^8|s] integral
228 1284 : prefac = sqrt_pi3/16.0_dp/sqrt_zet**11
229 1284 : temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2
230 1284 : temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4
231 1284 : sr_int = prefac*temp
232 1284 : s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
233 : !** derivatives
234 1284 : k = sqrt_pi3*a**8/sqrt_zet**19
235 1284 : dsr_int(1) = 2520.0_dp*pfac + 3024.0_dp*pfac**2*rab2
236 : dsr_int(1) = dsr_int(1) + 864.0_dp*pfac**3*rab2**2 &
237 1284 : + 64.0_dp*pfac**4*rab2**3
238 1284 : dsr_int(1) = prefac*dsr_int(1)
239 1284 : dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2
240 1284 : dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2
241 1284 : dsr_int(2) = prefac*dsr_int(2)
242 1284 : dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2
243 1284 : dsr_int(3) = prefac*dsr_int(3)
244 3314 : DO ids = 2, nds
245 2030 : n = ids - 1
246 2030 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
247 2030 : dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
248 : dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
249 2030 : *exp_rab2*dsr_int(2)
250 : dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
251 2030 : *exp_rab2*dsr_int(3)
252 : dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) &
253 2030 : *k*exp_rab2
254 : s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
255 3314 : + dtemp(4) + dtemp(5)
256 : END DO
257 : CASE (5)
258 : ![s|r^10|s] integral
259 640 : prefac = sqrt_pi3/32.0_dp/sqrt_zet**13
260 640 : temp = 10395.0_dp + 34650.0_dp*pfac*rab2
261 640 : temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3
262 640 : temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5
263 640 : sr_int = prefac*temp
264 640 : s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
265 : !** derivatives
266 640 : k = sqrt_pi3*a**10/sqrt_zet**23
267 640 : dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2
268 640 : dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2
269 640 : dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3
270 640 : dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4
271 640 : dsr_int(1) = prefac*dsr_int(1)
272 640 : dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2
273 640 : dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2
274 640 : dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3
275 640 : dsr_int(2) = prefac*dsr_int(2)
276 640 : dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2
277 640 : dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2
278 640 : dsr_int(3) = prefac*dsr_int(3)
279 640 : dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2
280 640 : dsr_int(4) = prefac*dsr_int(4)
281 998 : DO ids = 2, nds
282 358 : n = ids - 1
283 358 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
284 358 : dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
285 : dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
286 358 : *exp_rab2*dsr_int(2)
287 : dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
288 358 : *exp_rab2*dsr_int(3)
289 : dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) &
290 358 : *exp_rab2*dsr_int(4)
291 : dtemp(6) = REAL(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) &
292 358 : *k*exp_rab2
293 : s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
294 998 : + dtemp(4) + dtemp(5) + dtemp(6)
295 : END DO
296 : CASE DEFAULT
297 : !*** general formula; factor 1.5-2 slower than explicit expressions
298 0 : prefac = exp_rab2/sqrt_zet**(2*il + 3)
299 0 : sr_int = 0.0_dp
300 0 : DO i = 0, il
301 0 : sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
302 : END DO
303 0 : s(ipgfa, jpgfb, kpgfb, il + 1, 1) = prefac*sr_int
304 6610516 : DO ids = 2, nds
305 0 : n = ids - 1
306 0 : nfac = 1
307 0 : dfsr_int = (-xhi)**n*sr_int
308 0 : DO j = 1, il
309 : temp = 0.0_dp
310 0 : DO i = j, il
311 0 : temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
312 : END DO
313 0 : nfac = nfac*(n - j + 1)
314 0 : dfsr_int = dfsr_int + temp*REAL(nfac, dp)/fac(j)*(-xhi)**(n - j)
315 : END DO
316 0 : s(ipgfa, jpgfb, kpgfb, il + 1, ids) = prefac*dfsr_int
317 : END DO
318 :
319 : END SELECT
320 :
321 : END DO
322 :
323 : END DO
324 : END DO
325 : END DO
326 :
327 71192 : DEALLOCATE (dtemp, dsr_int)
328 71192 : DEALLOCATE (coeff_srs)
329 :
330 71192 : END SUBROUTINE s_overlap_abb
331 :
332 : ! **************************************************************************************************
333 : !> \brief calculates the uncontracted, not normalized [0a|ra^(2m)|0b] two-center integral,
334 : !> where ra = r-Ra and Ra center of a
335 : !> \param la_max maximal l quantum number on a
336 : !> \param npgfa number of primitive Gaussian on a
337 : !> \param zeta set of exponents on a
338 : !> \param lb_max maximal l quantum number on b
339 : !> \param npgfb number of primitive Gaussian on a
340 : !> \param zetb set of exponents on a
341 : !> \param m exponent of the ra operator
342 : !> \param rab distance vector between a and b
343 : !> \param s ...
344 : !> \param calculate_forces ...
345 : ! **************************************************************************************************
346 4800 : SUBROUTINE s_ra2m_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, m, rab, s, calculate_forces)
347 :
348 : INTEGER, INTENT(IN) :: la_max, npgfa
349 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
350 : INTEGER, INTENT(IN) :: lb_max, npgfb
351 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
352 : INTEGER, INTENT(IN) :: m
353 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
354 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
355 : INTENT(INOUT) :: s
356 : LOGICAL, INTENT(IN) :: calculate_forces
357 :
358 : INTEGER :: i, ids, il, ipgfa, j, jpgfb, n, ndev, &
359 : nds, nfac
360 : REAL(KIND=dp) :: a, b, dfsr_int, exp_rab2, k, pfac, &
361 : prefac, rab2, sqrt_pi3, sqrt_zet, &
362 : sr_int, temp, xhi, zet
363 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dsr_int, dtemp
364 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: coeff_srs
365 :
366 : ! Distance of the centers a and b
367 4800 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
368 4800 : ndev = 0
369 4800 : IF (calculate_forces) ndev = 1
370 :
371 24000 : ALLOCATE (dtemp(m + 1), dsr_int(m + 1))
372 24000 : ALLOCATE (coeff_srs(m + 1, m + 1, m + 1))
373 4800 : CALL get_prefac_sabb(m, coeff_srs)
374 : !IF(m > 5) CALL get_prefac_sabb(m, coeff_srs)
375 4800 : sqrt_pi3 = SQRT(pi**3)
376 :
377 : ! Loops over all pairs of primitive Gaussian-type functions
378 9600 : DO ipgfa = 1, npgfa
379 14400 : DO jpgfb = 1, npgfb
380 :
381 : ! Calculate some prefactors
382 4800 : a = zeta(ipgfa)
383 4800 : b = zetb(jpgfb)
384 4800 : zet = a + b
385 4800 : xhi = a*b/zet
386 4800 : exp_rab2 = EXP(-xhi*rab2)
387 :
388 4800 : sqrt_zet = SQRT(zet)
389 4800 : pfac = b**2/zet
390 :
391 4800 : nds = la_max + lb_max + ndev + 1
392 28800 : DO il = 0, m
393 4800 : SELECT CASE (il)
394 : CASE (0)
395 : ! [s|s] integral
396 4800 : s(ipgfa, jpgfb, m - il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2
397 34080 : DO ids = 2, nds
398 29280 : n = ids - 1
399 34080 : s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, m - il + 1, 1)
400 : END DO
401 : CASE (1)
402 : ![s|r^2|s] integral
403 4800 : sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp
404 4800 : s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
405 4800 : k = sqrt_pi3*b**2/sqrt_zet**7
406 34080 : DO ids = 2, nds
407 29280 : n = ids - 1
408 : s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
409 34080 : + n*(-xhi)**(n - 1)*k*exp_rab2
410 : END DO
411 : CASE (2)
412 : ![s|r^4|s] integral
413 4800 : prefac = sqrt_pi3/4.0_dp/sqrt_zet**7
414 4800 : temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2
415 4800 : sr_int = prefac*temp
416 4800 : s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
417 : !** derivatives
418 4800 : k = sqrt_pi3*b**4/sqrt_zet**11
419 4800 : dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
420 34080 : DO ids = 2, nds
421 29280 : n = ids - 1
422 29280 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
423 29280 : dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
424 29280 : dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2
425 34080 : s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3)
426 : END DO
427 : CASE (3)
428 : ![s|r^6|s] integral
429 4800 : prefac = sqrt_pi3/8.0_dp/sqrt_zet**9
430 4800 : temp = 105.0_dp + 210.0_dp*pfac*rab2
431 4800 : temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3
432 4800 : sr_int = prefac*temp
433 4800 : s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
434 : !** derivatives
435 4800 : k = sqrt_pi3*b**6/sqrt_zet**15
436 : dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 &
437 4800 : + 24.0_dp*pfac**3*rab2**2)
438 4800 : dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2)
439 34080 : DO ids = 2, nds
440 29280 : n = ids - 1
441 29280 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
442 29280 : dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
443 : dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
444 29280 : *exp_rab2*dsr_int(2)
445 29280 : dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)*(-xhi)**(n - 3)*k*exp_rab2
446 : s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) &
447 34080 : + dtemp(3) + dtemp(4)
448 : END DO
449 : CASE (4)
450 : ![s|r^8|s] integral
451 0 : prefac = sqrt_pi3/16.0_dp/sqrt_zet**11
452 0 : temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2
453 0 : temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4
454 0 : sr_int = prefac*temp
455 0 : s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
456 : !** derivatives
457 0 : k = sqrt_pi3*b**8/sqrt_zet**19
458 0 : dsr_int(1) = 2520.0_dp*pfac + 3024.0_dp*pfac**2*rab2
459 : dsr_int(1) = dsr_int(1) + 864.0_dp*pfac**3*rab2**2 &
460 0 : + 64.0_dp*pfac**4*rab2**3
461 0 : dsr_int(1) = prefac*dsr_int(1)
462 0 : dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2
463 0 : dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2
464 0 : dsr_int(2) = prefac*dsr_int(2)
465 0 : dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2
466 0 : dsr_int(3) = prefac*dsr_int(3)
467 0 : DO ids = 2, nds
468 0 : n = ids - 1
469 0 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
470 0 : dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
471 : dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
472 0 : *exp_rab2*dsr_int(2)
473 : dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
474 0 : *exp_rab2*dsr_int(3)
475 : dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) &
476 0 : *k*exp_rab2
477 : s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
478 0 : + dtemp(4) + dtemp(5)
479 : END DO
480 : CASE (5)
481 : ![s|r^10|s] integral
482 0 : prefac = sqrt_pi3/32.0_dp/sqrt_zet**13
483 0 : temp = 10395.0_dp + 34650.0_dp*pfac*rab2
484 0 : temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3
485 0 : temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5
486 0 : sr_int = prefac*temp
487 0 : s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
488 : !** derivatives
489 0 : k = sqrt_pi3*b**10/sqrt_zet**23
490 0 : dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2
491 0 : dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2
492 0 : dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3
493 0 : dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4
494 0 : dsr_int(1) = prefac*dsr_int(1)
495 0 : dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2
496 0 : dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2
497 0 : dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3
498 0 : dsr_int(2) = prefac*dsr_int(2)
499 0 : dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2
500 0 : dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2
501 0 : dsr_int(3) = prefac*dsr_int(3)
502 0 : dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2
503 0 : dsr_int(4) = prefac*dsr_int(4)
504 0 : DO ids = 2, nds
505 0 : n = ids - 1
506 0 : dtemp(1) = (-xhi)**n*exp_rab2*sr_int
507 0 : dtemp(2) = REAL(n, dp)*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
508 : dtemp(3) = REAL(n**2 - n, dp)/2.0_dp*(-xhi)**(n - 2) &
509 0 : *exp_rab2*dsr_int(2)
510 : dtemp(4) = REAL(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
511 0 : *exp_rab2*dsr_int(3)
512 : dtemp(5) = REAL(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) &
513 0 : *exp_rab2*dsr_int(4)
514 : dtemp(6) = REAL(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) &
515 0 : *k*exp_rab2
516 : s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
517 0 : + dtemp(4) + dtemp(5) + dtemp(6)
518 : END DO
519 : CASE DEFAULT
520 0 : prefac = exp_rab2/sqrt_zet**(2*il + 3)
521 0 : sr_int = 0.0_dp
522 0 : DO i = 0, il
523 0 : sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
524 : END DO
525 0 : s(ipgfa, jpgfb, m - il + 1, 1) = prefac*sr_int
526 19200 : DO ids = 2, nds
527 0 : n = ids - 1
528 0 : nfac = 1
529 0 : dfsr_int = (-xhi)**n*sr_int
530 0 : DO j = 1, il
531 : temp = 0.0_dp
532 0 : DO i = j, il
533 0 : temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
534 : END DO
535 0 : nfac = nfac*(n - j + 1)
536 0 : dfsr_int = dfsr_int + temp*REAL(nfac, dp)/fac(j)*(-xhi)**(n - j)
537 : END DO
538 0 : s(ipgfa, jpgfb, m - il + 1, ids) = prefac*dfsr_int
539 : END DO
540 : END SELECT
541 : END DO
542 : END DO
543 : END DO
544 :
545 4800 : DEALLOCATE (coeff_srs)
546 4800 : DEALLOCATE (dtemp, dsr_int)
547 :
548 4800 : END SUBROUTINE s_ra2m_ab
549 :
550 : ! **************************************************************************************************
551 : !> \brief prefactor for the general formula to calculate the (0a|0b|0b) overlap integrals
552 : !> \param nl ...
553 : !> \param prefac ...
554 : ! **************************************************************************************************
555 4800 : SUBROUTINE get_prefac_sabb(nl, prefac)
556 : INTEGER, INTENT(IN) :: nl
557 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: prefac
558 :
559 : INTEGER :: il, j, k
560 : REAL(KIND=dp) :: sqrt_pi3, temp
561 :
562 4800 : sqrt_pi3 = SQRT(pi**3)
563 :
564 24000 : DO il = 0, nl
565 19200 : temp = dfac(2*il + 1)*sqrt_pi3*fac(il)/2.0_dp**il
566 72000 : DO j = 0, il
567 163200 : DO k = j, il
568 144000 : prefac(k + 1, j + 1, il + 1) = temp*2.0_dp**k/dfac(2*k + 1)/fac(il - k)/fac(k - j)
569 : END DO
570 : END DO
571 : END DO
572 :
573 4800 : END SUBROUTINE get_prefac_sabb
574 :
575 : ! **************************************************************************************************
576 : !> \brief calculates the uncontracted, not normalized [s|1/r12|s] two-center coulomb integral
577 : !> \param la_max maximal l quantum number on a
578 : !> \param npgfa number of primitive Gaussian on a
579 : !> \param zeta set of exponents on a
580 : !> \param lb_max maximal l quantum number on b
581 : !> \param npgfb number of primitive Gaussian on a
582 : !> \param zetb set of exponents on a
583 : !> \param omega parameter not needed, but given for the sake of the abstract interface
584 : !> \param rab distance vector between a and b
585 : !> \param v uncontracted coulomb integral of s functions
586 : !> \param calculate_forces ...
587 : ! **************************************************************************************************
588 17736992 : SUBROUTINE s_coulomb_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
589 :
590 : INTEGER, INTENT(IN) :: la_max, npgfa
591 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
592 : INTEGER, INTENT(IN) :: lb_max, npgfb
593 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
594 : REAL(KIND=dp), INTENT(IN) :: omega
595 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
596 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
597 : LOGICAL, INTENT(IN) :: calculate_forces
598 :
599 : INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
600 : REAL(KIND=dp) :: a, b, dummy, prefac, rab2, T, xhi, zet
601 17736992 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
602 :
603 17736992 : dummy = omega
604 :
605 : ! Distance of the centers a and b
606 17736992 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
607 17736992 : ndev = 0
608 17736992 : IF (calculate_forces) ndev = 1
609 17736992 : nmax = la_max + lb_max + ndev + 1
610 53210976 : ALLOCATE (f(0:nmax))
611 : ! Loops over all pairs of primitive Gaussian-type functions
612 35473984 : DO ipgfa = 1, npgfa
613 53210976 : DO jpgfb = 1, npgfb
614 :
615 : ! Calculate some prefactors
616 17736992 : a = zeta(ipgfa)
617 17736992 : b = zetb(jpgfb)
618 17736992 : zet = a + b
619 17736992 : xhi = a*b/zet
620 17736992 : prefac = 2.0_dp*SQRT(pi**5)/(a*b)/SQRT(zet)
621 17736992 : T = xhi*rab2
622 17736992 : CALL fgamma(nmax - 1, T, f)
623 :
624 84195840 : DO ids = 1, nmax
625 48721856 : n = ids - 1
626 66458848 : v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*f(n)
627 : END DO
628 :
629 : END DO
630 : END DO
631 17736992 : DEALLOCATE (f)
632 :
633 17736992 : END SUBROUTINE s_coulomb_ab
634 :
635 : ! **************************************************************************************************
636 : !> \brief calculates the uncontracted, not normalized [s|1/erf(omega*r12)/r12|s] two-center integral
637 : !> \param la_max maximal l quantum number on a
638 : !> \param npgfa number of primitive Gaussian on a
639 : !> \param zeta set of exponents on a
640 : !> \param lb_max maximal l quantum number on b
641 : !> \param npgfb number of primitive Gaussian on a
642 : !> \param zetb set of exponents on a
643 : !> \param omega parameter in the operator
644 : !> \param rab distance vector between a and b
645 : !> \param v uncontracted erf(r)/r integral of s functions
646 : !> \param calculate_forces ...
647 : ! **************************************************************************************************
648 4800 : SUBROUTINE s_verf_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
649 :
650 : INTEGER, INTENT(IN) :: la_max, npgfa
651 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
652 : INTEGER, INTENT(IN) :: lb_max, npgfb
653 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
654 : REAL(KIND=dp), INTENT(IN) :: omega
655 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
656 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
657 : LOGICAL, INTENT(IN) :: calculate_forces
658 :
659 : INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
660 : REAL(KIND=dp) :: a, Arg, b, comega, prefac, rab2, T, xhi, &
661 : zet
662 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
663 :
664 : ! Distance of the centers a and b
665 4800 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
666 4800 : ndev = 0
667 4800 : IF (calculate_forces) ndev = 1
668 4800 : nmax = la_max + lb_max + ndev + 1
669 14400 : ALLOCATE (f(0:nmax))
670 : ! Loops over all pairs of primitive Gaussian-type functions
671 9600 : DO ipgfa = 1, npgfa
672 14400 : DO jpgfb = 1, npgfb
673 :
674 : ! Calculate some prefactors
675 4800 : a = zeta(ipgfa)
676 4800 : b = zetb(jpgfb)
677 4800 : zet = a + b
678 4800 : xhi = a*b/zet
679 4800 : comega = omega**2/(omega**2 + xhi)
680 4800 : prefac = 2.0_dp*SQRT(pi**5)*SQRT(comega)/(a*b)/SQRT(zet)
681 4800 : T = xhi*rab2
682 4800 : Arg = comega*T
683 4800 : CALL fgamma(nmax - 1, Arg, f)
684 :
685 43680 : DO ids = 1, nmax
686 34080 : n = ids - 1
687 38880 : v(ipgfa, jpgfb, ids) = prefac*(-xhi*comega)**n*f(n)
688 : END DO
689 :
690 : END DO
691 : END DO
692 4800 : DEALLOCATE (f)
693 :
694 4800 : END SUBROUTINE s_verf_ab
695 :
696 : ! **************************************************************************************************
697 : !> \brief calculates the uncontracted, not normalized [s|1/erfc(omega*r12)/r12|s] two-center
698 : !> integral
699 : !> \param la_max maximal l quantum number on a
700 : !> \param npgfa number of primitive Gaussian on a
701 : !> \param zeta set of exponents on a
702 : !> \param lb_max maximal l quantum number on b
703 : !> \param npgfb number of primitive Gaussian on a
704 : !> \param zetb set of exponents on a
705 : !> \param omega parameter in the operator
706 : !> \param rab distance vector between a and b
707 : !> \param v uncontracted erf(r)/r integral of s functions
708 : !> \param calculate_forces ...
709 : ! **************************************************************************************************
710 4800 : SUBROUTINE s_verfc_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
711 :
712 : INTEGER, INTENT(IN) :: la_max, npgfa
713 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
714 : INTEGER, INTENT(IN) :: lb_max, npgfb
715 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
716 : REAL(KIND=dp), INTENT(IN) :: omega
717 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
718 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
719 : LOGICAL, INTENT(IN) :: calculate_forces
720 :
721 : INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
722 : REAL(KIND=dp) :: a, b, comega, comegaT, prefac, rab2, T, &
723 : xhi, zet
724 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf
725 :
726 : ! Distance of the centers a and b
727 4800 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
728 4800 : ndev = 0
729 4800 : IF (calculate_forces) ndev = 1
730 4800 : nmax = la_max + lb_max + ndev + 1
731 24000 : ALLOCATE (fv(0:nmax), fverf(0:nmax))
732 : ! Loops over all pairs of primitive Gaussian-type functions
733 9600 : DO ipgfa = 1, npgfa
734 14400 : DO jpgfb = 1, npgfb
735 :
736 : ! Calculate some prefactors
737 4800 : a = zeta(ipgfa)
738 4800 : b = zetb(jpgfb)
739 4800 : zet = a + b
740 4800 : xhi = a*b/zet
741 4800 : comega = omega**2/(omega**2 + xhi)
742 4800 : prefac = 2.0_dp*SQRT(pi**5)/(a*b)/SQRT(zet)
743 4800 : T = xhi*rab2
744 4800 : comegaT = comega*T
745 4800 : CALL fgamma(nmax - 1, T, fv)
746 4800 : CALL fgamma(nmax - 1, comegaT, fverf)
747 :
748 43680 : DO ids = 1, nmax
749 34080 : n = ids - 1
750 38880 : v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*(fv(n) - SQRT(comega)*comega**n*fverf(n))
751 : END DO
752 :
753 : END DO
754 : END DO
755 4800 : DEALLOCATE (fv, fverf)
756 :
757 4800 : END SUBROUTINE s_verfc_ab
758 :
759 : ! **************************************************************************************************
760 : !> \brief calculates the uncontracted, not normalized [s|exp(-omega*r12**2)/r12|s] two-center
761 : !> integral
762 : !> \param la_max maximal l quantum number on a
763 : !> \param npgfa number of primitive Gaussian on a
764 : !> \param zeta set of exponents on a
765 : !> \param lb_max maximal l quantum number on b
766 : !> \param npgfb number of primitive Gaussian on a
767 : !> \param zetb set of exponents on a
768 : !> \param omega parameter in the operator
769 : !> \param rab distance vector between a and b
770 : !> \param v uncontracted exp(-omega*r**2)/r integral of s functions
771 : !> \param calculate_forces ...
772 : ! **************************************************************************************************
773 4800 : SUBROUTINE s_vgauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
774 :
775 : INTEGER, INTENT(IN) :: la_max, npgfa
776 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
777 : INTEGER, INTENT(IN) :: lb_max, npgfb
778 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
779 : REAL(KIND=dp), INTENT(IN) :: omega
780 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
781 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
782 : LOGICAL, INTENT(IN) :: calculate_forces
783 :
784 : INTEGER :: ids, ipgfa, j, jpgfb, n, ndev, nmax
785 : REAL(KIND=dp) :: a, b, eta, etaT, expT, oeta, prefac, &
786 : rab2, T, xeta, xhi, zet
787 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
788 :
789 : ! Distance of the centers a and b
790 4800 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
791 4800 : ndev = 0
792 4800 : IF (calculate_forces) ndev = 1
793 4800 : nmax = la_max + lb_max + ndev + 1
794 14400 : ALLOCATE (f(0:nmax))
795 : ! Loops over all pairs of primitive Gaussian-type functions
796 134400 : v = 0.0_dp
797 9600 : DO ipgfa = 1, npgfa
798 14400 : DO jpgfb = 1, npgfb
799 :
800 : ! Calculate some prefactors
801 4800 : a = zeta(ipgfa)
802 4800 : b = zetb(jpgfb)
803 4800 : zet = a + b
804 4800 : xhi = a*b/zet
805 4800 : eta = xhi/(xhi + omega)
806 4800 : oeta = omega*eta
807 4800 : xeta = xhi*eta
808 4800 : T = xhi*rab2
809 4800 : expT = EXP(-omega/(omega + xhi)*T)
810 4800 : prefac = 2.0_dp*SQRT(pi**5/zet**3)/(xhi + omega)*expT
811 4800 : etaT = eta*T
812 4800 : CALL fgamma(nmax - 1, etaT, f)
813 :
814 43680 : DO ids = 1, nmax
815 34080 : n = ids - 1
816 184832 : DO j = 0, n
817 : v(ipgfa, jpgfb, ids) = v(ipgfa, jpgfb, ids) &
818 180032 : + prefac*fac(n)/fac(j)/fac(n - j)*(-oeta)**(n - j)*(-xeta)**j*f(j)
819 : END DO
820 : END DO
821 :
822 : END DO
823 : END DO
824 4800 : DEALLOCATE (f)
825 :
826 4800 : END SUBROUTINE s_vgauss_ab
827 :
828 : ! **************************************************************************************************
829 : !> \brief calculates the uncontracted, not normalized [s|exp(-omega*r12**2)|s] two-center
830 : !> integral
831 : !> \param la_max maximal l quantum number on a
832 : !> \param npgfa number of primitive Gaussian on a
833 : !> \param zeta set of exponents on a
834 : !> \param lb_max maximal l quantum number on b
835 : !> \param npgfb number of primitive Gaussian on a
836 : !> \param zetb set of exponents on a
837 : !> \param omega parameter in the operator
838 : !> \param rab distance vector between a and b
839 : !> \param v uncontracted exp(-omega*r**2) integral of s functions
840 : !> \param calculate_forces ...
841 : ! **************************************************************************************************
842 4800 : SUBROUTINE s_gauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
843 :
844 : INTEGER, INTENT(IN) :: la_max, npgfa
845 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
846 : INTEGER, INTENT(IN) :: lb_max, npgfb
847 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
848 : REAL(KIND=dp), INTENT(IN) :: omega
849 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
850 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
851 : LOGICAL, INTENT(IN) :: calculate_forces
852 :
853 : INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
854 : REAL(KIND=dp) :: a, b, eta, expT, oeta, prefac, rab2, T, &
855 : xhi, zet
856 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f
857 :
858 : ! Distance of the centers a and b
859 4800 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
860 4800 : ndev = 0
861 4800 : IF (calculate_forces) ndev = 1
862 4800 : nmax = la_max + lb_max + ndev + 1
863 14400 : ALLOCATE (f(0:nmax))
864 : ! Loops over all pairs of primitive Gaussian-type functions
865 9600 : DO ipgfa = 1, npgfa
866 14400 : DO jpgfb = 1, npgfb
867 :
868 : ! Calculate some prefactors
869 4800 : a = zeta(ipgfa)
870 4800 : b = zetb(jpgfb)
871 4800 : zet = a + b
872 4800 : xhi = a*b/zet
873 4800 : eta = xhi/(xhi + omega)
874 4800 : oeta = omega*eta
875 4800 : T = xhi*rab2
876 4800 : expT = EXP(-omega/(omega + xhi)*T)
877 4800 : prefac = pi**3/SQRT(zet**3)/SQRT((xhi + omega)**3)*expT
878 :
879 43680 : DO ids = 1, nmax
880 34080 : n = ids - 1
881 38880 : v(ipgfa, jpgfb, ids) = prefac*(-oeta)**n
882 : END DO
883 :
884 : END DO
885 : END DO
886 4800 : DEALLOCATE (f)
887 :
888 4800 : END SUBROUTINE s_gauss_ab
889 :
890 : ! **************************************************************************************************
891 : !> \brief Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives;
892 : !> this routine is more efficient for uncontracted basis sets (clow), e.g. for ri basis sets
893 : !> \param la set of l quantum numbers for a
894 : !> \param npgfa number of primitive Gaussian on a
895 : !> \param nshella number of shells for a
896 : !> \param scona_shg SHG contraction/normalization matrix for a
897 : !> \param lb set of l quantum numbers for b
898 : !> \param npgfb number of primitive Gaussian on b
899 : !> \param nshellb number of shells for b
900 : !> \param sconb_shg SHG contraction/normalization matrix for b
901 : !> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives
902 : !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
903 : !> \param calculate_forces ...
904 : ! **************************************************************************************************
905 100787 : SUBROUTINE contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, &
906 201574 : swork, swork_cont, calculate_forces)
907 :
908 : INTEGER, DIMENSION(:), INTENT(IN) :: la
909 : INTEGER, INTENT(IN) :: npgfa, nshella
910 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg
911 : INTEGER, DIMENSION(:), INTENT(IN) :: lb
912 : INTEGER, INTENT(IN) :: npgfb, nshellb
913 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg
914 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork
915 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont
916 : LOGICAL, INTENT(IN) :: calculate_forces
917 :
918 : INTEGER :: ids, ids_start, ipgfa, ishella, jpgfb, &
919 : jshellb, lai, lbj, ndev, nds
920 :
921 100787 : ndev = 0
922 100787 : IF (calculate_forces) ndev = 1
923 :
924 8777315 : swork_cont = 0.0_dp
925 352772 : DO ishella = 1, nshella
926 251985 : lai = la(ishella)
927 1009456 : DO jshellb = 1, nshellb
928 656684 : lbj = lb(jshellb)
929 656684 : nds = lai + lbj + 1
930 656684 : ids_start = nds - MIN(lai, lbj)
931 1571977 : DO ipgfa = 1, npgfa
932 2010984 : DO jpgfb = 1, npgfb
933 2745677 : DO ids = ids_start, nds + ndev
934 : swork_cont(ids, ishella, jshellb) = swork_cont(ids, ishella, jshellb) &
935 : + scona_shg(ipgfa, ishella) &
936 : *sconb_shg(jpgfb, jshellb) &
937 2082369 : *swork(ipgfa, jpgfb, ids)
938 : END DO
939 : END DO
940 : END DO
941 : END DO
942 : END DO
943 :
944 100787 : END SUBROUTINE contract_sint_ab_clow
945 :
946 : ! **************************************************************************************************
947 : !> \brief Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives;
948 : !> this routine is more efficient for highly contracted basis sets (chigh)
949 : !> \param npgfa number of primitive Gaussian on a
950 : !> \param nshella number of shells for a
951 : !> \param scona SHG contraction/normalization matrix for a
952 : !> \param npgfb number of primitive Gaussian on b
953 : !> \param nshellb number of shells for b
954 : !> \param sconb SHG contraction/normalization matrix for b
955 : !> \param nds maximal derivative of [s|O(r12)|s] with respect to rab2
956 : !> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives
957 : !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
958 : ! **************************************************************************************************
959 53328033 : SUBROUTINE contract_sint_ab_chigh(npgfa, nshella, scona, &
960 17776011 : npgfb, nshellb, sconb, &
961 17776011 : nds, swork, swork_cont)
962 :
963 : INTEGER, INTENT(IN) :: npgfa, nshella
964 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona
965 : INTEGER, INTENT(IN) :: npgfb, nshellb
966 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb
967 : INTEGER, INTENT(IN) :: nds
968 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: swork
969 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont
970 :
971 17776011 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc
972 :
973 107488974 : swork_cont = 0.0_dp
974 88880055 : ALLOCATE (work_pc(npgfb, nds, nshella))
975 :
976 : CALL dgemm("T", "N", npgfb*nds, nshella, npgfa, 1._dp, swork, npgfa, scona, npgfa, &
977 17808785 : 0.0_dp, work_pc, npgfb*nds)
978 : CALL dgemm("T", "N", nds*nshella, nshellb, npgfb, 1._dp, work_pc, npgfb, sconb, npgfb, &
979 194901801 : 0.0_dp, swork_cont, nds*nshella)
980 :
981 17776011 : DEALLOCATE (work_pc)
982 :
983 17776011 : END SUBROUTINE contract_sint_ab_chigh
984 :
985 : ! **************************************************************************************************
986 : !> \brief Contraction, normalization and combinatiorial combination of the [0a|(r-Ra)^(2m)|0b]
987 : !> integrals and their scalar derivatives
988 : !> \param npgfa number of primitive Gaussian on a
989 : !> \param nshella number of shells for a
990 : !> \param scon_ra2m contraction matrix on a containg the combinatorial factors
991 : !> \param npgfb number of primitive Gaussian on b
992 : !> \param nshellb number of shells for b
993 : !> \param sconb SHG contraction/normalization matrix for b
994 : !> \param swork matrix storing the uncontracted and unnormalized s-integrals and derivatives
995 : !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
996 : !> \param m exponent in operator (r-Ra)^(2m)
997 : !> \param nds_max maximal derivative with respect to rab2
998 : ! **************************************************************************************************
999 4800 : SUBROUTINE contract_s_ra2m_ab(npgfa, nshella, scon_ra2m, npgfb, nshellb, sconb, swork, swork_cont, &
1000 : m, nds_max)
1001 :
1002 : INTEGER, INTENT(IN) :: npgfa, nshella
1003 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_ra2m
1004 : INTEGER, INTENT(IN) :: npgfb, nshellb
1005 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb
1006 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
1007 : INTENT(INOUT) :: swork
1008 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: swork_cont
1009 : INTEGER, INTENT(IN) :: m, nds_max
1010 :
1011 : INTEGER :: i, my_m
1012 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc
1013 :
1014 4800 : my_m = m + 1
1015 24000 : ALLOCATE (work_pc(npgfb, nds_max, nshella))
1016 255264 : work_pc = 0.0_dp
1017 565344 : swork_cont = 0.0_dp
1018 24000 : DO i = 1, my_m
1019 : CALL dgemm("T", "N", npgfb*nds_max, nshella, npgfa, 1.0_dp, swork(1:npgfa, 1:npgfb, i, 1:nds_max), npgfa, &
1020 1095360 : scon_ra2m(1:npgfa, i, 1:nshella), npgfa, 1.0_dp, work_pc, npgfb*nds_max)
1021 : END DO
1022 : CALL dgemm("T", "N", nds_max*nshella, nshellb, npgfb, 1.0_dp, work_pc, npgfb, sconb, npgfb, 0.0_dp, &
1023 499392 : swork_cont, nds_max*nshella)
1024 :
1025 4800 : DEALLOCATE (work_pc)
1026 :
1027 4800 : END SUBROUTINE contract_s_ra2m_ab
1028 :
1029 : ! **************************************************************************************************
1030 : !> \brief Contraction and normalization of the (abb) overlap
1031 : !> \param la set of l quantum numbers on a
1032 : !> \param npgfa number of primitive Gaussians functions on a; orbital basis
1033 : !> \param nshella number of shells for a; orbital basis
1034 : !> \param scona_shg SHG contraction/normalization matrix for a; orbital basis
1035 : !> \param lb set of l quantum numbers on b; orbital basis
1036 : !> \param npgfb number of primitive Gaussians functions on b; orbital basis
1037 : !> \param nshellb number of shells for b; orbital basis
1038 : !> \param lcb set of l quantum numbers on b; ri basis
1039 : !> \param npgfcb number of primitive Gaussians functions on b; ri basis
1040 : !> \param nshellcb number of shells for b; ri basis
1041 : !> \param orbb_index index for orbital basis at B for sconb_mix
1042 : !> \param rib_index index for ri basis at B for sconb_mix
1043 : !> \param sconb_mix mixed contraction matrix for orbital and ri basis at B
1044 : !> \param nl_max related to the parameter m in (a|rb^(2m)|b)
1045 : !> \param nds_max derivative with respect to rab**2
1046 : !> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives
1047 : !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
1048 : !> \param calculate_forces ...
1049 : ! **************************************************************************************************
1050 4004 : SUBROUTINE contract_s_overlap_abb(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, &
1051 8008 : lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, &
1052 4004 : nl_max, nds_max, swork, swork_cont, calculate_forces)
1053 :
1054 : INTEGER, DIMENSION(:), INTENT(IN) :: la
1055 : INTEGER, INTENT(IN) :: npgfa, nshella
1056 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scona_shg
1057 : INTEGER, DIMENSION(:), INTENT(IN) :: lb
1058 : INTEGER, INTENT(IN) :: npgfb, nshellb
1059 : INTEGER, DIMENSION(:), INTENT(IN) :: lcb
1060 : INTEGER, INTENT(IN) :: npgfcb, nshellcb
1061 : INTEGER, DIMENSION(:, :), INTENT(IN) :: orbb_index, rib_index
1062 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix
1063 : INTEGER, INTENT(IN) :: nl_max, nds_max
1064 : REAL(KIND=dp), DIMENSION(:, :, :, :, :), &
1065 : INTENT(IN) :: swork
1066 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
1067 : INTENT(INOUT) :: swork_cont
1068 : LOGICAL, INTENT(IN) :: calculate_forces
1069 :
1070 : INTEGER :: forb, fri, ids, ids_start, iil, il, &
1071 : ishella, jpgfb, jshellb, kpgfb, &
1072 : kshellb, lai, lbb, lbj, lbk, ndev, &
1073 : nds, nl
1074 : REAL(KIND=dp), ALLOCATABLE, &
1075 4004 : DIMENSION(:, :, :, :, :) :: work_ppc
1076 :
1077 4004 : ndev = 0
1078 4004 : IF (calculate_forces) ndev = 1
1079 :
1080 28028 : ALLOCATE (work_ppc(npgfb, npgfcb, nl_max, nds_max, nshella))
1081 2382679 : work_ppc = 0.0_dp
1082 : CALL dgemm("T", "N", npgfb*npgfcb*nl_max*nds_max, nshella, npgfa, 1._dp, swork, npgfa, &
1083 2747209 : scona_shg, npgfa, 0.0_dp, work_ppc, npgfb*npgfcb*nl_max*nds_max)
1084 6816765 : swork_cont = 0.0_dp
1085 17594 : DO kshellb = 1, nshellcb
1086 13590 : lbk = lcb(kshellb)
1087 69296 : DO jshellb = 1, nshellb
1088 51702 : lbj = lb(jshellb)
1089 51702 : nl = INT((lbj + lbk)/2)
1090 51702 : IF (lbj == 0 .OR. lbk == 0) nl = 0
1091 227859 : DO ishella = 1, nshella
1092 162567 : lai = la(ishella)
1093 480531 : DO il = 0, nl
1094 266262 : lbb = lbj + lbk - 2*il
1095 266262 : nds = lai + lbb + 1
1096 266262 : ids_start = nds - MIN(lai, lbb)
1097 2238579 : DO jpgfb = 1, npgfb
1098 1809750 : forb = orbb_index(jpgfb, jshellb)
1099 3988158 : DO kpgfb = 1, npgfcb
1100 1912146 : fri = rib_index(kpgfb, kshellb)
1101 6522309 : DO ids = ids_start, nds + ndev
1102 8857283 : DO iil = 0, il
1103 : swork_cont(ids, il, ishella, jshellb, kshellb) = &
1104 : swork_cont(ids, il, ishella, jshellb, kshellb) &
1105 6945137 : + sconb_mix(iil + 1, fri, forb, il + 1)*work_ppc(jpgfb, kpgfb, il - iil + 1, ids, ishella)
1106 : END DO
1107 : END DO
1108 : END DO
1109 : END DO
1110 : END DO
1111 : END DO
1112 : END DO
1113 : END DO
1114 4004 : DEALLOCATE (work_ppc)
1115 :
1116 4004 : END SUBROUTINE contract_s_overlap_abb
1117 :
1118 : ! **************************************************************************************************
1119 : !> \brief Contraction and normalization of the (aba) overlap
1120 : !> \param la set of l quantum numbers on a; orbital basis
1121 : !> \param npgfa number of primitive Gaussians functions on a; orbital basis
1122 : !> \param nshella number of shells for a; orbital basis
1123 : !> \param lb set of l quantum numbers on b; orbital basis
1124 : !> \param npgfb number of primitive Gaussians functions on b; orbital basis
1125 : !> \param nshellb number of shells for b; orbital basis
1126 : !> \param sconb_shg SHG contraction/normalization matrix for b; orbital basis
1127 : !> \param lca set of l quantum numbers on a; ri basis
1128 : !> \param npgfca number of primitive Gaussians functions on a; ri basis
1129 : !> \param nshellca number of shells for a; ri basis
1130 : !> \param orba_index index for orbital basis at A for scona_mix
1131 : !> \param ria_index index for ri basis at A for scona_mix
1132 : !> \param scona_mix mixed contraction matrix for orbital and ri basis at A
1133 : !> \param nl_max related to the parameter m in (a|ra^(2m)|b)
1134 : !> \param nds_max maximal derivative with respect to rab**2
1135 : !> \param swork matrix of storing the uncontracted and unnormalized s-integrals and derivatives
1136 : !> \param swork_cont matrix storing the contracted and normalized s-integrals and derivatives
1137 : !> \param calculate_forces ...
1138 : ! **************************************************************************************************
1139 67188 : SUBROUTINE contract_s_overlap_aba(la, npgfa, nshella, lb, npgfb, nshellb, sconb_shg, &
1140 134376 : lca, npgfca, nshellca, orba_index, ria_index, scona_mix, &
1141 67188 : nl_max, nds_max, swork, swork_cont, calculate_forces)
1142 :
1143 : INTEGER, DIMENSION(:), INTENT(IN) :: la
1144 : INTEGER, INTENT(IN) :: npgfa, nshella
1145 : INTEGER, DIMENSION(:), INTENT(IN) :: lb
1146 : INTEGER, INTENT(IN) :: npgfb, nshellb
1147 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sconb_shg
1148 : INTEGER, DIMENSION(:), INTENT(IN) :: lca
1149 : INTEGER, INTENT(IN) :: npgfca, nshellca
1150 : INTEGER, DIMENSION(:, :), INTENT(IN) :: orba_index, ria_index
1151 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix
1152 : INTEGER, INTENT(IN) :: nl_max, nds_max
1153 : REAL(KIND=dp), DIMENSION(:, :, :, :, :), &
1154 : INTENT(IN) :: swork
1155 : REAL(KIND=dp), DIMENSION(:, 0:, :, :, :), &
1156 : INTENT(INOUT) :: swork_cont
1157 : LOGICAL, INTENT(IN) :: calculate_forces
1158 :
1159 : INTEGER :: forb, fri, ids, ids_start, iil, il, &
1160 : ipgfa, ishella, jshellb, kpgfa, &
1161 : kshella, laa, lai, lak, lbj, ndev, &
1162 : nds, nl
1163 : REAL(KIND=dp), ALLOCATABLE, &
1164 67188 : DIMENSION(:, :, :, :, :) :: work_ppc
1165 :
1166 67188 : ndev = 0
1167 67188 : IF (calculate_forces) ndev = 1
1168 :
1169 470316 : ALLOCATE (work_ppc(npgfa, npgfca, nl_max, nds_max, nshellb))
1170 53116313 : work_ppc = 0.0_dp
1171 : CALL dgemm("T", "N", npgfa*npgfca*nl_max*nds_max, nshellb, npgfb, 1._dp, swork, npgfb, &
1172 12730845 : sconb_shg, npgfb, 0.0_dp, work_ppc, npgfa*npgfca*nl_max*nds_max)
1173 136140716 : swork_cont = 0.0_dp
1174 237666 : DO kshella = 1, nshellca
1175 170478 : lak = lca(kshella)
1176 1022027 : DO jshellb = 1, nshellb
1177 784361 : lbj = lb(jshellb)
1178 4730610 : DO ishella = 1, nshella
1179 3775771 : lai = la(ishella)
1180 3775771 : nl = INT((lai + lak)/2)
1181 3775771 : IF (lai == 0 .OR. lak == 0) nl = 0
1182 10386951 : DO il = 0, nl
1183 5826819 : laa = lai + lak - 2*il
1184 5826819 : nds = laa + lbj + 1
1185 5826819 : ids_start = nds - MIN(laa, lbj)
1186 15489187 : DO kpgfa = 1, npgfca
1187 5886597 : fri = ria_index(kpgfa, kshella)
1188 42069363 : DO ipgfa = 1, npgfa
1189 30355947 : forb = orba_index(ipgfa, ishella)
1190 95756521 : DO ids = ids_start, nds + ndev
1191 172445835 : DO iil = 0, il
1192 : swork_cont(ids, il, ishella, jshellb, kshella) = &
1193 : swork_cont(ids, il, ishella, jshellb, kshella) &
1194 142089888 : + scona_mix(iil + 1, fri, forb, il + 1)*work_ppc(ipgfa, kpgfa, il - iil + 1, ids, jshellb)
1195 : END DO
1196 : END DO
1197 : END DO
1198 : END DO
1199 : END DO
1200 : END DO
1201 : END DO
1202 : END DO
1203 :
1204 67188 : DEALLOCATE (work_ppc)
1205 67188 : END SUBROUTINE contract_s_overlap_aba
1206 :
1207 : END MODULE s_contract_shg
|
