(git:0de0cc2)
ai_eri_debug.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 Calculation of Coulomb integrals over Cartesian Gaussian-type functions
10 !> (electron repulsion integrals, ERIs).
11 !> \par Literature
12 !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 !> \par History
14 !> none
15 !> \author J. Hutter (07.2009)
16 ! **************************************************************************************************
18 
19  USE gamma, ONLY: fgamma_ref
20  USE kinds, ONLY: dp
21  USE mathconstants, ONLY: pi
22 #include "../base/base_uses.f90"
23 
24  IMPLICIT NONE
25 
26  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_eri_debug'
27 
28  INTEGER, PARAMETER :: lmax = 5
29 
30  REAL(dp) :: xa, xb, xc, xd
31  REAL(dp), DIMENSION(3) :: a, b, c, d
32  REAL(dp), DIMENSION(3) :: p, q, w
33  REAL(dp) :: xsi, eta, rho, t
34 
35  REAL(dp), DIMENSION(0:4*lmax) :: fm, i0m
36 
37  PRIVATE
38 
39 CONTAINS
40 
41 ! **************************************************************************************************
42 !> \brief Calculation of the primitive two-center Coulomb integrals over
43 !> Cartesian Gaussian-type functions.
44 !> \param ya ...
45 !> \param yb ...
46 !> \param yc ...
47 !> \param yd ...
48 !> \param rA ...
49 !> \param rB ...
50 !> \param rC ...
51 !> \param rD ...
52 !> \date 07.2009
53 !> \author J. Hutter
54 !> \version 1.0
55 ! **************************************************************************************************
56  SUBROUTINE init_os(ya, yb, yc, yd, rA, rB, rC, rD)
57  REAL(dp) :: ya, yb, yc, yd
58  REAL(dp), DIMENSION(3) :: ra, rb, rc, rd
59 
60  REAL(dp) :: eab, ecd, kab, kcd
61 
62  xa = ya
63  xb = yb
64  xc = yc
65  xd = yd
66  a = ra
67  b = rb
68  c = rc
69  d = rd
70 
71  xsi = xa + xb
72  eta = xc + xd
73 
74  p = (xa*a + xb*b)/xsi
75  q = (xc*c + xd*d)/eta
76  w = (xsi*p + eta*q)/(xsi + eta)
77 
78  rho = xsi*eta/(xsi + eta)
79 
80  t = rho*sum((p - q)**2)
81 
82  fm = fgamma_ref(4*lmax, t)
83 
84  eab = -xa*xb/xsi*sum((a - b)**2)
85  kab = sqrt(2._dp)*pi**1.25_dp/xsi*exp(eab)
86 
87  ecd = -xc*xd/eta*sum((c - d)**2)
88  kcd = sqrt(2._dp)*pi**1.25_dp/eta*exp(ecd)
89 
90  i0m = kab*kcd/sqrt(xsi + eta)*fm
91 
92  END SUBROUTINE init_os
93 
94 ! **************************************************************************************************
95 
96 ! **************************************************************************************************
97 !> \brief ...
98 !> \param an ...
99 !> \param bn ...
100 !> \param cn ...
101 !> \param dn ...
102 !> \param mi ...
103 !> \return ...
104 ! **************************************************************************************************
105  RECURSIVE FUNCTION os(an, bn, cn, dn, mi) RESULT(IABCD)
106  INTEGER, DIMENSION(3) :: an, bn, cn, dn
107  INTEGER, OPTIONAL :: mi
108  REAL(dp) :: iabcd
109 
110  INTEGER, DIMENSION(3), PARAMETER :: i1 = (/1, 0, 0/), i2 = (/0, 1, 0/), &
111  i3 = (/0, 0, 1/)
112 
113  INTEGER :: m
114 
115  m = 0
116  IF (PRESENT(mi)) m = mi
117 
118  iabcd = 0._dp
119  IF (any(an < 0)) RETURN
120  IF (any(bn < 0)) RETURN
121  IF (any(cn < 0)) RETURN
122  IF (any(dn < 0)) RETURN
123 
124  IF (sum(an + bn + cn + dn) == 0) THEN
125  iabcd = i0m(m)
126  RETURN
127  END IF
128 
129  IF (dn(1) > 0) THEN
130  iabcd = os(an, bn, cn + i1, dn - i1) - (d(1) - c(1))*os(an, bn, cn, dn - i1)
131  ELSEIF (dn(2) > 0) THEN
132  iabcd = os(an, bn, cn + i2, dn - i2) - (d(2) - c(2))*os(an, bn, cn, dn - i2)
133  ELSEIF (dn(3) > 0) THEN
134  iabcd = os(an, bn, cn + i3, dn - i3) - (d(3) - c(3))*os(an, bn, cn, dn - i3)
135  ELSE
136  IF (bn(1) > 0) THEN
137  iabcd = os(an + i1, bn - i1, cn, dn) - (b(1) - a(1))*os(an, bn - i1, cn, dn)
138  ELSEIF (bn(2) > 0) THEN
139  iabcd = os(an + i2, bn - i2, cn, dn) - (b(2) - a(2))*os(an, bn - i2, cn, dn)
140  ELSEIF (bn(3) > 0) THEN
141  iabcd = os(an + i3, bn - i3, cn, dn) - (b(3) - a(3))*os(an, bn - i3, cn, dn)
142  ELSE
143  IF (cn(1) > 0) THEN
144  iabcd = ((q(1) - c(1)) + xsi/eta*(p(1) - a(1)))*os(an, bn, cn - i1, dn) + &
145  0.5_dp*an(1)/eta*os(an - i1, bn, cn - i1, dn) + &
146  0.5_dp*(cn(1) - 1)/eta*os(an, bn, cn - i1 - i1, dn) - &
147  xsi/eta*os(an + i1, bn, cn - i1, dn)
148  ELSEIF (cn(2) > 0) THEN
149  iabcd = ((q(2) - c(2)) + xsi/eta*(p(2) - a(2)))*os(an, bn, cn - i2, dn) + &
150  0.5_dp*an(2)/eta*os(an - i2, bn, cn - i2, dn) + &
151  0.5_dp*(cn(2) - 1)/eta*os(an, bn, cn - i2 - i2, dn) - &
152  xsi/eta*os(an + i2, bn, cn - i2, dn)
153  ELSEIF (cn(3) > 0) THEN
154  iabcd = ((q(3) - c(3)) + xsi/eta*(p(3) - a(3)))*os(an, bn, cn - i3, dn) + &
155  0.5_dp*an(3)/eta*os(an - i3, bn, cn - i3, dn) + &
156  0.5_dp*(cn(3) - 1)/eta*os(an, bn, cn - i3 - i3, dn) - &
157  xsi/eta*os(an + i3, bn, cn - i3, dn)
158  ELSE
159  IF (an(1) > 0) THEN
160  iabcd = (p(1) - a(1))*os(an - i1, bn, cn, dn, m) + &
161  (w(1) - p(1))*os(an - i1, bn, cn, dn, m + 1) + &
162  0.5_dp*(an(1) - 1)/xsi*os(an - i1 - i1, bn, cn, dn, m) - &
163  0.5_dp*(an(1) - 1)/xsi*rho/xsi*os(an - i1 - i1, bn, cn, dn, m + 1)
164  ELSEIF (an(2) > 0) THEN
165  iabcd = (p(2) - a(2))*os(an - i2, bn, cn, dn, m) + &
166  (w(2) - p(2))*os(an - i2, bn, cn, dn, m + 1) + &
167  0.5_dp*(an(2) - 1)/xsi*os(an - i2 - i2, bn, cn, dn, m) - &
168  0.5_dp*(an(2) - 1)/xsi*rho/xsi*os(an - i2 - i2, bn, cn, dn, m + 1)
169  ELSEIF (an(3) > 0) THEN
170  iabcd = (p(3) - a(3))*os(an - i3, bn, cn, dn, m) + &
171  (w(3) - p(3))*os(an - i3, bn, cn, dn, m + 1) + &
172  0.5_dp*(an(3) - 1)/xsi*os(an - i3 - i3, bn, cn, dn, m) - &
173  0.5_dp*(an(3) - 1)/xsi*rho/xsi*os(an - i3 - i3, bn, cn, dn, m + 1)
174  ELSE
175  cpabort("I(0000)")
176  END IF
177  END IF
178  END IF
179  END IF
180 
181  END FUNCTION os
182 
183 ! **************************************************************************************************
184 
185 END MODULE ai_eri_debug
Calculation of Coulomb integrals over Cartesian Gaussian-type functions (electron repulsion integrals...
Definition: ai_eri_debug.F:17
real(dp) xa
Definition: ai_eri_debug.F:30
real(dp) xb
Definition: ai_eri_debug.F:30
real(dp), dimension(0:4 *lmax) i0m
Definition: ai_eri_debug.F:35
real(dp), dimension(3) w
Definition: ai_eri_debug.F:32
real(dp) xsi
Definition: ai_eri_debug.F:33
real(dp), dimension(3) c
Definition: ai_eri_debug.F:31
real(dp) t
Definition: ai_eri_debug.F:33
real(dp), dimension(3) a
Definition: ai_eri_debug.F:31
real(dp), dimension(3) d
Definition: ai_eri_debug.F:31
real(dp) rho
Definition: ai_eri_debug.F:33
real(dp), dimension(3) q
Definition: ai_eri_debug.F:32
real(dp) xd
Definition: ai_eri_debug.F:30
real(dp), dimension(3) b
Definition: ai_eri_debug.F:31
integer, parameter lmax
Definition: ai_eri_debug.F:28
real(dp), dimension(0:4 *lmax) fm
Definition: ai_eri_debug.F:35
real(dp) eta
Definition: ai_eri_debug.F:33
real(dp), dimension(3) p
Definition: ai_eri_debug.F:32
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
real(kind=dp) function, dimension(0:nmax), public fgamma_ref(nmax, t)
Calculation of the incomplete Gamma function F_n(t) using a spherical Bessel function expansion....
Definition: gamma.F:423
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
Exchange and Correlation functional calculations.
Definition: xc.F:17