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