28 #include "../base/base_uses.f90"
34 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
's_contract_shg'
57 SUBROUTINE s_overlap_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, rab, s, calculate_forces)
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
67 INTEGER :: ids, ipgfa, jpgfb, ndev
68 REAL(kind=
dp) :: a, b, rab2, xhi, zet
71 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
73 IF (calculate_forces) ndev = 1
87 s(ipgfa, jpgfb, 1) = (
pi/zet)**(1.5_dp)*exp(-xhi*rab2)
89 DO ids = 2, la_max + lb_max + ndev + 1
90 s(ipgfa, jpgfb, ids) = -xhi*s(ipgfa, jpgfb, ids - 1)
114 SUBROUTINE s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, &
115 rab, s, calculate_forces)
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(:, :, :, :, :), &
126 LOGICAL,
INTENT(IN) :: calculate_forces
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
137 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
139 IF (calculate_forces) ndev = 1
141 lbb_max = lb_max + lcb_max
143 IF (lb_max == 0 .OR. lcb_max == 0) nl = 0
144 lmax = la_max + lbb_max
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)
150 sqrt_pi3 = sqrt(
pi**3)
159 b = zetb(jpgfb) + zetcb(kpgfb)
163 exp_rab2 = exp(-xhi*rab2)
169 nds = lmax - 2*il + ndev + 1
173 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = (
pi/zet)**(1.5_dp)*exp_rab2
176 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, kpgfb, il + 1, 1)
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
185 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
186 + n*(-xhi)**(n - 1)*k*exp_rab2
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
193 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
195 k = sqrt_pi3*a**4/sqrt_zet**11
196 dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
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)
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
210 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
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)
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) &
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)
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
232 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
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)
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) &
250 dtemp(4) = real(n*(n - 1)*(n - 2),
dp)/6.0_dp*(-xhi)**(n - 3) &
252 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3),
dp)*(-xhi)**(n - 4) &
254 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
255 + dtemp(4) + dtemp(5)
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
264 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = exp_rab2*sr_int
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)
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) &
287 dtemp(4) = real(n*(n - 1)*(n - 2),
dp)/6.0_dp*(-xhi)**(n - 3) &
289 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3),
dp)/24.0_dp*(-xhi)**(n - 4) &
291 dtemp(6) = real(n*(n - 1)*(n - 2)*(n - 3)*(n - 4),
dp)*(-xhi)**(n - 5) &
293 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
294 + dtemp(4) + dtemp(5) + dtemp(6)
298 prefac = exp_rab2/sqrt_zet**(2*il + 3)
301 sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
303 s(ipgfa, jpgfb, kpgfb, il + 1, 1) = prefac*sr_int
307 dfsr_int = (-xhi)**n*sr_int
311 temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
313 nfac = nfac*(n - j + 1)
314 dfsr_int = dfsr_int + temp*real(nfac,
dp)/
fac(j)*(-xhi)**(n - j)
316 s(ipgfa, jpgfb, kpgfb, il + 1, ids) = prefac*dfsr_int
327 DEALLOCATE (dtemp, dsr_int)
328 DEALLOCATE (coeff_srs)
346 SUBROUTINE s_ra2m_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, m, rab, s, calculate_forces)
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(:, :, :, :), &
356 LOGICAL,
INTENT(IN) :: calculate_forces
358 INTEGER :: i, ids, il, ipgfa, j, jpgfb, n, ndev, &
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
367 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
369 IF (calculate_forces) ndev = 1
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)
375 sqrt_pi3 = sqrt(
pi**3)
386 exp_rab2 = exp(-xhi*rab2)
391 nds = la_max + lb_max + ndev + 1
396 s(ipgfa, jpgfb, m - il + 1, 1) = (
pi/zet)**(1.5_dp)*exp_rab2
399 s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*s(ipgfa, jpgfb, m - il + 1, 1)
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
408 s(ipgfa, jpgfb, m - il + 1, ids) = (-xhi)**n*exp_rab2*sr_int &
409 + n*(-xhi)**(n - 1)*k*exp_rab2
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
416 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
418 k = sqrt_pi3*b**4/sqrt_zet**11
419 dsr_int(1) = prefac*(20.0_dp*pfac + 8.0_dp*pfac**2*rab2)
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)
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
433 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
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)
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) &
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)
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
455 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
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)
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) &
473 dtemp(4) = real(n*(n - 1)*(n - 2),
dp)/6.0_dp*(-xhi)**(n - 3) &
475 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3),
dp)*(-xhi)**(n - 4) &
477 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
478 + dtemp(4) + dtemp(5)
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
487 s(ipgfa, jpgfb, m - il + 1, 1) = exp_rab2*sr_int
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)
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) &
510 dtemp(4) = real(n*(n - 1)*(n - 2),
dp)/6.0_dp*(-xhi)**(n - 3) &
512 dtemp(5) = real(n*(n - 1)*(n - 2)*(n - 3),
dp)/24.0_dp*(-xhi)**(n - 4) &
514 dtemp(6) = real(n*(n - 1)*(n - 2)*(n - 3)*(n - 4),
dp)*(-xhi)**(n - 5) &
516 s(ipgfa, jpgfb, m - il + 1, ids) = dtemp(1) + dtemp(2) + dtemp(3) &
517 + dtemp(4) + dtemp(5) + dtemp(6)
520 prefac = exp_rab2/sqrt_zet**(2*il + 3)
523 sr_int = sr_int + coeff_srs(i + 1, 1, il + 1)*(pfac)**i*rab2**i
525 s(ipgfa, jpgfb, m - il + 1, 1) = prefac*sr_int
529 dfsr_int = (-xhi)**n*sr_int
533 temp = temp + coeff_srs(i + 1, j + 1, il + 1)*(pfac)**i*rab2**(i - j)
535 nfac = nfac*(n - j + 1)
536 dfsr_int = dfsr_int + temp*real(nfac,
dp)/
fac(j)*(-xhi)**(n - j)
538 s(ipgfa, jpgfb, m - il + 1, ids) = prefac*dfsr_int
545 DEALLOCATE (coeff_srs)
546 DEALLOCATE (dtemp, dsr_int)
555 SUBROUTINE get_prefac_sabb(nl, prefac)
556 INTEGER,
INTENT(IN) :: nl
557 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: prefac
560 REAL(kind=
dp) :: sqrt_pi3, temp
562 sqrt_pi3 = sqrt(
pi**3)
565 temp =
dfac(2*il + 1)*sqrt_pi3*
fac(il)/2.0_dp**il
568 prefac(k + 1, j + 1, il + 1) = temp*2.0_dp**k/
dfac(2*k + 1)/
fac(il - k)/
fac(k - j)
573 END SUBROUTINE get_prefac_sabb
588 SUBROUTINE s_coulomb_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
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
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
606 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
608 IF (calculate_forces) ndev = 1
609 nmax = la_max + lb_max + ndev + 1
620 prefac = 2.0_dp*sqrt(
pi**5)/(a*b)/sqrt(zet)
622 CALL fgamma(nmax - 1, t, f)
626 v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*f(n)
648 SUBROUTINE s_verf_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
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
659 INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
660 REAL(kind=
dp) :: a, arg, b, comega, prefac, rab2, t, xhi, &
662 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: f
665 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
667 IF (calculate_forces) ndev = 1
668 nmax = la_max + lb_max + ndev + 1
679 comega = omega**2/(omega**2 + xhi)
680 prefac = 2.0_dp*sqrt(
pi**5)*sqrt(comega)/(a*b)/sqrt(zet)
683 CALL fgamma(nmax - 1, arg, f)
687 v(ipgfa, jpgfb, ids) = prefac*(-xhi*comega)**n*f(n)
710 SUBROUTINE s_verfc_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
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
721 INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
722 REAL(kind=
dp) :: a, b, comega, comegat, prefac, rab2, t, &
724 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: fv, fverf
727 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
729 IF (calculate_forces) ndev = 1
730 nmax = la_max + lb_max + ndev + 1
731 ALLOCATE (fv(0:nmax), fverf(0:nmax))
741 comega = omega**2/(omega**2 + xhi)
742 prefac = 2.0_dp*sqrt(
pi**5)/(a*b)/sqrt(zet)
745 CALL fgamma(nmax - 1, t, fv)
746 CALL fgamma(nmax - 1, comegat, fverf)
750 v(ipgfa, jpgfb, ids) = prefac*(-xhi)**n*(fv(n) - sqrt(comega)*comega**n*fverf(n))
755 DEALLOCATE (fv, fverf)
773 SUBROUTINE s_vgauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
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
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
790 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
792 IF (calculate_forces) ndev = 1
793 nmax = la_max + lb_max + ndev + 1
805 eta = xhi/(xhi + omega)
809 expt = exp(-omega/(omega + xhi)*t)
810 prefac = 2.0_dp*sqrt(
pi**5/zet**3)/(xhi + omega)*expt
812 CALL fgamma(nmax - 1, etat, f)
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)
842 SUBROUTINE s_gauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
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
853 INTEGER :: ids, ipgfa, jpgfb, n, ndev, nmax
854 REAL(kind=
dp) :: a, b, eta, expt, oeta, prefac, rab2, t, &
856 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: f
859 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
861 IF (calculate_forces) ndev = 1
862 nmax = la_max + lb_max + ndev + 1
873 eta = xhi/(xhi + omega)
876 expt = exp(-omega/(omega + xhi)*t)
877 prefac =
pi**3/sqrt(zet**3)/sqrt((xhi + omega)**3)*expt
881 v(ipgfa, jpgfb, ids) = prefac*(-oeta)**n
906 swork, swork_cont, calculate_forces)
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
918 INTEGER :: ids, ids_start, ipgfa, ishella, jpgfb, &
919 jshellb, lai, lbj, ndev, nds
922 IF (calculate_forces) ndev = 1
925 DO ishella = 1, nshella
927 DO jshellb = 1, nshellb
930 ids_start = nds - min(lai, lbj)
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)
960 npgfb, nshellb, sconb, &
961 nds, swork, swork_cont)
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
971 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: work_pc
974 ALLOCATE (work_pc(npgfb, nds, nshella))
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)
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
1012 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: work_pc
1015 ALLOCATE (work_pc(npgfb, nds_max, nshella))
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)
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)
1025 DEALLOCATE (work_pc)
1051 lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, &
1052 nl_max, nds_max, swork, swork_cont, calculate_forces)
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(:, :, :, :, :), &
1066 REAL(kind=
dp),
DIMENSION(:, 0:, :, :, :), &
1067 INTENT(INOUT) :: swork_cont
1068 LOGICAL,
INTENT(IN) :: calculate_forces
1070 INTEGER :: forb, fri, ids, ids_start, iil, il, &
1071 ishella, jpgfb, jshellb, kpgfb, &
1072 kshellb, lai, lbb, lbj, lbk, ndev, &
1074 REAL(kind=
dp),
ALLOCATABLE, &
1075 DIMENSION(:, :, :, :, :) :: work_ppc
1078 IF (calculate_forces) ndev = 1
1080 ALLOCATE (work_ppc(npgfb, npgfcb, nl_max, nds_max, nshella))
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)
1085 DO kshellb = 1, nshellcb
1087 DO jshellb = 1, nshellb
1089 nl = int((lbj + lbk)/2)
1090 IF (lbj == 0 .OR. lbk == 0) nl = 0
1091 DO ishella = 1, nshella
1094 lbb = lbj + lbk - 2*il
1096 ids_start = nds - min(lai, lbb)
1098 forb = orbb_index(jpgfb, jshellb)
1099 DO kpgfb = 1, npgfcb
1100 fri = rib_index(kpgfb, kshellb)
1101 DO ids = ids_start, nds + ndev
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)
1114 DEALLOCATE (work_ppc)
1140 lca, npgfca, nshellca, orba_index, ria_index, scona_mix, &
1141 nl_max, nds_max, swork, swork_cont, calculate_forces)
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(:, :, :, :, :), &
1155 REAL(kind=
dp),
DIMENSION(:, 0:, :, :, :), &
1156 INTENT(INOUT) :: swork_cont
1157 LOGICAL,
INTENT(IN) :: calculate_forces
1159 INTEGER :: forb, fri, ids, ids_start, iil, il, &
1160 ipgfa, ishella, jshellb, kpgfa, &
1161 kshella, laa, lai, lak, lbj, ndev, &
1163 REAL(kind=
dp),
ALLOCATABLE, &
1164 DIMENSION(:, :, :, :, :) :: work_ppc
1167 IF (calculate_forces) ndev = 1
1169 ALLOCATE (work_ppc(npgfa, npgfca, nl_max, nds_max, nshellb))
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)
1174 DO kshella = 1, nshellca
1176 DO jshellb = 1, nshellb
1178 DO ishella = 1, nshella
1180 nl = int((lai + lak)/2)
1181 IF (lai == 0 .OR. lak == 0) nl = 0
1183 laa = lai + lak - 2*il
1185 ids_start = nds - min(laa, lbj)
1186 DO kpgfa = 1, npgfca
1187 fri = ria_index(kpgfa, kshella)
1189 forb = orba_index(ipgfa, ishella)
1190 DO ids = ids_start, nds + ndev
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)
1204 DEALLOCATE (work_ppc)
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.
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Routines for calculating the s-integrals and their scalar derivatives with respect to rab2 over solid...
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.
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
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
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]
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
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...
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...
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
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.
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
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...
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
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 ...