Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : !> \par Parameters
16 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
17 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
18 : !> - cx,cy,cz : Angular momentum index numbers of orbital c.
19 : !> - coset : Cartesian orbital set pointer.
20 : !> - dab : Distance between the atomic centers a and b.
21 : !> - dac : Distance between the atomic centers a and c.
22 : !> - dbc : Distance between the atomic centers b and c.
23 : !> - gccc : Prefactor of the primitive Gaussian function c.
24 : !> - l{a,b,c} : Angular momentum quantum number of shell a, b or c.
25 : !> - l{a,b,c}_max: Maximum angular momentum quantum number of shell a, b or c.
26 : !> - l{a,b,c}_min: Minimum angular momentum quantum number of shell a, b or c.
27 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
28 : !> - npgf{a,b} : Degree of contraction of shell a or b.
29 : !> - rab : Distance vector between the atomic centers a and b.
30 : !> - rab2 : Square of the distance between the atomic centers a and b.
31 : !> - rac : Distance vector between the atomic centers a and c.
32 : !> - rac2 : Square of the distance between the atomic centers a and c.
33 : !> - rbc : Distance vector between the atomic centers b and c.
34 : !> - rbc2 : Square of the distance between the atomic centers b and c.
35 : !> - rpgf{a,b,c} : Radius of the primitive Gaussian-type function a, b or c.
36 : !> - zet{a,b,c} : Exponents of the Gaussian-type functions a, b or c.
37 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
38 : !> - zetw : Reciprocal of the sum of the exponents of orbital a, b and c.
39 : !> \author Matthias Krack (22.08.2000)
40 : ! **************************************************************************************************
41 : MODULE ai_coulomb
42 :
43 : USE ai_operators_r12, ONLY: operator2_recurrence
44 : USE ai_os_rr, ONLY: os_rr_coulomb3
45 : USE gamma, ONLY: fgamma => fgamma_0
46 : USE kinds, ONLY: dp
47 : USE mathconstants, ONLY: pi
48 : USE orbital_pointers, ONLY: ncoset
49 : #include "../base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_coulomb'
53 : PRIVATE
54 :
55 : ! *** Public subroutines ***
56 :
57 : PUBLIC :: coulomb2, coulomb3
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Calculation of the primitive two-center Coulomb integrals over
63 : !> Cartesian Gaussian-type functions.
64 : !> \param la_max ...
65 : !> \param npgfa ...
66 : !> \param zeta ...
67 : !> \param rpgfa ...
68 : !> \param la_min ...
69 : !> \param lc_max ...
70 : !> \param npgfc ...
71 : !> \param zetc ...
72 : !> \param rpgfc ...
73 : !> \param lc_min ...
74 : !> \param rac ...
75 : !> \param rac2 ...
76 : !> \param vac ...
77 : !> \param v ...
78 : !> \param f ...
79 : !> \param screening optional primitive-pair screening switch
80 : !> \date 05.12.2000
81 : !> \author Matthias Krack
82 : !> \version 1.0
83 : ! **************************************************************************************************
84 3108 : SUBROUTINE coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lc_max, npgfc, zetc, rpgfc, lc_min, &
85 4662 : rac, rac2, vac, v, f, screening)
86 : INTEGER, INTENT(IN) :: la_max, npgfa
87 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
88 : INTEGER, INTENT(IN) :: la_min, lc_max, npgfc
89 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
90 : INTEGER, INTENT(IN) :: lc_min
91 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
92 : REAL(KIND=dp), INTENT(IN) :: rac2
93 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
94 : REAL(KIND=dp), DIMENSION(:, :, :) :: v
95 : REAL(KIND=dp), DIMENSION(0:) :: f
96 : LOGICAL, INTENT(IN), OPTIONAL :: screening
97 :
98 : INTEGER :: i, ipgf, j, jpgf, n, na, nap, nc, ncp, &
99 : nmax
100 : LOGICAL :: do_screening
101 : REAL(KIND=dp) :: dac, f0, rho, t, zetp, zetq, zetw
102 :
103 1554 : do_screening = .TRUE.
104 1554 : IF (PRESENT(screening)) do_screening = screening
105 :
106 5324566 : v = 0.0_dp
107 :
108 1554 : nmax = la_max + lc_max + 1
109 :
110 1554 : dac = SQRT(rac2)
111 :
112 1554 : na = 0
113 1554 : nap = 0
114 3558 : DO ipgf = 1, npgfa
115 :
116 2004 : nc = 0
117 2004 : ncp = 0
118 :
119 5808 : DO jpgf = 1, npgfc
120 :
121 3804 : IF (do_screening .AND. rpgfa(ipgf) + rpgfc(jpgf) < dac) THEN
122 0 : DO j = nc + ncoset(lc_min - 1) + 1, nc + ncoset(lc_max)
123 0 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
124 0 : vac(i, j) = 0.0_dp
125 : END DO
126 : END DO
127 0 : nc = nc + ncoset(lc_max)
128 0 : CYCLE
129 : END IF
130 :
131 3804 : zetp = 1.0_dp/zeta(ipgf)
132 3804 : zetq = 1.0_dp/zetc(jpgf)
133 3804 : zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
134 :
135 3804 : rho = zeta(ipgf)*zetc(jpgf)*zetw
136 :
137 3804 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
138 :
139 3804 : t = rho*rac2
140 :
141 3804 : CALL fgamma(nmax - 1, t, f)
142 :
143 13367 : DO n = 1, nmax
144 13367 : v(1, 1, n) = f0*f(n - 1)
145 : END DO
146 :
147 : CALL operator2_recurrence(la_max, la_min, lc_max, lc_min, zeta(ipgf), zetc(jpgf), &
148 5808 : zetp, zetq, zetw, rho, rac, vac, v, na, nc, nap, ncp)
149 :
150 : END DO
151 :
152 2004 : na = na + ncoset(la_max)
153 3558 : nap = nap + ncoset(la_max)
154 :
155 : END DO
156 :
157 1554 : END SUBROUTINE coulomb2
158 : ! **************************************************************************************************
159 : !> \brief Calculation of the primitive three-center Coulomb integrals over
160 : !> Cartesian Gaussian-type functions (electron repulsion integrals,
161 : !> ERIs).
162 : !> \param la_max ...
163 : !> \param npgfa ...
164 : !> \param zeta ...
165 : !> \param rpgfa ...
166 : !> \param la_min ...
167 : !> \param lb_max ...
168 : !> \param npgfb ...
169 : !> \param zetb ...
170 : !> \param rpgfb ...
171 : !> \param lb_min ...
172 : !> \param lc_max ...
173 : !> \param zetc ...
174 : !> \param rpgfc ...
175 : !> \param lc_min ...
176 : !> \param gccc ...
177 : !> \param rab ...
178 : !> \param rab2 ...
179 : !> \param rac ...
180 : !> \param rac2 ...
181 : !> \param rbc2 ...
182 : !> \param vabc ...
183 : !> \param int_abc ...
184 : !> \param v ...
185 : !> \param f ...
186 : !> \param maxder ...
187 : !> \param vabc_plus ...
188 : !> \date 06.11.2000
189 : !> \author Matthias Krack
190 : !> \version 1.0
191 : ! **************************************************************************************************
192 1620 : SUBROUTINE coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
193 3240 : lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, &
194 4860 : v, f, maxder, vabc_plus)
195 : INTEGER, INTENT(IN) :: la_max, npgfa
196 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
197 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
198 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
199 : INTEGER, INTENT(IN) :: lb_min, lc_max
200 : REAL(KIND=dp), INTENT(IN) :: zetc, rpgfc
201 : INTEGER, INTENT(IN) :: lc_min
202 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: gccc
203 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
204 : REAL(KIND=dp), INTENT(IN) :: rab2
205 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
206 : REAL(KIND=dp), INTENT(IN) :: rac2, rbc2
207 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vabc
208 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: int_abc
209 : REAL(KIND=dp), DIMENSION(:, :, :, :) :: v
210 : REAL(KIND=dp), DIMENSION(0:) :: f
211 : INTEGER, INTENT(IN), OPTIONAL :: maxder
212 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: vabc_plus
213 :
214 : INTEGER :: i, ipgf, j, jpgf, k, kk, maxder_local, &
215 : n, na, nap, nb, nmax
216 : REAL(dp) :: dab, dac, dbc, f0, f1, f4, rap(3), &
217 : rbp(3), rcp(3), t, zetp, zetq, zetw
218 :
219 1374540 : v = 0.0_dp
220 :
221 1620 : maxder_local = 0
222 1620 : IF (PRESENT(maxder)) THEN
223 0 : maxder_local = maxder
224 : END IF
225 :
226 1620 : nmax = la_max + lb_max + lc_max + 1
227 :
228 : ! *** Calculate the distances of the centers a, b and c ***
229 :
230 1620 : dab = SQRT(rab2)
231 1620 : dac = SQRT(rac2)
232 1620 : dbc = SQRT(rbc2)
233 :
234 : ! *** Initialize integrals array
235 168570 : int_abc = 0.0_dp
236 :
237 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
238 :
239 1620 : na = 0
240 1620 : nap = 0
241 :
242 4320 : DO ipgf = 1, npgfa
243 :
244 : ! *** Screening ***
245 2700 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
246 0 : na = na + ncoset(la_max - maxder_local)
247 0 : nap = nap + ncoset(la_max)
248 0 : CYCLE
249 : END IF
250 :
251 2700 : nb = 0
252 :
253 7200 : DO jpgf = 1, npgfb
254 :
255 : ! *** Screening ***
256 : IF ( &
257 4500 : (rpgfb(jpgf) + rpgfc < dbc) .OR. &
258 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
259 0 : nb = nb + ncoset(lb_max)
260 0 : CYCLE
261 : END IF
262 :
263 : ! *** Calculate some prefactors ***
264 :
265 4500 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
266 4500 : zetq = 1.0_dp/zetc
267 4500 : zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
268 :
269 4500 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
270 4500 : f1 = zetb(jpgf)*zetp
271 4500 : f4 = -zetc*zetw
272 :
273 4500 : f0 = f0*EXP(-zeta(ipgf)*f1*rab2)
274 :
275 18000 : rap(:) = f1*rab(:)
276 18000 : rcp(:) = rap(:) - rac(:)
277 :
278 : ! *** Calculate the incomplete Gamma function ***
279 :
280 4500 : t = -f4*(rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3))/zetp
281 :
282 4500 : CALL fgamma(nmax - 1, t, f)
283 :
284 : ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
285 :
286 18300 : DO n = 1, nmax
287 18300 : v(1, 1, 1, n) = f0*f(n - 1)
288 : END DO
289 :
290 18000 : rbp = rap - rab
291 4500 : CALL os_rr_coulomb3(rap, la_max, rbp, lb_max, rcp, lc_max, zetp, zetq, zetw, v)
292 :
293 : ! *** Add the contribution of the current pair ***
294 : ! *** of primitive Gaussian-type functions ***
295 :
296 20250 : DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
297 15750 : kk = k - ncoset(lc_min - 1)
298 58050 : DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
299 152145 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
300 98595 : vabc(na + i, nb + j) = vabc(na + i, nb + j) + gccc(kk)*v(i, j, k, 1)
301 136395 : int_abc(na + i, nb + j, kk) = v(i, j, k, 1)
302 : END DO
303 : END DO
304 : END DO
305 :
306 4500 : IF (PRESENT(maxder)) THEN
307 0 : DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
308 0 : kk = k - ncoset(lc_min - 1)
309 0 : DO j = 1, ncoset(lb_max)
310 0 : DO i = 1, ncoset(la_max)
311 0 : vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) + gccc(kk)*v(i, j, k, 1)
312 : END DO
313 : END DO
314 : END DO
315 : END IF
316 :
317 7200 : nb = nb + ncoset(lb_max)
318 :
319 : END DO
320 :
321 2700 : na = na + ncoset(la_max - maxder_local)
322 4320 : nap = nap + ncoset(la_max)
323 :
324 : END DO
325 :
326 1620 : END SUBROUTINE coulomb3
327 :
328 : END MODULE ai_coulomb
|
