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 integrals over Cartesian Gaussian-type functions for [a|(r-Ra)^(2m)|b]
10 : !> Ra is the position of center a
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 : !> - cx,cy,cz : Angular momentum index numbers of orbital c.
18 : !> - coset : Cartesian orbital set pointer.
19 : !> - dac : Distance between the atomic centers a and c.
20 : !> - l{a,c} : Angular momentum quantum number of shell a or c.
21 : !> - l{a,c}_max : Maximum angular momentum quantum number of shell a or c.
22 : !> - l{a,c}_min : Minimum angular momentum quantum number of shell a or c.
23 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
24 : !> - npgf{a,c} : Degree of contraction of shell a or c.
25 : !> - rac : Distance vector between the atomic centers a and c.
26 : !> - rac2 : Square of the distance between the atomic centers a and c.
27 : !> - zet{a,c} : Exponents of the Gaussian-type functions a or c.
28 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
29 : !> - zetw : Reciprocal of the sum of the exponents of orbital a and c.
30 : !> \author Dorothea Golze (08.2016)
31 : ! **************************************************************************************************
32 : MODULE ai_operator_ra2m
33 :
34 : USE ai_os_rr, ONLY: os_rr_ovlp
35 : USE kinds, ONLY: dp
36 : USE mathconstants, ONLY: fac,&
37 : pi
38 : USE orbital_pointers, ONLY: coset,&
39 : ncoset
40 : #include "../base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_operator_ra2m'
44 :
45 : PRIVATE
46 :
47 : ! *** Public subroutines ***
48 :
49 : PUBLIC :: operator_ra2m
50 :
51 : ! **************************************************************************************************
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Calculation of the primitive two-center [a|(r-Ra)^(2m)|b] integrals over Cartesian
57 : !> Gaussian-type functions; operator is here (r-Ra)^(2m)
58 : !> \param la_max ...
59 : !> \param la_min ...
60 : !> \param npgfa ...
61 : !> \param zeta ...
62 : !> \param lb_max ...
63 : !> \param lb_min ...
64 : !> \param npgfb ...
65 : !> \param zetb ...
66 : !> \param m exponent in (r-Ra)^(2m) operator
67 : !> \param rab ...
68 : !> \param sab ...
69 : !> \param dsab ...
70 : !> \param calculate_forces ...
71 : ! **************************************************************************************************
72 14400 : SUBROUTINE operator_ra2m(la_max, la_min, npgfa, zeta, &
73 4800 : lb_max, lb_min, npgfb, zetb, &
74 4800 : m, rab, sab, dsab, calculate_forces)
75 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
76 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
77 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
78 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
79 : INTEGER, INTENT(IN) :: m
80 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
81 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
82 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dsab
83 : LOGICAL, INTENT(IN) :: calculate_forces
84 :
85 : CHARACTER(len=*), PARAMETER :: routineN = 'operator_ra2m'
86 :
87 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, &
88 : handle, i, ia, ib, ipgf, j, jpgf, k, &
89 : la, lb, ldrr, lma, lmb, ma, mb, na, &
90 : nb, ofa, ofb
91 : REAL(KIND=dp) :: a, b, dumx, dumy, dumz, f0, prefac, &
92 : rab2, tab, xhi, zet
93 4800 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr
94 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
95 :
96 4800 : CALL timeset(routineN, handle)
97 :
98 6052800 : sab = 0.0_dp
99 18163200 : IF (calculate_forces) dsab = 0.0_dp
100 :
101 : ! Distance of the centers a and b
102 :
103 4800 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
104 : tab = SQRT(rab2)
105 :
106 : ! Maximum l for auxiliary integrals
107 4800 : lma = la_max + 2*m
108 4800 : lmb = lb_max
109 4800 : IF (calculate_forces) lma = lma + 1
110 4800 : ldrr = MAX(lma, lmb) + 1
111 :
112 : ! Allocate space for auxiliary integrals
113 24000 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
114 :
115 : ! Number of integrals, check size of arrays
116 4800 : ofa = ncoset(la_min - 1)
117 4800 : ofb = ncoset(lb_min - 1)
118 4800 : na = ncoset(la_max) - ofa
119 4800 : nb = ncoset(lb_max) - ofb
120 4800 : CPASSERT((SIZE(sab, 1) >= na*npgfa))
121 4800 : CPASSERT((SIZE(sab, 2) >= nb*npgfb))
122 :
123 : ! Loops over all pairs of primitive Gaussian-type functions
124 4800 : ma = 0
125 9600 : DO ipgf = 1, npgfa
126 : mb = 0
127 9600 : DO jpgf = 1, npgfb
128 :
129 : ! Calculate some prefactors
130 4800 : a = zeta(ipgf)
131 4800 : b = zetb(jpgf)
132 4800 : zet = a + b
133 4800 : xhi = a*b/zet
134 19200 : rap = b*rab/zet
135 19200 : rbp = -a*rab/zet
136 :
137 : ! [s|s] integral
138 4800 : f0 = (pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
139 :
140 : ! Calculate the recurrence relation
141 4800 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
142 :
143 23040 : DO lb = lb_min, lb_max
144 72320 : DO bx = 0, lb
145 176320 : DO by = 0, lb - bx
146 108800 : bz = lb - bx - by
147 108800 : cob = coset(bx, by, bz) - ofb
148 108800 : ib = mb + cob
149 517120 : DO la = la_min, la_max
150 1294720 : DO ax = 0, la
151 2763520 : DO ay = 0, la - ax
152 1577600 : az = la - ax - ay
153 1577600 : coa = coset(ax, ay, az) - ofa
154 1577600 : ia = ma + coa
155 8714880 : DO i = 0, m
156 33129600 : DO j = 0, m
157 132518400 : DO k = 0, m
158 100966400 : IF (i + j + k /= m) CYCLE
159 15776000 : prefac = fac(m)/fac(i)/fac(j)/fac(k)
160 : sab(ia, ib) = sab(ia, ib) + prefac*f0 &
161 15776000 : *rr(ax + 2*i, bx, 1)*rr(ay + 2*j, by, 2)*rr(az + 2*k, bz, 3)
162 41017600 : IF (calculate_forces) THEN
163 : ! (da|b) = 2*a*(a+1|b) - N(a)*(a-1|b)
164 : ! dx
165 15776000 : dumx = 2.0_dp*a*rr(ax + 2*i + 1, bx, 1)
166 15776000 : IF (ax + 2*i > 0) dumx = dumx - REAL(ax + 2*i, dp)*rr(ax + 2*i - 1, bx, 1)
167 15776000 : dsab(ia, ib, 1) = dsab(ia, ib, 1) + prefac*f0*dumx*rr(ay + 2*j, by, 2)*rr(az + 2*k, bz, 3)
168 : ! dy
169 15776000 : dumy = 2.0_dp*a*rr(ay + 2*j + 1, by, 2)
170 15776000 : IF (ay + 2*j > 0) dumy = dumy - REAL(ay + 2*j, dp)*rr(ay + 2*j - 1, by, 2)
171 15776000 : dsab(ia, ib, 2) = dsab(ia, ib, 2) + prefac*f0*rr(ax + 2*i, bx, 1)*dumy*rr(az + 2*k, bz, 3)
172 : ! dz
173 15776000 : dumz = 2.0_dp*a*rr(az + 2*k + 1, bz, 3)
174 15776000 : IF (az + 2*k > 0) dumz = dumz - REAL(az + 2*k, dp)*rr(az + 2*k - 1, bz, 3)
175 15776000 : dsab(ia, ib, 3) = dsab(ia, ib, 3) + prefac*f0*rr(ax + 2*i, bx, 1)*rr(ay + 2*j, by, 2)*dumz
176 : END IF
177 : END DO
178 : END DO
179 : END DO
180 : !
181 : END DO
182 : END DO
183 : END DO !la
184 : END DO
185 : END DO
186 : END DO !lb
187 :
188 9600 : mb = mb + nb
189 : END DO
190 9600 : ma = ma + na
191 : END DO
192 :
193 4800 : DEALLOCATE (rr)
194 :
195 4800 : CALL timestop(handle)
196 :
197 4800 : END SUBROUTINE operator_ra2m
198 :
199 : END MODULE ai_operator_ra2m
|
