(git:ccc2433)
ai_coulomb.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Calculation of Coulomb integrals over Cartesian Gaussian-type functions
10 !> (electron repulsion integrals, ERIs).
11 !> \par Literature
12 !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 !> \par History
14 !> none
15 !> \par Parameters
16 !> - ax,ay,az : Angular momentum index numbers of orbital a.
17 !> - bx,by,bz : Angular momentum index numbers of orbital b.
18 !> - cx,cy,cz : Angular momentum index numbers of orbital c.
19 !> - coset : Cartesian orbital set pointer.
20 !> - dab : Distance between the atomic centers a and b.
21 !> - dac : Distance between the atomic centers a and c.
22 !> - dbc : Distance between the atomic centers b and c.
23 !> - gccc : Prefactor of the primitive Gaussian function c.
24 !> - l{a,b,c} : Angular momentum quantum number of shell a, b or c.
25 !> - l{a,b,c}_max: Maximum angular momentum quantum number of shell a, b or c.
26 !> - l{a,b,c}_min: Minimum angular momentum quantum number of shell a, b or c.
27 !> - ncoset : Number of orbitals in a Cartesian orbital set.
28 !> - npgf{a,b} : Degree of contraction of shell a or b.
29 !> - rab : Distance vector between the atomic centers a and b.
30 !> - rab2 : Square of the distance between the atomic centers a and b.
31 !> - rac : Distance vector between the atomic centers a and c.
32 !> - rac2 : Square of the distance between the atomic centers a and c.
33 !> - rbc : Distance vector between the atomic centers b and c.
34 !> - rbc2 : Square of the distance between the atomic centers b and c.
35 !> - rpgf{a,b,c} : Radius of the primitive Gaussian-type function a, b or c.
36 !> - zet{a,b,c} : Exponents of the Gaussian-type functions a, b or c.
37 !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
38 !> - zetw : Reciprocal of the sum of the exponents of orbital a, b and c.
39 !> \author Matthias Krack (22.08.2000)
40 ! **************************************************************************************************
41 MODULE ai_coulomb
42 
43  USE gamma, ONLY: fgamma => fgamma_0
44  USE kinds, ONLY: dp
45  USE mathconstants, ONLY: pi
46  USE orbital_pointers, ONLY: coset,&
47  ncoset
48 #include "../base/base_uses.f90"
49 
50  IMPLICIT NONE
51  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_coulomb'
52  PRIVATE
53 
54  ! *** Public subroutines ***
55 
56  PUBLIC :: coulomb2, coulomb2_new, coulomb3
57 
58 CONTAINS
59 
60 ! **************************************************************************************************
61 !> \brief Calculation of the primitive two-center Coulomb integrals over
62 !> Cartesian Gaussian-type functions.
63 !> \param la_max ...
64 !> \param npgfa ...
65 !> \param zeta ...
66 !> \param rpgfa ...
67 !> \param la_min ...
68 !> \param lc_max ...
69 !> \param npgfc ...
70 !> \param zetc ...
71 !> \param rpgfc ...
72 !> \param lc_min ...
73 !> \param rac ...
74 !> \param rac2 ...
75 !> \param vac ...
76 !> \param v ...
77 !> \param f ...
78 !> \param maxder ...
79 !> \param vac_plus ...
80 !> \date 05.12.2000
81 !> \author Matthias Krack
82 !> \version 1.0
83 ! **************************************************************************************************
84  SUBROUTINE coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lc_max, npgfc, zetc, rpgfc, lc_min, &
85  rac, rac2, vac, v, f, maxder, vac_plus)
86  INTEGER, INTENT(IN) :: la_max, npgfa
87  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
88  INTEGER, INTENT(IN) :: la_min, lc_max, npgfc
89  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
90  INTEGER, INTENT(IN) :: lc_min
91  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
92  REAL(kind=dp), INTENT(IN) :: rac2
93  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
94  REAL(kind=dp), DIMENSION(:, :, :) :: v
95  REAL(kind=dp), DIMENSION(0:) :: f
96  INTEGER, INTENT(IN), OPTIONAL :: maxder
97  REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vac_plus
98 
99  INTEGER :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
100  cy, cz, i, ipgf, j, jpgf, la, lc, &
101  maxder_local, n, na, nap, nc, nmax
102  REAL(kind=dp) :: dac, f0, f1, f2, f3, f4, f5, f6, fcx, &
103  fcy, fcz, rho, t, zetp, zetq, zetw
104  REAL(kind=dp), DIMENSION(3) :: raw, rcw
105 
106  v = 0.0_dp
107 
108  maxder_local = 0
109  IF (PRESENT(maxder)) THEN
110  maxder_local = maxder
111  vac_plus = 0.0_dp
112  END IF
113 
114  nmax = la_max + lc_max + 1
115 
116  ! *** Calculate the distance of the centers a and c ***
117 
118  dac = sqrt(rac2)
119 
120  ! *** Loop over all pairs of primitive Gaussian-type functions ***
121 
122  na = 0
123  nap = 0
124 
125  DO ipgf = 1, npgfa
126 
127  nc = 0
128 
129  DO jpgf = 1, npgfc
130 
131  ! *** Screening ***
132 
133  IF (rpgfa(ipgf) + rpgfc(jpgf) < dac) THEN
134  DO j = nc + ncoset(lc_min - 1) + 1, nc + ncoset(lc_max)
135  DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max - maxder_local)
136  vac(i, j) = 0.0_dp
137  END DO
138  END DO
139  nc = nc + ncoset(lc_max)
140  cycle
141  END IF
142 
143  ! *** Calculate some prefactors ***
144 
145  zetp = 1.0_dp/zeta(ipgf)
146  zetq = 1.0_dp/zetc(jpgf)
147  zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
148 
149  rho = zeta(ipgf)*zetc(jpgf)*zetw
150 
151  f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
152 
153  ! *** Calculate the incomplete Gamma function ***
154 
155  t = rho*rac2
156 
157  CALL fgamma(nmax - 1, t, f)
158 
159  ! *** Calculate the basic two-center Coulomb integrals [s||s]{n} ***
160 
161  DO n = 1, nmax
162  v(1, 1, n) = f0*f(n - 1)
163  END DO
164 
165  ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
166 
167  IF (lc_max > 0) THEN
168 
169  f1 = 0.5_dp*zetq
170  f2 = -rho*zetq
171 
172  rcw(:) = -zeta(ipgf)*zetw*rac(:)
173 
174  ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
175 
176  DO n = 1, nmax - 1
177  v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
178  v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
179  v(1, 4, n) = rcw(3)*v(1, 1, n + 1)
180  END DO
181 
182  ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + ***
183  ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + ***
184  ! ** f2*[s||c-2i]{n+1} ***
185 
186  DO lc = 2, lc_max
187 
188  DO n = 1, nmax - lc
189 
190  ! **** Increase the angular momentum component z of c ***
191 
192  v(1, coset(0, 0, lc), n) = &
193  rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + &
194  f1*real(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + &
195  f2*v(1, coset(0, 0, lc - 2), n + 1))
196 
197  ! *** Increase the angular momentum component y of c ***
198 
199  cz = lc - 1
200  v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
201 
202  DO cy = 2, lc
203  cz = lc - cy
204  v(1, coset(0, cy, cz), n) = &
205  rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + &
206  f1*real(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + &
207  f2*v(1, coset(0, cy - 2, cz), n + 1))
208  END DO
209 
210  ! *** Increase the angular momentum component x of c ***
211 
212  DO cy = 0, lc - 1
213  cz = lc - 1 - cy
214  v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
215  END DO
216 
217  DO cx = 2, lc
218  f6 = f1*real(cx - 1, dp)
219  DO cy = 0, lc - cx
220  cz = lc - cx - cy
221  v(1, coset(cx, cy, cz), n) = &
222  rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + &
223  f6*(v(1, coset(cx - 2, cy, cz), n) + &
224  f2*v(1, coset(cx - 2, cy, cz), n + 1))
225  END DO
226  END DO
227 
228  END DO
229 
230  END DO
231 
232  END IF
233 
234  ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
235 
236  IF (la_max > 0) THEN
237 
238  f3 = 0.5_dp*zetp
239  f4 = -rho*zetp
240  f5 = 0.5_dp*zetw
241 
242  raw(:) = zetc(jpgf)*zetw*rac(:)
243 
244  ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
245 
246  DO n = 1, nmax - 1
247  v(2, 1, n) = raw(1)*v(1, 1, n + 1)
248  v(3, 1, n) = raw(2)*v(1, 1, n + 1)
249  v(4, 1, n) = raw(3)*v(1, 1, n + 1)
250  END DO
251 
252  ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + ***
253  ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + ***
254  ! *** f4*[a-2i||s]{n+1}) ***
255 
256  DO la = 2, la_max
257 
258  DO n = 1, nmax - la
259 
260  ! *** Increase the angular momentum component z of a ***
261 
262  v(coset(0, 0, la), 1, n) = &
263  raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + &
264  f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + &
265  f4*v(coset(0, 0, la - 2), 1, n + 1))
266 
267  ! *** Increase the angular momentum component y of a ***
268 
269  az = la - 1
270  v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
271 
272  DO ay = 2, la
273  az = la - ay
274  v(coset(0, ay, az), 1, n) = &
275  raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + &
276  f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + &
277  f4*v(coset(0, ay - 2, az), 1, n + 1))
278  END DO
279 
280  ! *** Increase the angular momentum component x of a ***
281 
282  DO ay = 0, la - 1
283  az = la - 1 - ay
284  v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
285  END DO
286 
287  DO ax = 2, la
288  f6 = f3*real(ax - 1, dp)
289  DO ay = 0, la - ax
290  az = la - ax - ay
291  v(coset(ax, ay, az), 1, n) = &
292  raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + &
293  f6*(v(coset(ax - 2, ay, az), 1, n) + &
294  f4*v(coset(ax - 2, ay, az), 1, n + 1))
295  END DO
296  END DO
297 
298  END DO
299 
300  END DO
301 
302  DO lc = 1, lc_max
303 
304  DO cx = 0, lc
305  DO cy = 0, lc - cx
306  cz = lc - cx - cy
307 
308  coc = coset(cx, cy, cz)
309  cocx = coset(max(0, cx - 1), cy, cz)
310  cocy = coset(cx, max(0, cy - 1), cz)
311  cocz = coset(cx, cy, max(0, cz - 1))
312 
313  fcx = f5*real(cx, dp)
314  fcy = f5*real(cy, dp)
315  fcz = f5*real(cz, dp)
316 
317  ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
318  ! *** f5*Ni(c)*[s||c-1i]{n+1} ***
319 
320  DO n = 1, nmax - 1 - lc
321  v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
322  v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
323  v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1)
324  END DO
325 
326  ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + ***
327  ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + ***
328  ! *** f4*[a-2i||c]{n+1}) + ***
329  ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} ***
330 
331  DO la = 2, la_max
332 
333  DO n = 1, nmax - la - lc
334 
335  ! *** Increase the angular momentum component z of a ***
336 
337  v(coset(0, 0, la), coc, n) = &
338  raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + &
339  f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + &
340  f4*v(coset(0, 0, la - 2), coc, n + 1)) + &
341  fcz*v(coset(0, 0, la - 1), cocz, n + 1)
342 
343  ! *** Increase the angular momentum component y of a ***
344 
345  az = la - 1
346  v(coset(0, 1, az), coc, n) = &
347  raw(2)*v(coset(0, 0, az), coc, n + 1) + &
348  fcy*v(coset(0, 0, az), cocy, n + 1)
349 
350  DO ay = 2, la
351  az = la - ay
352  v(coset(0, ay, az), coc, n) = &
353  raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + &
354  f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + &
355  f4*v(coset(0, ay - 2, az), coc, n + 1)) + &
356  fcy*v(coset(0, ay - 1, az), cocy, n + 1)
357  END DO
358 
359  ! *** Increase the angular momentum component x of a ***
360 
361  DO ay = 0, la - 1
362  az = la - 1 - ay
363  v(coset(1, ay, az), coc, n) = &
364  raw(1)*v(coset(0, ay, az), coc, n + 1) + &
365  fcx*v(coset(0, ay, az), cocx, n + 1)
366  END DO
367 
368  DO ax = 2, la
369  f6 = f3*real(ax - 1, dp)
370  DO ay = 0, la - ax
371  az = la - ax - ay
372  v(coset(ax, ay, az), coc, n) = &
373  raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + &
374  f6*(v(coset(ax - 2, ay, az), coc, n) + &
375  f4*v(coset(ax - 2, ay, az), coc, n + 1)) + &
376  fcx*v(coset(ax - 1, ay, az), cocx, n + 1)
377  END DO
378  END DO
379 
380  END DO
381 
382  END DO
383 
384  END DO
385  END DO
386 
387  END DO
388 
389  END IF
390 
391  DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max)
392  DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
393  vac(na + i, nc + j) = v(i, j, 1)
394  END DO
395  END DO
396 
397  IF (PRESENT(maxder)) THEN
398  DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max)
399  DO i = 1, ncoset(la_max)
400  vac_plus(nap + i, nc + j) = v(i, j, 1)
401  END DO
402  END DO
403  END IF
404 
405  nc = nc + ncoset(lc_max)
406 
407  END DO
408 
409  na = na + ncoset(la_max - maxder_local)
410  nap = nap + ncoset(la_max)
411 
412  END DO
413 
414  END SUBROUTINE coulomb2
415 ! **************************************************************************************************
416 !> \brief Calculation of the primitive two-center Coulomb integrals over
417 !> Cartesian Gaussian-type functions. Same as coulomb2 treats derivative
418 !> different (now vac_plus is symmetric)
419 !> \param la_max ...
420 !> \param npgfa ...
421 !> \param zeta ...
422 !> \param la_min ...
423 !> \param lc_max ...
424 !> \param npgfc ...
425 !> \param zetc ...
426 !> \param lc_min ...
427 !> \param rac ...
428 !> \param rac2 ...
429 !> \param vac ...
430 !> \param v ...
431 !> \param f ...
432 !> \param maxder ...
433 !> \param vac_plus ...
434 !> \date 6.02.2008
435 !> \author CJM
436 !> \version 1.0
437 ! **************************************************************************************************
438 
439  SUBROUTINE coulomb2_new(la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, &
440  rac, rac2, vac, v, f, maxder, vac_plus)
441  INTEGER, INTENT(IN) :: la_max, npgfa
442  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
443  INTEGER, INTENT(IN) :: la_min, lc_max, npgfc
444  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetc
445  INTEGER, INTENT(IN) :: lc_min
446  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
447  REAL(kind=dp), INTENT(IN) :: rac2
448  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
449  REAL(kind=dp), DIMENSION(:, :, :) :: v
450  REAL(kind=dp), DIMENSION(0:) :: f
451  INTEGER, INTENT(IN), OPTIONAL :: maxder
452  REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vac_plus
453 
454  INTEGER :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
455  cy, cz, i, ipgf, j, jpgf, la, lc, &
456  maxder_local, n, na, nap, nc, ncp, nmax
457  REAL(kind=dp) :: dac, f0, f1, f2, f3, f4, f5, f6, fcx, &
458  fcy, fcz, rho, t, zetp, zetq, zetw
459  REAL(kind=dp), DIMENSION(3) :: raw, rcw
460 
461  v = 0.0_dp
462 
463  maxder_local = 0
464  IF (PRESENT(maxder)) THEN
465  maxder_local = maxder
466  vac_plus = 0.0_dp
467  END IF
468 
469  nmax = la_max + lc_max + 1
470 
471  ! *** Calculate the distance of the centers a and c ***
472 
473  dac = sqrt(rac2)
474 
475  ! *** Loop over all pairs of primitive Gaussian-type functions ***
476 
477  na = 0
478  nap = 0
479 
480  DO ipgf = 1, npgfa
481 
482  nc = 0
483  ncp = 0
484 
485  DO jpgf = 1, npgfc
486 
487  ! *** Calculate some prefactors ***
488 
489  zetp = 1.0_dp/zeta(ipgf)
490  zetq = 1.0_dp/zetc(jpgf)
491  zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
492 
493  rho = zeta(ipgf)*zetc(jpgf)*zetw
494 
495  f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
496 
497  ! *** Calculate the incomplete Gamma function ***
498 
499  t = rho*rac2
500 
501  CALL fgamma(nmax - 1, t, f)
502 
503  ! *** Calculate the basic two-center Coulomb integrals [s||s]{n} ***
504 
505  DO n = 1, nmax
506  v(1, 1, n) = f0*f(n - 1)
507  END DO
508 
509  ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
510 
511  IF (lc_max > 0) THEN
512 
513  f1 = 0.5_dp*zetq
514  f2 = -rho*zetq
515 
516  rcw(:) = -zeta(ipgf)*zetw*rac(:)
517 
518  ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
519 
520  DO n = 1, nmax - 1
521  v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
522  v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
523  v(1, 4, n) = rcw(3)*v(1, 1, n + 1)
524  END DO
525 
526  ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + ***
527  ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + ***
528  ! ** f2*[s||c-2i]{n+1} ***
529 
530  DO lc = 2, lc_max
531 
532  DO n = 1, nmax - lc
533 
534  ! **** Increase the angular momentum component z of c ***
535 
536  v(1, coset(0, 0, lc), n) = &
537  rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + &
538  f1*real(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + &
539  f2*v(1, coset(0, 0, lc - 2), n + 1))
540 
541  ! *** Increase the angular momentum component y of c ***
542 
543  cz = lc - 1
544  v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
545 
546  DO cy = 2, lc
547  cz = lc - cy
548  v(1, coset(0, cy, cz), n) = &
549  rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + &
550  f1*real(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + &
551  f2*v(1, coset(0, cy - 2, cz), n + 1))
552  END DO
553 
554  ! *** Increase the angular momentum component x of c ***
555 
556  DO cy = 0, lc - 1
557  cz = lc - 1 - cy
558  v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
559  END DO
560 
561  DO cx = 2, lc
562  f6 = f1*real(cx - 1, dp)
563  DO cy = 0, lc - cx
564  cz = lc - cx - cy
565  v(1, coset(cx, cy, cz), n) = &
566  rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + &
567  f6*(v(1, coset(cx - 2, cy, cz), n) + &
568  f2*v(1, coset(cx - 2, cy, cz), n + 1))
569  END DO
570  END DO
571 
572  END DO
573 
574  END DO
575 
576  END IF
577 
578  ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
579 
580  IF (la_max > 0) THEN
581 
582  f3 = 0.5_dp*zetp
583  f4 = -rho*zetp
584  f5 = 0.5_dp*zetw
585 
586  raw(:) = zetc(jpgf)*zetw*rac(:)
587 
588  ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
589 
590  DO n = 1, nmax - 1
591  v(2, 1, n) = raw(1)*v(1, 1, n + 1)
592  v(3, 1, n) = raw(2)*v(1, 1, n + 1)
593  v(4, 1, n) = raw(3)*v(1, 1, n + 1)
594  END DO
595 
596  ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + ***
597  ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + ***
598  ! *** f4*[a-2i||s]{n+1}) ***
599 
600  DO la = 2, la_max
601 
602  DO n = 1, nmax - la
603 
604  ! *** Increase the angular momentum component z of a ***
605 
606  v(coset(0, 0, la), 1, n) = &
607  raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + &
608  f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + &
609  f4*v(coset(0, 0, la - 2), 1, n + 1))
610 
611  ! *** Increase the angular momentum component y of a ***
612 
613  az = la - 1
614  v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
615 
616  DO ay = 2, la
617  az = la - ay
618  v(coset(0, ay, az), 1, n) = &
619  raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + &
620  f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + &
621  f4*v(coset(0, ay - 2, az), 1, n + 1))
622  END DO
623 
624  ! *** Increase the angular momentum component x of a ***
625 
626  DO ay = 0, la - 1
627  az = la - 1 - ay
628  v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
629  END DO
630 
631  DO ax = 2, la
632  f6 = f3*real(ax - 1, dp)
633  DO ay = 0, la - ax
634  az = la - ax - ay
635  v(coset(ax, ay, az), 1, n) = &
636  raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + &
637  f6*(v(coset(ax - 2, ay, az), 1, n) + &
638  f4*v(coset(ax - 2, ay, az), 1, n + 1))
639  END DO
640  END DO
641 
642  END DO
643 
644  END DO
645 
646  DO lc = 1, lc_max
647 
648  DO cx = 0, lc
649  DO cy = 0, lc - cx
650  cz = lc - cx - cy
651 
652  coc = coset(cx, cy, cz)
653  cocx = coset(max(0, cx - 1), cy, cz)
654  cocy = coset(cx, max(0, cy - 1), cz)
655  cocz = coset(cx, cy, max(0, cz - 1))
656 
657  fcx = f5*real(cx, dp)
658  fcy = f5*real(cy, dp)
659  fcz = f5*real(cz, dp)
660 
661  ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
662  ! *** f5*Ni(c)*[s||c-1i]{n+1} ***
663 
664  DO n = 1, nmax - 1 - lc
665  v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
666  v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
667  v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1)
668  END DO
669 
670  ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + ***
671  ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + ***
672  ! *** f4*[a-2i||c]{n+1}) + ***
673  ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} ***
674 
675  DO la = 2, la_max
676 
677  DO n = 1, nmax - la - lc
678 
679  ! *** Increase the angular momentum component z of a ***
680 
681  v(coset(0, 0, la), coc, n) = &
682  raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + &
683  f3*real(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + &
684  f4*v(coset(0, 0, la - 2), coc, n + 1)) + &
685  fcz*v(coset(0, 0, la - 1), cocz, n + 1)
686 
687  ! *** Increase the angular momentum component y of a ***
688 
689  az = la - 1
690  v(coset(0, 1, az), coc, n) = &
691  raw(2)*v(coset(0, 0, az), coc, n + 1) + &
692  fcy*v(coset(0, 0, az), cocy, n + 1)
693 
694  DO ay = 2, la
695  az = la - ay
696  v(coset(0, ay, az), coc, n) = &
697  raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + &
698  f3*real(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + &
699  f4*v(coset(0, ay - 2, az), coc, n + 1)) + &
700  fcy*v(coset(0, ay - 1, az), cocy, n + 1)
701  END DO
702 
703  ! *** Increase the angular momentum component x of a ***
704 
705  DO ay = 0, la - 1
706  az = la - 1 - ay
707  v(coset(1, ay, az), coc, n) = &
708  raw(1)*v(coset(0, ay, az), coc, n + 1) + &
709  fcx*v(coset(0, ay, az), cocx, n + 1)
710  END DO
711 
712  DO ax = 2, la
713  f6 = f3*real(ax - 1, dp)
714  DO ay = 0, la - ax
715  az = la - ax - ay
716  v(coset(ax, ay, az), coc, n) = &
717  raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + &
718  f6*(v(coset(ax - 2, ay, az), coc, n) + &
719  f4*v(coset(ax - 2, ay, az), coc, n + 1)) + &
720  fcx*v(coset(ax - 1, ay, az), cocx, n + 1)
721  END DO
722  END DO
723 
724  END DO
725 
726  END DO
727 
728  END DO
729  END DO
730 
731  END DO
732 
733  END IF
734 
735  DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max - maxder_local)
736  DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
737  vac(na + i, nc + j) = v(i, j, 1)
738  END DO
739  END DO
740 
741  IF (PRESENT(maxder)) THEN
742  DO j = 1, ncoset(lc_max)
743  DO i = 1, ncoset(la_max)
744  vac_plus(nap + i, ncp + j) = v(i, j, 1)
745  END DO
746  END DO
747  END IF
748 
749  nc = nc + ncoset(lc_max - maxder_local)
750  ncp = ncp + ncoset(lc_max)
751 
752  END DO
753 
754  na = na + ncoset(la_max - maxder_local)
755  nap = nap + ncoset(la_max)
756 
757  END DO
758 
759  END SUBROUTINE coulomb2_new
760 
761 ! **************************************************************************************************
762 !> \brief Calculation of the primitive three-center Coulomb integrals over
763 !> Cartesian Gaussian-type functions (electron repulsion integrals,
764 !> ERIs).
765 !> \param la_max ...
766 !> \param npgfa ...
767 !> \param zeta ...
768 !> \param rpgfa ...
769 !> \param la_min ...
770 !> \param lb_max ...
771 !> \param npgfb ...
772 !> \param zetb ...
773 !> \param rpgfb ...
774 !> \param lb_min ...
775 !> \param lc_max ...
776 !> \param zetc ...
777 !> \param rpgfc ...
778 !> \param lc_min ...
779 !> \param gccc ...
780 !> \param rab ...
781 !> \param rab2 ...
782 !> \param rac ...
783 !> \param rac2 ...
784 !> \param rbc2 ...
785 !> \param vabc ...
786 !> \param int_abc ...
787 !> \param v ...
788 !> \param f ...
789 !> \param maxder ...
790 !> \param vabc_plus ...
791 !> \date 06.11.2000
792 !> \author Matthias Krack
793 !> \version 1.0
794 ! **************************************************************************************************
795  SUBROUTINE coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
796  lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, &
797  v, f, maxder, vabc_plus)
798  INTEGER, INTENT(IN) :: la_max, npgfa
799  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
800  INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
801  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
802  INTEGER, INTENT(IN) :: lb_min, lc_max
803  REAL(kind=dp), INTENT(IN) :: zetc, rpgfc
804  INTEGER, INTENT(IN) :: lc_min
805  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: gccc
806  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
807  REAL(kind=dp), INTENT(IN) :: rab2
808  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rac
809  REAL(kind=dp), INTENT(IN) :: rac2, rbc2
810  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vabc
811  REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT) :: int_abc
812  REAL(kind=dp), DIMENSION(:, :, :, :) :: v
813  REAL(kind=dp), DIMENSION(0:) :: f
814  INTEGER, INTENT(IN), OPTIONAL :: maxder
815  REAL(kind=dp), DIMENSION(:, :), OPTIONAL :: vabc_plus
816 
817  INTEGER :: ax, ay, az, bx, by, bz, coc, cocx, cocy, &
818  cocz, cx, cy, cz, i, ipgf, j, jpgf, k, &
819  kk, la, la_start, lb, lc, &
820  maxder_local, n, na, nap, nb, nmax
821  REAL(kind=dp) :: dab, dac, dbc, f0, f1, f2, f3, f4, f5, &
822  f6, f7, fcx, fcy, fcz, fx, fy, fz, t, &
823  zetp, zetq, zetw
824  REAL(kind=dp), DIMENSION(3) :: rap, rbp, rcp, rcw, rpw
825 
826  v = 0.0_dp
827 
828  maxder_local = 0
829  IF (PRESENT(maxder)) THEN
830  maxder_local = maxder
831  END IF
832 
833  nmax = la_max + lb_max + lc_max + 1
834 
835  ! *** Calculate the distances of the centers a, b and c ***
836 
837  dab = sqrt(rab2)
838  dac = sqrt(rac2)
839  dbc = sqrt(rbc2)
840 
841  ! *** Initialize integrals array
842  int_abc = 0.0_dp
843 
844  ! *** Loop over all pairs of primitive Gaussian-type functions ***
845 
846  na = 0
847  nap = 0
848 
849  DO ipgf = 1, npgfa
850 
851  ! *** Screening ***
852  IF (rpgfa(ipgf) + rpgfc < dac) THEN
853  na = na + ncoset(la_max - maxder_local)
854  nap = nap + ncoset(la_max)
855  cycle
856  END IF
857 
858  nb = 0
859 
860  DO jpgf = 1, npgfb
861 
862  ! *** Screening ***
863  IF ( &
864  (rpgfb(jpgf) + rpgfc < dbc) .OR. &
865  (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
866  nb = nb + ncoset(lb_max)
867  cycle
868  END IF
869 
870  ! *** Calculate some prefactors ***
871 
872  zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
873  zetq = 1.0_dp/zetc
874  zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
875 
876  f0 = 2.0_dp*sqrt(pi**5*zetw)*zetp*zetq
877  f1 = zetb(jpgf)*zetp
878  f2 = 0.5_dp*zetp
879  f4 = -zetc*zetw
880 
881  f0 = f0*exp(-zeta(ipgf)*f1*rab2)
882 
883  rap(:) = f1*rab(:)
884  rcp(:) = rap(:) - rac(:)
885  rpw(:) = f4*rcp(:)
886 
887  ! *** Calculate the incomplete Gamma function ***
888 
889  t = -f4*(rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3))/zetp
890 
891  CALL fgamma(nmax - 1, t, f)
892 
893  ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
894 
895  DO n = 1, nmax
896  v(1, 1, 1, n) = f0*f(n - 1)
897  END DO
898 
899  ! *** Recurrence steps: [ss||s] -> [as||s] ***
900 
901  IF (la_max > 0) THEN
902 
903  ! *** Vertical recurrence steps: [ss||s] -> [as||s] ***
904 
905  ! *** [ps||s]{n} = (Pi - Ai)*[ss||s]{n} + ***
906  ! *** (Wi - Pi)*[ss||s]{n+1} (i = x,y,z) ***
907 
908  DO n = 1, nmax - 1
909  v(2, 1, 1, n) = rap(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
910  v(3, 1, 1, n) = rap(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
911  v(4, 1, 1, n) = rap(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
912  END DO
913 
914  ! *** [as||s]{n} = (Pi - Ai)*[(a-1i)s||s]{n} + ***
915  ! *** (Wi - Pi)*[(a-1i)s||s]{n+1} + ***
916  ! *** f2*Ni(a-1i)*( [(a-2i)s||s]{n} + ***
917  ! *** f4*[(a-2i)s||s]{n+1}) ***
918 
919  DO la = 2, la_max
920 
921  DO n = 1, nmax - la
922 
923  ! *** Increase the angular momentum component z of a ***
924 
925  v(coset(0, 0, la), 1, 1, n) = &
926  rap(3)*v(coset(0, 0, la - 1), 1, 1, n) + &
927  rpw(3)*v(coset(0, 0, la - 1), 1, 1, n + 1) + &
928  f2*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, 1, n) + &
929  f4*v(coset(0, 0, la - 2), 1, 1, n + 1))
930 
931  ! *** Increase the angular momentum component y of a ***
932 
933  az = la - 1
934  v(coset(0, 1, az), 1, 1, n) = &
935  rap(2)*v(coset(0, 0, az), 1, 1, n) + &
936  rpw(2)*v(coset(0, 0, az), 1, 1, n + 1)
937 
938  DO ay = 2, la
939  az = la - ay
940  v(coset(0, ay, az), 1, 1, n) = &
941  rap(2)*v(coset(0, ay - 1, az), 1, 1, n) + &
942  rpw(2)*v(coset(0, ay - 1, az), 1, 1, n + 1) + &
943  f2*real(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, 1, n) + &
944  f4*v(coset(0, ay - 2, az), 1, 1, n + 1))
945  END DO
946 
947  ! *** Increase the angular momentum component x of a ***
948 
949  DO ay = 0, la - 1
950  az = la - 1 - ay
951  v(coset(1, ay, az), 1, 1, n) = &
952  rap(1)*v(coset(0, ay, az), 1, 1, n) + &
953  rpw(1)*v(coset(0, ay, az), 1, 1, n + 1)
954  END DO
955 
956  DO ax = 2, la
957  f3 = f2*real(ax - 1, dp)
958  DO ay = 0, la - ax
959  az = la - ax - ay
960  v(coset(ax, ay, az), 1, 1, n) = &
961  rap(1)*v(coset(ax - 1, ay, az), 1, 1, n) + &
962  rpw(1)*v(coset(ax - 1, ay, az), 1, 1, n + 1) + &
963  f3*(v(coset(ax - 2, ay, az), 1, 1, n) + &
964  f4*v(coset(ax - 2, ay, az), 1, 1, n + 1))
965  END DO
966  END DO
967 
968  END DO
969 
970  END DO
971 
972  ! *** Recurrence steps: [as||s] -> [ab||s] ***
973 
974  IF (lb_max > 0) THEN
975 
976  ! *** Horizontal recurrence steps ***
977 
978  rbp(:) = rap(:) - rab(:)
979 
980  ! *** [ap||s]{n} = [(a+1i)s||s]{n} - (Bi - Ai)*[as||s]{n} ***
981 
982  la_start = max(0, la_min - 1)
983 
984  DO la = la_start, la_max - 1
985  DO n = 1, nmax - la - 1
986  DO ax = 0, la
987  DO ay = 0, la - ax
988  az = la - ax - ay
989  v(coset(ax, ay, az), 2, 1, n) = &
990  v(coset(ax + 1, ay, az), 1, 1, n) - &
991  rab(1)*v(coset(ax, ay, az), 1, 1, n)
992  v(coset(ax, ay, az), 3, 1, n) = &
993  v(coset(ax, ay + 1, az), 1, 1, n) - &
994  rab(2)*v(coset(ax, ay, az), 1, 1, n)
995  v(coset(ax, ay, az), 4, 1, n) = &
996  v(coset(ax, ay, az + 1), 1, 1, n) - &
997  rab(3)*v(coset(ax, ay, az), 1, 1, n)
998  END DO
999  END DO
1000  END DO
1001  END DO
1002 
1003  ! *** Vertical recurrence step ***
1004 
1005  ! *** [ap||s]{n} = (Pi - Bi)*[as||s]{n} + ***
1006  ! *** (Wi - Pi)*[as||s]{n+1} + ***
1007  ! *** f2*Ni(a)*( [(a-1i)s||s]{n} + ***
1008  ! *** f4*[(a-1i)s||s]{n+1}) ***
1009 
1010  DO n = 1, nmax - la_max - 1
1011  DO ax = 0, la_max
1012  fx = f2*real(ax, dp)
1013  DO ay = 0, la_max - ax
1014  fy = f2*real(ay, dp)
1015  az = la_max - ax - ay
1016  fz = f2*real(az, dp)
1017 
1018  IF (ax == 0) THEN
1019  v(coset(ax, ay, az), 2, 1, n) = &
1020  rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
1021  rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1)
1022  ELSE
1023  v(coset(ax, ay, az), 2, 1, n) = &
1024  rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
1025  rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1) + &
1026  fx*(v(coset(ax - 1, ay, az), 1, 1, n) + &
1027  f4*v(coset(ax - 1, ay, az), 1, 1, n + 1))
1028  END IF
1029 
1030  IF (ay == 0) THEN
1031  v(coset(ax, ay, az), 3, 1, n) = &
1032  rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
1033  rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1)
1034  ELSE
1035  v(coset(ax, ay, az), 3, 1, n) = &
1036  rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
1037  rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1) + &
1038  fy*(v(coset(ax, ay - 1, az), 1, 1, n) + &
1039  f4*v(coset(ax, ay - 1, az), 1, 1, n + 1))
1040  END IF
1041 
1042  IF (az == 0) THEN
1043  v(coset(ax, ay, az), 4, 1, n) = &
1044  rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
1045  rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1)
1046  ELSE
1047  v(coset(ax, ay, az), 4, 1, n) = &
1048  rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
1049  rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1) + &
1050  fz*(v(coset(ax, ay, az - 1), 1, 1, n) + &
1051  f4*v(coset(ax, ay, az - 1), 1, 1, n + 1))
1052  END IF
1053 
1054  END DO
1055  END DO
1056  END DO
1057 
1058  ! *** Recurrence steps: [ap||s] -> [ab||s] ***
1059 
1060  DO lb = 2, lb_max
1061 
1062  ! *** Horizontal recurrence steps ***
1063 
1064  ! *** [ab||s]{n} = [(a+1i)(b-1i)||s]{n} - ***
1065  ! *** (Bi - Ai)*[a(b-1i)||s]{n} ***
1066 
1067  la_start = max(0, la_min - 1)
1068 
1069  DO la = la_start, la_max - 1
1070  DO n = 1, nmax - la - lb
1071  DO ax = 0, la
1072  DO ay = 0, la - ax
1073  az = la - ax - ay
1074 
1075  ! *** Shift of angular momentum component z from a to b ***
1076 
1077  v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
1078  v(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 1, n) - &
1079  rab(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n)
1080 
1081  ! *** Shift of angular momentum component y from a to b ***
1082 
1083  DO by = 1, lb
1084  bz = lb - by
1085  v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
1086  v(coset(ax, ay + 1, az), coset(0, by - 1, bz), 1, n) - &
1087  rab(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n)
1088  END DO
1089 
1090  ! *** Shift of angular momentum component x from a to b ***
1091 
1092  DO bx = 1, lb
1093  DO by = 0, lb - bx
1094  bz = lb - bx - by
1095  v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
1096  v(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 1, n) - &
1097  rab(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n)
1098  END DO
1099  END DO
1100 
1101  END DO
1102  END DO
1103  END DO
1104  END DO
1105 
1106  ! *** Vertical recurrence step ***
1107 
1108  ! *** [ab||s]{n} = (Pi - Bi)*[a(b-1i)||s]{n} + ***
1109  ! *** (Wi - Pi)*[a(b-1i)||s]{n+1} + ***
1110  ! *** f2*Ni(a)*( [(a-1i)(b-1i)||s]{n} + ***
1111  ! *** f4*[(a-1i)(b-1i)||s]{n+1}) ***
1112  ! *** f2*Ni(b-1i)*( [a(b-2i)||s]{n} + ***
1113  ! *** f4*[a(b-2i)||s]{n+1}) ***
1114 
1115  DO n = 1, nmax - la_max - lb
1116  DO ax = 0, la_max
1117  fx = f2*real(ax, dp)
1118  DO ay = 0, la_max - ax
1119  fy = f2*real(ay, dp)
1120  az = la_max - ax - ay
1121  fz = f2*real(az, dp)
1122 
1123  ! *** Shift of angular momentum component z from a to b ***
1124 
1125  f3 = f2*real(lb - 1, dp)
1126 
1127  IF (az == 0) THEN
1128  v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
1129  rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
1130  rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
1131  f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
1132  f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
1133  ELSE
1134  v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
1135  rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
1136  rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
1137  fz*(v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n) + &
1138  f4*v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n + 1)) + &
1139  f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
1140  f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
1141  END IF
1142 
1143  ! *** Shift of angular momentum component y from a to b ***
1144 
1145  IF (ay == 0) THEN
1146  bz = lb - 1
1147  v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
1148  rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
1149  rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1)
1150  DO by = 2, lb
1151  bz = lb - by
1152  f3 = f2*real(by - 1, dp)
1153  v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
1154  rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
1155  rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
1156  f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
1157  f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
1158  END DO
1159  ELSE
1160  bz = lb - 1
1161  v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
1162  rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
1163  rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1) + &
1164  fy*(v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n) + &
1165  f4*v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n + 1))
1166  DO by = 2, lb
1167  bz = lb - by
1168  f3 = f2*real(by - 1, dp)
1169  v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
1170  rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
1171  rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
1172  fy*(v(coset(ax, ay - 1, az), coset(0, by - 1, bz), 1, n) + &
1173  f4*v(coset(ax, ay - 1, az), &
1174  coset(0, by - 1, bz), 1, n + 1)) + &
1175  f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
1176  f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
1177  END DO
1178  END IF
1179 
1180  ! *** Shift of angular momentum component x from a to b ***
1181 
1182  IF (ax == 0) THEN
1183  DO by = 0, lb - 1
1184  bz = lb - 1 - by
1185  v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
1186  rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
1187  rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1)
1188  END DO
1189  DO bx = 2, lb
1190  f3 = f2*real(bx - 1, dp)
1191  DO by = 0, lb - bx
1192  bz = lb - bx - by
1193  v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
1194  rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
1195  rpw(1)*v(coset(ax, ay, az), &
1196  coset(bx - 1, by, bz), 1, n + 1) + &
1197  f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
1198  f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
1199  END DO
1200  END DO
1201  ELSE
1202  DO by = 0, lb - 1
1203  bz = lb - 1 - by
1204  v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
1205  rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
1206  rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1) + &
1207  fx*(v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n) + &
1208  f4*v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n + 1))
1209  END DO
1210  DO bx = 2, lb
1211  f3 = f2*real(bx - 1, dp)
1212  DO by = 0, lb - bx
1213  bz = lb - bx - by
1214  v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
1215  rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
1216  rpw(1)*v(coset(ax, ay, az), &
1217  coset(bx - 1, by, bz), 1, n + 1) + &
1218  fx*(v(coset(ax - 1, ay, az), &
1219  coset(bx - 1, by, bz), 1, n) + &
1220  f4*v(coset(ax - 1, ay, az), &
1221  coset(bx - 1, by, bz), 1, n + 1)) + &
1222  f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
1223  f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
1224  END DO
1225  END DO
1226  END IF
1227 
1228  END DO
1229  END DO
1230  END DO
1231 
1232  END DO
1233 
1234  END IF
1235 
1236  ELSE
1237 
1238  IF (lb_max > 0) THEN
1239 
1240  ! *** Vertical recurrence steps: [ss||s] -> [sb||s] ***
1241 
1242  rbp(:) = rap(:) - rab(:)
1243 
1244  ! *** [sp||s]{n} = (Pi - Bi)*[ss||s]{n} + ***
1245  ! *** (Wi - Pi)*[ss||s]{n+1} ***
1246 
1247  DO n = 1, nmax - 1
1248  v(1, 2, 1, n) = rbp(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
1249  v(1, 3, 1, n) = rbp(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
1250  v(1, 4, 1, n) = rbp(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
1251  END DO
1252 
1253  ! *** [sb||s]{n} = (Pi - Bi)*[s(b-1i)||s]{n} + ***
1254  ! *** (Wi - Pi)*[s(b-1i)||s]{n+1} + ***
1255  ! *** f2*Ni(b-1i)*( [s(b-2i)||s]{n} + ***
1256  ! *** f4*[s(b-2i)||s]{n+1}) ***
1257 
1258  DO lb = 2, lb_max
1259 
1260  DO n = 1, nmax - lb
1261 
1262  ! *** Increase the angular momentum component z of b ***
1263 
1264  v(1, coset(0, 0, lb), 1, n) = &
1265  rbp(3)*v(1, coset(0, 0, lb - 1), 1, n) + &
1266  rpw(3)*v(1, coset(0, 0, lb - 1), 1, n + 1) + &
1267  f2*real(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), 1, n) + &
1268  f4*v(1, coset(0, 0, lb - 2), 1, n + 1))
1269 
1270  ! *** Increase the angular momentum component y of b ***
1271 
1272  bz = lb - 1
1273  v(1, coset(0, 1, bz), 1, n) = &
1274  rbp(2)*v(1, coset(0, 0, bz), 1, n) + &
1275  rpw(2)*v(1, coset(0, 0, bz), 1, n + 1)
1276 
1277  DO by = 2, lb
1278  bz = lb - by
1279  v(1, coset(0, by, bz), 1, n) = &
1280  rbp(2)*v(1, coset(0, by - 1, bz), 1, n) + &
1281  rpw(2)*v(1, coset(0, by - 1, bz), 1, n + 1) + &
1282  f2*real(by - 1, dp)*(v(1, coset(0, by - 2, bz), 1, n) + &
1283  f4*v(1, coset(0, by - 2, bz), 1, n + 1))
1284  END DO
1285 
1286  ! *** Increase the angular momentum component x of b ***
1287 
1288  DO by = 0, lb - 1
1289  bz = lb - 1 - by
1290  v(1, coset(1, by, bz), 1, n) = &
1291  rbp(1)*v(1, coset(0, by, bz), 1, n) + &
1292  rpw(1)*v(1, coset(0, by, bz), 1, n + 1)
1293  END DO
1294 
1295  DO bx = 2, lb
1296  f3 = f2*real(bx - 1, dp)
1297  DO by = 0, lb - bx
1298  bz = lb - bx - by
1299  v(1, coset(bx, by, bz), 1, n) = &
1300  rbp(1)*v(1, coset(bx - 1, by, bz), 1, n) + &
1301  rpw(1)*v(1, coset(bx - 1, by, bz), 1, n + 1) + &
1302  f3*(v(1, coset(bx - 2, by, bz), 1, n) + &
1303  f4*v(1, coset(bx - 2, by, bz), 1, n + 1))
1304  END DO
1305  END DO
1306 
1307  END DO
1308 
1309  END DO
1310 
1311  END IF
1312 
1313  END IF
1314 
1315  ! *** Recurrence steps: [ab||s] -> [ab||c] ***
1316 
1317  IF (lc_max > 0) THEN
1318 
1319  ! *** Vertical recurrence steps: [ss||s] -> [ss||c] ***
1320 
1321  f5 = -zetw/zetp
1322  f6 = 0.5_dp*zetw
1323  f7 = 0.5_dp*zetq
1324 
1325  rcw(:) = rcp(:) + rpw(:)
1326 
1327  ! *** [ss||p]{n} = (Wi - Ci)*[ss||s]{n+1} (i = x,y,z) ***
1328 
1329  DO n = 1, nmax - 1
1330  v(1, 1, 2, n) = rcw(1)*v(1, 1, 1, n + 1)
1331  v(1, 1, 3, n) = rcw(2)*v(1, 1, 1, n + 1)
1332  v(1, 1, 4, n) = rcw(3)*v(1, 1, 1, n + 1)
1333  END DO
1334 
1335  ! *** [ss||c]{n} = (Wi - Ci)*[ss||c-1i]{n+1} + ***
1336  ! *** f7*Ni(c-1i)*[ss||c-2i]{n} + ***
1337  ! *** f5*[ss||c-2i]{n+1} ***
1338 
1339  DO lc = 2, lc_max
1340 
1341  DO n = 1, nmax - lc
1342 
1343  ! *** Increase the angular momentum component z of c ***
1344 
1345  v(1, 1, coset(0, 0, lc), n) = &
1346  rcw(3)*v(1, 1, coset(0, 0, lc - 1), n + 1) + &
1347  f7*real(lc - 1, dp)*(v(1, 1, coset(0, 0, lc - 2), n) + &
1348  f5*v(1, 1, coset(0, 0, lc - 2), n + 1))
1349 
1350  ! *** Increase the angular momentum component y of c ***
1351 
1352  cz = lc - 1
1353  v(1, 1, coset(0, 1, cz), n) = rcw(2)*v(1, 1, coset(0, 0, cz), n + 1)
1354 
1355  DO cy = 2, lc
1356  cz = lc - cy
1357  v(1, 1, coset(0, cy, cz), n) = &
1358  rcw(2)*v(1, 1, coset(0, cy - 1, cz), n + 1) + &
1359  f7*real(cy - 1, dp)*(v(1, 1, coset(0, cy - 2, cz), n) + &
1360  f5*v(1, 1, coset(0, cy - 2, cz), n + 1))
1361  END DO
1362 
1363  ! *** Increase the angular momentum component x of c ***
1364 
1365  DO cy = 0, lc - 1
1366  cz = lc - 1 - cy
1367  v(1, 1, coset(1, cy, cz), n) = rcw(1)*v(1, 1, coset(0, cy, cz), n + 1)
1368  END DO
1369 
1370  DO cx = 2, lc
1371  DO cy = 0, lc - cx
1372  cz = lc - cx - cy
1373  v(1, 1, coset(cx, cy, cz), n) = &
1374  rcw(1)*v(1, 1, coset(cx - 1, cy, cz), n + 1) + &
1375  f7*real(cx - 1, dp)*(v(1, 1, coset(cx - 2, cy, cz), n) + &
1376  f5*v(1, 1, coset(cx - 2, cy, cz), n + 1))
1377  END DO
1378  END DO
1379 
1380  END DO
1381 
1382  END DO
1383 
1384  ! *** Recurrence steps: [ss||c] -> [ab||c] ***
1385 
1386  DO lc = 1, lc_max
1387 
1388  DO cx = 0, lc
1389  DO cy = 0, lc - cx
1390  cz = lc - cx - cy
1391 
1392  coc = coset(cx, cy, cz)
1393  cocx = coset(max(0, cx - 1), cy, cz)
1394  cocy = coset(cx, max(0, cy - 1), cz)
1395  cocz = coset(cx, cy, max(0, cz - 1))
1396 
1397  fcx = f6*real(cx, dp)
1398  fcy = f6*real(cy, dp)
1399  fcz = f6*real(cz, dp)
1400 
1401  ! *** Recurrence steps: [ss||c] -> [as||c] ***
1402 
1403  IF (la_max > 0) THEN
1404 
1405  ! *** Vertical recurrence steps: [ss||c] -> [as||c] ***
1406 
1407  ! *** [ps||c]{n} = (Pi - Ai)*[ss||c]{n} + ***
1408  ! *** (Wi - Pi)*[ss||c]{n+1} + ***
1409  ! *** f6*Ni(c)*[ss||c-1i]{n+1} (i = x,y,z) ***
1410 
1411  DO n = 1, nmax - 1 - lc
1412  v(2, 1, coc, n) = rap(1)*v(1, 1, coc, n) + &
1413  rpw(1)*v(1, 1, coc, n + 1) + &
1414  fcx*v(1, 1, cocx, n + 1)
1415  v(3, 1, coc, n) = rap(2)*v(1, 1, coc, n) + &
1416  rpw(2)*v(1, 1, coc, n + 1) + &
1417  fcy*v(1, 1, cocy, n + 1)
1418  v(4, 1, coc, n) = rap(3)*v(1, 1, coc, n) + &
1419  rpw(3)*v(1, 1, coc, n + 1) + &
1420  fcz*v(1, 1, cocz, n + 1)
1421  END DO
1422 
1423  ! *** [as||c]{n} = (Pi - Ai)*[(a-1i)s||c]{n} + ***
1424  ! *** (Wi - Pi)*[(a-1i)s||c]{n+1} + ***
1425  ! *** f2*Ni(a-1i)*( [(a-2i)s||c]{n} + ***
1426  ! *** f4*[(a-2i)s||c]{n+1}) + ***
1427  ! *** f6*Ni(c)*[(a-1i)s||c-1i]{n+1} ***
1428 
1429  DO la = 2, la_max
1430 
1431  DO n = 1, nmax - la - lc
1432 
1433  ! *** Increase the angular momentum component z of a ***
1434 
1435  v(coset(0, 0, la), 1, coc, n) = &
1436  rap(3)*v(coset(0, 0, la - 1), 1, coc, n) + &
1437  rpw(3)*v(coset(0, 0, la - 1), 1, coc, n + 1) + &
1438  f2*real(la - 1, dp)*(v(coset(0, 0, la - 2), 1, coc, n) + &
1439  f4*v(coset(0, 0, la - 2), 1, coc, n + 1)) + &
1440  fcz*v(coset(0, 0, la - 1), 1, cocz, n + 1)
1441 
1442  ! *** Increase the angular momentum component y of a ***
1443 
1444  az = la - 1
1445  v(coset(0, 1, az), 1, coc, n) = &
1446  rap(2)*v(coset(0, 0, az), 1, coc, n) + &
1447  rpw(2)*v(coset(0, 0, az), 1, coc, n + 1) + &
1448  fcy*v(coset(0, 0, az), 1, cocy, n + 1)
1449 
1450  DO ay = 2, la
1451  f3 = f2*real(ay - 1, dp)
1452  az = la - ay
1453  v(coset(0, ay, az), 1, coc, n) = &
1454  rap(2)*v(coset(0, ay - 1, az), 1, coc, n) + &
1455  rpw(2)*v(coset(0, ay - 1, az), 1, coc, n + 1) + &
1456  f3*(v(coset(0, ay - 2, az), 1, coc, n) + &
1457  f4*v(coset(0, ay - 2, az), 1, coc, n + 1)) + &
1458  fcy*v(coset(0, ay - 1, az), 1, cocy, n + 1)
1459  END DO
1460 
1461  ! *** Increase the angular momentum component x of a ***
1462 
1463  DO ay = 0, la - 1
1464  az = la - 1 - ay
1465  v(coset(1, ay, az), 1, coc, n) = &
1466  rap(1)*v(coset(0, ay, az), 1, coc, n) + &
1467  rpw(1)*v(coset(0, ay, az), 1, coc, n + 1) + &
1468  fcx*v(coset(0, ay, az), 1, cocx, n + 1)
1469  END DO
1470 
1471  DO ax = 2, la
1472  f3 = f2*real(ax - 1, dp)
1473  DO ay = 0, la - ax
1474  az = la - ax - ay
1475  v(coset(ax, ay, az), 1, coc, n) = &
1476  rap(1)*v(coset(ax - 1, ay, az), 1, coc, n) + &
1477  rpw(1)*v(coset(ax - 1, ay, az), 1, coc, n + 1) + &
1478  f3*(v(coset(ax - 2, ay, az), 1, coc, n) + &
1479  f4*v(coset(ax - 2, ay, az), 1, coc, n + 1)) + &
1480  fcx*v(coset(ax - 1, ay, az), 1, cocx, n + 1)
1481  END DO
1482  END DO
1483 
1484  END DO
1485 
1486  END DO
1487 
1488  ! *** Recurrence steps: [as||c] -> [ab||c] ***
1489 
1490  IF (lb_max > 0) THEN
1491 
1492  ! *** Horizontal recurrence steps ***
1493 
1494  ! *** [ap||c]{n} = [(a+1i)s||c]{n} - (Bi - Ai)*[as||c]{n} ***
1495 
1496  la_start = max(0, la_min - 1)
1497 
1498  DO la = la_start, la_max - 1
1499  DO n = 1, nmax - la - 1 - lc
1500  DO ax = 0, la
1501  DO ay = 0, la - ax
1502  az = la - ax - ay
1503  v(coset(ax, ay, az), 2, coc, n) = &
1504  v(coset(ax + 1, ay, az), 1, coc, n) - &
1505  rab(1)*v(coset(ax, ay, az), 1, coc, n)
1506  v(coset(ax, ay, az), 3, coc, n) = &
1507  v(coset(ax, ay + 1, az), 1, coc, n) - &
1508  rab(2)*v(coset(ax, ay, az), 1, coc, n)
1509  v(coset(ax, ay, az), 4, coc, n) = &
1510  v(coset(ax, ay, az + 1), 1, coc, n) - &
1511  rab(3)*v(coset(ax, ay, az), 1, coc, n)
1512  END DO
1513  END DO
1514  END DO
1515  END DO
1516 
1517  ! *** Vertical recurrence step ***
1518 
1519  ! *** [ap||c]{n} = (Pi - Bi)*[as||c]{n} + ***
1520  ! *** (Wi - Pi)*[as||c]{n+1} + ***
1521  ! *** f2*Ni(a)*( [(a-1i)s||c]{n} + ***
1522  ! *** f4*[(a-1i)s||c]{n+1}) + ***
1523  ! *** f6*Ni(c)*[(as||c-1i]{n+1}) ***
1524 
1525  DO n = 1, nmax - la_max - 1 - lc
1526  DO ax = 0, la_max
1527  fx = f2*real(ax, dp)
1528  DO ay = 0, la_max - ax
1529  fy = f2*real(ay, dp)
1530  az = la_max - ax - ay
1531  fz = f2*real(az, dp)
1532 
1533  IF (ax == 0) THEN
1534  v(coset(ax, ay, az), 2, coc, n) = &
1535  rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
1536  rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1537  fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
1538  ELSE
1539  v(coset(ax, ay, az), 2, coc, n) = &
1540  rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
1541  rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1542  fx*(v(coset(ax - 1, ay, az), 1, coc, n) + &
1543  f4*v(coset(ax - 1, ay, az), 1, coc, n + 1)) + &
1544  fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
1545  END IF
1546 
1547  IF (ay == 0) THEN
1548  v(coset(ax, ay, az), 3, coc, n) = &
1549  rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
1550  rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1551  fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
1552  ELSE
1553  v(coset(ax, ay, az), 3, coc, n) = &
1554  rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
1555  rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1556  fy*(v(coset(ax, ay - 1, az), 1, coc, n) + &
1557  f4*v(coset(ax, ay - 1, az), 1, coc, n + 1)) + &
1558  fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
1559  END IF
1560 
1561  IF (az == 0) THEN
1562  v(coset(ax, ay, az), 4, coc, n) = &
1563  rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
1564  rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1565  fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
1566  ELSE
1567  v(coset(ax, ay, az), 4, coc, n) = &
1568  rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
1569  rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1570  fz*(v(coset(ax, ay, az - 1), 1, coc, n) + &
1571  f4*v(coset(ax, ay, az - 1), 1, coc, n + 1)) + &
1572  fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
1573  END IF
1574 
1575  END DO
1576  END DO
1577  END DO
1578 
1579  ! *** Recurrence steps: [ap||c] -> [ab||c] ***
1580 
1581  DO lb = 2, lb_max
1582 
1583  ! *** Horizontal recurrence steps ***
1584 
1585  ! *** [ab||c]{n} = [(a+1i)(b-1i)||c]{n} - ***
1586  ! *** (Bi - Ai)*[a(b-1i)||c]{n} ***
1587 
1588  la_start = max(0, la_min - 1)
1589 
1590  DO la = la_start, la_max - 1
1591  DO n = 1, nmax - la - lb - lc
1592  DO ax = 0, la
1593  DO ay = 0, la - ax
1594  az = la - ax - ay
1595 
1596  ! *** Shift of angular momentum component z ***
1597 
1598  v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1599  v(coset(ax, ay, az + 1), &
1600  coset(0, 0, lb - 1), coc, n) - &
1601  rab(3)*v(coset(ax, ay, az), &
1602  coset(0, 0, lb - 1), coc, n)
1603 
1604  ! *** Shift of angular momentum component y ***
1605 
1606  DO by = 1, lb
1607  bz = lb - by
1608  v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1609  v(coset(ax, ay + 1, az), &
1610  coset(0, by - 1, bz), coc, n) - &
1611  rab(2)*v(coset(ax, ay, az), &
1612  coset(0, by - 1, bz), coc, n)
1613  END DO
1614 
1615  ! *** Shift of angular momentum component x ***
1616 
1617  DO bx = 1, lb
1618  DO by = 0, lb - bx
1619  bz = lb - bx - by
1620  v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1621  v(coset(ax + 1, ay, az), &
1622  coset(bx - 1, by, bz), coc, n) - &
1623  rab(1)*v(coset(ax, ay, az), &
1624  coset(bx - 1, by, bz), coc, n)
1625  END DO
1626  END DO
1627 
1628  END DO
1629  END DO
1630  END DO
1631  END DO
1632 
1633  ! *** Vertical recurrence step ***
1634 
1635  ! *** [ab||c]{n} = (Pi - Bi)*[a(b-1i)||c]{n} + ***
1636  ! *** (Wi - Pi)*[a(b-1i)||c]{n+1} + ***
1637  ! *** f2*Ni(a)*( [(a-1i)(b-1i)||c]{n} + ***
1638  ! *** f4*[(a-1i)(b-1i)||c]{n+1}) ***
1639  ! *** f2*Ni(b-1i)*( [a(b-2i)||c]{n} + ***
1640  ! *** f4*[a(b-2i)||c]{n+1}) + ***
1641  ! *** f6*Ni(c)*[a(b-1i)||c-1i]{n+1}) ***
1642 
1643  DO n = 1, nmax - la_max - lb - lc
1644  DO ax = 0, la_max
1645  fx = f2*real(ax, dp)
1646  DO ay = 0, la_max - ax
1647  fy = f2*real(ay, dp)
1648  az = la_max - ax - ay
1649  fz = f2*real(az, dp)
1650 
1651  ! *** Shift of angular momentum component z from a to b ***
1652 
1653  f3 = f2*real(lb - 1, dp)
1654 
1655  IF (az == 0) THEN
1656  v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1657  rbp(3)*v(coset(ax, ay, az), &
1658  coset(0, 0, lb - 1), coc, n) + &
1659  rpw(3)*v(coset(ax, ay, az), &
1660  coset(0, 0, lb - 1), coc, n + 1) + &
1661  f3*(v(coset(ax, ay, az), &
1662  coset(0, 0, lb - 2), coc, n) + &
1663  f4*v(coset(ax, ay, az), &
1664  coset(0, 0, lb - 2), coc, n + 1)) + &
1665  fcz*v(coset(ax, ay, az), &
1666  coset(0, 0, lb - 1), cocz, n + 1)
1667  ELSE
1668  v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1669  rbp(3)*v(coset(ax, ay, az), &
1670  coset(0, 0, lb - 1), coc, n) + &
1671  rpw(3)*v(coset(ax, ay, az), &
1672  coset(0, 0, lb - 1), coc, n + 1) + &
1673  fz*(v(coset(ax, ay, az - 1), &
1674  coset(0, 0, lb - 1), coc, n) + &
1675  f4*v(coset(ax, ay, az - 1), &
1676  coset(0, 0, lb - 1), coc, n + 1)) + &
1677  f3*(v(coset(ax, ay, az), &
1678  coset(0, 0, lb - 2), coc, n) + &
1679  f4*v(coset(ax, ay, az), &
1680  coset(0, 0, lb - 2), coc, n + 1)) + &
1681  fcz*v(coset(ax, ay, az), &
1682  coset(0, 0, lb - 1), cocz, n + 1)
1683  END IF
1684 
1685  ! *** Shift of angular momentum component y from a to b ***
1686 
1687  IF (ay == 0) THEN
1688  bz = lb - 1
1689  v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1690  rbp(2)*v(coset(ax, ay, az), &
1691  coset(0, 0, bz), coc, n) + &
1692  rpw(2)*v(coset(ax, ay, az), &
1693  coset(0, 0, bz), coc, n + 1) + &
1694  fcy*v(coset(ax, ay, az), &
1695  coset(0, 0, bz), cocy, n + 1)
1696  DO by = 2, lb
1697  bz = lb - by
1698  f3 = f2*real(by - 1, dp)
1699  v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1700  rbp(2)*v(coset(ax, ay, az), &
1701  coset(0, by - 1, bz), coc, n) + &
1702  rpw(2)*v(coset(ax, ay, az), &
1703  coset(0, by - 1, bz), coc, n + 1) + &
1704  f3*(v(coset(ax, ay, az), &
1705  coset(0, by - 2, bz), coc, n) + &
1706  f4*v(coset(ax, ay, az), &
1707  coset(0, by - 2, bz), coc, n + 1)) + &
1708  fcy*v(coset(ax, ay, az), &
1709  coset(0, by - 1, bz), cocy, n + 1)
1710  END DO
1711  ELSE
1712  bz = lb - 1
1713  v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1714  rbp(2)*v(coset(ax, ay, az), &
1715  coset(0, 0, bz), coc, n) + &
1716  rpw(2)*v(coset(ax, ay, az), &
1717  coset(0, 0, bz), coc, n + 1) + &
1718  fy*(v(coset(ax, ay - 1, az), &
1719  coset(0, 0, bz), coc, n) + &
1720  f4*v(coset(ax, ay - 1, az), &
1721  coset(0, 0, bz), coc, n + 1)) + &
1722  fcy*v(coset(ax, ay, az), &
1723  coset(0, 0, bz), cocy, n + 1)
1724  DO by = 2, lb
1725  bz = lb - by
1726  f3 = f2*real(by - 1, dp)
1727  v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1728  rbp(2)*v(coset(ax, ay, az), &
1729  coset(0, by - 1, bz), coc, n) + &
1730  rpw(2)*v(coset(ax, ay, az), &
1731  coset(0, by - 1, bz), coc, n + 1) + &
1732  fy*(v(coset(ax, ay - 1, az), &
1733  coset(0, by - 1, bz), coc, n) + &
1734  f4*v(coset(ax, ay - 1, az), &
1735  coset(0, by - 1, bz), coc, n + 1)) + &
1736  f3*(v(coset(ax, ay, az), &
1737  coset(0, by - 2, bz), coc, n) + &
1738  f4*v(coset(ax, ay, az), &
1739  coset(0, by - 2, bz), coc, n + 1)) + &
1740  fcy*v(coset(ax, ay, az), &
1741  coset(0, by - 1, bz), cocy, n + 1)
1742  END DO
1743  END IF
1744 
1745  ! *** Shift of angular momentum component x from a to b ***
1746 
1747  IF (ax == 0) THEN
1748  DO by = 0, lb - 1
1749  bz = lb - 1 - by
1750  v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1751  rbp(1)*v(coset(ax, ay, az), &
1752  coset(0, by, bz), coc, n) + &
1753  rpw(1)*v(coset(ax, ay, az), &
1754  coset(0, by, bz), coc, n + 1) + &
1755  fcx*v(coset(ax, ay, az), &
1756  coset(0, by, bz), cocx, n + 1)
1757  END DO
1758  DO bx = 2, lb
1759  f3 = f2*real(bx - 1, dp)
1760  DO by = 0, lb - bx
1761  bz = lb - bx - by
1762  v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1763  rbp(1)*v(coset(ax, ay, az), &
1764  coset(bx - 1, by, bz), coc, n) + &
1765  rpw(1)*v(coset(ax, ay, az), &
1766  coset(bx - 1, by, bz), coc, n + 1) + &
1767  f3*(v(coset(ax, ay, az), &
1768  coset(bx - 2, by, bz), coc, n) + &
1769  f4*v(coset(ax, ay, az), &
1770  coset(bx - 2, by, bz), coc, n + 1)) + &
1771  fcx*v(coset(ax, ay, az), &
1772  coset(bx - 1, by, bz), cocx, n + 1)
1773  END DO
1774  END DO
1775  ELSE
1776  DO by = 0, lb - 1
1777  bz = lb - 1 - by
1778  v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1779  rbp(1)*v(coset(ax, ay, az), &
1780  coset(0, by, bz), coc, n) + &
1781  rpw(1)*v(coset(ax, ay, az), &
1782  coset(0, by, bz), coc, n + 1) + &
1783  fx*(v(coset(ax - 1, ay, az), &
1784  coset(0, by, bz), coc, n) + &
1785  f4*v(coset(ax - 1, ay, az), &
1786  coset(0, by, bz), coc, n + 1)) + &
1787  fcx*v(coset(ax, ay, az), &
1788  coset(0, by, bz), cocx, n + 1)
1789  END DO
1790  DO bx = 2, lb
1791  f3 = f2*real(bx - 1, dp)
1792  DO by = 0, lb - bx
1793  bz = lb - bx - by
1794  v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1795  rbp(1)*v(coset(ax, ay, az), &
1796  coset(bx - 1, by, bz), coc, n) + &
1797  rpw(1)*v(coset(ax, ay, az), &
1798  coset(bx - 1, by, bz), coc, n + 1) + &
1799  fx*(v(coset(ax - 1, ay, az), &
1800  coset(bx - 1, by, bz), coc, n) + &
1801  f4*v(coset(ax - 1, ay, az), &
1802  coset(bx - 1, by, bz), coc, n + 1)) + &
1803  f3*(v(coset(ax, ay, az), &
1804  coset(bx - 2, by, bz), coc, n) + &
1805  f4*v(coset(ax, ay, az), &
1806  coset(bx - 2, by, bz), coc, n + 1)) + &
1807  fcx*v(coset(ax, ay, az), &
1808  coset(bx - 1, by, bz), cocx, n + 1)
1809  END DO
1810  END DO
1811  END IF
1812 
1813  END DO
1814  END DO
1815  END DO
1816 
1817  END DO
1818  END IF
1819 
1820  ELSE
1821 
1822  IF (lb_max > 0) THEN
1823 
1824  ! *** Vertical recurrence steps: [ss||c] -> [sb||c] ***
1825 
1826  ! *** [sp||c]{n} = (Pi - Bi)*[ss||c]{n} + ***
1827  ! *** (Wi - Pi)*[ss||c]{n+1} + ***
1828  ! *** f6*Ni(c)**[ss||c-1i]{n+1} ***
1829 
1830  DO n = 1, nmax - 1 - lc
1831  v(1, 2, coc, n) = rbp(1)*v(1, 1, coc, n) + &
1832  rpw(1)*v(1, 1, coc, n + 1) + &
1833  fcx*v(1, 1, cocx, n + 1)
1834  v(1, 3, coc, n) = rbp(2)*v(1, 1, coc, n) + &
1835  rpw(2)*v(1, 1, coc, n + 1) + &
1836  fcy*v(1, 1, cocy, n + 1)
1837  v(1, 4, coc, n) = rbp(3)*v(1, 1, coc, n) + &
1838  rpw(3)*v(1, 1, coc, n + 1) + &
1839  fcz*v(1, 1, cocz, n + 1)
1840  END DO
1841 
1842  ! *** [sb||c]{n} = (Pi - Bi)*[s(b-1i)||c]{n} + ***
1843  ! *** (Wi - Pi)*[s(b-1i)||c]{n+1} + ***
1844  ! *** f2*Ni(b-1i)*( [s(b-2i)||c]{n} + ***
1845  ! *** f4*[s(b-2i)||c]{n+1}) + ***
1846  ! *** f6*Ni(c)**[s(b-1i)||c-1i]{n+1} ***
1847 
1848  DO lb = 2, lb_max
1849 
1850  DO n = 1, nmax - lb - lc
1851 
1852  ! *** Increase the angular momentum component z of b ***
1853 
1854  v(1, coset(0, 0, lb), coc, n) = &
1855  rbp(3)*v(1, coset(0, 0, lb - 1), coc, n) + &
1856  rpw(3)*v(1, coset(0, 0, lb - 1), coc, n + 1) + &
1857  f2*real(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), coc, n) + &
1858  f4*v(1, coset(0, 0, lb - 2), coc, n + 1)) + &
1859  fcz*v(1, coset(0, 0, lb - 1), cocz, n + 1)
1860 
1861  ! *** Increase the angular momentum component y of b ***
1862 
1863  bz = lb - 1
1864  v(1, coset(0, 1, bz), coc, n) = &
1865  rbp(2)*v(1, coset(0, 0, bz), coc, n) + &
1866  rpw(2)*v(1, coset(0, 0, bz), coc, n + 1) + &
1867  fcy*v(1, coset(0, 0, bz), cocy, n + 1)
1868 
1869  DO by = 2, lb
1870  f3 = f2*real(by - 1, dp)
1871  bz = lb - by
1872  v(1, coset(0, by, bz), coc, n) = &
1873  rbp(2)*v(1, coset(0, by - 1, bz), coc, n) + &
1874  rpw(2)*v(1, coset(0, by - 1, bz), coc, n + 1) + &
1875  f3*(v(1, coset(0, by - 2, bz), coc, n) + &
1876  f4*v(1, coset(0, by - 2, bz), coc, n + 1)) + &
1877  fcy*v(1, coset(0, by - 1, bz), cocy, n + 1)
1878  END DO
1879 
1880  ! *** Increase the angular momentum component x of b ***
1881 
1882  DO by = 0, lb - 1
1883  bz = lb - 1 - by
1884  v(1, coset(1, by, bz), coc, n) = &
1885  rbp(1)*v(1, coset(0, by, bz), coc, n) + &
1886  rpw(1)*v(1, coset(0, by, bz), coc, n + 1) + &
1887  fcx*v(1, coset(0, by, bz), cocx, n + 1)
1888  END DO
1889 
1890  DO bx = 2, lb
1891  f3 = f2*real(bx - 1, dp)
1892  DO by = 0, lb - bx
1893  bz = lb - bx - by
1894  v(1, coset(bx, by, bz), coc, n) = &
1895  rbp(1)*v(1, coset(bx - 1, by, bz), coc, n) + &
1896  rpw(1)*v(1, coset(bx - 1, by, bz), coc, n + 1) + &
1897  f3*(v(1, coset(bx - 2, by, bz), coc, n) + &
1898  f4*v(1, coset(bx - 2, by, bz), coc, n + 1)) + &
1899  fcx*v(1, coset(bx - 1, by, bz), cocx, n + 1)
1900  END DO
1901  END DO
1902 
1903  END DO
1904 
1905  END DO
1906 
1907  END IF
1908 
1909  END IF
1910 
1911  END DO
1912  END DO
1913 
1914  END DO
1915 
1916  END IF
1917 
1918  ! *** Add the contribution of the current pair ***
1919  ! *** of primitive Gaussian-type functions ***
1920 
1921  DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1922  kk = k - ncoset(lc_min - 1)
1923  DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
1924  DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
1925  vabc(na + i, nb + j) = vabc(na + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1926  int_abc(na + i, nb + j, kk) = v(i, j, k, 1)
1927  END DO
1928  END DO
1929  END DO
1930 
1931  IF (PRESENT(maxder)) THEN
1932  DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1933  kk = k - ncoset(lc_min - 1)
1934  DO j = 1, ncoset(lb_max)
1935  DO i = 1, ncoset(la_max)
1936  vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1937  END DO
1938  END DO
1939  END DO
1940  END IF
1941 
1942  nb = nb + ncoset(lb_max)
1943 
1944  END DO
1945 
1946  na = na + ncoset(la_max - maxder_local)
1947  nap = nap + ncoset(la_max)
1948 
1949  END DO
1950 
1951  END SUBROUTINE coulomb3
1952 
1953 END MODULE ai_coulomb
ai_coulomb
Calculation of Coulomb integrals over Cartesian Gaussian-type functions (electron repulsion integrals...
Definition: ai_coulomb.F:41
ai_coulomb::coulomb2
subroutine, public coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lc_max, npgfc, zetc, rpgfc, lc_min, rac, rac2, vac, v, f, maxder, vac_plus)
Calculation of the primitive two-center Coulomb integrals over Cartesian Gaussian-type functions.
Definition: ai_coulomb.F:86
ai_coulomb::coulomb3
subroutine, public coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, v, f, maxder, vabc_plus)
Calculation of the primitive three-center Coulomb integrals over Cartesian Gaussian-type functions (e...
Definition: ai_coulomb.F:798
ai_coulomb::coulomb2_new
subroutine, public coulomb2_new(la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, rac, rac2, vac, v, f, maxder, vac_plus)
Calculation of the primitive two-center Coulomb integrals over Cartesian Gaussian-type functions....
Definition: ai_coulomb.F:441
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition: gamma.F:154
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition: orbital_pointers.F:45