Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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, 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 screening optional primitive-pair screening switch
79 : !> \date 05.12.2000
80 : !> \author Matthias Krack
81 : !> \version 1.0
82 : ! **************************************************************************************************
83 3108 : SUBROUTINE coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lc_max, npgfc, zetc, rpgfc, lc_min, &
84 4662 : rac, rac2, vac, v, f, screening)
85 : INTEGER, INTENT(IN) :: la_max, npgfa
86 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
87 : INTEGER, INTENT(IN) :: la_min, lc_max, npgfc
88 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
89 : INTEGER, INTENT(IN) :: lc_min
90 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
91 : REAL(KIND=dp), INTENT(IN) :: rac2
92 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vac
93 : REAL(KIND=dp), DIMENSION(:, :, :) :: v
94 : REAL(KIND=dp), DIMENSION(0:) :: f
95 : LOGICAL, INTENT(IN), OPTIONAL :: screening
96 :
97 : INTEGER :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
98 : cy, cz, i, ipgf, j, jpgf, la, lc, n, &
99 : na, nc, nmax
100 : LOGICAL :: do_screening
101 : REAL(KIND=dp) :: dac, f0, f1, f2, f3, f4, f5, f6, fcx, &
102 : fcy, fcz, rho, t, zetp, zetq, zetw
103 : REAL(KIND=dp), DIMENSION(3) :: raw, rcw
104 :
105 1554 : do_screening = .TRUE.
106 1554 : IF (PRESENT(screening)) do_screening = screening
107 :
108 5324566 : v = 0.0_dp
109 :
110 1554 : nmax = la_max + lc_max + 1
111 :
112 1554 : dac = SQRT(rac2)
113 :
114 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
115 :
116 1554 : na = 0
117 3558 : DO ipgf = 1, npgfa
118 :
119 2004 : nc = 0
120 :
121 5808 : DO jpgf = 1, npgfc
122 :
123 3804 : IF (do_screening .AND. rpgfa(ipgf) + rpgfc(jpgf) < dac) THEN
124 0 : DO j = nc + ncoset(lc_min - 1) + 1, nc + ncoset(lc_max)
125 0 : DO i = na + ncoset(la_min - 1) + 1, na + ncoset(la_max)
126 0 : vac(i, j) = 0.0_dp
127 : END DO
128 : END DO
129 0 : nc = nc + ncoset(lc_max)
130 0 : CYCLE
131 : END IF
132 :
133 : ! *** Calculate some prefactors ***
134 :
135 3804 : zetp = 1.0_dp/zeta(ipgf)
136 3804 : zetq = 1.0_dp/zetc(jpgf)
137 3804 : zetw = 1.0_dp/(zeta(ipgf) + zetc(jpgf))
138 :
139 3804 : rho = zeta(ipgf)*zetc(jpgf)*zetw
140 :
141 3804 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
142 :
143 : ! *** Calculate the incomplete Gamma function ***
144 :
145 3804 : t = rho*rac2
146 :
147 3804 : CALL fgamma(nmax - 1, t, f)
148 :
149 : ! *** Calculate the basic two-center Coulomb integrals [s||s]{n} ***
150 :
151 13367 : DO n = 1, nmax
152 13367 : v(1, 1, n) = f0*f(n - 1)
153 : END DO
154 :
155 : ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
156 :
157 3804 : IF (lc_max > 0) THEN
158 :
159 1607 : f1 = 0.5_dp*zetq
160 1607 : f2 = -rho*zetq
161 :
162 6428 : rcw(:) = -zeta(ipgf)*zetw*rac(:)
163 :
164 : ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1} (i = x,y,z) ***
165 :
166 6853 : DO n = 1, nmax - 1
167 5246 : v(1, 2, n) = rcw(1)*v(1, 1, n + 1)
168 5246 : v(1, 3, n) = rcw(2)*v(1, 1, n + 1)
169 6853 : v(1, 4, n) = rcw(3)*v(1, 1, n + 1)
170 : END DO
171 :
172 : ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} + ***
173 : ! ** f1*Ni(c-1i)*( [s||c-2i]{n} + ***
174 : ! ** f2*[s||c-2i]{n+1} ***
175 :
176 2899 : DO lc = 2, lc_max
177 :
178 8745 : DO n = 1, nmax - lc
179 :
180 : ! **** Increase the angular momentum component z of c ***
181 :
182 : v(1, coset(0, 0, lc), n) = &
183 : rcw(3)*v(1, coset(0, 0, lc - 1), n + 1) + &
184 : f1*REAL(lc - 1, dp)*(v(1, coset(0, 0, lc - 2), n) + &
185 5846 : f2*v(1, coset(0, 0, lc - 2), n + 1))
186 :
187 : ! *** Increase the angular momentum component y of c ***
188 :
189 5846 : cz = lc - 1
190 5846 : v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n + 1)
191 :
192 17072 : DO cy = 2, lc
193 11226 : cz = lc - cy
194 : v(1, coset(0, cy, cz), n) = &
195 : rcw(2)*v(1, coset(0, cy - 1, cz), n + 1) + &
196 : f1*REAL(cy - 1, dp)*(v(1, coset(0, cy - 2, cz), n) + &
197 17072 : f2*v(1, coset(0, cy - 2, cz), n + 1))
198 : END DO
199 :
200 : ! *** Increase the angular momentum component x of c ***
201 :
202 22918 : DO cy = 0, lc - 1
203 17072 : cz = lc - 1 - cy
204 22918 : v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n + 1)
205 : END DO
206 :
207 18364 : DO cx = 2, lc
208 11226 : f6 = f1*REAL(cx - 1, dp)
209 37198 : DO cy = 0, lc - cx
210 20126 : cz = lc - cx - cy
211 : v(1, coset(cx, cy, cz), n) = &
212 : rcw(1)*v(1, coset(cx - 1, cy, cz), n + 1) + &
213 : f6*(v(1, coset(cx - 2, cy, cz), n) + &
214 31352 : f2*v(1, coset(cx - 2, cy, cz), n + 1))
215 : END DO
216 : END DO
217 :
218 : END DO
219 :
220 : END DO
221 :
222 : END IF
223 :
224 : ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
225 :
226 3804 : IF (la_max > 0) THEN
227 :
228 1592 : f3 = 0.5_dp*zetp
229 1592 : f4 = -rho*zetp
230 1592 : f5 = 0.5_dp*zetw
231 :
232 6368 : raw(:) = zetc(jpgf)*zetw*rac(:)
233 :
234 : ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1} (i = x,y,z) ***
235 :
236 6808 : DO n = 1, nmax - 1
237 5216 : v(2, 1, n) = raw(1)*v(1, 1, n + 1)
238 5216 : v(3, 1, n) = raw(2)*v(1, 1, n + 1)
239 6808 : v(4, 1, n) = raw(3)*v(1, 1, n + 1)
240 : END DO
241 :
242 : ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} + ***
243 : ! *** f3*Ni(a-1i)*( [a-2i||s]{n} + ***
244 : ! *** f4*[a-2i||s]{n+1}) ***
245 :
246 2860 : DO la = 2, la_max
247 :
248 8675 : DO n = 1, nmax - la
249 :
250 : ! *** Increase the angular momentum component z of a ***
251 :
252 : v(coset(0, 0, la), 1, n) = &
253 : raw(3)*v(coset(0, 0, la - 1), 1, n + 1) + &
254 : f3*REAL(la - 1, dp)*(v(coset(0, 0, la - 2), 1, n) + &
255 5815 : f4*v(coset(0, 0, la - 2), 1, n + 1))
256 :
257 : ! *** Increase the angular momentum component y of a ***
258 :
259 5815 : az = la - 1
260 5815 : v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n + 1)
261 :
262 17016 : DO ay = 2, la
263 11201 : az = la - ay
264 : v(coset(0, ay, az), 1, n) = &
265 : raw(2)*v(coset(0, ay - 1, az), 1, n + 1) + &
266 : f3*REAL(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, n) + &
267 17016 : f4*v(coset(0, ay - 2, az), 1, n + 1))
268 : END DO
269 :
270 : ! *** Increase the angular momentum component x of a ***
271 :
272 22831 : DO ay = 0, la - 1
273 17016 : az = la - 1 - ay
274 22831 : v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n + 1)
275 : END DO
276 :
277 18284 : DO ax = 2, la
278 11201 : f6 = f3*REAL(ax - 1, dp)
279 37123 : DO ay = 0, la - ax
280 20107 : az = la - ax - ay
281 : v(coset(ax, ay, az), 1, n) = &
282 : raw(1)*v(coset(ax - 1, ay, az), 1, n + 1) + &
283 : f6*(v(coset(ax - 2, ay, az), 1, n) + &
284 31308 : f4*v(coset(ax - 2, ay, az), 1, n + 1))
285 : END DO
286 : END DO
287 :
288 : END DO
289 :
290 : END DO
291 :
292 3948 : DO lc = 1, lc_max
293 :
294 10600 : DO cx = 0, lc
295 23248 : DO cy = 0, lc - cx
296 14240 : cz = lc - cx - cy
297 :
298 14240 : coc = coset(cx, cy, cz)
299 14240 : cocx = coset(MAX(0, cx - 1), cy, cz)
300 14240 : cocy = coset(cx, MAX(0, cy - 1), cz)
301 14240 : cocz = coset(cx, cy, MAX(0, cz - 1))
302 :
303 14240 : fcx = f5*REAL(cx, dp)
304 14240 : fcy = f5*REAL(cy, dp)
305 14240 : fcz = f5*REAL(cz, dp)
306 :
307 : ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
308 : ! *** f5*Ni(c)*[s||c-1i]{n+1} ***
309 :
310 72773 : DO n = 1, nmax - 1 - lc
311 58533 : v(2, coc, n) = raw(1)*v(1, coc, n + 1) + fcx*v(1, cocx, n + 1)
312 58533 : v(3, coc, n) = raw(2)*v(1, coc, n + 1) + fcy*v(1, cocy, n + 1)
313 72773 : v(4, coc, n) = raw(3)*v(1, coc, n + 1) + fcz*v(1, cocz, n + 1)
314 : END DO
315 :
316 : ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} + ***
317 : ! *** f3*Ni(a-1i)*( [a-2i||c]{n} + ***
318 : ! *** f4*[a-2i||c]{n+1}) + ***
319 : ! *** f5*Ni(c)*[a-1i||c-1i]{n+1} ***
320 :
321 54605 : DO la = 2, la_max
322 :
323 164268 : DO n = 1, nmax - la - lc
324 :
325 : ! *** Increase the angular momentum component z of a ***
326 :
327 : v(coset(0, 0, la), coc, n) = &
328 : raw(3)*v(coset(0, 0, la - 1), coc, n + 1) + &
329 : f3*REAL(la - 1, dp)*(v(coset(0, 0, la - 2), coc, n) + &
330 : f4*v(coset(0, 0, la - 2), coc, n + 1)) + &
331 116315 : fcz*v(coset(0, 0, la - 1), cocz, n + 1)
332 :
333 : ! *** Increase the angular momentum component y of a ***
334 :
335 116315 : az = la - 1
336 : v(coset(0, 1, az), coc, n) = &
337 : raw(2)*v(coset(0, 0, az), coc, n + 1) + &
338 116315 : fcy*v(coset(0, 0, az), cocy, n + 1)
339 :
340 372244 : DO ay = 2, la
341 255929 : az = la - ay
342 : v(coset(0, ay, az), coc, n) = &
343 : raw(2)*v(coset(0, ay - 1, az), coc, n + 1) + &
344 : f3*REAL(ay - 1, dp)*(v(coset(0, ay - 2, az), coc, n) + &
345 : f4*v(coset(0, ay - 2, az), coc, n + 1)) + &
346 372244 : fcy*v(coset(0, ay - 1, az), cocy, n + 1)
347 : END DO
348 :
349 : ! *** Increase the angular momentum component x of a ***
350 :
351 488559 : DO ay = 0, la - 1
352 372244 : az = la - 1 - ay
353 : v(coset(1, ay, az), coc, n) = &
354 : raw(1)*v(coset(0, ay, az), coc, n + 1) + &
355 488559 : fcx*v(coset(0, ay, az), cocx, n + 1)
356 : END DO
357 :
358 405957 : DO ax = 2, la
359 255929 : f6 = f3*REAL(ax - 1, dp)
360 867307 : DO ay = 0, la - ax
361 495063 : az = la - ax - ay
362 : v(coset(ax, ay, az), coc, n) = &
363 : raw(1)*v(coset(ax - 1, ay, az), coc, n + 1) + &
364 : f6*(v(coset(ax - 2, ay, az), coc, n) + &
365 : f4*v(coset(ax - 2, ay, az), coc, n + 1)) + &
366 750992 : fcx*v(coset(ax - 1, ay, az), cocx, n + 1)
367 : END DO
368 : END DO
369 :
370 : END DO
371 :
372 : END DO
373 :
374 : END DO
375 : END DO
376 :
377 : END DO
378 :
379 : END IF
380 :
381 15568 : DO j = ncoset(lc_min - 1) + 1, ncoset(lc_max)
382 103964 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max)
383 100160 : vac(na + i, nc + j) = v(i, j, 1)
384 : END DO
385 : END DO
386 :
387 5808 : nc = nc + ncoset(lc_max)
388 :
389 : END DO
390 :
391 3558 : na = na + ncoset(la_max)
392 :
393 : END DO
394 :
395 1554 : END SUBROUTINE coulomb2
396 : ! **************************************************************************************************
397 : !> \brief Calculation of the primitive three-center Coulomb integrals over
398 : !> Cartesian Gaussian-type functions (electron repulsion integrals,
399 : !> ERIs).
400 : !> \param la_max ...
401 : !> \param npgfa ...
402 : !> \param zeta ...
403 : !> \param rpgfa ...
404 : !> \param la_min ...
405 : !> \param lb_max ...
406 : !> \param npgfb ...
407 : !> \param zetb ...
408 : !> \param rpgfb ...
409 : !> \param lb_min ...
410 : !> \param lc_max ...
411 : !> \param zetc ...
412 : !> \param rpgfc ...
413 : !> \param lc_min ...
414 : !> \param gccc ...
415 : !> \param rab ...
416 : !> \param rab2 ...
417 : !> \param rac ...
418 : !> \param rac2 ...
419 : !> \param rbc2 ...
420 : !> \param vabc ...
421 : !> \param int_abc ...
422 : !> \param v ...
423 : !> \param f ...
424 : !> \param maxder ...
425 : !> \param vabc_plus ...
426 : !> \date 06.11.2000
427 : !> \author Matthias Krack
428 : !> \version 1.0
429 : ! **************************************************************************************************
430 1620 : SUBROUTINE coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
431 3240 : lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, &
432 4860 : v, f, maxder, vabc_plus)
433 : INTEGER, INTENT(IN) :: la_max, npgfa
434 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
435 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
436 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
437 : INTEGER, INTENT(IN) :: lb_min, lc_max
438 : REAL(KIND=dp), INTENT(IN) :: zetc, rpgfc
439 : INTEGER, INTENT(IN) :: lc_min
440 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: gccc
441 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
442 : REAL(KIND=dp), INTENT(IN) :: rab2
443 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
444 : REAL(KIND=dp), INTENT(IN) :: rac2, rbc2
445 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vabc
446 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: int_abc
447 : REAL(KIND=dp), DIMENSION(:, :, :, :) :: v
448 : REAL(KIND=dp), DIMENSION(0:) :: f
449 : INTEGER, INTENT(IN), OPTIONAL :: maxder
450 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: vabc_plus
451 :
452 : INTEGER :: ax, ay, az, bx, by, bz, coc, cocx, cocy, &
453 : cocz, cx, cy, cz, i, ipgf, j, jpgf, k, &
454 : kk, la, la_start, lb, lc, &
455 : maxder_local, n, na, nap, nb, nmax
456 : REAL(KIND=dp) :: dab, dac, dbc, f0, f1, f2, f3, f4, f5, &
457 : f6, f7, fcx, fcy, fcz, fx, fy, fz, t, &
458 : zetp, zetq, zetw
459 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp, rcp, rcw, rpw
460 :
461 1374540 : v = 0.0_dp
462 :
463 1620 : maxder_local = 0
464 1620 : IF (PRESENT(maxder)) THEN
465 0 : maxder_local = maxder
466 : END IF
467 :
468 1620 : nmax = la_max + lb_max + lc_max + 1
469 :
470 : ! *** Calculate the distances of the centers a, b and c ***
471 :
472 1620 : dab = SQRT(rab2)
473 1620 : dac = SQRT(rac2)
474 1620 : dbc = SQRT(rbc2)
475 :
476 : ! *** Initialize integrals array
477 168570 : int_abc = 0.0_dp
478 :
479 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
480 :
481 1620 : na = 0
482 1620 : nap = 0
483 :
484 4320 : DO ipgf = 1, npgfa
485 :
486 : ! *** Screening ***
487 2700 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
488 0 : na = na + ncoset(la_max - maxder_local)
489 0 : nap = nap + ncoset(la_max)
490 0 : CYCLE
491 : END IF
492 :
493 2700 : nb = 0
494 :
495 7200 : DO jpgf = 1, npgfb
496 :
497 : ! *** Screening ***
498 : IF ( &
499 4500 : (rpgfb(jpgf) + rpgfc < dbc) .OR. &
500 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
501 0 : nb = nb + ncoset(lb_max)
502 0 : CYCLE
503 : END IF
504 :
505 : ! *** Calculate some prefactors ***
506 :
507 4500 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
508 4500 : zetq = 1.0_dp/zetc
509 4500 : zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
510 :
511 4500 : f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
512 4500 : f1 = zetb(jpgf)*zetp
513 4500 : f2 = 0.5_dp*zetp
514 4500 : f4 = -zetc*zetw
515 :
516 4500 : f0 = f0*EXP(-zeta(ipgf)*f1*rab2)
517 :
518 18000 : rap(:) = f1*rab(:)
519 18000 : rcp(:) = rap(:) - rac(:)
520 18000 : rpw(:) = f4*rcp(:)
521 :
522 : ! *** Calculate the incomplete Gamma function ***
523 :
524 4500 : t = -f4*(rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3))/zetp
525 :
526 4500 : CALL fgamma(nmax - 1, t, f)
527 :
528 : ! *** Calculate the basic three-center Coulomb integrals [ss||s]{n} ***
529 :
530 18300 : DO n = 1, nmax
531 18300 : v(1, 1, 1, n) = f0*f(n - 1)
532 : END DO
533 :
534 : ! *** Recurrence steps: [ss||s] -> [as||s] ***
535 :
536 4500 : IF (la_max > 0) THEN
537 :
538 : ! *** Vertical recurrence steps: [ss||s] -> [as||s] ***
539 :
540 : ! *** [ps||s]{n} = (Pi - Ai)*[ss||s]{n} + ***
541 : ! *** (Wi - Pi)*[ss||s]{n+1} (i = x,y,z) ***
542 :
543 7920 : DO n = 1, nmax - 1
544 5820 : v(2, 1, 1, n) = rap(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
545 5820 : v(3, 1, 1, n) = rap(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
546 7920 : v(4, 1, 1, n) = rap(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
547 : END DO
548 :
549 : ! *** [as||s]{n} = (Pi - Ai)*[(a-1i)s||s]{n} + ***
550 : ! *** (Wi - Pi)*[(a-1i)s||s]{n+1} + ***
551 : ! *** f2*Ni(a-1i)*( [(a-2i)s||s]{n} + ***
552 : ! *** f4*[(a-2i)s||s]{n+1}) ***
553 :
554 2400 : DO la = 2, la_max
555 :
556 3210 : DO n = 1, nmax - la
557 :
558 : ! *** Increase the angular momentum component z of a ***
559 :
560 : v(coset(0, 0, la), 1, 1, n) = &
561 : rap(3)*v(coset(0, 0, la - 1), 1, 1, n) + &
562 : rpw(3)*v(coset(0, 0, la - 1), 1, 1, n + 1) + &
563 : f2*REAL(la - 1, dp)*(v(coset(0, 0, la - 2), 1, 1, n) + &
564 810 : f4*v(coset(0, 0, la - 2), 1, 1, n + 1))
565 :
566 : ! *** Increase the angular momentum component y of a ***
567 :
568 810 : az = la - 1
569 : v(coset(0, 1, az), 1, 1, n) = &
570 : rap(2)*v(coset(0, 0, az), 1, 1, n) + &
571 810 : rpw(2)*v(coset(0, 0, az), 1, 1, n + 1)
572 :
573 1620 : DO ay = 2, la
574 810 : az = la - ay
575 : v(coset(0, ay, az), 1, 1, n) = &
576 : rap(2)*v(coset(0, ay - 1, az), 1, 1, n) + &
577 : rpw(2)*v(coset(0, ay - 1, az), 1, 1, n + 1) + &
578 : f2*REAL(ay - 1, dp)*(v(coset(0, ay - 2, az), 1, 1, n) + &
579 1620 : f4*v(coset(0, ay - 2, az), 1, 1, n + 1))
580 : END DO
581 :
582 : ! *** Increase the angular momentum component x of a ***
583 :
584 2430 : DO ay = 0, la - 1
585 1620 : az = la - 1 - ay
586 : v(coset(1, ay, az), 1, 1, n) = &
587 : rap(1)*v(coset(0, ay, az), 1, 1, n) + &
588 2430 : rpw(1)*v(coset(0, ay, az), 1, 1, n + 1)
589 : END DO
590 :
591 1920 : DO ax = 2, la
592 810 : f3 = f2*REAL(ax - 1, dp)
593 2430 : DO ay = 0, la - ax
594 810 : az = la - ax - ay
595 : v(coset(ax, ay, az), 1, 1, n) = &
596 : rap(1)*v(coset(ax - 1, ay, az), 1, 1, n) + &
597 : rpw(1)*v(coset(ax - 1, ay, az), 1, 1, n + 1) + &
598 : f3*(v(coset(ax - 2, ay, az), 1, 1, n) + &
599 1620 : f4*v(coset(ax - 2, ay, az), 1, 1, n + 1))
600 : END DO
601 : END DO
602 :
603 : END DO
604 :
605 : END DO
606 :
607 : ! *** Recurrence steps: [as||s] -> [ab||s] ***
608 :
609 2100 : IF (lb_max > 0) THEN
610 :
611 : ! *** Horizontal recurrence steps ***
612 :
613 4560 : rbp(:) = rap(:) - rab(:)
614 :
615 : ! *** [ap||s]{n} = [(a+1i)s||s]{n} - (Bi - Ai)*[as||s]{n} ***
616 :
617 1140 : la_start = MAX(0, la_min - 1)
618 :
619 2280 : DO la = la_start, la_max - 1
620 5880 : DO n = 1, nmax - la - 1
621 8910 : DO ax = 0, la
622 12510 : DO ay = 0, la - ax
623 4740 : az = la - ax - ay
624 : v(coset(ax, ay, az), 2, 1, n) = &
625 : v(coset(ax + 1, ay, az), 1, 1, n) - &
626 4740 : rab(1)*v(coset(ax, ay, az), 1, 1, n)
627 : v(coset(ax, ay, az), 3, 1, n) = &
628 : v(coset(ax, ay + 1, az), 1, 1, n) - &
629 4740 : rab(2)*v(coset(ax, ay, az), 1, 1, n)
630 : v(coset(ax, ay, az), 4, 1, n) = &
631 : v(coset(ax, ay, az + 1), 1, 1, n) - &
632 8910 : rab(3)*v(coset(ax, ay, az), 1, 1, n)
633 : END DO
634 : END DO
635 : END DO
636 : END DO
637 :
638 : ! *** Vertical recurrence step ***
639 :
640 : ! *** [ap||s]{n} = (Pi - Bi)*[as||s]{n} + ***
641 : ! *** (Wi - Pi)*[as||s]{n+1} + ***
642 : ! *** f2*Ni(a)*( [(a-1i)s||s]{n} + ***
643 : ! *** f4*[(a-1i)s||s]{n+1}) ***
644 :
645 3600 : DO n = 1, nmax - la_max - 1
646 8910 : DO ax = 0, la_max
647 5310 : fx = f2*REAL(ax, dp)
648 16320 : DO ay = 0, la_max - ax
649 8550 : fy = f2*REAL(ay, dp)
650 8550 : az = la_max - ax - ay
651 8550 : fz = f2*REAL(az, dp)
652 :
653 8550 : IF (ax == 0) THEN
654 : v(coset(ax, ay, az), 2, 1, n) = &
655 : rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
656 5310 : rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1)
657 : ELSE
658 : v(coset(ax, ay, az), 2, 1, n) = &
659 : rbp(1)*v(coset(ax, ay, az), 1, 1, n) + &
660 : rpw(1)*v(coset(ax, ay, az), 1, 1, n + 1) + &
661 : fx*(v(coset(ax - 1, ay, az), 1, 1, n) + &
662 3240 : f4*v(coset(ax - 1, ay, az), 1, 1, n + 1))
663 : END IF
664 :
665 8550 : IF (ay == 0) THEN
666 : v(coset(ax, ay, az), 3, 1, n) = &
667 : rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
668 5310 : rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1)
669 : ELSE
670 : v(coset(ax, ay, az), 3, 1, n) = &
671 : rbp(2)*v(coset(ax, ay, az), 1, 1, n) + &
672 : rpw(2)*v(coset(ax, ay, az), 1, 1, n + 1) + &
673 : fy*(v(coset(ax, ay - 1, az), 1, 1, n) + &
674 3240 : f4*v(coset(ax, ay - 1, az), 1, 1, n + 1))
675 : END IF
676 :
677 13860 : IF (az == 0) THEN
678 : v(coset(ax, ay, az), 4, 1, n) = &
679 : rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
680 5310 : rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1)
681 : ELSE
682 : v(coset(ax, ay, az), 4, 1, n) = &
683 : rbp(3)*v(coset(ax, ay, az), 1, 1, n) + &
684 : rpw(3)*v(coset(ax, ay, az), 1, 1, n + 1) + &
685 : fz*(v(coset(ax, ay, az - 1), 1, 1, n) + &
686 3240 : f4*v(coset(ax, ay, az - 1), 1, 1, n + 1))
687 : END IF
688 :
689 : END DO
690 : END DO
691 : END DO
692 :
693 : ! *** Recurrence steps: [ap||s] -> [ab||s] ***
694 :
695 1320 : DO lb = 2, lb_max
696 :
697 : ! *** Horizontal recurrence steps ***
698 :
699 : ! *** [ab||s]{n} = [(a+1i)(b-1i)||s]{n} - ***
700 : ! *** (Bi - Ai)*[a(b-1i)||s]{n} ***
701 :
702 360 : la_start = MAX(0, la_min - 1)
703 :
704 360 : DO la = la_start, la_max - 1
705 900 : DO n = 1, nmax - la - lb
706 1350 : DO ax = 0, la
707 1890 : DO ay = 0, la - ax
708 720 : az = la - ax - ay
709 :
710 : ! *** Shift of angular momentum component z from a to b ***
711 :
712 : v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
713 : v(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 1, n) - &
714 720 : rab(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n)
715 :
716 : ! *** Shift of angular momentum component y from a to b ***
717 :
718 2160 : DO by = 1, lb
719 1440 : bz = lb - by
720 : v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
721 : v(coset(ax, ay + 1, az), coset(0, by - 1, bz), 1, n) - &
722 2160 : rab(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n)
723 : END DO
724 :
725 : ! *** Shift of angular momentum component x from a to b ***
726 :
727 2790 : DO bx = 1, lb
728 4320 : DO by = 0, lb - bx
729 2160 : bz = lb - bx - by
730 : v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
731 : v(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 1, n) - &
732 3600 : rab(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n)
733 : END DO
734 : END DO
735 :
736 : END DO
737 : END DO
738 : END DO
739 : END DO
740 :
741 : ! *** Vertical recurrence step ***
742 :
743 : ! *** [ab||s]{n} = (Pi - Bi)*[a(b-1i)||s]{n} + ***
744 : ! *** (Wi - Pi)*[a(b-1i)||s]{n+1} + ***
745 : ! *** f2*Ni(a)*( [(a-1i)(b-1i)||s]{n} + ***
746 : ! *** f4*[(a-1i)(b-1i)||s]{n+1}) ***
747 : ! *** f2*Ni(b-1i)*( [a(b-2i)||s]{n} + ***
748 : ! *** f4*[a(b-2i)||s]{n+1}) ***
749 :
750 1680 : DO n = 1, nmax - la_max - lb
751 1320 : DO ax = 0, la_max
752 780 : fx = f2*REAL(ax, dp)
753 2400 : DO ay = 0, la_max - ax
754 1260 : fy = f2*REAL(ay, dp)
755 1260 : az = la_max - ax - ay
756 1260 : fz = f2*REAL(az, dp)
757 :
758 : ! *** Shift of angular momentum component z from a to b ***
759 :
760 1260 : f3 = f2*REAL(lb - 1, dp)
761 :
762 1260 : IF (az == 0) THEN
763 : v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
764 : rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
765 : rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
766 : f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
767 780 : f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
768 : ELSE
769 : v(coset(ax, ay, az), coset(0, 0, lb), 1, n) = &
770 : rbp(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n) + &
771 : rpw(3)*v(coset(ax, ay, az), coset(0, 0, lb - 1), 1, n + 1) + &
772 : fz*(v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n) + &
773 : f4*v(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1, n + 1)) + &
774 : f3*(v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n) + &
775 480 : f4*v(coset(ax, ay, az), coset(0, 0, lb - 2), 1, n + 1))
776 : END IF
777 :
778 : ! *** Shift of angular momentum component y from a to b ***
779 :
780 1260 : IF (ay == 0) THEN
781 780 : bz = lb - 1
782 : v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
783 : rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
784 780 : rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1)
785 1560 : DO by = 2, lb
786 780 : bz = lb - by
787 780 : f3 = f2*REAL(by - 1, dp)
788 : v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
789 : rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
790 : rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
791 : f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
792 1560 : f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
793 : END DO
794 : ELSE
795 480 : bz = lb - 1
796 : v(coset(ax, ay, az), coset(0, 1, bz), 1, n) = &
797 : rbp(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n) + &
798 : rpw(2)*v(coset(ax, ay, az), coset(0, 0, bz), 1, n + 1) + &
799 : fy*(v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n) + &
800 480 : f4*v(coset(ax, ay - 1, az), coset(0, 0, bz), 1, n + 1))
801 960 : DO by = 2, lb
802 480 : bz = lb - by
803 480 : f3 = f2*REAL(by - 1, dp)
804 : v(coset(ax, ay, az), coset(0, by, bz), 1, n) = &
805 : rbp(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n) + &
806 : rpw(2)*v(coset(ax, ay, az), coset(0, by - 1, bz), 1, n + 1) + &
807 : fy*(v(coset(ax, ay - 1, az), coset(0, by - 1, bz), 1, n) + &
808 : f4*v(coset(ax, ay - 1, az), &
809 : coset(0, by - 1, bz), 1, n + 1)) + &
810 : f3*(v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n) + &
811 960 : f4*v(coset(ax, ay, az), coset(0, by - 2, bz), 1, n + 1))
812 : END DO
813 : END IF
814 :
815 : ! *** Shift of angular momentum component x from a to b ***
816 :
817 2040 : IF (ax == 0) THEN
818 2340 : DO by = 0, lb - 1
819 1560 : bz = lb - 1 - by
820 : v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
821 : rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
822 2340 : rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1)
823 : END DO
824 1560 : DO bx = 2, lb
825 780 : f3 = f2*REAL(bx - 1, dp)
826 2340 : DO by = 0, lb - bx
827 780 : bz = lb - bx - by
828 : v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
829 : rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
830 : rpw(1)*v(coset(ax, ay, az), &
831 : coset(bx - 1, by, bz), 1, n + 1) + &
832 : f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
833 1560 : f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
834 : END DO
835 : END DO
836 : ELSE
837 1440 : DO by = 0, lb - 1
838 960 : bz = lb - 1 - by
839 : v(coset(ax, ay, az), coset(1, by, bz), 1, n) = &
840 : rbp(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n) + &
841 : rpw(1)*v(coset(ax, ay, az), coset(0, by, bz), 1, n + 1) + &
842 : fx*(v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n) + &
843 1440 : f4*v(coset(ax - 1, ay, az), coset(0, by, bz), 1, n + 1))
844 : END DO
845 960 : DO bx = 2, lb
846 480 : f3 = f2*REAL(bx - 1, dp)
847 1440 : DO by = 0, lb - bx
848 480 : bz = lb - bx - by
849 : v(coset(ax, ay, az), coset(bx, by, bz), 1, n) = &
850 : rbp(1)*v(coset(ax, ay, az), coset(bx - 1, by, bz), 1, n) + &
851 : rpw(1)*v(coset(ax, ay, az), &
852 : coset(bx - 1, by, bz), 1, n + 1) + &
853 : fx*(v(coset(ax - 1, ay, az), &
854 : coset(bx - 1, by, bz), 1, n) + &
855 : f4*v(coset(ax - 1, ay, az), &
856 : coset(bx - 1, by, bz), 1, n + 1)) + &
857 : f3*(v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n) + &
858 960 : f4*v(coset(ax, ay, az), coset(bx - 2, by, bz), 1, n + 1))
859 : END DO
860 : END DO
861 : END IF
862 :
863 : END DO
864 : END DO
865 : END DO
866 :
867 : END DO
868 :
869 : END IF
870 :
871 : ELSE
872 :
873 2400 : IF (lb_max > 0) THEN
874 :
875 : ! *** Vertical recurrence steps: [ss||s] -> [sb||s] ***
876 :
877 3840 : rbp(:) = rap(:) - rab(:)
878 :
879 : ! *** [sp||s]{n} = (Pi - Bi)*[ss||s]{n} + ***
880 : ! *** (Wi - Pi)*[ss||s]{n+1} ***
881 :
882 3000 : DO n = 1, nmax - 1
883 2040 : v(1, 2, 1, n) = rbp(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
884 2040 : v(1, 3, 1, n) = rbp(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
885 3000 : v(1, 4, 1, n) = rbp(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
886 : END DO
887 :
888 : ! *** [sb||s]{n} = (Pi - Bi)*[s(b-1i)||s]{n} + ***
889 : ! *** (Wi - Pi)*[s(b-1i)||s]{n+1} + ***
890 : ! *** f2*Ni(b-1i)*( [s(b-2i)||s]{n} + ***
891 : ! *** f4*[s(b-2i)||s]{n+1}) ***
892 :
893 1080 : DO lb = 2, lb_max
894 :
895 1320 : DO n = 1, nmax - lb
896 :
897 : ! *** Increase the angular momentum component z of b ***
898 :
899 : v(1, coset(0, 0, lb), 1, n) = &
900 : rbp(3)*v(1, coset(0, 0, lb - 1), 1, n) + &
901 : rpw(3)*v(1, coset(0, 0, lb - 1), 1, n + 1) + &
902 : f2*REAL(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), 1, n) + &
903 240 : f4*v(1, coset(0, 0, lb - 2), 1, n + 1))
904 :
905 : ! *** Increase the angular momentum component y of b ***
906 :
907 240 : bz = lb - 1
908 : v(1, coset(0, 1, bz), 1, n) = &
909 : rbp(2)*v(1, coset(0, 0, bz), 1, n) + &
910 240 : rpw(2)*v(1, coset(0, 0, bz), 1, n + 1)
911 :
912 480 : DO by = 2, lb
913 240 : bz = lb - by
914 : v(1, coset(0, by, bz), 1, n) = &
915 : rbp(2)*v(1, coset(0, by - 1, bz), 1, n) + &
916 : rpw(2)*v(1, coset(0, by - 1, bz), 1, n + 1) + &
917 : f2*REAL(by - 1, dp)*(v(1, coset(0, by - 2, bz), 1, n) + &
918 480 : f4*v(1, coset(0, by - 2, bz), 1, n + 1))
919 : END DO
920 :
921 : ! *** Increase the angular momentum component x of b ***
922 :
923 720 : DO by = 0, lb - 1
924 480 : bz = lb - 1 - by
925 : v(1, coset(1, by, bz), 1, n) = &
926 : rbp(1)*v(1, coset(0, by, bz), 1, n) + &
927 720 : rpw(1)*v(1, coset(0, by, bz), 1, n + 1)
928 : END DO
929 :
930 600 : DO bx = 2, lb
931 240 : f3 = f2*REAL(bx - 1, dp)
932 720 : DO by = 0, lb - bx
933 240 : bz = lb - bx - by
934 : v(1, coset(bx, by, bz), 1, n) = &
935 : rbp(1)*v(1, coset(bx - 1, by, bz), 1, n) + &
936 : rpw(1)*v(1, coset(bx - 1, by, bz), 1, n + 1) + &
937 : f3*(v(1, coset(bx - 2, by, bz), 1, n) + &
938 480 : f4*v(1, coset(bx - 2, by, bz), 1, n + 1))
939 : END DO
940 : END DO
941 :
942 : END DO
943 :
944 : END DO
945 :
946 : END IF
947 :
948 : END IF
949 :
950 : ! *** Recurrence steps: [ab||s] -> [ab||c] ***
951 :
952 4500 : IF (lc_max > 0) THEN
953 :
954 : ! *** Vertical recurrence steps: [ss||s] -> [ss||c] ***
955 :
956 2700 : f5 = -zetw/zetp
957 2700 : f6 = 0.5_dp*zetw
958 2700 : f7 = 0.5_dp*zetq
959 :
960 10800 : rcw(:) = rcp(:) + rpw(:)
961 :
962 : ! *** [ss||p]{n} = (Wi - Ci)*[ss||s]{n+1} (i = x,y,z) ***
963 :
964 10080 : DO n = 1, nmax - 1
965 7380 : v(1, 1, 2, n) = rcw(1)*v(1, 1, 1, n + 1)
966 7380 : v(1, 1, 3, n) = rcw(2)*v(1, 1, 1, n + 1)
967 10080 : v(1, 1, 4, n) = rcw(3)*v(1, 1, 1, n + 1)
968 : END DO
969 :
970 : ! *** [ss||c]{n} = (Wi - Ci)*[ss||c-1i]{n+1} + ***
971 : ! *** f7*Ni(c-1i)*[ss||c-2i]{n} + ***
972 : ! *** f5*[ss||c-2i]{n+1} ***
973 :
974 4500 : DO lc = 2, lc_max
975 :
976 8670 : DO n = 1, nmax - lc
977 :
978 : ! *** Increase the angular momentum component z of c ***
979 :
980 : v(1, 1, coset(0, 0, lc), n) = &
981 : rcw(3)*v(1, 1, coset(0, 0, lc - 1), n + 1) + &
982 : f7*REAL(lc - 1, dp)*(v(1, 1, coset(0, 0, lc - 2), n) + &
983 4170 : f5*v(1, 1, coset(0, 0, lc - 2), n + 1))
984 :
985 : ! *** Increase the angular momentum component y of c ***
986 :
987 4170 : cz = lc - 1
988 4170 : v(1, 1, coset(0, 1, cz), n) = rcw(2)*v(1, 1, coset(0, 0, cz), n + 1)
989 :
990 9270 : DO cy = 2, lc
991 5100 : cz = lc - cy
992 : v(1, 1, coset(0, cy, cz), n) = &
993 : rcw(2)*v(1, 1, coset(0, cy - 1, cz), n + 1) + &
994 : f7*REAL(cy - 1, dp)*(v(1, 1, coset(0, cy - 2, cz), n) + &
995 9270 : f5*v(1, 1, coset(0, cy - 2, cz), n + 1))
996 : END DO
997 :
998 : ! *** Increase the angular momentum component x of c ***
999 :
1000 13440 : DO cy = 0, lc - 1
1001 9270 : cz = lc - 1 - cy
1002 13440 : v(1, 1, coset(1, cy, cz), n) = rcw(1)*v(1, 1, coset(0, cy, cz), n + 1)
1003 : END DO
1004 :
1005 11070 : DO cx = 2, lc
1006 15300 : DO cy = 0, lc - cx
1007 6030 : cz = lc - cx - cy
1008 : v(1, 1, coset(cx, cy, cz), n) = &
1009 : rcw(1)*v(1, 1, coset(cx - 1, cy, cz), n + 1) + &
1010 : f7*REAL(cx - 1, dp)*(v(1, 1, coset(cx - 2, cy, cz), n) + &
1011 11130 : f5*v(1, 1, coset(cx - 2, cy, cz), n + 1))
1012 : END DO
1013 : END DO
1014 :
1015 : END DO
1016 :
1017 : END DO
1018 :
1019 : ! *** Recurrence steps: [ss||c] -> [ab||c] ***
1020 :
1021 7200 : DO lc = 1, lc_max
1022 :
1023 18450 : DO cx = 0, lc
1024 36450 : DO cy = 0, lc - cx
1025 20700 : cz = lc - cx - cy
1026 :
1027 20700 : coc = coset(cx, cy, cz)
1028 20700 : cocx = coset(MAX(0, cx - 1), cy, cz)
1029 20700 : cocy = coset(cx, MAX(0, cy - 1), cz)
1030 20700 : cocz = coset(cx, cy, MAX(0, cz - 1))
1031 :
1032 20700 : fcx = f6*REAL(cx, dp)
1033 20700 : fcy = f6*REAL(cy, dp)
1034 20700 : fcz = f6*REAL(cz, dp)
1035 :
1036 : ! *** Recurrence steps: [ss||c] -> [as||c] ***
1037 :
1038 31950 : IF (la_max > 0) THEN
1039 :
1040 : ! *** Vertical recurrence steps: [ss||c] -> [as||c] ***
1041 :
1042 : ! *** [ps||c]{n} = (Pi - Ai)*[ss||c]{n} + ***
1043 : ! *** (Wi - Pi)*[ss||c]{n+1} + ***
1044 : ! *** f6*Ni(c)*[ss||c-1i]{n+1} (i = x,y,z) ***
1045 :
1046 30552 : DO n = 1, nmax - 1 - lc
1047 : v(2, 1, coc, n) = rap(1)*v(1, 1, coc, n) + &
1048 : rpw(1)*v(1, 1, coc, n + 1) + &
1049 20892 : fcx*v(1, 1, cocx, n + 1)
1050 : v(3, 1, coc, n) = rap(2)*v(1, 1, coc, n) + &
1051 : rpw(2)*v(1, 1, coc, n + 1) + &
1052 20892 : fcy*v(1, 1, cocy, n + 1)
1053 : v(4, 1, coc, n) = rap(3)*v(1, 1, coc, n) + &
1054 : rpw(3)*v(1, 1, coc, n + 1) + &
1055 30552 : fcz*v(1, 1, cocz, n + 1)
1056 : END DO
1057 :
1058 : ! *** [as||c]{n} = (Pi - Ai)*[(a-1i)s||c]{n} + ***
1059 : ! *** (Wi - Pi)*[(a-1i)s||c]{n+1} + ***
1060 : ! *** f2*Ni(a-1i)*( [(a-2i)s||c]{n} + ***
1061 : ! *** f4*[(a-2i)s||c]{n+1}) + ***
1062 : ! *** f6*Ni(c)*[(a-1i)s||c-1i]{n+1} ***
1063 :
1064 11040 : DO la = 2, la_max
1065 :
1066 13926 : DO n = 1, nmax - la - lc
1067 :
1068 : ! *** Increase the angular momentum component z of a ***
1069 :
1070 : v(coset(0, 0, la), 1, coc, n) = &
1071 : rap(3)*v(coset(0, 0, la - 1), 1, coc, n) + &
1072 : rpw(3)*v(coset(0, 0, la - 1), 1, coc, n + 1) + &
1073 : f2*REAL(la - 1, dp)*(v(coset(0, 0, la - 2), 1, coc, n) + &
1074 : f4*v(coset(0, 0, la - 2), 1, coc, n + 1)) + &
1075 2886 : fcz*v(coset(0, 0, la - 1), 1, cocz, n + 1)
1076 :
1077 : ! *** Increase the angular momentum component y of a ***
1078 :
1079 2886 : az = la - 1
1080 : v(coset(0, 1, az), 1, coc, n) = &
1081 : rap(2)*v(coset(0, 0, az), 1, coc, n) + &
1082 : rpw(2)*v(coset(0, 0, az), 1, coc, n + 1) + &
1083 2886 : fcy*v(coset(0, 0, az), 1, cocy, n + 1)
1084 :
1085 5772 : DO ay = 2, la
1086 2886 : f3 = f2*REAL(ay - 1, dp)
1087 2886 : az = la - ay
1088 : v(coset(0, ay, az), 1, coc, n) = &
1089 : rap(2)*v(coset(0, ay - 1, az), 1, coc, n) + &
1090 : rpw(2)*v(coset(0, ay - 1, az), 1, coc, n + 1) + &
1091 : f3*(v(coset(0, ay - 2, az), 1, coc, n) + &
1092 : f4*v(coset(0, ay - 2, az), 1, coc, n + 1)) + &
1093 5772 : fcy*v(coset(0, ay - 1, az), 1, cocy, n + 1)
1094 : END DO
1095 :
1096 : ! *** Increase the angular momentum component x of a ***
1097 :
1098 8658 : DO ay = 0, la - 1
1099 5772 : az = la - 1 - ay
1100 : v(coset(1, ay, az), 1, coc, n) = &
1101 : rap(1)*v(coset(0, ay, az), 1, coc, n) + &
1102 : rpw(1)*v(coset(0, ay, az), 1, coc, n + 1) + &
1103 8658 : fcx*v(coset(0, ay, az), 1, cocx, n + 1)
1104 : END DO
1105 :
1106 7152 : DO ax = 2, la
1107 2886 : f3 = f2*REAL(ax - 1, dp)
1108 8658 : DO ay = 0, la - ax
1109 2886 : az = la - ax - ay
1110 : v(coset(ax, ay, az), 1, coc, n) = &
1111 : rap(1)*v(coset(ax - 1, ay, az), 1, coc, n) + &
1112 : rpw(1)*v(coset(ax - 1, ay, az), 1, coc, n + 1) + &
1113 : f3*(v(coset(ax - 2, ay, az), 1, coc, n) + &
1114 : f4*v(coset(ax - 2, ay, az), 1, coc, n + 1)) + &
1115 5772 : fcx*v(coset(ax - 1, ay, az), 1, cocx, n + 1)
1116 : END DO
1117 : END DO
1118 :
1119 : END DO
1120 :
1121 : END DO
1122 :
1123 : ! *** Recurrence steps: [as||c] -> [ab||c] ***
1124 :
1125 9660 : IF (lb_max > 0) THEN
1126 :
1127 : ! *** Horizontal recurrence steps ***
1128 :
1129 : ! *** [ap||c]{n} = [(a+1i)s||c]{n} - (Bi - Ai)*[as||c]{n} ***
1130 :
1131 5244 : la_start = MAX(0, la_min - 1)
1132 :
1133 10488 : DO la = la_start, la_max - 1
1134 23856 : DO n = 1, nmax - la - 1 - lc
1135 34098 : DO ax = 0, la
1136 46458 : DO ay = 0, la - ax
1137 17604 : az = la - ax - ay
1138 : v(coset(ax, ay, az), 2, coc, n) = &
1139 : v(coset(ax + 1, ay, az), 1, coc, n) - &
1140 17604 : rab(1)*v(coset(ax, ay, az), 1, coc, n)
1141 : v(coset(ax, ay, az), 3, coc, n) = &
1142 : v(coset(ax, ay + 1, az), 1, coc, n) - &
1143 17604 : rab(2)*v(coset(ax, ay, az), 1, coc, n)
1144 : v(coset(ax, ay, az), 4, coc, n) = &
1145 : v(coset(ax, ay, az + 1), 1, coc, n) - &
1146 33090 : rab(3)*v(coset(ax, ay, az), 1, coc, n)
1147 : END DO
1148 : END DO
1149 : END DO
1150 : END DO
1151 :
1152 : ! *** Vertical recurrence step ***
1153 :
1154 : ! *** [ap||c]{n} = (Pi - Bi)*[as||c]{n} + ***
1155 : ! *** (Wi - Pi)*[as||c]{n+1} + ***
1156 : ! *** f2*Ni(a)*( [(a-1i)s||c]{n} + ***
1157 : ! *** f4*[(a-1i)s||c]{n+1}) + ***
1158 : ! *** f6*Ni(c)*[(as||c-1i]{n+1}) ***
1159 :
1160 13368 : DO n = 1, nmax - la_max - 1 - lc
1161 30906 : DO ax = 0, la_max
1162 17538 : fx = f2*REAL(ax, dp)
1163 53904 : DO ay = 0, la_max - ax
1164 28242 : fy = f2*REAL(ay, dp)
1165 28242 : az = la_max - ax - ay
1166 28242 : fz = f2*REAL(az, dp)
1167 :
1168 28242 : IF (ax == 0) THEN
1169 : v(coset(ax, ay, az), 2, coc, n) = &
1170 : rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
1171 : rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1172 17538 : fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
1173 : ELSE
1174 : v(coset(ax, ay, az), 2, coc, n) = &
1175 : rbp(1)*v(coset(ax, ay, az), 1, coc, n) + &
1176 : rpw(1)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1177 : fx*(v(coset(ax - 1, ay, az), 1, coc, n) + &
1178 : f4*v(coset(ax - 1, ay, az), 1, coc, n + 1)) + &
1179 10704 : fcx*v(coset(ax, ay, az), 1, cocx, n + 1)
1180 : END IF
1181 :
1182 28242 : IF (ay == 0) THEN
1183 : v(coset(ax, ay, az), 3, coc, n) = &
1184 : rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
1185 : rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1186 17538 : fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
1187 : ELSE
1188 : v(coset(ax, ay, az), 3, coc, n) = &
1189 : rbp(2)*v(coset(ax, ay, az), 1, coc, n) + &
1190 : rpw(2)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1191 : fy*(v(coset(ax, ay - 1, az), 1, coc, n) + &
1192 : f4*v(coset(ax, ay - 1, az), 1, coc, n + 1)) + &
1193 10704 : fcy*v(coset(ax, ay, az), 1, cocy, n + 1)
1194 : END IF
1195 :
1196 45780 : IF (az == 0) THEN
1197 : v(coset(ax, ay, az), 4, coc, n) = &
1198 : rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
1199 : rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1200 17538 : fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
1201 : ELSE
1202 : v(coset(ax, ay, az), 4, coc, n) = &
1203 : rbp(3)*v(coset(ax, ay, az), 1, coc, n) + &
1204 : rpw(3)*v(coset(ax, ay, az), 1, coc, n + 1) + &
1205 : fz*(v(coset(ax, ay, az - 1), 1, coc, n) + &
1206 : f4*v(coset(ax, ay, az - 1), 1, coc, n + 1)) + &
1207 10704 : fcz*v(coset(ax, ay, az), 1, cocz, n + 1)
1208 : END IF
1209 :
1210 : END DO
1211 : END DO
1212 : END DO
1213 :
1214 : ! *** Recurrence steps: [ap||c] -> [ab||c] ***
1215 :
1216 6072 : DO lb = 2, lb_max
1217 :
1218 : ! *** Horizontal recurrence steps ***
1219 :
1220 : ! *** [ab||c]{n} = [(a+1i)(b-1i)||c]{n} - ***
1221 : ! *** (Bi - Ai)*[a(b-1i)||c]{n} ***
1222 :
1223 1656 : la_start = MAX(0, la_min - 1)
1224 :
1225 1656 : DO la = la_start, la_max - 1
1226 3636 : DO n = 1, nmax - la - lb - lc
1227 5118 : DO ax = 0, la
1228 6930 : DO ay = 0, la - ax
1229 2640 : az = la - ax - ay
1230 :
1231 : ! *** Shift of angular momentum component z ***
1232 :
1233 : v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1234 : v(coset(ax, ay, az + 1), &
1235 : coset(0, 0, lb - 1), coc, n) - &
1236 : rab(3)*v(coset(ax, ay, az), &
1237 2640 : coset(0, 0, lb - 1), coc, n)
1238 :
1239 : ! *** Shift of angular momentum component y ***
1240 :
1241 7920 : DO by = 1, lb
1242 5280 : bz = lb - by
1243 : v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1244 : v(coset(ax, ay + 1, az), &
1245 : coset(0, by - 1, bz), coc, n) - &
1246 : rab(2)*v(coset(ax, ay, az), &
1247 7920 : coset(0, by - 1, bz), coc, n)
1248 : END DO
1249 :
1250 : ! *** Shift of angular momentum component x ***
1251 :
1252 10230 : DO bx = 1, lb
1253 15840 : DO by = 0, lb - bx
1254 7920 : bz = lb - bx - by
1255 : v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1256 : v(coset(ax + 1, ay, az), &
1257 : coset(bx - 1, by, bz), coc, n) - &
1258 : rab(1)*v(coset(ax, ay, az), &
1259 13200 : coset(bx - 1, by, bz), coc, n)
1260 : END DO
1261 : END DO
1262 :
1263 : END DO
1264 : END DO
1265 : END DO
1266 : END DO
1267 :
1268 : ! *** Vertical recurrence step ***
1269 :
1270 : ! *** [ab||c]{n} = (Pi - Bi)*[a(b-1i)||c]{n} + ***
1271 : ! *** (Wi - Pi)*[a(b-1i)||c]{n+1} + ***
1272 : ! *** f2*Ni(a)*( [(a-1i)(b-1i)||c]{n} + ***
1273 : ! *** f4*[(a-1i)(b-1i)||c]{n+1}) ***
1274 : ! *** f2*Ni(b-1i)*( [a(b-2i)||c]{n} + ***
1275 : ! *** f4*[a(b-2i)||c]{n+1}) + ***
1276 : ! *** f6*Ni(c)*[a(b-1i)||c-1i]{n+1}) ***
1277 :
1278 7224 : DO n = 1, nmax - la_max - lb - lc
1279 4476 : DO ax = 0, la_max
1280 2496 : fx = f2*REAL(ax, dp)
1281 7680 : DO ay = 0, la_max - ax
1282 4032 : fy = f2*REAL(ay, dp)
1283 4032 : az = la_max - ax - ay
1284 4032 : fz = f2*REAL(az, dp)
1285 :
1286 : ! *** Shift of angular momentum component z from a to b ***
1287 :
1288 4032 : f3 = f2*REAL(lb - 1, dp)
1289 :
1290 4032 : IF (az == 0) THEN
1291 : v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1292 : rbp(3)*v(coset(ax, ay, az), &
1293 : coset(0, 0, lb - 1), coc, n) + &
1294 : rpw(3)*v(coset(ax, ay, az), &
1295 : coset(0, 0, lb - 1), coc, n + 1) + &
1296 : f3*(v(coset(ax, ay, az), &
1297 : coset(0, 0, lb - 2), coc, n) + &
1298 : f4*v(coset(ax, ay, az), &
1299 : coset(0, 0, lb - 2), coc, n + 1)) + &
1300 : fcz*v(coset(ax, ay, az), &
1301 2496 : coset(0, 0, lb - 1), cocz, n + 1)
1302 : ELSE
1303 : v(coset(ax, ay, az), coset(0, 0, lb), coc, n) = &
1304 : rbp(3)*v(coset(ax, ay, az), &
1305 : coset(0, 0, lb - 1), coc, n) + &
1306 : rpw(3)*v(coset(ax, ay, az), &
1307 : coset(0, 0, lb - 1), coc, n + 1) + &
1308 : fz*(v(coset(ax, ay, az - 1), &
1309 : coset(0, 0, lb - 1), coc, n) + &
1310 : f4*v(coset(ax, ay, az - 1), &
1311 : coset(0, 0, lb - 1), coc, n + 1)) + &
1312 : f3*(v(coset(ax, ay, az), &
1313 : coset(0, 0, lb - 2), coc, n) + &
1314 : f4*v(coset(ax, ay, az), &
1315 : coset(0, 0, lb - 2), coc, n + 1)) + &
1316 : fcz*v(coset(ax, ay, az), &
1317 1536 : coset(0, 0, lb - 1), cocz, n + 1)
1318 : END IF
1319 :
1320 : ! *** Shift of angular momentum component y from a to b ***
1321 :
1322 4032 : IF (ay == 0) THEN
1323 2496 : bz = lb - 1
1324 : v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1325 : rbp(2)*v(coset(ax, ay, az), &
1326 : coset(0, 0, bz), coc, n) + &
1327 : rpw(2)*v(coset(ax, ay, az), &
1328 : coset(0, 0, bz), coc, n + 1) + &
1329 : fcy*v(coset(ax, ay, az), &
1330 2496 : coset(0, 0, bz), cocy, n + 1)
1331 4992 : DO by = 2, lb
1332 2496 : bz = lb - by
1333 2496 : f3 = f2*REAL(by - 1, dp)
1334 : v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1335 : rbp(2)*v(coset(ax, ay, az), &
1336 : coset(0, by - 1, bz), coc, n) + &
1337 : rpw(2)*v(coset(ax, ay, az), &
1338 : coset(0, by - 1, bz), coc, n + 1) + &
1339 : f3*(v(coset(ax, ay, az), &
1340 : coset(0, by - 2, bz), coc, n) + &
1341 : f4*v(coset(ax, ay, az), &
1342 : coset(0, by - 2, bz), coc, n + 1)) + &
1343 : fcy*v(coset(ax, ay, az), &
1344 4992 : coset(0, by - 1, bz), cocy, n + 1)
1345 : END DO
1346 : ELSE
1347 1536 : bz = lb - 1
1348 : v(coset(ax, ay, az), coset(0, 1, bz), coc, n) = &
1349 : rbp(2)*v(coset(ax, ay, az), &
1350 : coset(0, 0, bz), coc, n) + &
1351 : rpw(2)*v(coset(ax, ay, az), &
1352 : coset(0, 0, bz), coc, n + 1) + &
1353 : fy*(v(coset(ax, ay - 1, az), &
1354 : coset(0, 0, bz), coc, n) + &
1355 : f4*v(coset(ax, ay - 1, az), &
1356 : coset(0, 0, bz), coc, n + 1)) + &
1357 : fcy*v(coset(ax, ay, az), &
1358 1536 : coset(0, 0, bz), cocy, n + 1)
1359 3072 : DO by = 2, lb
1360 1536 : bz = lb - by
1361 1536 : f3 = f2*REAL(by - 1, dp)
1362 : v(coset(ax, ay, az), coset(0, by, bz), coc, n) = &
1363 : rbp(2)*v(coset(ax, ay, az), &
1364 : coset(0, by - 1, bz), coc, n) + &
1365 : rpw(2)*v(coset(ax, ay, az), &
1366 : coset(0, by - 1, bz), coc, n + 1) + &
1367 : fy*(v(coset(ax, ay - 1, az), &
1368 : coset(0, by - 1, bz), coc, n) + &
1369 : f4*v(coset(ax, ay - 1, az), &
1370 : coset(0, by - 1, bz), coc, n + 1)) + &
1371 : f3*(v(coset(ax, ay, az), &
1372 : coset(0, by - 2, bz), coc, n) + &
1373 : f4*v(coset(ax, ay, az), &
1374 : coset(0, by - 2, bz), coc, n + 1)) + &
1375 : fcy*v(coset(ax, ay, az), &
1376 3072 : coset(0, by - 1, bz), cocy, n + 1)
1377 : END DO
1378 : END IF
1379 :
1380 : ! *** Shift of angular momentum component x from a to b ***
1381 :
1382 6528 : IF (ax == 0) THEN
1383 7488 : DO by = 0, lb - 1
1384 4992 : bz = lb - 1 - by
1385 : v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1386 : rbp(1)*v(coset(ax, ay, az), &
1387 : coset(0, by, bz), coc, n) + &
1388 : rpw(1)*v(coset(ax, ay, az), &
1389 : coset(0, by, bz), coc, n + 1) + &
1390 : fcx*v(coset(ax, ay, az), &
1391 7488 : coset(0, by, bz), cocx, n + 1)
1392 : END DO
1393 4992 : DO bx = 2, lb
1394 2496 : f3 = f2*REAL(bx - 1, dp)
1395 7488 : DO by = 0, lb - bx
1396 2496 : bz = lb - bx - by
1397 : v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1398 : rbp(1)*v(coset(ax, ay, az), &
1399 : coset(bx - 1, by, bz), coc, n) + &
1400 : rpw(1)*v(coset(ax, ay, az), &
1401 : coset(bx - 1, by, bz), coc, n + 1) + &
1402 : f3*(v(coset(ax, ay, az), &
1403 : coset(bx - 2, by, bz), coc, n) + &
1404 : f4*v(coset(ax, ay, az), &
1405 : coset(bx - 2, by, bz), coc, n + 1)) + &
1406 : fcx*v(coset(ax, ay, az), &
1407 4992 : coset(bx - 1, by, bz), cocx, n + 1)
1408 : END DO
1409 : END DO
1410 : ELSE
1411 4608 : DO by = 0, lb - 1
1412 3072 : bz = lb - 1 - by
1413 : v(coset(ax, ay, az), coset(1, by, bz), coc, n) = &
1414 : rbp(1)*v(coset(ax, ay, az), &
1415 : coset(0, by, bz), coc, n) + &
1416 : rpw(1)*v(coset(ax, ay, az), &
1417 : coset(0, by, bz), coc, n + 1) + &
1418 : fx*(v(coset(ax - 1, ay, az), &
1419 : coset(0, by, bz), coc, n) + &
1420 : f4*v(coset(ax - 1, ay, az), &
1421 : coset(0, by, bz), coc, n + 1)) + &
1422 : fcx*v(coset(ax, ay, az), &
1423 4608 : coset(0, by, bz), cocx, n + 1)
1424 : END DO
1425 3072 : DO bx = 2, lb
1426 1536 : f3 = f2*REAL(bx - 1, dp)
1427 4608 : DO by = 0, lb - bx
1428 1536 : bz = lb - bx - by
1429 : v(coset(ax, ay, az), coset(bx, by, bz), coc, n) = &
1430 : rbp(1)*v(coset(ax, ay, az), &
1431 : coset(bx - 1, by, bz), coc, n) + &
1432 : rpw(1)*v(coset(ax, ay, az), &
1433 : coset(bx - 1, by, bz), coc, n + 1) + &
1434 : fx*(v(coset(ax - 1, ay, az), &
1435 : coset(bx - 1, by, bz), coc, n) + &
1436 : f4*v(coset(ax - 1, ay, az), &
1437 : coset(bx - 1, by, bz), coc, n + 1)) + &
1438 : f3*(v(coset(ax, ay, az), &
1439 : coset(bx - 2, by, bz), coc, n) + &
1440 : f4*v(coset(ax, ay, az), &
1441 : coset(bx - 2, by, bz), coc, n + 1)) + &
1442 : fcx*v(coset(ax, ay, az), &
1443 3072 : coset(bx - 1, by, bz), cocx, n + 1)
1444 : END DO
1445 : END DO
1446 : END IF
1447 :
1448 : END DO
1449 : END DO
1450 : END DO
1451 :
1452 : END DO
1453 : END IF
1454 :
1455 : ELSE
1456 :
1457 11040 : IF (lb_max > 0) THEN
1458 :
1459 : ! *** Vertical recurrence steps: [ss||c] -> [sb||c] ***
1460 :
1461 : ! *** [sp||c]{n} = (Pi - Bi)*[ss||c]{n} + ***
1462 : ! *** (Wi - Pi)*[ss||c]{n+1} + ***
1463 : ! *** f6*Ni(c)**[ss||c-1i]{n+1} ***
1464 :
1465 11112 : DO n = 1, nmax - 1 - lc
1466 : v(1, 2, coc, n) = rbp(1)*v(1, 1, coc, n) + &
1467 : rpw(1)*v(1, 1, coc, n + 1) + &
1468 6696 : fcx*v(1, 1, cocx, n + 1)
1469 : v(1, 3, coc, n) = rbp(2)*v(1, 1, coc, n) + &
1470 : rpw(2)*v(1, 1, coc, n + 1) + &
1471 6696 : fcy*v(1, 1, cocy, n + 1)
1472 : v(1, 4, coc, n) = rbp(3)*v(1, 1, coc, n) + &
1473 : rpw(3)*v(1, 1, coc, n + 1) + &
1474 11112 : fcz*v(1, 1, cocz, n + 1)
1475 : END DO
1476 :
1477 : ! *** [sb||c]{n} = (Pi - Bi)*[s(b-1i)||c]{n} + ***
1478 : ! *** (Wi - Pi)*[s(b-1i)||c]{n+1} + ***
1479 : ! *** f2*Ni(b-1i)*( [s(b-2i)||c]{n} + ***
1480 : ! *** f4*[s(b-2i)||c]{n+1}) + ***
1481 : ! *** f6*Ni(c)**[s(b-1i)||c-1i]{n+1} ***
1482 :
1483 4968 : DO lb = 2, lb_max
1484 :
1485 5736 : DO n = 1, nmax - lb - lc
1486 :
1487 : ! *** Increase the angular momentum component z of b ***
1488 :
1489 : v(1, coset(0, 0, lb), coc, n) = &
1490 : rbp(3)*v(1, coset(0, 0, lb - 1), coc, n) + &
1491 : rpw(3)*v(1, coset(0, 0, lb - 1), coc, n + 1) + &
1492 : f2*REAL(lb - 1, dp)*(v(1, coset(0, 0, lb - 2), coc, n) + &
1493 : f4*v(1, coset(0, 0, lb - 2), coc, n + 1)) + &
1494 768 : fcz*v(1, coset(0, 0, lb - 1), cocz, n + 1)
1495 :
1496 : ! *** Increase the angular momentum component y of b ***
1497 :
1498 768 : bz = lb - 1
1499 : v(1, coset(0, 1, bz), coc, n) = &
1500 : rbp(2)*v(1, coset(0, 0, bz), coc, n) + &
1501 : rpw(2)*v(1, coset(0, 0, bz), coc, n + 1) + &
1502 768 : fcy*v(1, coset(0, 0, bz), cocy, n + 1)
1503 :
1504 1536 : DO by = 2, lb
1505 768 : f3 = f2*REAL(by - 1, dp)
1506 768 : bz = lb - by
1507 : v(1, coset(0, by, bz), coc, n) = &
1508 : rbp(2)*v(1, coset(0, by - 1, bz), coc, n) + &
1509 : rpw(2)*v(1, coset(0, by - 1, bz), coc, n + 1) + &
1510 : f3*(v(1, coset(0, by - 2, bz), coc, n) + &
1511 : f4*v(1, coset(0, by - 2, bz), coc, n + 1)) + &
1512 1536 : fcy*v(1, coset(0, by - 1, bz), cocy, n + 1)
1513 : END DO
1514 :
1515 : ! *** Increase the angular momentum component x of b ***
1516 :
1517 2304 : DO by = 0, lb - 1
1518 1536 : bz = lb - 1 - by
1519 : v(1, coset(1, by, bz), coc, n) = &
1520 : rbp(1)*v(1, coset(0, by, bz), coc, n) + &
1521 : rpw(1)*v(1, coset(0, by, bz), coc, n + 1) + &
1522 2304 : fcx*v(1, coset(0, by, bz), cocx, n + 1)
1523 : END DO
1524 :
1525 2088 : DO bx = 2, lb
1526 768 : f3 = f2*REAL(bx - 1, dp)
1527 2304 : DO by = 0, lb - bx
1528 768 : bz = lb - bx - by
1529 : v(1, coset(bx, by, bz), coc, n) = &
1530 : rbp(1)*v(1, coset(bx - 1, by, bz), coc, n) + &
1531 : rpw(1)*v(1, coset(bx - 1, by, bz), coc, n + 1) + &
1532 : f3*(v(1, coset(bx - 2, by, bz), coc, n) + &
1533 : f4*v(1, coset(bx - 2, by, bz), coc, n + 1)) + &
1534 1536 : fcx*v(1, coset(bx - 1, by, bz), cocx, n + 1)
1535 : END DO
1536 : END DO
1537 :
1538 : END DO
1539 :
1540 : END DO
1541 :
1542 : END IF
1543 :
1544 : END IF
1545 :
1546 : END DO
1547 : END DO
1548 :
1549 : END DO
1550 :
1551 : END IF
1552 :
1553 : ! *** Add the contribution of the current pair ***
1554 : ! *** of primitive Gaussian-type functions ***
1555 :
1556 20250 : DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1557 15750 : kk = k - ncoset(lc_min - 1)
1558 58050 : DO j = ncoset(lb_min - 1) + 1, ncoset(lb_max)
1559 152145 : DO i = ncoset(la_min - 1) + 1, ncoset(la_max - maxder_local)
1560 98595 : vabc(na + i, nb + j) = vabc(na + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1561 136395 : int_abc(na + i, nb + j, kk) = v(i, j, k, 1)
1562 : END DO
1563 : END DO
1564 : END DO
1565 :
1566 4500 : IF (PRESENT(maxder)) THEN
1567 0 : DO k = ncoset(lc_min - 1) + 1, ncoset(lc_max)
1568 0 : kk = k - ncoset(lc_min - 1)
1569 0 : DO j = 1, ncoset(lb_max)
1570 0 : DO i = 1, ncoset(la_max)
1571 0 : vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1572 : END DO
1573 : END DO
1574 : END DO
1575 : END IF
1576 :
1577 7200 : nb = nb + ncoset(lb_max)
1578 :
1579 : END DO
1580 :
1581 2700 : na = na + ncoset(la_max - maxder_local)
1582 4320 : nap = nap + ncoset(la_max)
1583 :
1584 : END DO
1585 :
1586 1620 : END SUBROUTINE coulomb3
1587 :
1588 : END MODULE ai_coulomb
|
