Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : ! **************************************************************************************************
17 : MODULE ai_eri_debug
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 0 : 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 0 : xa = ya
63 0 : xb = yb
64 0 : xc = yc
65 0 : xd = yd
66 0 : A = rA
67 0 : B = rB
68 0 : C = rC
69 0 : D = rD
70 :
71 0 : xsi = xa + xb
72 0 : eta = xc + xd
73 :
74 0 : P = (xa*A + xb*B)/xsi
75 0 : Q = (xc*C + xd*D)/eta
76 0 : W = (xsi*P + eta*Q)/(xsi + eta)
77 :
78 0 : rho = xsi*eta/(xsi + eta)
79 :
80 0 : T = rho*SUM((P - Q)**2)
81 :
82 0 : fm = fgamma_ref(4*lmax, T)
83 :
84 0 : eab = -xa*xb/xsi*SUM((A - B)**2)
85 0 : kab = SQRT(2._dp)*pi**1.25_dp/xsi*EXP(eab)
86 :
87 0 : ecd = -xc*xd/eta*SUM((C - D)**2)
88 0 : kcd = SQRT(2._dp)*pi**1.25_dp/eta*EXP(ecd)
89 :
90 0 : I0M = kab*kcd/SQRT(xsi + eta)*fm
91 :
92 0 : 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 0 : 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 0 : m = 0
116 0 : IF (PRESENT(mi)) m = mi
117 :
118 0 : IABCD = 0._dp
119 0 : IF (ANY(an < 0)) RETURN
120 0 : IF (ANY(bn < 0)) RETURN
121 0 : IF (ANY(cn < 0)) RETURN
122 0 : IF (ANY(dn < 0)) RETURN
123 :
124 0 : IF (SUM(an + bn + cn + dn) == 0) THEN
125 0 : IABCD = I0M(m)
126 0 : RETURN
127 : END IF
128 :
129 0 : IF (dn(1) > 0) THEN
130 0 : IABCD = os(an, bn, cn + i1, dn - i1) - (D(1) - C(1))*os(an, bn, cn, dn - i1)
131 0 : ELSEIF (dn(2) > 0) THEN
132 0 : IABCD = os(an, bn, cn + i2, dn - i2) - (D(2) - C(2))*os(an, bn, cn, dn - i2)
133 0 : ELSEIF (dn(3) > 0) THEN
134 0 : IABCD = os(an, bn, cn + i3, dn - i3) - (D(3) - C(3))*os(an, bn, cn, dn - i3)
135 : ELSE
136 0 : IF (bn(1) > 0) THEN
137 0 : IABCD = os(an + i1, bn - i1, cn, dn) - (B(1) - A(1))*os(an, bn - i1, cn, dn)
138 0 : ELSEIF (bn(2) > 0) THEN
139 0 : IABCD = os(an + i2, bn - i2, cn, dn) - (B(2) - A(2))*os(an, bn - i2, cn, dn)
140 0 : ELSEIF (bn(3) > 0) THEN
141 0 : IABCD = os(an + i3, bn - i3, cn, dn) - (B(3) - A(3))*os(an, bn - i3, cn, dn)
142 : ELSE
143 0 : 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 0 : xsi/eta*os(an + i1, bn, cn - i1, dn)
148 0 : 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 0 : xsi/eta*os(an + i2, bn, cn - i2, dn)
153 0 : 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 0 : xsi/eta*os(an + i3, bn, cn - i3, dn)
158 : ELSE
159 0 : 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 : 0.5_dp*(an(1) - 1)/xsi*rho/xsi*os(an - i1 - i1, bn, cn, dn, m + 1)
164 0 : 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 : 0.5_dp*(an(2) - 1)/xsi*rho/xsi*os(an - i2 - i2, bn, cn, dn, m + 1)
169 0 : 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 : 0.5_dp*(an(3) - 1)/xsi*rho/xsi*os(an - i3 - i3, bn, cn, dn, m + 1)
174 : ELSE
175 0 : 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
|
