Line data Source code
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 general three-center integrals over Cartesian
10 : !> Gaussian-type functions and a spherical operator centered at position C
11 : !>
12 : !> <a|V(local)|b> = <a|F(|r-C|)|b>
13 : !> \par Literature
14 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
15 : !> \par History
16 : !> - Based in part on code by MK
17 : !> \par Parameters
18 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
19 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
20 : !> - coset : Cartesian orbital set pointer.
21 : !> - dab : Distance between the atomic centers a and b.
22 : !> - dac : Distance between the atomic centers a and c.
23 : !> - dbc : Distance between the atomic centers b and c.
24 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
25 : !> - l{a,b}_max : Maximum angular momentum quantum number of shell a or b.
26 : !> - ncoset : Number of Cartesian orbitals up to l.
27 : !> - rab : Distance vector between the atomic centers a and b.
28 : !> - rac : Distance vector between the atomic centers a and c.
29 : !> - rbc : Distance vector between the atomic centers b and c.
30 : !> - rpgf{a,b,c}: Radius of the primitive Gaussian-type function a or b.
31 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
32 : !> \author jhu (05.2011)
33 : ! **************************************************************************************************
34 : MODULE ai_oneelectron
35 :
36 : USE kinds, ONLY: dp
37 : USE orbital_pointers, ONLY: coset,&
38 : ncoset
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_oneelectron'
46 :
47 : ! *** Public subroutines ***
48 :
49 : PUBLIC :: os_3center, os_2center
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief Calculation of three-center integrals <a|c|b> over
55 : !> Cartesian Gaussian functions and a spherical potential
56 : !>
57 : !> \param la_max_set ...
58 : !> \param la_min_set ...
59 : !> \param npgfa ...
60 : !> \param rpgfa ...
61 : !> \param zeta ...
62 : !> \param lb_max_set ...
63 : !> \param lb_min_set ...
64 : !> \param npgfb ...
65 : !> \param rpgfb ...
66 : !> \param zetb ...
67 : !> \param auxint ...
68 : !> \param rpgfc ...
69 : !> \param rab ...
70 : !> \param dab ...
71 : !> \param rac ...
72 : !> \param dac ...
73 : !> \param rbc ...
74 : !> \param dbc ...
75 : !> \param vab ...
76 : !> \param s ...
77 : !> \param pab ...
78 : !> \param force_a ...
79 : !> \param force_b ...
80 : !> \param fs ...
81 : !> \param vab2 The derivative of the 3-center integrals according to the weighting factors.
82 : !> \param vab2_work ...
83 : !> \param deltaR DIMENSION(3, natoms), weighting factors of the derivatives for each atom and direction
84 : !> \param iatom ...
85 : !> \param jatom ...
86 : !> \param katom ...
87 : !> \date May 2011
88 : !> \author Juerg Hutter
89 : !> \version 1.0
90 : !> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
91 : ! **************************************************************************************************
92 34492497 : SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
93 34492497 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
94 68984994 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
95 68984994 : vab2, vab2_work, deltaR, iatom, jatom, katom)
96 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
97 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
98 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
99 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
100 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: auxint
101 : REAL(KIND=dp), INTENT(IN) :: rpgfc
102 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
103 : REAL(KIND=dp), INTENT(IN) :: dab
104 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
105 : REAL(KIND=dp), INTENT(IN) :: dac
106 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
107 : REAL(KIND=dp), INTENT(IN) :: dbc
108 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
109 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
110 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
111 : OPTIONAL :: pab
112 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
113 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
114 : OPTIONAL :: fs, vab2, vab2_work
115 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
116 : OPTIONAL :: deltaR
117 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
118 :
119 : INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, cdb, cdbx, cdby, cdbz, coa, coamx, &
120 : coamy, coamz, coapx, coapy, coapz, cob, cobmx, cobmy, cobmz, cobpx, cobpy, cobpz, da, &
121 : da_max, dax, day, daz, db, db_max, dbx, dby, dbz, i, ia, iap, iax, iay, iaz, ib, ibm, &
122 : ibx, iby, ibz, idir, ii(3), iim(3), ij, ipgf, ir, ir1, ir2, irm(3), irr(3), irx, iry, &
123 : irz, ix, ixx(1), j, jj(3), jjp(3), jpgf, la, la_max, la_min, lb, lb_max, lb_min, llr, m, &
124 : ma, mb, mmax, na, nb
125 34492497 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iiap
126 : LOGICAL :: calculate_force_a, calculate_force_b
127 : REAL(KIND=dp) :: aai, abx, fax, fay, faz, fbx, fby, fbz, &
128 : ftz, orho, rho, s1, s2
129 : REAL(KIND=dp), DIMENSION(3) :: pai, pbi, pci
130 :
131 34492497 : IF (PRESENT(pab)) THEN
132 8774386 : CPASSERT(PRESENT(fs))
133 8774386 : IF (PRESENT(force_a)) THEN
134 : calculate_force_a = .TRUE.
135 : ELSE
136 0 : calculate_force_a = .FALSE.
137 : END IF
138 8774386 : IF (PRESENT(force_b)) THEN
139 : calculate_force_b = .TRUE.
140 : ELSE
141 0 : calculate_force_b = .FALSE.
142 : END IF
143 : ELSE
144 : calculate_force_a = .FALSE.
145 : calculate_force_b = .FALSE.
146 : END IF
147 :
148 8774386 : IF (calculate_force_a) THEN
149 8774386 : da_max = 1
150 8774386 : force_a = 0.0_dp
151 : ELSE
152 : da_max = 0
153 : END IF
154 :
155 34492497 : IF (calculate_force_b) THEN
156 8774386 : db_max = 1
157 8774386 : force_b = 0.0_dp
158 : ELSE
159 : db_max = 0
160 : END IF
161 :
162 34492497 : IF (PRESENT(vab2)) THEN
163 870 : da_max = 1
164 870 : db_max = 1
165 : END IF
166 :
167 34492497 : la_max = la_max_set + da_max
168 34492497 : la_min = MAX(0, la_min_set - da_max)
169 :
170 34492497 : lb_max = lb_max_set + db_max
171 34492497 : lb_min = MAX(0, lb_min_set - db_max)
172 :
173 34492497 : mmax = la_max + lb_max
174 :
175 : ! precalculate indices for horizontal recursion
176 103477491 : ALLOCATE (iiap(ncoset(mmax), 3))
177 156625738 : DO ma = 0, mmax
178 470015277 : DO iax = 0, ma
179 1114691448 : DO iay = 0, ma - iax
180 679168668 : iaz = ma - iax - iay
181 679168668 : ia = coset(iax, iay, iaz)
182 679168668 : jj(1) = iax; jj(2) = iay; jj(3) = iaz
183 679168668 : jjp = jj
184 679168668 : jjp(1) = jjp(1) + 1
185 679168668 : iap = coset(jjp(1), jjp(2), jjp(3))
186 679168668 : iiap(ia, 1) = iap
187 679168668 : jjp = jj
188 679168668 : jjp(2) = jjp(2) + 1
189 679168668 : iap = coset(jjp(1), jjp(2), jjp(3))
190 679168668 : iiap(ia, 2) = iap
191 679168668 : jjp = jj
192 679168668 : jjp(3) = jjp(3) + 1
193 679168668 : iap = coset(jjp(1), jjp(2), jjp(3))
194 992558207 : iiap(ia, 3) = iap
195 : END DO
196 : END DO
197 : END DO
198 :
199 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
200 :
201 34492497 : na = 0
202 :
203 152671131 : DO ipgf = 1, npgfa
204 :
205 : ! *** Screening ***
206 118178634 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
207 74021846 : na = na + ncoset(la_max_set)
208 74021846 : CYCLE
209 : END IF
210 :
211 44156788 : nb = 0
212 :
213 201455297 : DO jpgf = 1, npgfb
214 :
215 : ! *** Screening ***
216 157298509 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
217 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
218 99682740 : nb = nb + ncoset(lb_max_set)
219 99682740 : CYCLE
220 : END IF
221 :
222 : ! *** Calculate some prefactors ***
223 57615769 : rho = zeta(ipgf) + zetb(jpgf)
224 230463076 : pai(:) = zetb(jpgf)/rho*rab(:)
225 230463076 : pbi(:) = -zeta(ipgf)/rho*rab(:)
226 230463076 : pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
227 57615769 : orho = 0.5_dp/rho
228 :
229 57615769 : ij = (ipgf - 1)*npgfb + jpgf
230 266992174 : s(1, 1, 1:mmax + 1) = auxint(0:mmax, ij)
231 :
232 57615769 : IF (la_max > 0) THEN
233 : ! *** Recurrence steps: [s|c|s] -> [a|c|s] ***
234 : ! *** [a|c|s](m) = (Pi - Ai)*[a-1i|c|s](m) - ***
235 : ! *** (Pi - Ci)*[a-1i|c|s](m+1)) + ***
236 : ! *** Ni(a-1i)/2(a+b)*[a-2i|c|s](m) - ***
237 : ! *** Ni(a-1i)/2(a+b)*[a-2i|c|s](m+1) ***
238 199838400 : DO llr = 1, mmax
239 199838400 : IF (llr == 1) THEN
240 199838400 : DO m = 0, mmax - llr
241 149414599 : s1 = s(1, 1, m + 1)
242 149414599 : s2 = s(1, 1, m + 2)
243 149414599 : s(2, 1, m + 1) = pai(1)*s1 - pci(1)*s2 ! [px|o|s]
244 149414599 : s(3, 1, m + 1) = pai(2)*s1 - pci(2)*s2 ! [py|o|s]
245 199838400 : s(4, 1, m + 1) = pai(3)*s1 - pci(3)*s2 ! [pz|o|s]
246 : END DO
247 98990798 : ELSE IF (llr == 2) THEN
248 147539029 : DO m = 0, mmax - llr
249 98990798 : s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
250 98990798 : s(5, 1, m + 1) = pai(1)*s(2, 1, m + 1) - pci(1)*s(2, 1, m + 2) + orho*s1 ! [dx2|o|s]
251 98990798 : s(6, 1, m + 1) = pai(1)*s(3, 1, m + 1) - pci(1)*s(3, 1, m + 2) ! [dxy|o|s]
252 98990798 : s(7, 1, m + 1) = pai(1)*s(4, 1, m + 1) - pci(1)*s(4, 1, m + 2) ! [dxz|o|s]
253 98990798 : s(8, 1, m + 1) = pai(2)*s(3, 1, m + 1) - pci(2)*s(3, 1, m + 2) + orho*s1 ! [dy2|o|s]
254 98990798 : s(9, 1, m + 1) = pai(2)*s(4, 1, m + 1) - pci(2)*s(4, 1, m + 2) ! [dyz|o|s]
255 147539029 : s(10, 1, m + 1) = pai(3)*s(4, 1, m + 1) - pci(3)*s(4, 1, m + 2) + orho*s1 ! [dz2|o|s]
256 : END DO
257 50442567 : ELSE IF (llr == 3) THEN
258 75863243 : DO m = 0, mmax - llr
259 : s(11, 1, m + 1) = pai(1)*s(5, 1, m + 1) - pci(1)*s(5, 1, m + 2) & ! [fx3 |o|s]
260 50442567 : + 2._dp*orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
261 : s(12, 1, m + 1) = pai(1)*s(6, 1, m + 1) - pci(1)*s(6, 1, m + 2) & ! [fx2y|o|s]
262 50442567 : + orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
263 : s(13, 1, m + 1) = pai(1)*s(7, 1, m + 1) - pci(1)*s(7, 1, m + 2) & ! [fx2z|o|s]
264 50442567 : + orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
265 : s(14, 1, m + 1) = pai(2)*s(6, 1, m + 1) - pci(2)*s(6, 1, m + 2) & ! [fxy2|o|s]
266 50442567 : + orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
267 50442567 : s(15, 1, m + 1) = pai(1)*s(9, 1, m + 1) - pci(1)*s(9, 1, m + 2) ! [fxyz|o|s]
268 : s(16, 1, m + 1) = pai(3)*s(7, 1, m + 1) - pci(3)*s(7, 1, m + 2) & ! [fxz2|o|s]
269 50442567 : + orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
270 : s(17, 1, m + 1) = pai(2)*s(8, 1, m + 1) - pci(2)*s(8, 1, m + 2) & ! [fy3 |o|s]
271 50442567 : + 2._dp*orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
272 : s(18, 1, m + 1) = pai(2)*s(9, 1, m + 1) - pci(2)*s(9, 1, m + 2) & ! [fy2z|o|s]
273 50442567 : + orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
274 : s(19, 1, m + 1) = pai(3)*s(9, 1, m + 1) - pci(3)*s(9, 1, m + 2) & ! [fyz2|o|s]
275 50442567 : + orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
276 : s(20, 1, m + 1) = pai(3)*s(10, 1, m + 1) - pci(3)*s(10, 1, m + 2) & ! [fz3 |o|s]
277 75863243 : + 2._dp*orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
278 : END DO
279 25021891 : ELSE IF (llr == 4) THEN
280 42958353 : DO m = 0, mmax - llr
281 : s(21, 1, m + 1) = pai(1)*s(11, 1, m + 1) - pci(1)*s(11, 1, m + 2) & ! [gx4 |s|s]
282 25021891 : + 3._dp*orho*(s(5, 1, m + 1) - s(5, 1, m + 2))
283 : s(22, 1, m + 1) = pai(1)*s(12, 1, m + 1) - pci(1)*s(12, 1, m + 2) & ! [gx3y |s|s]
284 25021891 : + 2._dp*orho*(s(6, 1, m + 1) - s(6, 1, m + 2))
285 : s(23, 1, m + 1) = pai(1)*s(13, 1, m + 1) - pci(1)*s(13, 1, m + 2) & ! [gx3z |s|s]
286 25021891 : + 2._dp*orho*(s(7, 1, m + 1) - s(7, 1, m + 2))
287 : s(24, 1, m + 1) = pai(1)*s(14, 1, m + 1) - pci(1)*s(14, 1, m + 2) & ! [gx2y2|s|s]
288 25021891 : + orho*(s(8, 1, m + 1) - s(8, 1, m + 2))
289 : s(25, 1, m + 1) = pai(1)*s(15, 1, m + 1) - pci(1)*s(15, 1, m + 2) & ! [gx2yz|s|s]
290 25021891 : + orho*(s(9, 1, m + 1) - s(9, 1, m + 2))
291 : s(26, 1, m + 1) = pai(1)*s(16, 1, m + 1) - pci(1)*s(16, 1, m + 2) & ! [gx2z2|s|s]
292 25021891 : + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
293 25021891 : s(27, 1, m + 1) = pai(1)*s(17, 1, m + 1) - pci(1)*s(17, 1, m + 2) ! [gxy3 |s|s]
294 25021891 : s(28, 1, m + 1) = pai(1)*s(18, 1, m + 1) - pci(1)*s(18, 1, m + 2) ! [gxy2z|s|s]
295 25021891 : s(29, 1, m + 1) = pai(1)*s(19, 1, m + 1) - pci(1)*s(19, 1, m + 2) ! [gxyz2|s|s]
296 25021891 : s(30, 1, m + 1) = pai(1)*s(20, 1, m + 1) - pci(1)*s(20, 1, m + 2) ! [gxz3 |s|s]
297 : s(31, 1, m + 1) = pai(2)*s(17, 1, m + 1) - pci(2)*s(17, 1, m + 2) & ! [gy4 |s|s]
298 25021891 : + 3._dp*orho*(s(8, 1, m + 1) - s(8, 1, m + 2))
299 : s(32, 1, m + 1) = pai(2)*s(18, 1, m + 1) - pci(2)*s(18, 1, m + 2) & ! [gy3z |s|s]
300 25021891 : + 2._dp*orho*(s(9, 1, m + 1) - s(9, 1, m + 2))
301 : s(33, 1, m + 1) = pai(2)*s(19, 1, m + 1) - pci(2)*s(19, 1, m + 2) & ! [gy2z2|s|s]
302 25021891 : + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
303 25021891 : s(34, 1, m + 1) = pai(2)*s(20, 1, m + 1) - pci(2)*s(20, 1, m + 2) ! [gyz3 |s|s]
304 : s(35, 1, m + 1) = pai(3)*s(20, 1, m + 1) - pci(3)*s(20, 1, m + 2) & ! [gz4 |s|s]
305 42958353 : + 3._dp*orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
306 : END DO
307 : ELSE
308 52331726 : DO irx = 0, llr
309 220848278 : DO iry = 0, llr - irx
310 168516552 : irz = llr - irx - iry
311 168516552 : irr(1) = irx; irr(2) = iry; irr(3) = irz
312 842582760 : ixx = MAXLOC(irr)
313 168516552 : ix = ixx(1)
314 168516552 : ir = coset(irx, iry, irz)
315 168516552 : irm = irr
316 168516552 : irm(ix) = irm(ix) - 1
317 168516552 : aai = REAL(MAX(irm(ix), 0), dp)*orho
318 168516552 : ir1 = coset(irm(1), irm(2), irm(3))
319 168516552 : irm(ix) = irm(ix) - 1
320 168516552 : ir2 = coset(irm(1), irm(2), irm(3))
321 443939716 : DO m = 0, mmax - llr
322 : s(ir, 1, m + 1) = pai(ix)*s(ir1, 1, m + 1) - pci(ix)*s(ir1, 1, m + 2) &
323 398693419 : + aai*(s(ir2, 1, m + 1) - s(ir2, 1, m + 2))
324 : END DO
325 : END DO
326 : END DO
327 : END IF
328 : END DO
329 :
330 : ! *** Horizontal recurrence steps ***
331 : ! *** [a|c|b+1i] = [a+1i|c|b] + (Ai - Bi)*[a|c|b] ***
332 :
333 124000630 : DO mb = 1, lb_max
334 300146053 : DO ibx = 0, mb
335 561324615 : DO iby = 0, mb - ibx
336 311602363 : ibz = mb - ibx - iby
337 311602363 : ib = coset(ibx, iby, ibz)
338 311602363 : ii(1) = ibx; ii(2) = iby; ii(3) = ibz
339 1558011815 : ixx = MAXLOC(ii)
340 311602363 : ix = ixx(1)
341 311602363 : abx = -rab(ix)
342 311602363 : iim = ii
343 311602363 : iim(ix) = iim(ix) - 1
344 311602363 : ibm = coset(iim(1), iim(2), iim(3))
345 4778594840 : DO ia = 1, ncoset(mmax - mb)
346 4290847054 : iap = iiap(ia, ix)
347 4602449417 : s(ia, ib, 1) = s(iap, ibm, 1) + abx*s(ia, ibm, 1)
348 : END DO
349 : END DO
350 : END DO
351 : END DO
352 :
353 7191968 : ELSE IF (lb_max > 0) THEN
354 :
355 : ! *** Recurrence steps: [s|c|s] -> [s|c|b] ***
356 : ! *** [s|c|b](m) = (Pi - Bi)*[s|c|b-1i](m) - ***
357 : ! *** (Pi - Ci)*[s|c|b-1i](m+1)) + ***
358 : ! *** Ni(b-1i)/2(a+b)*[s|c|b-2i](m) - ***
359 : ! *** Ni(b-1i)/2(a+b)*[s|c|b-2i](m+1) ***
360 4464835 : DO llr = 1, lb_max
361 4464835 : IF (llr == 1) THEN
362 4464835 : DO m = 0, lb_max - llr
363 2346037 : s1 = s(1, 1, m + 1)
364 2346037 : s2 = s(1, 1, m + 2)
365 2346037 : s(1, 2, m + 1) = pbi(1)*s1 - pci(1)*s2 ! [px|o|s]
366 2346037 : s(1, 3, m + 1) = pbi(2)*s1 - pci(2)*s2 ! [py|o|s]
367 4464835 : s(1, 4, m + 1) = pbi(3)*s1 - pci(3)*s2 ! [pz|o|s]
368 : END DO
369 227239 : ELSE IF (llr == 2) THEN
370 442850 : DO m = 0, lb_max - llr
371 227239 : s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
372 227239 : s(1, 5, m + 1) = pbi(1)*s(1, 2, m + 1) - pci(1)*s(1, 2, m + 2) + orho*s1 ! [dx2|o|s]
373 227239 : s(1, 6, m + 1) = pbi(1)*s(1, 3, m + 1) - pci(1)*s(1, 3, m + 2) ! [dxy|o|s]
374 227239 : s(1, 7, m + 1) = pbi(1)*s(1, 4, m + 1) - pci(1)*s(1, 4, m + 2) ! [dxz|o|s]
375 227239 : s(1, 8, m + 1) = pbi(2)*s(1, 3, m + 1) - pci(2)*s(1, 3, m + 2) + orho*s1 ! [dy2|o|s]
376 227239 : s(1, 9, m + 1) = pbi(2)*s(1, 4, m + 1) - pci(2)*s(1, 4, m + 2) ! [dyz|o|s]
377 442850 : s(1, 10, m + 1) = pbi(3)*s(1, 4, m + 1) - pci(3)*s(1, 4, m + 2) + orho*s1 ! [dz2|o|s]
378 : END DO
379 11628 : ELSE IF (llr == 3) THEN
380 23174 : DO m = 0, lb_max - llr
381 : s(1, 11, m + 1) = pbi(1)*s(1, 5, m + 1) - pci(1)*s(1, 5, m + 2) & ! [fx3 |o|s]
382 11628 : + 2._dp*orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
383 : s(1, 12, m + 1) = pbi(1)*s(1, 6, m + 1) - pci(1)*s(1, 6, m + 2) & ! [fx2y|o|s]
384 11628 : + orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
385 : s(1, 13, m + 1) = pbi(1)*s(1, 7, m + 1) - pci(1)*s(1, 7, m + 2) & ! [fx2z|o|s]
386 11628 : + orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
387 : s(1, 14, m + 1) = pbi(2)*s(1, 6, m + 1) - pci(2)*s(1, 6, m + 2) & ! [fxy2|o|s]
388 11628 : + orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
389 11628 : s(1, 15, m + 1) = pbi(1)*s(1, 9, m + 1) - pci(1)*s(1, 9, m + 2) ! [fxyz|o|s]
390 : s(1, 16, m + 1) = pbi(3)*s(1, 7, m + 1) - pci(3)*s(1, 7, m + 2) & ! [fxz2|o|s]
391 11628 : + orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
392 : s(1, 17, m + 1) = pbi(2)*s(1, 8, m + 1) - pci(2)*s(1, 8, m + 2) & ! [fy3 |o|s]
393 11628 : + 2._dp*orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
394 : s(1, 18, m + 1) = pbi(2)*s(1, 9, m + 1) - pci(2)*s(1, 9, m + 2) & ! [fy2z|o|s]
395 11628 : + orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
396 : s(1, 19, m + 1) = pbi(3)*s(1, 9, m + 1) - pci(3)*s(1, 9, m + 2) & ! [fyz2|o|s]
397 11628 : + orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
398 : s(1, 20, m + 1) = pbi(3)*s(1, 10, m + 1) - pci(3)*s(1, 10, m + 2) & ! [fz3 |o|s]
399 23174 : + 2._dp*orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
400 : END DO
401 : ELSE
402 492 : DO irx = 0, llr
403 1722 : DO iry = 0, llr - irx
404 1230 : irz = llr - irx - iry
405 1230 : irr(1) = irx; irr(2) = iry; irr(3) = irz
406 6150 : ixx = MAXLOC(irr)
407 1230 : ix = ixx(1)
408 1230 : ir = coset(irx, iry, irz)
409 1230 : irm = irr
410 1230 : irm(ix) = irm(ix) - 1
411 1230 : aai = REAL(MAX(irm(ix), 0), dp)
412 1230 : ir1 = coset(irm(1), irm(2), irm(3))
413 1230 : irm(ix) = irm(ix) - 1
414 1230 : ir2 = coset(irm(1), irm(2), irm(3))
415 2870 : DO m = 0, lb_max - llr
416 : s(1, ir, m + 1) = pbi(ix)*s(1, ir1, m + 1) - pci(ix)*s(1, ir1, m + 2) &
417 2460 : + aai*orho*(s(1, ir2, m + 1) - s(1, ir2, m + 2))
418 : END DO
419 : END DO
420 : END DO
421 : END IF
422 : END DO
423 :
424 : END IF
425 :
426 : ! *** Store the primitive three-center overlap integrals ***
427 306821050 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
428 1693890490 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
429 1636274721 : vab(na + i, nb + j) = vab(na + i, nb + j) + s(i, j, 1)
430 : END DO
431 : END DO
432 :
433 : ! *** Calculate the requested derivatives with respect ***
434 : ! *** to the nuclear coordinates of the atomic center a ***
435 :
436 73578595 : DO da = 0, da_max - 1
437 15962826 : ftz = 2.0_dp*zeta(ipgf)
438 89541421 : DO dax = 0, da
439 47888478 : DO day = 0, da - dax
440 15962826 : daz = da - dax - day
441 15962826 : cda = coset(dax, day, daz)
442 15962826 : cdax = coset(dax + 1, day, daz)
443 15962826 : cday = coset(dax, day + 1, daz)
444 15962826 : cdaz = coset(dax, day, daz + 1)
445 :
446 : ! *** [da/dAi|c|b] = 2*zeta*[a+1i|c|b] - Ni(a)[a-1i|c|b] ***
447 :
448 62326157 : DO la = la_min_set, la_max - da - 1
449 96468147 : DO ax = 0, la
450 50104816 : fax = REAL(ax, dp)
451 154004881 : DO ay = 0, la - ax
452 73499560 : fay = REAL(ay, dp)
453 73499560 : az = la - ax - ay
454 73499560 : faz = REAL(az, dp)
455 73499560 : coa = coset(ax, ay, az)
456 73499560 : coamx = coset(ax - 1, ay, az)
457 73499560 : coamy = coset(ax, ay - 1, az)
458 73499560 : coamz = coset(ax, ay, az - 1)
459 73499560 : coapx = coset(ax + 1, ay, az)
460 73499560 : coapy = coset(ax, ay + 1, az)
461 73499560 : coapz = coset(ax, ay, az + 1)
462 557843474 : DO cob = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
463 434239098 : fs(coa, cob, cdax) = ftz*s(coapx, cob, cda) - fax*s(coamx, cob, cda)
464 434239098 : fs(coa, cob, cday) = ftz*s(coapy, cob, cda) - fay*s(coamy, cob, cda)
465 507738658 : fs(coa, cob, cdaz) = ftz*s(coapz, cob, cda) - faz*s(coamz, cob, cda)
466 : END DO
467 : END DO
468 : END DO
469 : END DO
470 : END DO
471 :
472 : END DO
473 : END DO
474 :
475 : ! DFPT for APTs
476 57615769 : IF (PRESENT(vab2_work)) THEN
477 28512 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
478 72090 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
479 43578 : vab2_work(na + i, nb + j, 1) = vab2_work(na + i, nb + j, 1) + fs(i, j, 2)
480 43578 : vab2_work(na + i, nb + j, 2) = vab2_work(na + i, nb + j, 2) + fs(i, j, 3)
481 62676 : vab2_work(na + i, nb + j, 3) = vab2_work(na + i, nb + j, 3) + fs(i, j, 4)
482 : END DO
483 : END DO
484 : END IF
485 :
486 : ! *** Calculate the force contribution for the atomic center a ***
487 :
488 57615769 : IF (calculate_force_a) THEN
489 89585554 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
490 523781074 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
491 434195520 : force_a(1) = force_a(1) + pab(na + i, nb + j)*fs(i, j, 2)
492 434195520 : force_a(2) = force_a(2) + pab(na + i, nb + j)*fs(i, j, 3)
493 507827662 : force_a(3) = force_a(3) + pab(na + i, nb + j)*fs(i, j, 4)
494 : END DO
495 : END DO
496 : END IF
497 :
498 : ! *** Calculate the requested derivatives with respect ***
499 : ! *** to the nuclear coordinates of the atomic center b ***
500 :
501 73578595 : DO db = 0, db_max - 1
502 15962826 : ftz = 2.0_dp*zetb(jpgf)
503 89541421 : DO dbx = 0, db
504 47888478 : DO dby = 0, db - dbx
505 15962826 : dbz = db - dbx - dby
506 15962826 : cdb = coset(dbx, dby, dbz)
507 15962826 : cdbx = coset(dbx + 1, dby, dbz)
508 15962826 : cdby = coset(dbx, dby + 1, dbz)
509 15962826 : cdbz = coset(dbx, dby, dbz + 1)
510 :
511 : ! *** [a|c|db/dBi] = 2*zetb*[a|c|b+1i] - Ni(b)[a|c|b-1i] ***
512 :
513 62374369 : DO lb = lb_min_set, lb_max - db - 1
514 96615181 : DO bx = 0, lb
515 50203638 : fbx = REAL(bx, dp)
516 154303595 : DO by = 0, lb - bx
517 73651240 : fby = REAL(by, dp)
518 73651240 : bz = lb - bx - by
519 73651240 : fbz = REAL(bz, dp)
520 73651240 : cob = coset(bx, by, bz)
521 73651240 : cobmx = coset(bx - 1, by, bz)
522 73651240 : cobmy = coset(bx, by - 1, bz)
523 73651240 : cobmz = coset(bx, by, bz - 1)
524 73651240 : cobpx = coset(bx + 1, by, bz)
525 73651240 : cobpy = coset(bx, by + 1, bz)
526 73651240 : cobpz = coset(bx, by, bz + 1)
527 558093976 : DO coa = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
528 434239098 : fs(coa, cob, cdbx) = ftz*s(coa, cobpx, cdb) - fbx*s(coa, cobmx, cdb)
529 434239098 : fs(coa, cob, cdby) = ftz*s(coa, cobpy, cdb) - fby*s(coa, cobmy, cdb)
530 507890338 : fs(coa, cob, cdbz) = ftz*s(coa, cobpz, cdb) - fbz*s(coa, cobmz, cdb)
531 : END DO
532 : END DO
533 : END DO
534 : END DO
535 :
536 : END DO
537 : END DO
538 : END DO
539 :
540 : ! DFPT for APTs
541 57615769 : IF (PRESENT(vab2_work)) THEN
542 28512 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
543 72090 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
544 43578 : vab2_work(na + i, nb + j, 4) = vab2_work(na + i, nb + j, 4) + fs(i, j, 2)
545 43578 : vab2_work(na + i, nb + j, 5) = vab2_work(na + i, nb + j, 5) + fs(i, j, 3)
546 62676 : vab2_work(na + i, nb + j, 6) = vab2_work(na + i, nb + j, 6) + fs(i, j, 4)
547 : END DO
548 : END DO
549 :
550 37656 : DO idir = 1, 3
551 94950 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
552 216270 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
553 : vab2(na + i, nb + j, idir) = vab2(na + i, nb + j, idir) &
554 : + vab2_work(na + i, nb + j, idir)*deltaR(idir, iatom) &
555 : - vab2_work(na + i, nb + j, idir)*deltaR(idir, katom) &
556 : + vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, jatom) &
557 188028 : - vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, katom)
558 : END DO
559 : END DO
560 : END DO
561 : END IF
562 :
563 : ! *** Calculate the force contribution for the atomic center b ***
564 :
565 57615769 : IF (calculate_force_b) THEN
566 89585554 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
567 523781074 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
568 434195520 : force_b(1) = force_b(1) + pab(na + i, nb + j)*fs(i, j, 2)
569 434195520 : force_b(2) = force_b(2) + pab(na + i, nb + j)*fs(i, j, 3)
570 507827662 : force_b(3) = force_b(3) + pab(na + i, nb + j)*fs(i, j, 4)
571 : END DO
572 : END DO
573 : END IF
574 :
575 101772557 : nb = nb + ncoset(lb_max_set)
576 :
577 : END DO
578 :
579 78649285 : na = na + ncoset(la_max_set)
580 :
581 : END DO
582 :
583 34492497 : DEALLOCATE (iiap)
584 :
585 34492497 : END SUBROUTINE os_3center
586 : ! **************************************************************************************************
587 : !> \brief Calculation of two-center integrals <a|c> over
588 : !> Cartesian Gaussian functions and a spherical potential
589 : !>
590 : !> \param la_max_set ...
591 : !> \param la_min_set ...
592 : !> \param npgfa ...
593 : !> \param rpgfa ...
594 : !> \param zeta ...
595 : !> \param auxint ...
596 : !> \param rpgfc ...
597 : !> \param rac ...
598 : !> \param dac ...
599 : !> \param va ...
600 : !> \param dva ...
601 : !> \date December 2017
602 : !> \author Juerg Hutter
603 : !> \version 1.0
604 : ! **************************************************************************************************
605 328 : SUBROUTINE os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
606 656 : auxint, rpgfc, rac, dac, va, dva)
607 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
608 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
609 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: auxint
610 : REAL(KIND=dp), INTENT(IN) :: rpgfc
611 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
612 : REAL(KIND=dp), INTENT(IN) :: dac
613 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: va
614 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
615 : OPTIONAL :: dva
616 :
617 : INTEGER :: ax, ay, az, coa, coamx, coamy, coamz, coapx, coapy, coapz, da_max, i, ipgf, ir, &
618 : ir1, ir2, irm(3), irr(3), irx, iry, irz, ix, ixx(1), la, la_max, la_min, llr, m, mmax, na
619 : REAL(KIND=dp) :: aai, fax, fay, faz, ftz, orho, s1
620 328 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: s
621 :
622 328 : IF (PRESENT(dva)) THEN
623 : da_max = 1
624 : ELSE
625 164 : da_max = 0
626 : END IF
627 :
628 328 : la_max = la_max_set + da_max
629 328 : la_min = MAX(0, la_min_set - da_max)
630 :
631 328 : mmax = la_max
632 :
633 1312 : ALLOCATE (s(ncoset(mmax), mmax + 1))
634 328 : na = 0
635 656 : DO ipgf = 1, npgfa
636 328 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
637 0 : na = na + ncoset(la_max_set)
638 0 : CYCLE
639 : END IF
640 1260 : s(1, 1:mmax + 1) = auxint(0:mmax, ipgf)
641 328 : IF (la_max > 0) THEN
642 : ! Recurrence steps: [s|c] -> [a|c]
643 : ! [a|c](m) = (Ci - Ai)*[a-1i|c](m+1) +
644 : ! Ni(a-1i)/2a*[a-2i|c](m) -
645 : ! Ni(a-1i)/2a*[a-2i|c](m+1)
646 : !
647 :
648 268 : orho = 0.5_dp/zeta(ipgf)
649 :
650 872 : DO llr = 1, mmax
651 872 : IF (llr == 1) THEN
652 872 : DO m = 0, mmax - llr
653 604 : s1 = s(1, m + 2)
654 604 : s(2, m + 1) = -rac(1)*s1 ! [px|o]
655 604 : s(3, m + 1) = -rac(2)*s1 ! [py|o]
656 872 : s(4, m + 1) = -rac(3)*s1 ! [pz|o]
657 : END DO
658 336 : ELSE IF (llr == 2) THEN
659 544 : DO m = 0, mmax - llr
660 336 : s1 = s(1, m + 1) - s(1, m + 2)
661 336 : s(5, m + 1) = -rac(1)*s(2, m + 2) + orho*s1 ! [dx2|o]
662 336 : s(6, m + 1) = -rac(1)*s(3, m + 2) ! [dxy|o]
663 336 : s(7, m + 1) = -rac(1)*s(4, m + 2) ! [dxz|o]
664 336 : s(8, m + 1) = -rac(2)*s(3, m + 2) + orho*s1 ! [dy2|o]
665 336 : s(9, m + 1) = -rac(2)*s(4, m + 2) ! [dyz|o]
666 544 : s(10, m + 1) = -rac(3)*s(4, m + 2) + orho*s1 ! [dz2|o]
667 : END DO
668 128 : ELSE IF (llr == 3) THEN
669 244 : DO m = 0, mmax - llr
670 128 : s(11, m + 1) = -rac(1)*s(5, m + 2) + 2._dp*orho*(s(2, m + 1) - s(2, m + 2)) ! [fx3 |o]
671 128 : s(12, m + 1) = -rac(1)*s(6, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fx2y|o]
672 128 : s(13, m + 1) = -rac(1)*s(7, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fx2z|o]
673 128 : s(14, m + 1) = -rac(2)*s(6, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxy2|o]
674 128 : s(15, m + 1) = -rac(1)*s(9, m + 2) ! [fxyz|o]
675 128 : s(16, m + 1) = -rac(3)*s(7, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxz2|o]
676 128 : s(17, m + 1) = -rac(2)*s(8, m + 2) + 2._dp*orho*(s(3, m + 1) - s(3, m + 2)) ! [fy3 |o]
677 128 : s(18, m + 1) = -rac(2)*s(9, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fy2z|o]
678 128 : s(19, m + 1) = -rac(3)*s(9, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fyz2|o]
679 244 : s(20, m + 1) = -rac(3)*s(10, m + 2) + 2._dp*orho*(s(4, m + 1) - s(4, m + 2)) ! [fz3 |o]
680 : END DO
681 12 : ELSE IF (llr == 4) THEN
682 24 : DO m = 0, mmax - llr
683 12 : s(21, m + 1) = -rac(1)*s(11, m + 2) + 3._dp*orho*(s(5, m + 1) - s(5, m + 2)) ! [gx4 |s]
684 12 : s(22, m + 1) = -rac(1)*s(12, m + 2) + 2._dp*orho*(s(6, m + 1) - s(6, m + 2)) ! [gx3y |s]
685 12 : s(23, m + 1) = -rac(1)*s(13, m + 2) + 2._dp*orho*(s(7, m + 1) - s(7, m + 2)) ! [gx3z |s]
686 12 : s(24, m + 1) = -rac(1)*s(14, m + 2) + orho*(s(8, m + 1) - s(8, m + 2)) ! [gx2y2|s]
687 12 : s(25, m + 1) = -rac(1)*s(15, m + 2) + orho*(s(9, m + 1) - s(9, m + 2)) ! [gx2yz|s]
688 12 : s(26, m + 1) = -rac(1)*s(16, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gx2z2|s]
689 12 : s(27, m + 1) = -rac(1)*s(17, m + 2) ! [gxy3 |s]
690 12 : s(28, m + 1) = -rac(1)*s(18, m + 2) ! [gxy2z|s]
691 12 : s(29, m + 1) = -rac(1)*s(19, m + 2) ! [gxyz2|s]
692 12 : s(30, m + 1) = -rac(1)*s(20, m + 2) ! [gxz3 |s]
693 12 : s(31, m + 1) = -rac(2)*s(17, m + 2) + 3._dp*orho*(s(8, m + 1) - s(8, m + 2)) ! [gy4 |s]
694 12 : s(32, m + 1) = -rac(2)*s(18, m + 2) + 2._dp*orho*(s(9, m + 1) - s(9, m + 2)) ! [gy3z |s]
695 12 : s(33, m + 1) = -rac(2)*s(19, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gy2z2|s]
696 12 : s(34, m + 1) = -rac(2)*s(20, m + 2) ! [gyz3 |s]
697 24 : s(35, m + 1) = -rac(3)*s(20, m + 2) + 3._dp*orho*(s(10, m + 1) - s(10, m + 2)) ! [gz4 |s]
698 : END DO
699 : ELSE
700 0 : DO irx = 0, llr
701 0 : DO iry = 0, llr - irx
702 0 : irz = llr - irx - iry
703 0 : irr(1) = irx; irr(2) = iry; irr(3) = irz
704 0 : ixx = MAXLOC(irr)
705 0 : ix = ixx(1)
706 0 : ir = coset(irx, iry, irz)
707 0 : irm = irr
708 0 : irm(ix) = irm(ix) - 1
709 0 : aai = REAL(MAX(irm(ix), 0), dp)*orho
710 0 : ir1 = coset(irm(1), irm(2), irm(3))
711 0 : irm(ix) = irm(ix) - 1
712 0 : ir2 = coset(irm(1), irm(2), irm(3))
713 0 : DO m = 0, mmax - llr
714 0 : s(ir, m + 1) = -rac(ix)*s(ir1, m + 2) + aai*(s(ir2, m + 1) - s(ir2, m + 2))
715 : END DO
716 : END DO
717 : END DO
718 : END IF
719 : END DO
720 :
721 : END IF
722 :
723 : ! Store the primitive three-center overlap integrals
724 2404 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
725 2404 : va(na + i) = va(na + i) + s(i, 1)
726 : END DO
727 :
728 : ! Calculate the requested derivatives with respect ***
729 : ! to the nuclear coordinates of the atomic center a ***
730 : ! [da/dAi|c] = 2*zeta*[a+1i|c] - Ni(a)[a-1i|c] ***
731 328 : IF (PRESENT(dva)) THEN
732 164 : ftz = 2.0_dp*zeta(ipgf)
733 450 : DO la = la_min_set, la_max_set
734 1048 : DO ax = 0, la
735 598 : fax = REAL(ax, dp)
736 1922 : DO ay = 0, la - ax
737 1038 : fay = REAL(ay, dp)
738 1038 : az = la - ax - ay
739 1038 : faz = REAL(az, dp)
740 1038 : coa = coset(ax, ay, az)
741 1038 : coamx = coset(ax - 1, ay, az)
742 1038 : coamy = coset(ax, ay - 1, az)
743 1038 : coamz = coset(ax, ay, az - 1)
744 1038 : coapx = coset(ax + 1, ay, az)
745 1038 : coapy = coset(ax, ay + 1, az)
746 1038 : coapz = coset(ax, ay, az + 1)
747 1038 : dva(na + coa, 1) = dva(na + coa, 1) + ftz*s(coapx, 1) - fax*s(coamx, 1)
748 1038 : dva(na + coa, 2) = dva(na + coa, 2) + ftz*s(coapy, 1) - fay*s(coamy, 1)
749 1636 : dva(na + coa, 3) = dva(na + coa, 3) + ftz*s(coapz, 1) - faz*s(coamz, 1)
750 : END DO
751 : END DO
752 : END DO
753 : END IF
754 :
755 656 : na = na + ncoset(la_max_set)
756 :
757 : END DO
758 :
759 328 : DEALLOCATE (s)
760 :
761 328 : END SUBROUTINE os_2center
762 : ! **************************************************************************************************
763 :
764 : END MODULE ai_oneelectron
|
