(git:1f285aa)
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 ! **************************************************************************************************
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  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 
323 END 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