(git:ccc2433)
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 ! **************************************************************************************************
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  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 
764 END 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_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
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
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