(git:374b731)
Loading...
Searching...
No Matches
ai_elec_field.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!> \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!> - coset : Cartesian orbital set pointer.
19!> - dab : Distance between the atomic centers a and b.
20!> - dac : Distance between the atomic centers a and c.
21!> - dbc : Distance between the atomic centers b and c.
22!> - l{a,b,c} : Angular momentum quantum number of shell a, b or c.
23!> - l{a,b,c}_max: Maximum angular momentum quantum number of shell a, b or c.
24!> - l{a,b,c}_min: Minimum angular momentum quantum number of shell a, b or c.
25!> - ncoset : Number of orbitals in a Cartesian orbital set.
26!> - npgf{a,b} : Degree of contraction of shell a or b.
27!> - rab : Distance vector between the atomic centers a and b.
28!> - rab2 : Square of the distance between the atomic centers a and b.
29!> - rac : Distance vector between the atomic centers a and c.
30!> - rac2 : Square of the distance between the atomic centers a and c.
31!> - rbc : Distance vector between the atomic centers b and c.
32!> - rbc2 : Square of the distance between the atomic centers b and c.
33!> - rpgf{a,b,c} : Radius of the primitive Gaussian-type function a, b or c.
34!> - zet{a,b,c} : Exponents of the Gaussian-type functions a, b or c.
35!> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
36!> \author VW
37! **************************************************************************************************
38MODULE ai_elec_field
39 USE ai_os_rr, ONLY: os_rr_coul
40 USE kinds, ONLY: dp
41 USE mathconstants, ONLY: pi
42 USE orbital_pointers, ONLY: coset,&
43 ncoset
44#include "../base/base_uses.f90"
45
46 IMPLICIT NONE
47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_elec_field'
48 PRIVATE
49
50 ! *** Public subroutines ***
51
52 PUBLIC :: efg
53
54CONTAINS
55
56! **************************************************************************************************
57!> \brief Calculation of the primitive electric field integrals over
58!> Cartesian Gaussian-type functions.
59!> \param la_max ...
60!> \param la_min ...
61!> \param npgfa ...
62!> \param rpgfa ...
63!> \param zeta ...
64!> \param lb_max ...
65!> \param lb_min ...
66!> \param npgfb ...
67!> \param rpgfb ...
68!> \param zetb ...
69!> \param rac ...
70!> \param rbc ...
71!> \param rab ...
72!> \param vab ...
73!> \param ldrr1 ...
74!> \param ldrr2 ...
75!> \param rr ...
76!> \date 02.03.2009
77!> \author VW
78!> \version 1.0
79! **************************************************************************************************
80 SUBROUTINE efg(la_max, la_min, npgfa, rpgfa, zeta, &
81 lb_max, lb_min, npgfb, rpgfb, zetb, &
82 rac, rbc, rab, vab, ldrr1, ldrr2, rr)
83 INTEGER, INTENT(IN) :: la_max, la_min, npgfa
84 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
85 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
86 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
87 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac, rbc, rab
88 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: vab
89 INTEGER, INTENT(IN) :: ldrr1, ldrr2
90 REAL(kind=dp), DIMENSION(0:ldrr1-1, ldrr2, *), &
91 INTENT(INOUT) :: rr
92
93 INTEGER :: ax, ay, az, bx, by, bz, coa, coam1x, coam1y, coam1z, coam2x, coam2y, coam2z, &
94 coamxpy, coamxpz, coamxy, coamxz, coamypx, coamypz, coamyz, coamzpx, coamzpy, coap1x, &
95 coap1y, coap1z, coap2x, coap2y, coap2z, coapxy, coapxz, coapyz, cob, cobm1x, cobm1y, &
96 cobm1z, cobm2x, cobm2y, cobm2z, cobmxpy, cobmxpz, cobmxy, cobmxz, cobmypx, cobmypz, &
97 cobmyz, cobmzpx, cobmzpy, cobp1x, cobp1y, cobp1z, cobp2x, cobp2y, cobp2z, cobpxy, cobpxz, &
98 cobpyz, i, ipgf, j, jpgf, la, lb, ma, mb, na, nb
99 REAL(kind=dp) :: dab, dum, dumxx, dumxy, dumxz, dumyy, &
100 dumyz, dumzz, f0, rab2, xhi, za, zb, &
101 zet
102 REAL(kind=dp), DIMENSION(3) :: rap, rbp, rcp
103
104! *** Calculate the distance of the centers a and c ***
105
106 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
107 dab = sqrt(rab2)
108
109 ! *** Loop over all pairs of primitive Gaussian-type functions ***
110
111 na = 0
112
113 DO ipgf = 1, npgfa
114
115 nb = 0
116
117 DO jpgf = 1, npgfb
118
119 ! *** Screening ***
120
121 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
122 DO j = nb + 1, nb + ncoset(lb_max)
123 DO i = na + 1, na + ncoset(la_max)
124 vab(i, j, 1) = 0.0_dp
125 vab(i, j, 2) = 0.0_dp
126 vab(i, j, 3) = 0.0_dp
127 vab(i, j, 4) = 0.0_dp
128 vab(i, j, 5) = 0.0_dp
129 vab(i, j, 6) = 0.0_dp
130 END DO
131 END DO
132 nb = nb + ncoset(lb_max)
133 cycle
134 END IF
135
136 ! *** Calculate some prefactors ***
137
138 za = zeta(ipgf)
139 zb = zetb(jpgf)
140 zet = za + zb
141 xhi = za*zb/zet
142 rap = zb*rab/zet
143 rbp = -za*rab/zet
144 rcp = -(za*rac + zb*rbc)/zet
145 f0 = 2.0_dp*sqrt(zet/pi)*(pi/zet)**(1.5_dp)*exp(-xhi*rab2)
146
147 ! *** Calculate the recurrence relation
148
149 CALL os_rr_coul(rap, la_max + 2, rbp, lb_max + 2, rcp, zet, ldrr1, ldrr2, rr)
150
151 ! *** Calculate the primitive electric field gradient integrals ***
152
153 DO lb = lb_min, lb_max
154 DO bx = 0, lb
155 DO by = 0, lb - bx
156 bz = lb - bx - by
157 cob = coset(bx, by, bz)
158 cobm1x = coset(max(bx - 1, 0), by, bz)
159 cobm1y = coset(bx, max(by - 1, 0), bz)
160 cobm1z = coset(bx, by, max(bz - 1, 0))
161 cobm2x = coset(max(bx - 2, 0), by, bz)
162 cobm2y = coset(bx, max(by - 2, 0), bz)
163 cobm2z = coset(bx, by, max(bz - 2, 0))
164 cobmxy = coset(max(bx - 1, 0), max(by - 1, 0), bz)
165 cobmxz = coset(max(bx - 1, 0), by, max(bz - 1, 0))
166 cobmyz = coset(bx, max(by - 1, 0), max(bz - 1, 0))
167 cobp1x = coset(bx + 1, by, bz)
168 cobp1y = coset(bx, by + 1, bz)
169 cobp1z = coset(bx, by, bz + 1)
170 cobp2x = coset(bx + 2, by, bz)
171 cobp2y = coset(bx, by + 2, bz)
172 cobp2z = coset(bx, by, bz + 2)
173 cobpxy = coset(bx + 1, by + 1, bz)
174 cobpxz = coset(bx + 1, by, bz + 1)
175 cobpyz = coset(bx, by + 1, bz + 1)
176 cobmxpy = coset(max(bx - 1, 0), by + 1, bz)
177 cobmxpz = coset(max(bx - 1, 0), by, bz + 1)
178 cobmypx = coset(bx + 1, max(by - 1, 0), bz)
179 cobmypz = coset(bx, max(by - 1, 0), bz + 1)
180 cobmzpx = coset(bx + 1, by, max(bz - 1, 0))
181 cobmzpy = coset(bx, by + 1, max(bz - 1, 0))
182 mb = nb + cob
183 DO la = la_min, la_max
184 DO ax = 0, la
185 DO ay = 0, la - ax
186 az = la - ax - ay
187 coa = coset(ax, ay, az)
188 coap1x = coset(ax + 1, ay, az)
189 coap1y = coset(ax, ay + 1, az)
190 coap1z = coset(ax, ay, az + 1)
191 coap2x = coset(ax + 2, ay, az)
192 coap2y = coset(ax, ay + 2, az)
193 coap2z = coset(ax, ay, az + 2)
194 coapxy = coset(ax + 1, ay + 1, az)
195 coapxz = coset(ax + 1, ay, az + 1)
196 coapyz = coset(ax, ay + 1, az + 1)
197 coam1x = coset(max(ax - 1, 0), ay, az)
198 coam1y = coset(ax, max(ay - 1, 0), az)
199 coam1z = coset(ax, ay, max(az - 1, 0))
200 coam2x = coset(max(ax - 2, 0), ay, az)
201 coam2y = coset(ax, max(ay - 2, 0), az)
202 coam2z = coset(ax, ay, max(az - 2, 0))
203 coamxy = coset(max(ax - 1, 0), max(ay - 1, 0), az)
204 coamxz = coset(max(ax - 1, 0), ay, max(az - 1, 0))
205 coamyz = coset(ax, max(ay - 1, 0), max(az - 1, 0))
206 coamxpy = coset(max(ax - 1, 0), ay + 1, az)
207 coamxpz = coset(max(ax - 1, 0), ay, az + 1)
208 coamypx = coset(ax + 1, max(ay - 1, 0), az)
209 coamypz = coset(ax, max(ay - 1, 0), az + 1)
210 coamzpx = coset(ax + 1, ay, max(az - 1, 0))
211 coamzpy = coset(ax, ay + 1, max(az - 1, 0))
212 ma = na + coa
213 !
214 ! (a|xx|b)
215 dum = 4.0_dp*(za**2*rr(0, coap2x, cob) + zb**2*rr(0, coa, cobp2x) &
216 & + 2.0_dp*za*zb*rr(0, coap1x, cobp1x)) &
217 - 2.0_dp*rr(0, coa, cob)*(za*real(2*ax + 1, dp) + zb*real(2*bx + 1, dp))
218 IF (ax .GT. 1) dum = dum + real(ax*(ax - 1), dp)*rr(0, coam2x, cob)
219 IF (bx .GT. 1) dum = dum + real(bx*(bx - 1), dp)*rr(0, coa, cobm2x)
220 IF (ax .GT. 0) dum = dum - 4.0_dp*zb*real(ax, dp)*rr(0, coam1x, cobp1x)
221 IF (bx .GT. 0) dum = dum - 4.0_dp*za*real(bx, dp)*rr(0, coap1x, cobm1x)
222 IF (ax .GT. 0 .AND. bx .GT. 0) dum = dum + 2.0_dp*real(ax*bx, dp)*rr(0, coam1x, cobm1x)
223 dumxx = f0*dum
224 !
225 ! (a|yy|b)
226 dum = 4.0_dp*(za**2*rr(0, coap2y, cob) + zb**2*rr(0, coa, cobp2y) &
227 & + 2.0_dp*za*zb*rr(0, coap1y, cobp1y)) &
228 - 2.0_dp*rr(0, coa, cob)*(za*real(2*ay + 1, dp) + zb*real(2*by + 1, dp))
229 IF (ay .GT. 1) dum = dum + real(ay*(ay - 1), dp)*rr(0, coam2y, cob)
230 IF (by .GT. 1) dum = dum + real(by*(by - 1), dp)*rr(0, coa, cobm2y)
231 IF (ay .GT. 0) dum = dum - 4.0_dp*zb*real(ay, dp)*rr(0, coam1y, cobp1y)
232 IF (by .GT. 0) dum = dum - 4.0_dp*za*real(by, dp)*rr(0, coap1y, cobm1y)
233 IF (ay .GT. 0 .AND. by .GT. 0) dum = dum + 2.0_dp*real(ay*by, dp)*rr(0, coam1y, cobm1y)
234 dumyy = f0*dum
235 !
236 ! (a|zz|b)
237 dum = 4.0_dp*(za**2*rr(0, coap2z, cob) + zb**2*rr(0, coa, cobp2z) &
238 & + 2.0_dp*za*zb*rr(0, coap1z, cobp1z)) &
239 - 2.0_dp*rr(0, coa, cob)*(za*real(2*az + 1, dp) + zb*real(2*bz + 1, dp))
240 IF (az .GT. 1) dum = dum + real(az*(az - 1), dp)*rr(0, coam2z, cob)
241 IF (bz .GT. 1) dum = dum + real(bz*(bz - 1), dp)*rr(0, coa, cobm2z)
242 IF (az .GT. 0) dum = dum - 4.0_dp*zb*real(az, dp)*rr(0, coam1z, cobp1z)
243 IF (bz .GT. 0) dum = dum - 4.0_dp*za*real(bz, dp)*rr(0, coap1z, cobm1z)
244 IF (az .GT. 0 .AND. bz .GT. 0) dum = dum + 2.0_dp*real(az*bz, dp)*rr(0, coam1z, cobm1z)
245 dumzz = f0*dum
246 !
247 ! (a|xy|b)
248 dum = 4.0_dp*(za**2*rr(0, coapxy, cob) + zb**2*rr(0, coa, cobpxy) &
249 & + za*zb*(rr(0, coap1x, cobp1y) + rr(0, coap1y, cobp1x)))
250 IF (ax .GT. 0) dum = dum - 2.0_dp*real(ax, dp)* &
251 & (za*rr(0, coamxpy, cob) + zb*rr(0, coam1x, cobp1y))
252 IF (ay .GT. 0) dum = dum - 2.0_dp*real(ay, dp)* &
253 & (za*rr(0, coamypx, cob) + zb*rr(0, coam1y, cobp1x))
254 IF (ax .GT. 0 .AND. ay .GT. 0) dum = dum + real(ax*ay, dp)*rr(0, coamxy, cob)
255 IF (bx .GT. 0) dum = dum - 2.0_dp*real(bx, dp)* &
256 & (zb*rr(0, coa, cobmxpy) + za*rr(0, coap1y, cobm1x))
257 IF (by .GT. 0) dum = dum - 2.0_dp*real(by, dp)* &
258 & (zb*rr(0, coa, cobmypx) + za*rr(0, coap1x, cobm1y))
259 IF (bx .GT. 0 .AND. by .GT. 0) dum = dum + real(bx*by, dp)*rr(0, coa, cobmxy)
260 IF (ax .GT. 0 .AND. by .GT. 0) dum = dum + real(ax*by, dp)*rr(0, coam1x, cobm1y)
261 IF (ay .GT. 0 .AND. bx .GT. 0) dum = dum + real(ay*bx, dp)*rr(0, coam1y, cobm1x)
262 dumxy = f0*dum
263 !
264 ! (a|xz|b)
265 dum = 4.0_dp*(za**2*rr(0, coapxz, cob) + zb**2*rr(0, coa, cobpxz) &
266 & + za*zb*(rr(0, coap1x, cobp1z) + rr(0, coap1z, cobp1x)))
267 IF (ax .GT. 0) dum = dum - 2.0_dp*real(ax, dp)* &
268 & (za*rr(0, coamxpz, cob) + zb*rr(0, coam1x, cobp1z))
269 IF (az .GT. 0) dum = dum - 2.0_dp*real(az, dp)* &
270 & (za*rr(0, coamzpx, cob) + zb*rr(0, coam1z, cobp1x))
271 IF (ax .GT. 0 .AND. az .GT. 0) dum = dum + real(ax*az, dp)*rr(0, coamxz, cob)
272 IF (bx .GT. 0) dum = dum - 2.0_dp*real(bx, dp)* &
273 & (zb*rr(0, coa, cobmxpz) + za*rr(0, coap1z, cobm1x))
274 IF (bz .GT. 0) dum = dum - 2.0_dp*real(bz, dp)* &
275 & (zb*rr(0, coa, cobmzpx) + za*rr(0, coap1x, cobm1z))
276 IF (bx .GT. 0 .AND. bz .GT. 0) dum = dum + real(bx*bz, dp)*rr(0, coa, cobmxz)
277 IF (ax .GT. 0 .AND. bz .GT. 0) dum = dum + real(ax*bz, dp)*rr(0, coam1x, cobm1z)
278 IF (az .GT. 0 .AND. bx .GT. 0) dum = dum + real(az*bx, dp)*rr(0, coam1z, cobm1x)
279 dumxz = f0*dum
280 !
281 ! (a|yz|b)
282 dum = 4.0_dp*(za**2*rr(0, coapyz, cob) + zb**2*rr(0, coa, cobpyz) &
283 & + za*zb*(rr(0, coap1y, cobp1z) + rr(0, coap1z, cobp1y)))
284 IF (ay .GT. 0) dum = dum - 2.0_dp*real(ay, dp)* &
285 & (za*rr(0, coamypz, cob) + zb*rr(0, coam1y, cobp1z))
286 IF (az .GT. 0) dum = dum - 2.0_dp*real(az, dp)* &
287 & (za*rr(0, coamzpy, cob) + zb*rr(0, coam1z, cobp1y))
288 IF (ay .GT. 0 .AND. az .GT. 0) dum = dum + real(ay*az, dp)*rr(0, coamyz, cob)
289 IF (by .GT. 0) dum = dum - 2.0_dp*real(by, dp)* &
290 & (zb*rr(0, coa, cobmypz) + za*rr(0, coap1z, cobm1y))
291 IF (bz .GT. 0) dum = dum - 2.0_dp*real(bz, dp)* &
292 & (zb*rr(0, coa, cobmzpy) + za*rr(0, coap1y, cobm1z))
293 IF (by .GT. 0 .AND. bz .GT. 0) dum = dum + real(by*bz, dp)*rr(0, coa, cobmyz)
294 IF (ay .GT. 0 .AND. bz .GT. 0) dum = dum + real(ay*bz, dp)*rr(0, coam1y, cobm1z)
295 IF (az .GT. 0 .AND. by .GT. 0) dum = dum + real(az*by, dp)*rr(0, coam1z, cobm1y)
296 dumyz = f0*dum
297 !
298 !
299 vab(ma, mb, 1) = (2.0_dp*dumxx - dumyy - dumzz)/3.0_dp !xx
300 vab(ma, mb, 2) = dumxy !xy
301 vab(ma, mb, 3) = dumxz !xz
302 vab(ma, mb, 4) = (2.0_dp*dumyy - dumzz - dumxx)/3.0_dp !yy
303 vab(ma, mb, 5) = dumyz !yz
304 vab(ma, mb, 6) = (2.0_dp*dumzz - dumxx - dumyy)/3.0_dp !zz
305 END DO
306 END DO
307 END DO !la
308
309 END DO
310 END DO
311 END DO !lb
312
313 nb = nb + ncoset(lb_max)
314
315 END DO
316
317 na = na + ncoset(la_max)
318
319 END DO
320
321 END SUBROUTINE efg
322
323END MODULE ai_elec_field
ai_elec_field
Calculation of Coulomb integrals over Cartesian Gaussian-type functions (electron repulsion integrals...
Definition ai_elec_field.F:38
ai_elec_field::efg
subroutine, public efg(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rac, rbc, rab, vab, ldrr1, ldrr2, rr)
Calculation of the primitive electric field integrals over Cartesian Gaussian-type functions.
Definition ai_elec_field.F:83
ai_os_rr
Definition ai_os_rr.F:8
ai_os_rr::os_rr_coul
subroutine, public os_rr_coul(rap, la_max, rbp, lb_max, rcp, zet, ldrr1, ldrr2, rr)
Calculation of the Obara-Saika recurrence relation for 1/r_C.
Definition ai_os_rr.F:118
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
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45