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 : !> \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 : ! **************************************************************************************************
38 : MODULE 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 :
54 : CONTAINS
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 9222 : SUBROUTINE efg(la_max, la_min, npgfa, rpgfa, zeta, &
81 18444 : lb_max, lb_min, npgfb, rpgfb, zetb, &
82 9222 : 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 9222 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
107 9222 : dab = SQRT(rab2)
108 :
109 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
110 :
111 9222 : na = 0
112 :
113 21050 : DO ipgf = 1, npgfa
114 :
115 11828 : nb = 0
116 :
117 27026 : DO jpgf = 1, npgfb
118 :
119 : ! *** Screening ***
120 :
121 15198 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
122 6888 : DO j = nb + 1, nb + ncoset(lb_max)
123 12801 : DO i = na + 1, na + ncoset(la_max)
124 5913 : vab(i, j, 1) = 0.0_dp
125 5913 : vab(i, j, 2) = 0.0_dp
126 5913 : vab(i, j, 3) = 0.0_dp
127 5913 : vab(i, j, 4) = 0.0_dp
128 5913 : vab(i, j, 5) = 0.0_dp
129 9999 : vab(i, j, 6) = 0.0_dp
130 : END DO
131 : END DO
132 2802 : nb = nb + ncoset(lb_max)
133 2802 : CYCLE
134 : END IF
135 :
136 : ! *** Calculate some prefactors ***
137 :
138 12396 : za = zeta(ipgf)
139 12396 : zb = zetb(jpgf)
140 12396 : zet = za + zb
141 12396 : xhi = za*zb/zet
142 49584 : rap = zb*rab/zet
143 49584 : rbp = -za*rab/zet
144 49584 : rcp = -(za*rac + zb*rbc)/zet
145 12396 : f0 = 2.0_dp*SQRT(zet/pi)*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
146 :
147 : ! *** Calculate the recurrence relation
148 :
149 12396 : 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 27668 : DO lb = lb_min, lb_max
154 45827 : DO bx = 0, lb
155 54477 : DO by = 0, lb - bx
156 21046 : bz = lb - bx - by
157 21046 : cob = coset(bx, by, bz)
158 21046 : cobm1x = coset(MAX(bx - 1, 0), by, bz)
159 21046 : cobm1y = coset(bx, MAX(by - 1, 0), bz)
160 21046 : cobm1z = coset(bx, by, MAX(bz - 1, 0))
161 21046 : cobm2x = coset(MAX(bx - 2, 0), by, bz)
162 21046 : cobm2y = coset(bx, MAX(by - 2, 0), bz)
163 21046 : cobm2z = coset(bx, by, MAX(bz - 2, 0))
164 21046 : cobmxy = coset(MAX(bx - 1, 0), MAX(by - 1, 0), bz)
165 21046 : cobmxz = coset(MAX(bx - 1, 0), by, MAX(bz - 1, 0))
166 21046 : cobmyz = coset(bx, MAX(by - 1, 0), MAX(bz - 1, 0))
167 21046 : cobp1x = coset(bx + 1, by, bz)
168 21046 : cobp1y = coset(bx, by + 1, bz)
169 21046 : cobp1z = coset(bx, by, bz + 1)
170 21046 : cobp2x = coset(bx + 2, by, bz)
171 21046 : cobp2y = coset(bx, by + 2, bz)
172 21046 : cobp2z = coset(bx, by, bz + 2)
173 21046 : cobpxy = coset(bx + 1, by + 1, bz)
174 21046 : cobpxz = coset(bx + 1, by, bz + 1)
175 21046 : cobpyz = coset(bx, by + 1, bz + 1)
176 21046 : cobmxpy = coset(MAX(bx - 1, 0), by + 1, bz)
177 21046 : cobmxpz = coset(MAX(bx - 1, 0), by, bz + 1)
178 21046 : cobmypx = coset(bx + 1, MAX(by - 1, 0), bz)
179 21046 : cobmypz = coset(bx, MAX(by - 1, 0), bz + 1)
180 21046 : cobmzpx = coset(bx + 1, by, MAX(bz - 1, 0))
181 21046 : cobmzpy = coset(bx, by + 1, MAX(bz - 1, 0))
182 21046 : mb = nb + cob
183 65949 : DO la = la_min, la_max
184 80245 : DO ax = 0, la
185 97365 : DO ay = 0, la - ax
186 38166 : az = la - ax - ay
187 38166 : coa = coset(ax, ay, az)
188 38166 : coap1x = coset(ax + 1, ay, az)
189 38166 : coap1y = coset(ax, ay + 1, az)
190 38166 : coap1z = coset(ax, ay, az + 1)
191 38166 : coap2x = coset(ax + 2, ay, az)
192 38166 : coap2y = coset(ax, ay + 2, az)
193 38166 : coap2z = coset(ax, ay, az + 2)
194 38166 : coapxy = coset(ax + 1, ay + 1, az)
195 38166 : coapxz = coset(ax + 1, ay, az + 1)
196 38166 : coapyz = coset(ax, ay + 1, az + 1)
197 38166 : coam1x = coset(MAX(ax - 1, 0), ay, az)
198 38166 : coam1y = coset(ax, MAX(ay - 1, 0), az)
199 38166 : coam1z = coset(ax, ay, MAX(az - 1, 0))
200 38166 : coam2x = coset(MAX(ax - 2, 0), ay, az)
201 38166 : coam2y = coset(ax, MAX(ay - 2, 0), az)
202 38166 : coam2z = coset(ax, ay, MAX(az - 2, 0))
203 38166 : coamxy = coset(MAX(ax - 1, 0), MAX(ay - 1, 0), az)
204 38166 : coamxz = coset(MAX(ax - 1, 0), ay, MAX(az - 1, 0))
205 38166 : coamyz = coset(ax, MAX(ay - 1, 0), MAX(az - 1, 0))
206 38166 : coamxpy = coset(MAX(ax - 1, 0), ay + 1, az)
207 38166 : coamxpz = coset(MAX(ax - 1, 0), ay, az + 1)
208 38166 : coamypx = coset(ax + 1, MAX(ay - 1, 0), az)
209 38166 : coamypz = coset(ax, MAX(ay - 1, 0), az + 1)
210 38166 : coamzpx = coset(ax + 1, ay, MAX(az - 1, 0))
211 38166 : coamzpy = coset(ax, ay + 1, MAX(az - 1, 0))
212 38166 : 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 38166 : - 2.0_dp*rr(0, coa, cob)*(za*REAL(2*ax + 1, dp) + zb*REAL(2*bx + 1, dp))
218 38166 : IF (ax > 1) dum = dum + REAL(ax*(ax - 1), dp)*rr(0, coam2x, cob)
219 38166 : IF (bx > 1) dum = dum + REAL(bx*(bx - 1), dp)*rr(0, coa, cobm2x)
220 38166 : IF (ax > 0) dum = dum - 4.0_dp*zb*REAL(ax, dp)*rr(0, coam1x, cobp1x)
221 38166 : IF (bx > 0) dum = dum - 4.0_dp*za*REAL(bx, dp)*rr(0, coap1x, cobm1x)
222 38166 : IF (ax > 0 .AND. bx > 0) dum = dum + 2.0_dp*REAL(ax*bx, dp)*rr(0, coam1x, cobm1x)
223 38166 : 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 38166 : - 2.0_dp*rr(0, coa, cob)*(za*REAL(2*ay + 1, dp) + zb*REAL(2*by + 1, dp))
229 38166 : IF (ay > 1) dum = dum + REAL(ay*(ay - 1), dp)*rr(0, coam2y, cob)
230 38166 : IF (by > 1) dum = dum + REAL(by*(by - 1), dp)*rr(0, coa, cobm2y)
231 38166 : IF (ay > 0) dum = dum - 4.0_dp*zb*REAL(ay, dp)*rr(0, coam1y, cobp1y)
232 38166 : IF (by > 0) dum = dum - 4.0_dp*za*REAL(by, dp)*rr(0, coap1y, cobm1y)
233 38166 : IF (ay > 0 .AND. by > 0) dum = dum + 2.0_dp*REAL(ay*by, dp)*rr(0, coam1y, cobm1y)
234 38166 : 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 38166 : - 2.0_dp*rr(0, coa, cob)*(za*REAL(2*az + 1, dp) + zb*REAL(2*bz + 1, dp))
240 38166 : IF (az > 1) dum = dum + REAL(az*(az - 1), dp)*rr(0, coam2z, cob)
241 38166 : IF (bz > 1) dum = dum + REAL(bz*(bz - 1), dp)*rr(0, coa, cobm2z)
242 38166 : IF (az > 0) dum = dum - 4.0_dp*zb*REAL(az, dp)*rr(0, coam1z, cobp1z)
243 38166 : IF (bz > 0) dum = dum - 4.0_dp*za*REAL(bz, dp)*rr(0, coap1z, cobm1z)
244 38166 : IF (az > 0 .AND. bz > 0) dum = dum + 2.0_dp*REAL(az*bz, dp)*rr(0, coam1z, cobm1z)
245 38166 : 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 38166 : & + za*zb*(rr(0, coap1x, cobp1y) + rr(0, coap1y, cobp1x)))
250 38166 : IF (ax > 0) dum = dum - 2.0_dp*REAL(ax, dp)* &
251 5711 : & (za*rr(0, coamxpy, cob) + zb*rr(0, coam1x, cobp1y))
252 38166 : IF (ay > 0) dum = dum - 2.0_dp*REAL(ay, dp)* &
253 5711 : & (za*rr(0, coamypx, cob) + zb*rr(0, coam1y, cobp1x))
254 38166 : IF (ax > 0 .AND. ay > 0) dum = dum + REAL(ax*ay, dp)*rr(0, coamxy, cob)
255 38166 : IF (bx > 0) dum = dum - 2.0_dp*REAL(bx, dp)* &
256 5898 : & (zb*rr(0, coa, cobmxpy) + za*rr(0, coap1y, cobm1x))
257 38166 : IF (by > 0) dum = dum - 2.0_dp*REAL(by, dp)* &
258 5898 : & (zb*rr(0, coa, cobmypx) + za*rr(0, coap1x, cobm1y))
259 38166 : IF (bx > 0 .AND. by > 0) dum = dum + REAL(bx*by, dp)*rr(0, coa, cobmxy)
260 38166 : IF (ax > 0 .AND. by > 0) dum = dum + REAL(ax*by, dp)*rr(0, coam1x, cobm1y)
261 38166 : IF (ay > 0 .AND. bx > 0) dum = dum + REAL(ay*bx, dp)*rr(0, coam1y, cobm1x)
262 38166 : 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 38166 : & + za*zb*(rr(0, coap1x, cobp1z) + rr(0, coap1z, cobp1x)))
267 38166 : IF (ax > 0) dum = dum - 2.0_dp*REAL(ax, dp)* &
268 5711 : & (za*rr(0, coamxpz, cob) + zb*rr(0, coam1x, cobp1z))
269 38166 : IF (az > 0) dum = dum - 2.0_dp*REAL(az, dp)* &
270 5711 : & (za*rr(0, coamzpx, cob) + zb*rr(0, coam1z, cobp1x))
271 38166 : IF (ax > 0 .AND. az > 0) dum = dum + REAL(ax*az, dp)*rr(0, coamxz, cob)
272 38166 : IF (bx > 0) dum = dum - 2.0_dp*REAL(bx, dp)* &
273 5898 : & (zb*rr(0, coa, cobmxpz) + za*rr(0, coap1z, cobm1x))
274 38166 : IF (bz > 0) dum = dum - 2.0_dp*REAL(bz, dp)* &
275 5898 : & (zb*rr(0, coa, cobmzpx) + za*rr(0, coap1x, cobm1z))
276 38166 : IF (bx > 0 .AND. bz > 0) dum = dum + REAL(bx*bz, dp)*rr(0, coa, cobmxz)
277 38166 : IF (ax > 0 .AND. bz > 0) dum = dum + REAL(ax*bz, dp)*rr(0, coam1x, cobm1z)
278 38166 : IF (az > 0 .AND. bx > 0) dum = dum + REAL(az*bx, dp)*rr(0, coam1z, cobm1x)
279 38166 : 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 38166 : & + za*zb*(rr(0, coap1y, cobp1z) + rr(0, coap1z, cobp1y)))
284 38166 : IF (ay > 0) dum = dum - 2.0_dp*REAL(ay, dp)* &
285 5711 : & (za*rr(0, coamypz, cob) + zb*rr(0, coam1y, cobp1z))
286 38166 : IF (az > 0) dum = dum - 2.0_dp*REAL(az, dp)* &
287 5711 : & (za*rr(0, coamzpy, cob) + zb*rr(0, coam1z, cobp1y))
288 38166 : IF (ay > 0 .AND. az > 0) dum = dum + REAL(ay*az, dp)*rr(0, coamyz, cob)
289 38166 : IF (by > 0) dum = dum - 2.0_dp*REAL(by, dp)* &
290 5898 : & (zb*rr(0, coa, cobmypz) + za*rr(0, coap1z, cobm1y))
291 38166 : IF (bz > 0) dum = dum - 2.0_dp*REAL(bz, dp)* &
292 5898 : & (zb*rr(0, coa, cobmzpy) + za*rr(0, coap1y, cobm1z))
293 38166 : IF (by > 0 .AND. bz > 0) dum = dum + REAL(by*bz, dp)*rr(0, coa, cobmyz)
294 38166 : IF (ay > 0 .AND. bz > 0) dum = dum + REAL(ay*bz, dp)*rr(0, coam1y, cobm1z)
295 38166 : IF (az > 0 .AND. by > 0) dum = dum + REAL(az*by, dp)*rr(0, coam1z, cobm1y)
296 38166 : dumyz = f0*dum
297 : !
298 : !
299 38166 : vab(ma, mb, 1) = (2.0_dp*dumxx - dumyy - dumzz)/3.0_dp !xx
300 38166 : vab(ma, mb, 2) = dumxy !xy
301 38166 : vab(ma, mb, 3) = dumxz !xz
302 38166 : vab(ma, mb, 4) = (2.0_dp*dumyy - dumzz - dumxx)/3.0_dp !yy
303 38166 : vab(ma, mb, 5) = dumyz !yz
304 70621 : 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 24224 : nb = nb + ncoset(lb_max)
314 :
315 : END DO
316 :
317 21050 : na = na + ncoset(la_max)
318 :
319 : END DO
320 :
321 9222 : END SUBROUTINE efg
322 :
323 : END MODULE ai_elec_field
|
