(git:374b731)
Loading...
Searching...
No Matches
s_contract_shg.F
Go to the documentation of this file.
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! **************************************************************************************************
22MODULE 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
43CONTAINS
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 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 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
72 ndev = 0
73 IF (calculate_forces) ndev = 1
74 ! Loops over all pairs of primitive Gaussian-type functions
75 DO ipgfa = 1, npgfa
76 DO jpgfb = 1, npgfb
77
78 ! Distance Screening !maybe later
79
80 ! Calculate some prefactors
81 a = zeta(ipgfa)
82 b = zetb(jpgfb)
83 zet = a + b
84 xhi = a*b/zet
85
86 ! [s|s] integral
87 s(ipgfa, jpgfb, 1) = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
88
89 DO ids = 2, la_max + lb_max + ndev + 1
90 s(ipgfa, jpgfb, ids) = -xhi*s(ipgfa, jpgfb, ids - 1)
91 END DO
92
93 END DO
94 END DO
95
96 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 SUBROUTINE s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, &
115 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 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 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
138 ndev = 0
139 IF (calculate_forces) ndev = 1
140
141 lbb_max = lb_max + lcb_max
142 nl = int(lbb_max/2)
143 IF (lb_max == 0 .OR. lcb_max == 0) nl = 0
144 lmax = la_max + lbb_max
145
146 ALLOCATE (dtemp(nl + 1), dsr_int(nl + 1))
147 ALLOCATE (coeff_srs(nl + 1, nl + 1, nl + 1))
148 IF (nl > 5) CALL get_prefac_sabb(nl, coeff_srs)
149
150 sqrt_pi3 = sqrt(pi**3)
151
152 ! Loops over all pairs of primitive Gaussian-type functions
153 DO kpgfb = 1, npgfcb
154 DO jpgfb = 1, npgfb
155 DO ipgfa = 1, npgfa
156
157 !Calculate some prefactors
158 a = zeta(ipgfa)
159 b = zetb(jpgfb) + zetcb(kpgfb)
160
161 zet = a + b
162 xhi = a*b/zet
163 exp_rab2 = exp(-xhi*rab2)
164
165 pfac = a**2/zet
166 sqrt_zet = sqrt(zet)
167
168 DO il = 0, nl
169 nds = lmax - 2*il + ndev + 1
170 SELECT CASE (il)
171 CASE (0)
172 ! [s|s] integral
173 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2
174 DO ids = 2, nds
175 n = ids - 1
176 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 sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp
181 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
182 k = sqrt_pi3*a**2/sqrt_zet**7
183 DO ids = 2, nds
184 n = ids - 1
185 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
186 + n*(-xhi)**(n - 1)*k*exp_rab2
187 END DO
188 CASE (2)
189 ![s|r^4|s] integral
190 prefac = sqrt_pi3/4.0_dp/sqrt_zet**7
191 temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2
192 sr_int = prefac*temp
193 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
194 !** derivatives
195 k = sqrt_pi3*a**4/sqrt_zet**11
196 dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
197 DO ids = 2, nds
198 n = ids - 1
199 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
200 dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
201 dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2
202 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 prefac = sqrt_pi3/8.0_dp/sqrt_zet**9
207 temp = 105.0_dp + 210.0_dp*pfac*rab2
208 temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3
209 sr_int = prefac*temp
210 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
211 !** derivatives
212 k = sqrt_pi3*a**6/sqrt_zet**15
213 dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 &
214 + 24.0_dp*pfac**3*rab2**2)
215 dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2)
216 DO ids = 2, nds
217 n = ids - 1
218 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
219 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 *exp_rab2*dsr_int(2)
222 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 + dtemp(3) + dtemp(4)
225 END DO
226 CASE (4)
227 ![s|r^8|s] integral
228 prefac = sqrt_pi3/16.0_dp/sqrt_zet**11
229 temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2
230 temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4
231 sr_int = prefac*temp
232 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
233 !** derivatives
234 k = sqrt_pi3*a**8/sqrt_zet**19
235 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 + 64.0_dp*pfac**4*rab2**3
238 dsr_int(1) = prefac*dsr_int(1)
239 dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2
240 dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2
241 dsr_int(2) = prefac*dsr_int(2)
242 dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2
243 dsr_int(3) = prefac*dsr_int(3)
244 DO ids = 2, nds
245 n = ids - 1
246 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
247 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 *exp_rab2*dsr_int(2)
250 dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
251 *exp_rab2*dsr_int(3)
252 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) &
253 *k*exp_rab2
254 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
255 + dtemp(4) + dtemp(5)
256 END DO
257 CASE (5)
258 ![s|r^10|s] integral
259 prefac = sqrt_pi3/32.0_dp/sqrt_zet**13
260 temp = 10395.0_dp + 34650.0_dp*pfac*rab2
261 temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3
262 temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5
263 sr_int = prefac*temp
264 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
265 !** derivatives
266 k = sqrt_pi3*a**10/sqrt_zet**23
267 dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2
268 dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2
269 dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3
270 dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4
271 dsr_int(1) = prefac*dsr_int(1)
272 dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2
273 dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2
274 dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3
275 dsr_int(2) = prefac*dsr_int(2)
276 dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2
277 dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2
278 dsr_int(3) = prefac*dsr_int(3)
279 dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2
280 dsr_int(4) = prefac*dsr_int(4)
281 DO ids = 2, nds
282 n = ids - 1
283 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
284 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 *exp_rab2*dsr_int(2)
287 dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
288 *exp_rab2*dsr_int(3)
289 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) &
290 *exp_rab2*dsr_int(4)
291 dtemp(6) = real(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) &
292 *k*exp_rab2
293 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
294 + dtemp(4) + dtemp(5) + dtemp(6)
295 END DO
296 CASE DEFAULT
297 !*** general formula; factor 1.5-2 slower than explicit expressions
298 prefac = exp_rab2/sqrt_zet**(2*il + 3)
299 sr_int = 0.0_dp
300 DO i = 0, il
301 sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
302 END DO
303 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = prefac*sr_int
304 DO ids = 2, nds
305 n = ids - 1
306 nfac = 1
307 dfsr_int = (-xhi)**n*sr_int
308 DO j = 1, il
309 temp = 0.0_dp
310 DO i = j, il
311 temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
312 END DO
313 nfac = nfac*(n - j + 1)
314 dfsr_int = dfsr_int + temp*real(nfac, dp)/fac(j)*(-xhi)**(n - j)
315 END DO
316 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 DEALLOCATE (dtemp, dsr_int)
328 DEALLOCATE (coeff_srs)
329
330 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 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 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 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
368 ndev = 0
369 IF (calculate_forces) ndev = 1
370
371 ALLOCATE (dtemp(m + 1), dsr_int(m + 1))
372 ALLOCATE (coeff_srs(m + 1, m + 1, m + 1))
373 CALL get_prefac_sabb(m, coeff_srs)
374 !IF(m > 5) CALL get_prefac_sabb(m, coeff_srs)
375 sqrt_pi3 = sqrt(pi**3)
376
377 ! Loops over all pairs of primitive Gaussian-type functions
378 DO ipgfa = 1, npgfa
379 DO jpgfb = 1, npgfb
380
381 ! Calculate some prefactors
382 a = zeta(ipgfa)
383 b = zetb(jpgfb)
384 zet = a + b
385 xhi = a*b/zet
386 exp_rab2 = exp(-xhi*rab2)
387
388 sqrt_zet = sqrt(zet)
389 pfac = b**2/zet
390
391 nds = la_max + lb_max + ndev + 1
392 DO il = 0, m
393 SELECT CASE (il)
394 CASE (0)
395 ! [s|s] integral
396 s(ipgfa, jpgfb, m - il + 1, 1) = (pi/zet)**(1.5_dp)*exp_rab2
397 DO ids = 2, nds
398 n = ids - 1
399 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 sr_int = sqrt_pi3/sqrt_zet**5*(3.0_dp + 2.0_dp*pfac*rab2)/2.0_dp
404 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
405 k = sqrt_pi3*b**2/sqrt_zet**7
406 DO ids = 2, nds
407 n = ids - 1
408 s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
409 + n*(-xhi)**(n - 1)*k*exp_rab2
410 END DO
411 CASE (2)
412 ![s|r^4|s] integral
413 prefac = sqrt_pi3/4.0_dp/sqrt_zet**7
414 temp = 15.0_dp + 20.0_dp*pfac*rab2 + 4.0_dp*(pfac*rab2)**2
415 sr_int = prefac*temp
416 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
417 !** derivatives
418 k = sqrt_pi3*b**4/sqrt_zet**11
419 dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
420 DO ids = 2, nds
421 n = ids - 1
422 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
423 dtemp(2) = n*(-xhi)**(n - 1)*exp_rab2*dsr_int(1)
424 dtemp(3) = (n**2 - n)*(-xhi)**(n - 2)*k*exp_rab2
425 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 prefac = sqrt_pi3/8.0_dp/sqrt_zet**9
430 temp = 105.0_dp + 210.0_dp*pfac*rab2
431 temp = temp + 84.0_dp*(pfac*rab2)**2 + 8.0_dp*(pfac*rab2)**3
432 sr_int = prefac*temp
433 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
434 !** derivatives
435 k = sqrt_pi3*b**6/sqrt_zet**15
436 dsr_int(1) = prefac*(210.0_dp*pfac + 168.0_dp*pfac**2*rab2 &
437 + 24.0_dp*pfac**3*rab2**2)
438 dsr_int(2) = prefac*(168.0_dp*pfac**2 + 48.0_dp*pfac**3*rab2)
439 DO ids = 2, nds
440 n = ids - 1
441 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
442 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 *exp_rab2*dsr_int(2)
445 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 + dtemp(3) + dtemp(4)
448 END DO
449 CASE (4)
450 ![s|r^8|s] integral
451 prefac = sqrt_pi3/16.0_dp/sqrt_zet**11
452 temp = 945.0_dp + 2520.0_dp*pfac*rab2 + 1512.0_dp*(pfac*rab2)**2
453 temp = temp + 288.0_dp*(pfac*rab2)**3 + 16.0_dp*(pfac*rab2)**4
454 sr_int = prefac*temp
455 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
456 !** derivatives
457 k = sqrt_pi3*b**8/sqrt_zet**19
458 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 + 64.0_dp*pfac**4*rab2**3
461 dsr_int(1) = prefac*dsr_int(1)
462 dsr_int(2) = 3024.0_dp*pfac**2 + 1728.0_dp*pfac**3*rab2
463 dsr_int(2) = dsr_int(2) + 192.0_dp*pfac**4*rab2**2
464 dsr_int(2) = prefac*dsr_int(2)
465 dsr_int(3) = 1728.0_dp*pfac**3 + 384.0_dp*pfac**4*rab2
466 dsr_int(3) = prefac*dsr_int(3)
467 DO ids = 2, nds
468 n = ids - 1
469 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
470 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 *exp_rab2*dsr_int(2)
473 dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
474 *exp_rab2*dsr_int(3)
475 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)*(-xhi)**(n - 4) &
476 *k*exp_rab2
477 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
478 + dtemp(4) + dtemp(5)
479 END DO
480 CASE (5)
481 ![s|r^10|s] integral
482 prefac = sqrt_pi3/32.0_dp/sqrt_zet**13
483 temp = 10395.0_dp + 34650.0_dp*pfac*rab2
484 temp = temp + 27720.0_dp*(pfac*rab2)**2 + 7920.0_dp*(pfac*rab2)**3
485 temp = temp + 880.0_dp*(pfac*rab2)**4 + 32.0_dp*(pfac*rab2)**5
486 sr_int = prefac*temp
487 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
488 !** derivatives
489 k = sqrt_pi3*b**10/sqrt_zet**23
490 dsr_int(1) = 34650.0_dp*pfac + 55440.0_dp*pfac**2*rab2
491 dsr_int(1) = dsr_int(1) + 23760.0_dp*pfac**3*rab2**2
492 dsr_int(1) = dsr_int(1) + 3520.0_dp*pfac**4*rab2**3
493 dsr_int(1) = dsr_int(1) + 160.0_dp*pfac**5*rab2**4
494 dsr_int(1) = prefac*dsr_int(1)
495 dsr_int(2) = 55440.0_dp*pfac**2 + 47520.0_dp*pfac**3*rab2
496 dsr_int(2) = dsr_int(2) + 10560.0_dp*pfac**4*rab2**2
497 dsr_int(2) = dsr_int(2) + 640.0_dp*pfac**5*rab2**3
498 dsr_int(2) = prefac*dsr_int(2)
499 dsr_int(3) = 47520.0_dp*pfac**3 + 21120.0_dp*pfac**4*rab2
500 dsr_int(3) = dsr_int(3) + 1920.0_dp*pfac**5*rab2**2
501 dsr_int(3) = prefac*dsr_int(3)
502 dsr_int(4) = 21120.0_dp*pfac**4 + 3840.0_dp*pfac**5*rab2
503 dsr_int(4) = prefac*dsr_int(4)
504 DO ids = 2, nds
505 n = ids - 1
506 dtemp(1) = (-xhi)**n*exp_rab2*sr_int
507 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 *exp_rab2*dsr_int(2)
510 dtemp(4) = real(n*(n - 1)*(n - 2), dp)/6.0_dp*(-xhi)**(n - 3) &
511 *exp_rab2*dsr_int(3)
512 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3), dp)/24.0_dp*(-xhi)**(n - 4) &
513 *exp_rab2*dsr_int(4)
514 dtemp(6) = real(n*(n - 1)*(n - 2)*(n - 3)*(n - 4), dp)*(-xhi)**(n - 5) &
515 *k*exp_rab2
516 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
517 + dtemp(4) + dtemp(5) + dtemp(6)
518 END DO
519 CASE DEFAULT
520 prefac = exp_rab2/sqrt_zet**(2*il + 3)
521 sr_int = 0.0_dp
522 DO i = 0, il
523 sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
524 END DO
525 s(ipgfa, jpgfb, m - il + 1, 1) = prefac*sr_int
526 DO ids = 2, nds
527 n = ids - 1
528 nfac = 1
529 dfsr_int = (-xhi)**n*sr_int
530 DO j = 1, il
531 temp = 0.0_dp
532 DO i = j, il
533 temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
534 END DO
535 nfac = nfac*(n - j + 1)
536 dfsr_int = dfsr_int + temp*real(nfac, dp)/fac(j)*(-xhi)**(n - j)
537 END DO
538 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 DEALLOCATE (coeff_srs)
546 DEALLOCATE (dtemp, dsr_int)
547
548 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 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 sqrt_pi3 = sqrt(pi**3)
563
564 DO il = 0, nl
565 temp = dfac(2*il + 1)*sqrt_pi3*fac(il)/2.0_dp**il
566 DO j = 0, il
567 DO k = j, il
568 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
602
603 dummy = omega
604
605 ! Distance of the centers a and b
606 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
607 ndev = 0
608 IF (calculate_forces) ndev = 1
609 nmax = la_max + lb_max + ndev + 1
610 ALLOCATE (f(0:nmax))
611 ! Loops over all pairs of primitive Gaussian-type functions
612 DO ipgfa = 1, npgfa
613 DO jpgfb = 1, npgfb
614
615 ! Calculate some prefactors
616 a = zeta(ipgfa)
617 b = zetb(jpgfb)
618 zet = a + b
619 xhi = a*b/zet
620 prefac = 2.0_dp*sqrt(pi**5)/(a*b)/sqrt(zet)
621 t = xhi*rab2
622 CALL fgamma(nmax - 1, t, f)
623
624 DO ids = 1, nmax
625 n = ids - 1
626 v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*f(n)
627 END DO
628
629 END DO
630 END DO
631 DEALLOCATE (f)
632
633 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
663
664 ! Distance of the centers a and b
665 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
666 ndev = 0
667 IF (calculate_forces) ndev = 1
668 nmax = la_max + lb_max + ndev + 1
669 ALLOCATE (f(0:nmax))
670 ! Loops over all pairs of primitive Gaussian-type functions
671 DO ipgfa = 1, npgfa
672 DO jpgfb = 1, npgfb
673
674 ! Calculate some prefactors
675 a = zeta(ipgfa)
676 b = zetb(jpgfb)
677 zet = a + b
678 xhi = a*b/zet
679 comega = omega**2/(omega**2 + xhi)
680 prefac = 2.0_dp*sqrt(pi**5)*sqrt(comega)/(a*b)/sqrt(zet)
681 t = xhi*rab2
682 arg = comega*t
683 CALL fgamma(nmax - 1, arg, f)
684
685 DO ids = 1, nmax
686 n = ids - 1
687 v(ipgfa, jpgfb, ids) = prefac*(-xhi*comega)**n*f(n)
688 END DO
689
690 END DO
691 END DO
692 DEALLOCATE (f)
693
694 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fv, fverf
725
726 ! Distance of the centers a and b
727 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
728 ndev = 0
729 IF (calculate_forces) ndev = 1
730 nmax = la_max + lb_max + ndev + 1
731 ALLOCATE (fv(0:nmax), fverf(0:nmax))
732 ! Loops over all pairs of primitive Gaussian-type functions
733 DO ipgfa = 1, npgfa
734 DO jpgfb = 1, npgfb
735
736 ! Calculate some prefactors
737 a = zeta(ipgfa)
738 b = zetb(jpgfb)
739 zet = a + b
740 xhi = a*b/zet
741 comega = omega**2/(omega**2 + xhi)
742 prefac = 2.0_dp*sqrt(pi**5)/(a*b)/sqrt(zet)
743 t = xhi*rab2
744 comegat = comega*t
745 CALL fgamma(nmax - 1, t, fv)
746 CALL fgamma(nmax - 1, comegat, fverf)
747
748 DO ids = 1, nmax
749 n = ids - 1
750 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 DEALLOCATE (fv, fverf)
756
757 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
788
789 ! Distance of the centers a and b
790 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
791 ndev = 0
792 IF (calculate_forces) ndev = 1
793 nmax = la_max + lb_max + ndev + 1
794 ALLOCATE (f(0:nmax))
795 ! Loops over all pairs of primitive Gaussian-type functions
796 v = 0.0_dp
797 DO ipgfa = 1, npgfa
798 DO jpgfb = 1, npgfb
799
800 ! Calculate some prefactors
801 a = zeta(ipgfa)
802 b = zetb(jpgfb)
803 zet = a + b
804 xhi = a*b/zet
805 eta = xhi/(xhi + omega)
806 oeta = omega*eta
807 xeta = xhi*eta
808 t = xhi*rab2
809 expt = exp(-omega/(omega + xhi)*t)
810 prefac = 2.0_dp*sqrt(pi**5/zet**3)/(xhi + omega)*expt
811 etat = eta*t
812 CALL fgamma(nmax - 1, etat, f)
813
814 DO ids = 1, nmax
815 n = ids - 1
816 DO j = 0, n
817 v(ipgfa, jpgfb, ids) = v(ipgfa, jpgfb, ids) &
818 + 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 DEALLOCATE (f)
825
826 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: f
857
858 ! Distance of the centers a and b
859 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
860 ndev = 0
861 IF (calculate_forces) ndev = 1
862 nmax = la_max + lb_max + ndev + 1
863 ALLOCATE (f(0:nmax))
864 ! Loops over all pairs of primitive Gaussian-type functions
865 DO ipgfa = 1, npgfa
866 DO jpgfb = 1, npgfb
867
868 ! Calculate some prefactors
869 a = zeta(ipgfa)
870 b = zetb(jpgfb)
871 zet = a + b
872 xhi = a*b/zet
873 eta = xhi/(xhi + omega)
874 oeta = omega*eta
875 t = xhi*rab2
876 expt = exp(-omega/(omega + xhi)*t)
877 prefac = pi**3/sqrt(zet**3)/sqrt((xhi + omega)**3)*expt
878
879 DO ids = 1, nmax
880 n = ids - 1
881 v(ipgfa, jpgfb, ids) = prefac*(-oeta)**n
882 END DO
883
884 END DO
885 END DO
886 DEALLOCATE (f)
887
888 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 SUBROUTINE contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, &
906 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 ndev = 0
922 IF (calculate_forces) ndev = 1
923
924 swork_cont = 0.0_dp
925 DO ishella = 1, nshella
926 lai = la(ishella)
927 DO jshellb = 1, nshellb
928 lbj = lb(jshellb)
929 nds = lai + lbj + 1
930 ids_start = nds - min(lai, lbj)
931 DO ipgfa = 1, npgfa
932 DO jpgfb = 1, npgfb
933 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 *swork(ipgfa, jpgfb, ids)
938 END DO
939 END DO
940 END DO
941 END DO
942 END DO
943
944 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 SUBROUTINE contract_sint_ab_chigh(npgfa, nshella, scona, &
960 npgfb, nshellb, sconb, &
961 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc
972
973 swork_cont = 0.0_dp
974 ALLOCATE (work_pc(npgfb, nds, nshella))
975
976 CALL dgemm("T", "N", npgfb*nds, nshella, npgfa, 1._dp, swork, npgfa, scona, npgfa, &
977 0.0_dp, work_pc, npgfb*nds)
978 CALL dgemm("T", "N", nds*nshella, nshellb, npgfb, 1._dp, work_pc, npgfb, sconb, npgfb, &
979 0.0_dp, swork_cont, nds*nshella)
980
981 DEALLOCATE (work_pc)
982
983 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work_pc
1013
1014 my_m = m + 1
1015 ALLOCATE (work_pc(npgfb, nds_max, nshella))
1016 work_pc = 0.0_dp
1017 swork_cont = 0.0_dp
1018 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 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 swork_cont, nds_max*nshella)
1024
1025 DEALLOCATE (work_pc)
1026
1027 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 SUBROUTINE contract_s_overlap_abb(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, &
1051 lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, &
1052 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 DIMENSION(:, :, :, :, :) :: work_ppc
1076
1077 ndev = 0
1078 IF (calculate_forces) ndev = 1
1079
1080 ALLOCATE (work_ppc(npgfb, npgfcb, nl_max, nds_max, nshella))
1081 work_ppc = 0.0_dp
1082 CALL dgemm("T", "N", npgfb*npgfcb*nl_max*nds_max, nshella, npgfa, 1._dp, swork, npgfa, &
1083 scona_shg, npgfa, 0.0_dp, work_ppc, npgfb*npgfcb*nl_max*nds_max)
1084 swork_cont = 0.0_dp
1085 DO kshellb = 1, nshellcb
1086 lbk = lcb(kshellb)
1087 DO jshellb = 1, nshellb
1088 lbj = lb(jshellb)
1089 nl = int((lbj + lbk)/2)
1090 IF (lbj == 0 .OR. lbk == 0) nl = 0
1091 DO ishella = 1, nshella
1092 lai = la(ishella)
1093 DO il = 0, nl
1094 lbb = lbj + lbk - 2*il
1095 nds = lai + lbb + 1
1096 ids_start = nds - min(lai, lbb)
1097 DO jpgfb = 1, npgfb
1098 forb = orbb_index(jpgfb, jshellb)
1099 DO kpgfb = 1, npgfcb
1100 fri = rib_index(kpgfb, kshellb)
1101 DO ids = ids_start, nds + ndev
1102 DO iil = 0, il
1103 swork_cont(ids, il, ishella, jshellb, kshellb) = &
1104 swork_cont(ids, il, ishella, jshellb, kshellb) &
1105 + 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 DEALLOCATE (work_ppc)
1115
1116 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 SUBROUTINE contract_s_overlap_aba(la, npgfa, nshella, lb, npgfb, nshellb, sconb_shg, &
1140 lca, npgfca, nshellca, orba_index, ria_index, scona_mix, &
1141 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 DIMENSION(:, :, :, :, :) :: work_ppc
1165
1166 ndev = 0
1167 IF (calculate_forces) ndev = 1
1168
1169 ALLOCATE (work_ppc(npgfa, npgfca, nl_max, nds_max, nshellb))
1170 work_ppc = 0.0_dp
1171 CALL dgemm("T", "N", npgfa*npgfca*nl_max*nds_max, nshellb, npgfb, 1._dp, swork, npgfb, &
1172 sconb_shg, npgfb, 0.0_dp, work_ppc, npgfa*npgfca*nl_max*nds_max)
1173 swork_cont = 0.0_dp
1174 DO kshella = 1, nshellca
1175 lak = lca(kshella)
1176 DO jshellb = 1, nshellb
1177 lbj = lb(jshellb)
1178 DO ishella = 1, nshella
1179 lai = la(ishella)
1180 nl = int((lai + lak)/2)
1181 IF (lai == 0 .OR. lak == 0) nl = 0
1182 DO il = 0, nl
1183 laa = lai + lak - 2*il
1184 nds = laa + lbj + 1
1185 ids_start = nds - min(laa, lbj)
1186 DO kpgfa = 1, npgfca
1187 fri = ria_index(kpgfa, kshella)
1188 DO ipgfa = 1, npgfa
1189 forb = orba_index(ipgfa, ishella)
1190 DO ids = ids_start, nds + ndev
1191 DO iil = 0, il
1192 swork_cont(ids, il, ishella, jshellb, kshella) = &
1193 swork_cont(ids, il, ishella, jshellb, kshella) &
1194 + 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 DEALLOCATE (work_ppc)
1205 END SUBROUTINE contract_s_overlap_aba
1206
1207END MODULE s_contract_shg
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:195
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
s_contract_shg
Routines for calculating the s-integrals and their scalar derivatives with respect to rab2 over solid...
Definition s_contract_shg.F:22
s_contract_shg::contract_s_overlap_aba
subroutine, public contract_s_overlap_aba(la, npgfa, nshella, lb, npgfb, nshellb, sconb_shg, lca, npgfca, nshellca, orba_index, ria_index, scona_mix, nl_max, nds_max, swork, swork_cont, calculate_forces)
Contraction and normalization of the (aba) overlap.
Definition s_contract_shg.F:1142
s_contract_shg::s_verf_ab
subroutine, public s_verf_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/erf(omega*r12)/r12|s] two-center integral
Definition s_contract_shg.F:649
s_contract_shg::s_gauss_ab
subroutine, public s_gauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|exp(-omega*r12**2)|s] two-center integral
Definition s_contract_shg.F:843
s_contract_shg::s_overlap_abb
subroutine, public s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, rab, s, calculate_forces)
calculates [s|ra^n|s] integrals for [aba] and the [s|rb^n|s] integrals for [abb]
Definition s_contract_shg.F:116
s_contract_shg::s_vgauss_ab
subroutine, public s_vgauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|exp(-omega*r12**2)/r12|s] two-center integral
Definition s_contract_shg.F:774
s_contract_shg::contract_s_ra2m_ab
subroutine, public contract_s_ra2m_ab(npgfa, nshella, scon_ra2m, npgfb, nshellb, sconb, swork, swork_cont, m, nds_max)
Contraction, normalization and combinatiorial combination of the [0a|(r-Ra)^(2m)|0b] integrals and th...
Definition s_contract_shg.F:1001
s_contract_shg::contract_sint_ab_chigh
subroutine, public contract_sint_ab_chigh(npgfa, nshella, scona, npgfb, nshellb, sconb, nds, swork, swork_cont)
Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; this routin...
Definition s_contract_shg.F:962
s_contract_shg::s_coulomb_ab
subroutine, public s_coulomb_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/r12|s] two-center coulomb integral
Definition s_contract_shg.F:589
s_contract_shg::contract_s_overlap_abb
subroutine, public contract_s_overlap_abb(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, nl_max, nds_max, swork, swork_cont, calculate_forces)
Contraction and normalization of the (abb) overlap.
Definition s_contract_shg.F:1053
s_contract_shg::s_verfc_ab
subroutine, public s_verfc_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/erfc(omega*r12)/r12|s] two-center integral
Definition s_contract_shg.F:711
s_contract_shg::contract_sint_ab_clow
subroutine, public contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, swork, swork_cont, calculate_forces)
Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; this routin...
Definition s_contract_shg.F:907
s_contract_shg::s_overlap_ab
subroutine, public s_overlap_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, rab, s, calculate_forces)
calculates the uncontracted, not normalized [s|s] overlap
Definition s_contract_shg.F:58
s_contract_shg::s_ra2m_ab
subroutine, public s_ra2m_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, m, rab, s, calculate_forces)
calculates the uncontracted, not normalized [0a|ra^(2m)|0b] two-center integral, where ra = r-Ra and ...
Definition s_contract_shg.F:347