(git:374b731)
Loading...
Searching...
No Matches
ai_oneelectron.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 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! **************************************************************************************************
34MODULE 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
51CONTAINS
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 SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
93 lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
94 rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
95 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 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 IF (PRESENT(pab)) THEN
132 cpassert(PRESENT(fs))
133 IF (PRESENT(force_a)) THEN
134 calculate_force_a = .true.
135 ELSE
136 calculate_force_a = .false.
137 END IF
138 IF (PRESENT(force_b)) THEN
139 calculate_force_b = .true.
140 ELSE
141 calculate_force_b = .false.
142 END IF
143 ELSE
144 calculate_force_a = .false.
145 calculate_force_b = .false.
146 END IF
147
148 IF (calculate_force_a) THEN
149 da_max = 1
150 force_a = 0.0_dp
151 ELSE
152 da_max = 0
153 END IF
154
155 IF (calculate_force_b) THEN
156 db_max = 1
157 force_b = 0.0_dp
158 ELSE
159 db_max = 0
160 END IF
161
162 IF (PRESENT(vab2)) THEN
163 da_max = 1
164 db_max = 1
165 END IF
166
167 la_max = la_max_set + da_max
168 la_min = max(0, la_min_set - da_max)
169
170 lb_max = lb_max_set + db_max
171 lb_min = max(0, lb_min_set - db_max)
172
173 mmax = la_max + lb_max
174
175 ! precalculate indices for horizontal recursion
176 ALLOCATE (iiap(ncoset(mmax), 3))
177 DO ma = 0, mmax
178 DO iax = 0, ma
179 DO iay = 0, ma - iax
180 iaz = ma - iax - iay
181 ia = coset(iax, iay, iaz)
182 jj(1) = iax; jj(2) = iay; jj(3) = iaz
183 jjp = jj
184 jjp(1) = jjp(1) + 1
185 iap = coset(jjp(1), jjp(2), jjp(3))
186 iiap(ia, 1) = iap
187 jjp = jj
188 jjp(2) = jjp(2) + 1
189 iap = coset(jjp(1), jjp(2), jjp(3))
190 iiap(ia, 2) = iap
191 jjp = jj
192 jjp(3) = jjp(3) + 1
193 iap = coset(jjp(1), jjp(2), jjp(3))
194 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 na = 0
202
203 DO ipgf = 1, npgfa
204
205! *** Screening ***
206 IF (rpgfa(ipgf) + rpgfc < dac) THEN
207 na = na + ncoset(la_max_set)
208 cycle
209 END IF
210
211 nb = 0
212
213 DO jpgf = 1, npgfb
214
215! *** Screening ***
216 IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
217 (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
218 nb = nb + ncoset(lb_max_set)
219 cycle
220 END IF
221
222! *** Calculate some prefactors ***
223 rho = zeta(ipgf) + zetb(jpgf)
224 pai(:) = zetb(jpgf)/rho*rab(:)
225 pbi(:) = -zeta(ipgf)/rho*rab(:)
226 pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
227 orho = 0.5_dp/rho
228
229 ij = (ipgf - 1)*npgfb + jpgf
230 s(1, 1, 1:mmax + 1) = auxint(0:mmax, ij)
231
232 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 DO llr = 1, mmax
239 IF (llr == 1) THEN
240 DO m = 0, mmax - llr
241 s1 = s(1, 1, m + 1)
242 s2 = s(1, 1, m + 2)
243 s(2, 1, m + 1) = pai(1)*s1 - pci(1)*s2 ! [px|o|s]
244 s(3, 1, m + 1) = pai(2)*s1 - pci(2)*s2 ! [py|o|s]
245 s(4, 1, m + 1) = pai(3)*s1 - pci(3)*s2 ! [pz|o|s]
246 END DO
247 ELSE IF (llr == 2) THEN
248 DO m = 0, mmax - llr
249 s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
250 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 s(6, 1, m + 1) = pai(1)*s(3, 1, m + 1) - pci(1)*s(3, 1, m + 2) ! [dxy|o|s]
252 s(7, 1, m + 1) = pai(1)*s(4, 1, m + 1) - pci(1)*s(4, 1, m + 2) ! [dxz|o|s]
253 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 s(9, 1, m + 1) = pai(2)*s(4, 1, m + 1) - pci(2)*s(4, 1, m + 2) ! [dyz|o|s]
255 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 ELSE IF (llr == 3) THEN
258 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 + 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 + 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 + 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 + orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
267 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 + 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 + 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 + 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 + 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 + 2._dp*orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
278 END DO
279 ELSE IF (llr == 4) THEN
280 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 + 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 + 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 + 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 + 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 + 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 + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
293 s(27, 1, m + 1) = pai(1)*s(17, 1, m + 1) - pci(1)*s(17, 1, m + 2) ! [gxy3 |s|s]
294 s(28, 1, m + 1) = pai(1)*s(18, 1, m + 1) - pci(1)*s(18, 1, m + 2) ! [gxy2z|s|s]
295 s(29, 1, m + 1) = pai(1)*s(19, 1, m + 1) - pci(1)*s(19, 1, m + 2) ! [gxyz2|s|s]
296 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 + 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 + 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 + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
303 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 + 3._dp*orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
306 END DO
307 ELSE
308 DO irx = 0, llr
309 DO iry = 0, llr - irx
310 irz = llr - irx - iry
311 irr(1) = irx; irr(2) = iry; irr(3) = irz
312 ixx = maxloc(irr)
313 ix = ixx(1)
314 ir = coset(irx, iry, irz)
315 irm = irr
316 irm(ix) = irm(ix) - 1
317 aai = real(max(irm(ix), 0), dp)*orho
318 ir1 = coset(irm(1), irm(2), irm(3))
319 irm(ix) = irm(ix) - 1
320 ir2 = coset(irm(1), irm(2), irm(3))
321 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 + 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 DO mb = 1, lb_max
334 DO ibx = 0, mb
335 DO iby = 0, mb - ibx
336 ibz = mb - ibx - iby
337 ib = coset(ibx, iby, ibz)
338 ii(1) = ibx; ii(2) = iby; ii(3) = ibz
339 ixx = maxloc(ii)
340 ix = ixx(1)
341 abx = -rab(ix)
342 iim = ii
343 iim(ix) = iim(ix) - 1
344 ibm = coset(iim(1), iim(2), iim(3))
345 DO ia = 1, ncoset(mmax - mb)
346 iap = iiap(ia, ix)
347 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 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 DO llr = 1, lb_max
361 IF (llr == 1) THEN
362 DO m = 0, lb_max - llr
363 s1 = s(1, 1, m + 1)
364 s2 = s(1, 1, m + 2)
365 s(1, 2, m + 1) = pbi(1)*s1 - pci(1)*s2 ! [px|o|s]
366 s(1, 3, m + 1) = pbi(2)*s1 - pci(2)*s2 ! [py|o|s]
367 s(1, 4, m + 1) = pbi(3)*s1 - pci(3)*s2 ! [pz|o|s]
368 END DO
369 ELSE IF (llr == 2) THEN
370 DO m = 0, lb_max - llr
371 s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
372 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 s(1, 6, m + 1) = pbi(1)*s(1, 3, m + 1) - pci(1)*s(1, 3, m + 2) ! [dxy|o|s]
374 s(1, 7, m + 1) = pbi(1)*s(1, 4, m + 1) - pci(1)*s(1, 4, m + 2) ! [dxz|o|s]
375 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 s(1, 9, m + 1) = pbi(2)*s(1, 4, m + 1) - pci(2)*s(1, 4, m + 2) ! [dyz|o|s]
377 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 ELSE IF (llr == 3) THEN
380 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 + 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 + 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 + 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 + orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
389 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 + 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 + 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 + 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 + 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 + 2._dp*orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
400 END DO
401 ELSE
402 DO irx = 0, llr
403 DO iry = 0, llr - irx
404 irz = llr - irx - iry
405 irr(1) = irx; irr(2) = iry; irr(3) = irz
406 ixx = maxloc(irr)
407 ix = ixx(1)
408 ir = coset(irx, iry, irz)
409 irm = irr
410 irm(ix) = irm(ix) - 1
411 aai = real(max(irm(ix), 0), dp)
412 ir1 = coset(irm(1), irm(2), irm(3))
413 irm(ix) = irm(ix) - 1
414 ir2 = coset(irm(1), irm(2), irm(3))
415 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 + 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 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
428 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
429 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 DO da = 0, da_max - 1
437 ftz = 2.0_dp*zeta(ipgf)
438 DO dax = 0, da
439 DO day = 0, da - dax
440 daz = da - dax - day
441 cda = coset(dax, day, daz)
442 cdax = coset(dax + 1, day, daz)
443 cday = coset(dax, day + 1, daz)
444 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 DO la = la_min_set, la_max - da - 1
449 DO ax = 0, la
450 fax = real(ax, dp)
451 DO ay = 0, la - ax
452 fay = real(ay, dp)
453 az = la - ax - ay
454 faz = real(az, dp)
455 coa = coset(ax, ay, az)
456 coamx = coset(ax - 1, ay, az)
457 coamy = coset(ax, ay - 1, az)
458 coamz = coset(ax, ay, az - 1)
459 coapx = coset(ax + 1, ay, az)
460 coapy = coset(ax, ay + 1, az)
461 coapz = coset(ax, ay, az + 1)
462 DO cob = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
463 fs(coa, cob, cdax) = ftz*s(coapx, cob, cda) - fax*s(coamx, cob, cda)
464 fs(coa, cob, cday) = ftz*s(coapy, cob, cda) - fay*s(coamy, cob, cda)
465 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 IF (PRESENT(vab2_work)) THEN
477 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
478 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
479 vab2_work(na + i, nb + j, 1) = vab2_work(na + i, nb + j, 1) + fs(i, j, 2)
480 vab2_work(na + i, nb + j, 2) = vab2_work(na + i, nb + j, 2) + fs(i, j, 3)
481 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 IF (calculate_force_a) THEN
489 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
490 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
491 force_a(1) = force_a(1) + pab(na + i, nb + j)*fs(i, j, 2)
492 force_a(2) = force_a(2) + pab(na + i, nb + j)*fs(i, j, 3)
493 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 DO db = 0, db_max - 1
502 ftz = 2.0_dp*zetb(jpgf)
503 DO dbx = 0, db
504 DO dby = 0, db - dbx
505 dbz = db - dbx - dby
506 cdb = coset(dbx, dby, dbz)
507 cdbx = coset(dbx + 1, dby, dbz)
508 cdby = coset(dbx, dby + 1, dbz)
509 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 DO lb = lb_min_set, lb_max - db - 1
514 DO bx = 0, lb
515 fbx = real(bx, dp)
516 DO by = 0, lb - bx
517 fby = real(by, dp)
518 bz = lb - bx - by
519 fbz = real(bz, dp)
520 cob = coset(bx, by, bz)
521 cobmx = coset(bx - 1, by, bz)
522 cobmy = coset(bx, by - 1, bz)
523 cobmz = coset(bx, by, bz - 1)
524 cobpx = coset(bx + 1, by, bz)
525 cobpy = coset(bx, by + 1, bz)
526 cobpz = coset(bx, by, bz + 1)
527 DO coa = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
528 fs(coa, cob, cdbx) = ftz*s(coa, cobpx, cdb) - fbx*s(coa, cobmx, cdb)
529 fs(coa, cob, cdby) = ftz*s(coa, cobpy, cdb) - fby*s(coa, cobmy, cdb)
530 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 IF (PRESENT(vab2_work)) THEN
542 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
543 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
544 vab2_work(na + i, nb + j, 4) = vab2_work(na + i, nb + j, 4) + fs(i, j, 2)
545 vab2_work(na + i, nb + j, 5) = vab2_work(na + i, nb + j, 5) + fs(i, j, 3)
546 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 DO idir = 1, 3
551 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
552 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 - 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 IF (calculate_force_b) THEN
566 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
567 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
568 force_b(1) = force_b(1) + pab(na + i, nb + j)*fs(i, j, 2)
569 force_b(2) = force_b(2) + pab(na + i, nb + j)*fs(i, j, 3)
570 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 nb = nb + ncoset(lb_max_set)
576
577 END DO
578
579 na = na + ncoset(la_max_set)
580
581 END DO
582
583 DEALLOCATE (iiap)
584
585 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 SUBROUTINE os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
606 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: s
621
622 IF (PRESENT(dva)) THEN
623 da_max = 1
624 ELSE
625 da_max = 0
626 END IF
627
628 la_max = la_max_set + da_max
629 la_min = max(0, la_min_set - da_max)
630
631 mmax = la_max
632
633 ALLOCATE (s(ncoset(mmax), mmax + 1))
634 na = 0
635 DO ipgf = 1, npgfa
636 IF (rpgfa(ipgf) + rpgfc < dac) THEN
637 na = na + ncoset(la_max_set)
638 cycle
639 END IF
640 s(1, 1:mmax + 1) = auxint(0:mmax, ipgf)
641 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 orho = 0.5_dp/zeta(ipgf)
649
650 DO llr = 1, mmax
651 IF (llr == 1) THEN
652 DO m = 0, mmax - llr
653 s1 = s(1, m + 2)
654 s(2, m + 1) = -rac(1)*s1 ! [px|o]
655 s(3, m + 1) = -rac(2)*s1 ! [py|o]
656 s(4, m + 1) = -rac(3)*s1 ! [pz|o]
657 END DO
658 ELSE IF (llr == 2) THEN
659 DO m = 0, mmax - llr
660 s1 = s(1, m + 1) - s(1, m + 2)
661 s(5, m + 1) = -rac(1)*s(2, m + 2) + orho*s1 ! [dx2|o]
662 s(6, m + 1) = -rac(1)*s(3, m + 2) ! [dxy|o]
663 s(7, m + 1) = -rac(1)*s(4, m + 2) ! [dxz|o]
664 s(8, m + 1) = -rac(2)*s(3, m + 2) + orho*s1 ! [dy2|o]
665 s(9, m + 1) = -rac(2)*s(4, m + 2) ! [dyz|o]
666 s(10, m + 1) = -rac(3)*s(4, m + 2) + orho*s1 ! [dz2|o]
667 END DO
668 ELSE IF (llr == 3) THEN
669 DO m = 0, mmax - llr
670 s(11, m + 1) = -rac(1)*s(5, m + 2) + 2._dp*orho*(s(2, m + 1) - s(2, m + 2)) ! [fx3 |o]
671 s(12, m + 1) = -rac(1)*s(6, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fx2y|o]
672 s(13, m + 1) = -rac(1)*s(7, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fx2z|o]
673 s(14, m + 1) = -rac(2)*s(6, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxy2|o]
674 s(15, m + 1) = -rac(1)*s(9, m + 2) ! [fxyz|o]
675 s(16, m + 1) = -rac(3)*s(7, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxz2|o]
676 s(17, m + 1) = -rac(2)*s(8, m + 2) + 2._dp*orho*(s(3, m + 1) - s(3, m + 2)) ! [fy3 |o]
677 s(18, m + 1) = -rac(2)*s(9, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fy2z|o]
678 s(19, m + 1) = -rac(3)*s(9, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fyz2|o]
679 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 ELSE IF (llr == 4) THEN
682 DO m = 0, mmax - llr
683 s(21, m + 1) = -rac(1)*s(11, m + 2) + 3._dp*orho*(s(5, m + 1) - s(5, m + 2)) ! [gx4 |s]
684 s(22, m + 1) = -rac(1)*s(12, m + 2) + 2._dp*orho*(s(6, m + 1) - s(6, m + 2)) ! [gx3y |s]
685 s(23, m + 1) = -rac(1)*s(13, m + 2) + 2._dp*orho*(s(7, m + 1) - s(7, m + 2)) ! [gx3z |s]
686 s(24, m + 1) = -rac(1)*s(14, m + 2) + orho*(s(8, m + 1) - s(8, m + 2)) ! [gx2y2|s]
687 s(25, m + 1) = -rac(1)*s(15, m + 2) + orho*(s(9, m + 1) - s(9, m + 2)) ! [gx2yz|s]
688 s(26, m + 1) = -rac(1)*s(16, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gx2z2|s]
689 s(27, m + 1) = -rac(1)*s(17, m + 2) ! [gxy3 |s]
690 s(28, m + 1) = -rac(1)*s(18, m + 2) ! [gxy2z|s]
691 s(29, m + 1) = -rac(1)*s(19, m + 2) ! [gxyz2|s]
692 s(30, m + 1) = -rac(1)*s(20, m + 2) ! [gxz3 |s]
693 s(31, m + 1) = -rac(2)*s(17, m + 2) + 3._dp*orho*(s(8, m + 1) - s(8, m + 2)) ! [gy4 |s]
694 s(32, m + 1) = -rac(2)*s(18, m + 2) + 2._dp*orho*(s(9, m + 1) - s(9, m + 2)) ! [gy3z |s]
695 s(33, m + 1) = -rac(2)*s(19, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gy2z2|s]
696 s(34, m + 1) = -rac(2)*s(20, m + 2) ! [gyz3 |s]
697 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 DO irx = 0, llr
701 DO iry = 0, llr - irx
702 irz = llr - irx - iry
703 irr(1) = irx; irr(2) = iry; irr(3) = irz
704 ixx = maxloc(irr)
705 ix = ixx(1)
706 ir = coset(irx, iry, irz)
707 irm = irr
708 irm(ix) = irm(ix) - 1
709 aai = real(max(irm(ix), 0), dp)*orho
710 ir1 = coset(irm(1), irm(2), irm(3))
711 irm(ix) = irm(ix) - 1
712 ir2 = coset(irm(1), irm(2), irm(3))
713 DO m = 0, mmax - llr
714 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 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
725 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 IF (PRESENT(dva)) THEN
732 ftz = 2.0_dp*zeta(ipgf)
733 DO la = la_min_set, la_max_set
734 DO ax = 0, la
735 fax = real(ax, dp)
736 DO ay = 0, la - ax
737 fay = real(ay, dp)
738 az = la - ax - ay
739 faz = real(az, dp)
740 coa = coset(ax, ay, az)
741 coamx = coset(ax - 1, ay, az)
742 coamy = coset(ax, ay - 1, az)
743 coamz = coset(ax, ay, az - 1)
744 coapx = coset(ax + 1, ay, az)
745 coapy = coset(ax, ay + 1, az)
746 coapz = coset(ax, ay, az + 1)
747 dva(na + coa, 1) = dva(na + coa, 1) + ftz*s(coapx, 1) - fax*s(coamx, 1)
748 dva(na + coa, 2) = dva(na + coa, 2) + ftz*s(coapy, 1) - fay*s(coamy, 1)
749 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 na = na + ncoset(la_max_set)
756
757 END DO
758
759 DEALLOCATE (s)
760
761 END SUBROUTINE os_2center
762! **************************************************************************************************
763
764END MODULE ai_oneelectron
ai_oneelectron
Calculation of general three-center integrals over Cartesian Gaussian-type functions and a spherical ...
Definition ai_oneelectron.F:34
ai_oneelectron::os_3center
subroutine, public os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, vab2, vab2_work, deltar, iatom, jatom, katom)
Calculation of three-center integrals <a|c|b> over Cartesian Gaussian functions and a spherical poten...
Definition ai_oneelectron.F:96
ai_oneelectron::os_2center
subroutine, public os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, auxint, rpgfc, rac, dac, va, dva)
Calculation of two-center integrals <a|c> over Cartesian Gaussian functions and a spherical potential...
Definition ai_oneelectron.F:607
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
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