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 the angular momentum integrals over
10 : !> Cartesian Gaussian-type functions.
11 : !> \par Literature
12 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 : !> \par History
14 : !> none
15 : !> \author JGH (20.02.2005)
16 : ! **************************************************************************************************
17 : MODULE ai_angmom
18 :
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: pi
21 : USE orbital_pointers, ONLY: coset,&
22 : ncoset
23 : #include "../base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 :
27 : PRIVATE
28 :
29 : ! *** Public subroutines ***
30 :
31 : PUBLIC :: angmom
32 :
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief ...
37 : !> \param la_max ...
38 : !> \param npgfa ...
39 : !> \param zeta ...
40 : !> \param rpgfa ...
41 : !> \param la_min ...
42 : !> \param lb_max ...
43 : !> \param npgfb ...
44 : !> \param zetb ...
45 : !> \param rpgfb ...
46 : !> \param rac ...
47 : !> \param rbc ...
48 : !> \param angab ...
49 : ! **************************************************************************************************
50 1057605 : SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, &
51 211521 : lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
52 :
53 : INTEGER, INTENT(IN) :: la_max, npgfa
54 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
55 : INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
56 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
57 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
58 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: angab
59 :
60 : INTEGER :: ax, ay, az, bx, by, bz, i, ipgf, j, &
61 : jpgf, la, la_start, lb, na, nb
62 : REAL(KIND=dp) :: dab, f0, f1, f2, f3, fx, fy, fz, rab2, &
63 : zetp
64 : REAL(KIND=dp), DIMENSION(3) :: ac1, ac2, ac3, bc1, bc2, bc3, rab, rap, &
65 : rbp
66 : REAL(KIND=dp), &
67 423042 : DIMENSION(ncoset(la_max), ncoset(lb_max)) :: s
68 : REAL(KIND=dp), &
69 211521 : DIMENSION(ncoset(la_max), ncoset(lb_max), 3) :: as
70 :
71 : ! ax,ay,az : Angular momentum index numbers of orbital a.
72 : ! bx,by,bz : Angular momentum index numbers of orbital b.
73 : ! coset : Cartesian orbital set pointer.
74 : ! dab : Distance between the atomic centers a and b.
75 : ! l{a,b} : Angular momentum quantum number of shell a or b.
76 : ! l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
77 : ! l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
78 : ! rac : Distance vector between the atomic center a and C.
79 : ! rbc : Distance vector between the atomic center b and C.
80 : ! rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
81 : ! angab : Shell set of angular momentum integrals.
82 : ! zet{a,b} : Exponents of the Gaussian-type functions a or b.
83 : ! zetp : Reciprocal of the sum of the exponents of orbital a and b.
84 : ! *** Calculate the distance between the centers a and b ***
85 :
86 846084 : rab = rbc - rac
87 211521 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
88 211521 : dab = SQRT(rab2)
89 :
90 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
91 117260616 : angab = 0.0_dp
92 3966428 : s = 0.0_dp
93 12110805 : as = 0.0_dp
94 :
95 211521 : na = 0
96 :
97 726666 : DO ipgf = 1, npgfa
98 :
99 515145 : nb = 0
100 :
101 1957149 : DO jpgf = 1, npgfb
102 :
103 22306382 : s = 0.0_dp
104 68361150 : as = 0.0_dp
105 :
106 : ! *** Screening ***
107 :
108 1442004 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109 3662863 : DO j = nb + 1, nb + ncoset(lb_max)
110 13557247 : DO i = na + 1, na + ncoset(la_max)
111 9894384 : angab(i, j, 1) = 0.0_dp
112 9894384 : angab(i, j, 2) = 0.0_dp
113 12658419 : angab(i, j, 3) = 0.0_dp
114 : END DO
115 : END DO
116 898828 : nb = nb + ncoset(lb_max)
117 898828 : CYCLE
118 : END IF
119 :
120 : ! *** Calculate some prefactors ***
121 :
122 543176 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
123 :
124 543176 : f0 = (pi*zetp)**1.5_dp
125 543176 : f1 = zetb(jpgf)*zetp
126 543176 : f2 = 0.5_dp*zetp
127 : !
128 543176 : bc1(1) = 0._dp
129 543176 : bc1(2) = -f1*rbc(3)
130 543176 : bc1(3) = f1*rbc(2)
131 : !
132 543176 : bc2(1) = f1*rbc(3)
133 543176 : bc2(2) = 0._dp
134 543176 : bc2(3) = -f1*rbc(1)
135 : !
136 543176 : bc3(1) = -f1*rbc(2)
137 543176 : bc3(2) = f1*rbc(1)
138 543176 : bc3(3) = 0._dp
139 : !
140 543176 : ac1(1) = 0._dp
141 543176 : ac1(2) = zeta(ipgf)*zetp*rac(3)
142 543176 : ac1(3) = -zeta(ipgf)*zetp*rac(2)
143 : !
144 543176 : ac2(1) = -zeta(ipgf)*zetp*rac(3)
145 543176 : ac2(2) = 0._dp
146 543176 : ac2(3) = zeta(ipgf)*zetp*rac(1)
147 : !
148 543176 : ac3(1) = zeta(ipgf)*zetp*rac(2)
149 543176 : ac3(2) = -zeta(ipgf)*zetp*rac(1)
150 543176 : ac3(3) = 0._dp
151 : !
152 : ! *** Calculate the basic two-center overlap integral [s|s] ***
153 : ! *** Calculate the basic two-center angular momentum integral [s|L|s] ***
154 :
155 543176 : s(1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2)
156 543176 : as(1, 1, 1) = 2._dp*zeta(ipgf)*f1*(rac(2)*rbc(3) - rac(3)*rbc(2))*s(1, 1)
157 543176 : as(1, 1, 2) = 2._dp*zeta(ipgf)*f1*(rac(3)*rbc(1) - rac(1)*rbc(3))*s(1, 1)
158 543176 : as(1, 1, 3) = 2._dp*zeta(ipgf)*f1*(rac(1)*rbc(2) - rac(2)*rbc(1))*s(1, 1)
159 :
160 : ! *** Recurrence steps: [s|L|s] -> [a|L|b] ***
161 :
162 543176 : IF (la_max > 0) THEN
163 :
164 : ! *** Vertical recurrence steps: [s|L|s] -> [a|L|s] ***
165 :
166 1313204 : rap(:) = f1*rab(:)
167 :
168 : ! *** [p|s] = (Pi - Ai)*[s|s] (i = x,y,z) ***
169 : ! *** [p|Ln|s] = (Pi - Ai)*[s|Ln|s]+xb/(xa+xb)*(1i x BC)*[s|s] (i = x,y,z) ***
170 :
171 328301 : s(2, 1) = rap(1)*s(1, 1)
172 328301 : s(3, 1) = rap(2)*s(1, 1)
173 328301 : s(4, 1) = rap(3)*s(1, 1)
174 328301 : as(2, 1, 1) = rap(1)*as(1, 1, 1) + bc1(1)*s(1, 1)
175 328301 : as(2, 1, 2) = rap(1)*as(1, 1, 2) + bc1(2)*s(1, 1)
176 328301 : as(2, 1, 3) = rap(1)*as(1, 1, 3) + bc1(3)*s(1, 1)
177 328301 : as(3, 1, 1) = rap(2)*as(1, 1, 1) + bc2(1)*s(1, 1)
178 328301 : as(3, 1, 2) = rap(2)*as(1, 1, 2) + bc2(2)*s(1, 1)
179 328301 : as(3, 1, 3) = rap(2)*as(1, 1, 3) + bc2(3)*s(1, 1)
180 328301 : as(4, 1, 1) = rap(3)*as(1, 1, 1) + bc3(1)*s(1, 1)
181 328301 : as(4, 1, 2) = rap(3)*as(1, 1, 2) + bc3(2)*s(1, 1)
182 328301 : as(4, 1, 3) = rap(3)*as(1, 1, 3) + bc3(3)*s(1, 1)
183 :
184 : ! *** [a|s] = (Pi - Ai)*[a-1i|s] + f2*Ni(a-1i)*[a-2i|s] ***
185 : ! *** [a|Ln|s] = (Pi - Ai)*[a-1i|Ln|s] + f2*Ni(a-1i)*[a-2i|Ln|s] ***
186 : ! *** + xb/(xa+xb)*{(1i x BC)}_n*[a-1i|s] ***
187 :
188 358668 : DO la = 2, la_max
189 :
190 : ! *** Increase the angular momentum component z of function a ***
191 :
192 : s(coset(0, 0, la), 1) = rap(3)*s(coset(0, 0, la - 1), 1) + &
193 30367 : f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1)
194 : as(coset(0, 0, la), 1, 1) = rap(3)*as(coset(0, 0, la - 1), 1, 1) + &
195 : f2*REAL(la - 1, dp)*as(coset(0, 0, la - 2), 1, 1) + &
196 30367 : bc3(1)*s(coset(0, 0, la - 1), 1)
197 : as(coset(0, 0, la), 1, 2) = rap(3)*as(coset(0, 0, la - 1), 1, 2) + &
198 : f2*REAL(la - 1, dp)*as(coset(0, 0, la - 2), 1, 2) + &
199 30367 : bc3(2)*s(coset(0, 0, la - 1), 1)
200 : as(coset(0, 0, la), 1, 3) = rap(3)*as(coset(0, 0, la - 1), 1, 3) + &
201 : f2*REAL(la - 1, dp)*as(coset(0, 0, la - 2), 1, 3) + &
202 30367 : bc3(3)*s(coset(0, 0, la - 1), 1)
203 :
204 : ! *** Increase the angular momentum component y of function a ***
205 :
206 30367 : az = la - 1
207 30367 : s(coset(0, 1, az), 1) = rap(2)*s(coset(0, 0, az), 1)
208 : as(coset(0, 1, az), 1, 1) = rap(2)*as(coset(0, 0, az), 1, 1) + &
209 30367 : bc2(1)*s(coset(0, 0, az), 1)
210 : as(coset(0, 1, az), 1, 2) = rap(2)*as(coset(0, 0, az), 1, 2) + &
211 30367 : bc2(2)*s(coset(0, 0, az), 1)
212 : as(coset(0, 1, az), 1, 3) = rap(2)*as(coset(0, 0, az), 1, 3) + &
213 30367 : bc2(3)*s(coset(0, 0, az), 1)
214 :
215 60734 : DO ay = 2, la
216 30367 : az = la - ay
217 : s(coset(0, ay, az), 1) = rap(2)*s(coset(0, ay - 1, az), 1) + &
218 30367 : f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1)
219 : as(coset(0, ay, az), 1, 1) = rap(2)*as(coset(0, ay - 1, az), 1, 1) + &
220 : f2*REAL(ay - 1, dp)*as(coset(0, ay - 2, az), 1, 1) + &
221 30367 : bc2(1)*s(coset(0, ay - 1, az), 1)
222 : as(coset(0, ay, az), 1, 2) = rap(2)*as(coset(0, ay - 1, az), 1, 2) + &
223 : f2*REAL(ay - 1, dp)*as(coset(0, ay - 2, az), 1, 2) + &
224 30367 : bc2(2)*s(coset(0, ay - 1, az), 1)
225 : as(coset(0, ay, az), 1, 3) = rap(2)*as(coset(0, ay - 1, az), 1, 3) + &
226 : f2*REAL(ay - 1, dp)*as(coset(0, ay - 2, az), 1, 3) + &
227 60734 : bc2(3)*s(coset(0, ay - 1, az), 1)
228 : END DO
229 :
230 : ! *** Increase the angular momentum component x of function a ***
231 :
232 91101 : DO ay = 0, la - 1
233 60734 : az = la - 1 - ay
234 60734 : s(coset(1, ay, az), 1) = rap(1)*s(coset(0, ay, az), 1)
235 : as(coset(1, ay, az), 1, 1) = rap(1)*as(coset(0, ay, az), 1, 1) + &
236 60734 : bc1(1)*s(coset(0, ay, az), 1)
237 : as(coset(1, ay, az), 1, 2) = rap(1)*as(coset(0, ay, az), 1, 2) + &
238 60734 : bc1(2)*s(coset(0, ay, az), 1)
239 : as(coset(1, ay, az), 1, 3) = rap(1)*as(coset(0, ay, az), 1, 3) + &
240 91101 : bc1(3)*s(coset(0, ay, az), 1)
241 : END DO
242 :
243 389035 : DO ax = 2, la
244 30367 : f3 = f2*REAL(ax - 1, dp)
245 91101 : DO ay = 0, la - ax
246 30367 : az = la - ax - ay
247 : s(coset(ax, ay, az), 1) = rap(1)*s(coset(ax - 1, ay, az), 1) + &
248 30367 : f3*s(coset(ax - 2, ay, az), 1)
249 : as(coset(ax, ay, az), 1, 1) = rap(1)*as(coset(ax - 1, ay, az), 1, 1) + &
250 : f3*as(coset(ax - 2, ay, az), 1, 1) + &
251 30367 : bc1(1)*s(coset(ax - 1, ay, az), 1)
252 : as(coset(ax, ay, az), 1, 2) = rap(1)*as(coset(ax - 1, ay, az), 1, 2) + &
253 : f3*as(coset(ax - 2, ay, az), 1, 2) + &
254 30367 : bc1(2)*s(coset(ax - 1, ay, az), 1)
255 : as(coset(ax, ay, az), 1, 3) = rap(1)*as(coset(ax - 1, ay, az), 1, 3) + &
256 : f3*as(coset(ax - 2, ay, az), 1, 3) + &
257 60734 : bc1(3)*s(coset(ax - 1, ay, az), 1)
258 : END DO
259 : END DO
260 :
261 : END DO
262 :
263 : ! *** Recurrence steps: [a|L|s] -> [a|L|b] ***
264 :
265 328301 : IF (lb_max > 0) THEN
266 :
267 1227546 : DO j = 2, ncoset(lb_max)
268 5657760 : DO i = 1, ncoset(la_max)
269 4430214 : s(i, j) = 0.0_dp
270 4430214 : as(i, j, 1) = 0.0_dp
271 4430214 : as(i, j, 2) = 0.0_dp
272 5400612 : as(i, j, 3) = 0.0_dp
273 : END DO
274 : END DO
275 :
276 : ! *** Horizontal recurrence steps ***
277 :
278 1028592 : rbp(:) = rap(:) - rab(:)
279 :
280 : ! *** [a|L|p] = [a+1i|Lm|s] - (Bi - Ai)*[a|Lm|s] ***
281 : ! *** + [a+1k|s] + (Ak - Ck)*[a|s] eps(i,m,k)
282 :
283 257148 : IF (lb_max == 1) THEN
284 228889 : la_start = la_min
285 : ELSE
286 28259 : la_start = MAX(0, la_min - 1)
287 : END IF
288 :
289 486954 : DO la = la_start, la_max - 1
290 719478 : DO ax = 0, la
291 697572 : DO ay = 0, la - ax
292 235242 : az = la - ax - ay
293 : s(coset(ax, ay, az), 2) = s(coset(ax + 1, ay, az), 1) - &
294 235242 : rab(1)*s(coset(ax, ay, az), 1)
295 : s(coset(ax, ay, az), 3) = s(coset(ax, ay + 1, az), 1) - &
296 235242 : rab(2)*s(coset(ax, ay, az), 1)
297 : s(coset(ax, ay, az), 4) = s(coset(ax, ay, az + 1), 1) - &
298 235242 : rab(3)*s(coset(ax, ay, az), 1)
299 : as(coset(ax, ay, az), 2, 1) = as(coset(ax + 1, ay, az), 1, 1) - &
300 235242 : rab(1)*as(coset(ax, ay, az), 1, 1)
301 : as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay + 1, az), 1, 1) - &
302 : rab(2)*as(coset(ax, ay, az), 1, 1) &
303 235242 : - s(coset(ax, ay, az + 1), 1) - rac(3)*s(coset(ax, ay, az), 1)
304 : as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az + 1), 1, 1) - &
305 : rab(3)*as(coset(ax, ay, az), 1, 1) &
306 235242 : + s(coset(ax, ay + 1, az), 1) + rac(2)*s(coset(ax, ay, az), 1)
307 : as(coset(ax, ay, az), 2, 2) = as(coset(ax + 1, ay, az), 1, 2) - &
308 : rab(1)*as(coset(ax, ay, az), 1, 2) &
309 235242 : + s(coset(ax, ay, az + 1), 1) + rac(3)*s(coset(ax, ay, az), 1)
310 : as(coset(ax, ay, az), 3, 2) = as(coset(ax, ay + 1, az), 1, 2) - &
311 235242 : rab(2)*as(coset(ax, ay, az), 1, 2)
312 : as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az + 1), 1, 2) - &
313 : rab(3)*as(coset(ax, ay, az), 1, 2) &
314 235242 : - s(coset(ax + 1, ay, az), 1) - rac(1)*s(coset(ax, ay, az), 1)
315 : as(coset(ax, ay, az), 2, 3) = as(coset(ax + 1, ay, az), 1, 3) - &
316 : rab(1)*as(coset(ax, ay, az), 1, 3) &
317 235242 : - s(coset(ax, ay + 1, az), 1) - rac(2)*s(coset(ax, ay, az), 1)
318 : as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay + 1, az), 1, 3) - &
319 : rab(2)*as(coset(ax, ay, az), 1, 3) &
320 235242 : + s(coset(ax + 1, ay, az), 1) + rac(1)*s(coset(ax, ay, az), 1)
321 : as(coset(ax, ay, az), 4, 3) = as(coset(ax, ay, az + 1), 1, 3) - &
322 467766 : rab(3)*as(coset(ax, ay, az), 1, 3)
323 : END DO
324 : END DO
325 : END DO
326 :
327 : ! *** Vertical recurrence step ***
328 :
329 : ! *** [a|p] = (Pi - Bi)*[a|s] + f2*Ni(a)*[a-1i|s] ***
330 : ! *** [a|L|p] = (Pi - Bi)*[a|L|s] + f2*Ni(a)*[a-1i|L|s] ***
331 : ! *** + xa/(xa+xb)*(AC x 1i)*[a|s] ***
332 : ! *** + 0.5*zetp*SUM_k Nk(a)*(1k x 1i)*[a-1k|s] ***
333 :
334 797197 : DO ax = 0, la_max
335 540049 : fx = f2*REAL(ax, dp)
336 1645900 : DO ay = 0, la_max - ax
337 848703 : fy = f2*REAL(ay, dp)
338 848703 : az = la_max - ax - ay
339 848703 : fz = f2*REAL(az, dp)
340 848703 : IF (ax == 0) THEN
341 540049 : s(coset(ax, ay, az), 2) = rbp(1)*s(coset(ax, ay, az), 1)
342 : as(coset(ax, ay, az), 2, 1) = rbp(1)*as(coset(ax, ay, az), 1, 1) + &
343 540049 : ac1(1)*s(coset(ax, ay, az), 1)
344 : as(coset(ax, ay, az), 2, 2) = rbp(1)*as(coset(ax, ay, az), 1, 2) + &
345 540049 : ac1(2)*s(coset(ax, ay, az), 1)
346 : as(coset(ax, ay, az), 2, 3) = rbp(1)*as(coset(ax, ay, az), 1, 3) + &
347 540049 : ac1(3)*s(coset(ax, ay, az), 1)
348 : ELSE
349 : s(coset(ax, ay, az), 2) = rbp(1)*s(coset(ax, ay, az), 1) + &
350 308654 : fx*s(coset(ax - 1, ay, az), 1)
351 : as(coset(ax, ay, az), 2, 1) = rbp(1)*as(coset(ax, ay, az), 1, 1) + &
352 : fx*as(coset(ax - 1, ay, az), 1, 1) + &
353 308654 : ac1(1)*s(coset(ax, ay, az), 1)
354 : as(coset(ax, ay, az), 2, 2) = rbp(1)*as(coset(ax, ay, az), 1, 2) + &
355 : fx*as(coset(ax - 1, ay, az), 1, 2) + &
356 308654 : ac1(2)*s(coset(ax, ay, az), 1)
357 : as(coset(ax, ay, az), 2, 3) = rbp(1)*as(coset(ax, ay, az), 1, 3) + &
358 : fx*as(coset(ax - 1, ay, az), 1, 3) + &
359 308654 : ac1(3)*s(coset(ax, ay, az), 1)
360 : END IF
361 848703 : IF (az > 0) as(coset(ax, ay, az), 2, 2) = as(coset(ax, ay, az), 2, 2) + &
362 308654 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
363 848703 : IF (ay > 0) as(coset(ax, ay, az), 2, 3) = as(coset(ax, ay, az), 2, 3) - &
364 308654 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
365 848703 : IF (ay == 0) THEN
366 540049 : s(coset(ax, ay, az), 3) = rbp(2)*s(coset(ax, ay, az), 1)
367 : as(coset(ax, ay, az), 3, 1) = rbp(2)*as(coset(ax, ay, az), 1, 1) + &
368 540049 : ac2(1)*s(coset(ax, ay, az), 1)
369 : as(coset(ax, ay, az), 3, 2) = rbp(2)*as(coset(ax, ay, az), 1, 2) + &
370 540049 : ac2(2)*s(coset(ax, ay, az), 1)
371 : as(coset(ax, ay, az), 3, 3) = rbp(2)*as(coset(ax, ay, az), 1, 3) + &
372 540049 : ac2(3)*s(coset(ax, ay, az), 1)
373 : ELSE
374 : s(coset(ax, ay, az), 3) = rbp(2)*s(coset(ax, ay, az), 1) + &
375 308654 : fy*s(coset(ax, ay - 1, az), 1)
376 : as(coset(ax, ay, az), 3, 1) = rbp(2)*as(coset(ax, ay, az), 1, 1) + &
377 : fy*as(coset(ax, ay - 1, az), 1, 1) + &
378 308654 : ac2(1)*s(coset(ax, ay, az), 1)
379 : as(coset(ax, ay, az), 3, 2) = rbp(2)*as(coset(ax, ay, az), 1, 2) + &
380 : fy*as(coset(ax, ay - 1, az), 1, 2) + &
381 308654 : ac2(2)*s(coset(ax, ay, az), 1)
382 : as(coset(ax, ay, az), 3, 3) = rbp(2)*as(coset(ax, ay, az), 1, 3) + &
383 : fy*as(coset(ax, ay - 1, az), 1, 3) + &
384 308654 : ac2(3)*s(coset(ax, ay, az), 1)
385 : END IF
386 848703 : IF (az > 0) as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay, az), 3, 1) - &
387 308654 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
388 848703 : IF (ax > 0) as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay, az), 3, 3) + &
389 308654 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1)
390 848703 : IF (az == 0) THEN
391 540049 : s(coset(ax, ay, az), 4) = rbp(3)*s(coset(ax, ay, az), 1)
392 : as(coset(ax, ay, az), 4, 1) = rbp(3)*as(coset(ax, ay, az), 1, 1) + &
393 540049 : ac3(1)*s(coset(ax, ay, az), 1)
394 : as(coset(ax, ay, az), 4, 2) = rbp(3)*as(coset(ax, ay, az), 1, 2) + &
395 540049 : ac3(2)*s(coset(ax, ay, az), 1)
396 : as(coset(ax, ay, az), 4, 3) = rbp(3)*as(coset(ax, ay, az), 1, 3) + &
397 540049 : ac3(3)*s(coset(ax, ay, az), 1)
398 : ELSE
399 : s(coset(ax, ay, az), 4) = rbp(3)*s(coset(ax, ay, az), 1) + &
400 308654 : fz*s(coset(ax, ay, az - 1), 1)
401 : as(coset(ax, ay, az), 4, 1) = rbp(3)*as(coset(ax, ay, az), 1, 1) + &
402 : fz*as(coset(ax, ay, az - 1), 1, 1) + &
403 308654 : ac3(1)*s(coset(ax, ay, az), 1)
404 : as(coset(ax, ay, az), 4, 2) = rbp(3)*as(coset(ax, ay, az), 1, 2) + &
405 : fz*as(coset(ax, ay, az - 1), 1, 2) + &
406 308654 : ac3(2)*s(coset(ax, ay, az), 1)
407 : as(coset(ax, ay, az), 4, 3) = rbp(3)*as(coset(ax, ay, az), 1, 3) + &
408 : fz*as(coset(ax, ay, az - 1), 1, 3) + &
409 308654 : ac3(3)*s(coset(ax, ay, az), 1)
410 : END IF
411 848703 : IF (ay > 0) as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az), 4, 1) + &
412 308654 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
413 848703 : IF (ax > 0) as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az), 4, 2) - &
414 848703 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1)
415 : END DO
416 : END DO
417 :
418 : ! *** Recurrence steps: [a|L|p] -> [a|L|b] ***
419 :
420 287927 : DO lb = 2, lb_max
421 :
422 : ! *** Horizontal recurrence steps ***
423 :
424 : ! *** [a|Lm|b] = [a+1i|Lm|b-1i] - (Bi - Ai)*[a|Lm|b-1i] ***
425 : ! *** + [a+1k|b-1i] + (Ak - Ck)*[a|b-1i] eps(i,m,k)
426 :
427 30779 : IF (lb == lb_max) THEN
428 28259 : la_start = la_min
429 : ELSE
430 2520 : la_start = MAX(0, la_min - 1)
431 : END IF
432 :
433 59220 : DO la = la_start, la_max - 1
434 88016 : DO ax = 0, la
435 86388 : DO ay = 0, la - ax
436 29151 : az = la - ax - ay
437 :
438 : ! *** Shift of angular momentum component z from a to b ***
439 :
440 : s(coset(ax, ay, az), coset(0, 0, lb)) = &
441 : s(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - &
442 29151 : rab(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
443 : as(coset(ax, ay, az), coset(0, 0, lb), 1) = &
444 : as(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 1) - &
445 : rab(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 1) &
446 : + s(coset(ax, ay + 1, az), coset(0, 0, lb - 1)) &
447 29151 : + rac(2)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
448 : as(coset(ax, ay, az), coset(0, 0, lb), 2) = &
449 : as(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 2) - &
450 : rab(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 2) &
451 : - s(coset(ax + 1, ay, az), coset(0, 0, lb - 1)) &
452 29151 : - rac(1)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
453 : as(coset(ax, ay, az), coset(0, 0, lb), 3) = &
454 : as(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 3) - &
455 29151 : rab(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 3)
456 :
457 : ! *** Shift of angular momentum component y from a to b ***
458 :
459 90399 : DO by = 1, lb
460 61248 : bz = lb - by
461 : s(coset(ax, ay, az), coset(0, by, bz)) = &
462 : s(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - &
463 61248 : rab(2)*s(coset(ax, ay, az), coset(0, by - 1, bz))
464 : as(coset(ax, ay, az), coset(0, by, bz), 1) = &
465 : as(coset(ax, ay + 1, az), coset(0, by - 1, bz), 1) - &
466 : rab(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 1) &
467 : - s(coset(ax, ay, az + 1), coset(0, by - 1, bz)) &
468 61248 : - rac(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
469 : as(coset(ax, ay, az), coset(0, by, bz), 2) = &
470 : as(coset(ax, ay + 1, az), coset(0, by - 1, bz), 2) - &
471 61248 : rab(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 2)
472 : as(coset(ax, ay, az), coset(0, by, bz), 3) = &
473 : as(coset(ax, ay + 1, az), coset(0, by - 1, bz), 3) - &
474 : rab(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 3) &
475 : + s(coset(ax + 1, ay, az), coset(0, by - 1, bz)) &
476 90399 : + rac(1)*s(coset(ax, ay, az), coset(0, by - 1, bz))
477 : END DO
478 :
479 : ! *** Shift of angular momentum component x from a to b ***
480 :
481 119195 : DO bx = 1, lb
482 187392 : DO by = 0, lb - bx
483 96993 : bz = lb - bx - by
484 : s(coset(ax, ay, az), coset(bx, by, bz)) = &
485 : s(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - &
486 96993 : rab(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
487 : as(coset(ax, ay, az), coset(bx, by, bz), 1) = &
488 : as(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 1) - &
489 96993 : rab(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 1)
490 : as(coset(ax, ay, az), coset(bx, by, bz), 2) = &
491 : as(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 2) - &
492 : rab(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 2) &
493 : + s(coset(ax, ay, az + 1), coset(bx - 1, by, bz)) &
494 96993 : + rac(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
495 : as(coset(ax, ay, az), coset(bx, by, bz), 3) = &
496 : as(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 3) - &
497 : rab(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 3) &
498 : - s(coset(ax, ay + 1, az), coset(bx - 1, by, bz)) &
499 158241 : - rac(2)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
500 : END DO
501 : END DO
502 :
503 : END DO
504 : END DO
505 : END DO
506 :
507 : ! *** Vertical recurrence step ***
508 :
509 : ! *** [a|b] = (Pi - Bi)*[a|b-1i] + f2*Ni(a)*[a-1i|b-1i] + ***
510 : ! *** f2*Ni(b-1i)*[a|b-2i] ***
511 : ! *** [a|L|b] = (Pi - Bi)*[a|L|b-1i] + f2*Ni(a)*[a-1i|L|b-1i] + ***
512 : ! *** f2*Ni(b-1i)*[a|L|b-2i] ***
513 : ! *** + xa/(xa+xb)*(AC x 1i)*[a|b-1i] ***
514 : ! *** + 0.5*zetp*SUM_k Nk(a)*(1k x 1i)*[a-1k|b-1i] ***
515 :
516 351848 : DO ax = 0, la_max
517 63921 : fx = f2*REAL(ax, dp)
518 194126 : DO ay = 0, la_max - ax
519 99426 : fy = f2*REAL(ay, dp)
520 99426 : az = la_max - ax - ay
521 99426 : fz = f2*REAL(az, dp)
522 :
523 : ! *** Increase the angular momentum component z of function b ***
524 :
525 99426 : f3 = f2*REAL(lb - 1, dp)
526 :
527 99426 : IF (az == 0) THEN
528 : s(coset(ax, ay, az), coset(0, 0, lb)) = &
529 : rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
530 63921 : f3*s(coset(ax, ay, az), coset(0, 0, lb - 2))
531 : as(coset(ax, ay, az), coset(0, 0, lb), 1) = &
532 : rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 1) + &
533 : f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 1) + &
534 63921 : ac3(1)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
535 : as(coset(ax, ay, az), coset(0, 0, lb), 2) = &
536 : rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 2) + &
537 : f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 2) + &
538 63921 : ac3(2)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
539 : as(coset(ax, ay, az), coset(0, 0, lb), 3) = &
540 : rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 3) + &
541 : f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 3) + &
542 63921 : ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
543 : ELSE
544 : s(coset(ax, ay, az), coset(0, 0, lb)) = &
545 : rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) + &
546 : fz*s(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + &
547 35505 : f3*s(coset(ax, ay, az), coset(0, 0, lb - 2))
548 : as(coset(ax, ay, az), coset(0, 0, lb), 1) = &
549 : rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 1) + &
550 : fz*as(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1) + &
551 : f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 1) + &
552 35505 : ac3(1)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
553 : as(coset(ax, ay, az), coset(0, 0, lb), 2) = &
554 : rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 2) + &
555 : fz*as(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 2) + &
556 : f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 2) + &
557 35505 : ac3(2)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
558 : as(coset(ax, ay, az), coset(0, 0, lb), 3) = &
559 : rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 3) + &
560 : fz*as(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 3) + &
561 : f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 3) + &
562 35505 : ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
563 : END IF
564 99426 : IF (ay > 0) as(coset(ax, ay, az), coset(0, 0, lb), 1) = &
565 : as(coset(ax, ay, az), coset(0, 0, lb), 1) + &
566 35505 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, 0, lb - 1))
567 99426 : IF (ax > 0) as(coset(ax, ay, az), coset(0, 0, lb), 2) = &
568 : as(coset(ax, ay, az), coset(0, 0, lb), 2) - &
569 35505 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, lb - 1))
570 :
571 : ! *** Increase the angular momentum component y of function b ***
572 :
573 99426 : IF (ay == 0) THEN
574 63921 : bz = lb - 1
575 : s(coset(ax, ay, az), coset(0, 1, bz)) = &
576 63921 : rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz))
577 : as(coset(ax, ay, az), coset(0, 1, bz), 1) = &
578 : rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 1) + &
579 63921 : ac2(1)*s(coset(ax, ay, az), coset(0, 0, bz))
580 : as(coset(ax, ay, az), coset(0, 1, bz), 2) = &
581 : rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 2) + &
582 63921 : ac2(2)*s(coset(ax, ay, az), coset(0, 0, bz))
583 : as(coset(ax, ay, az), coset(0, 1, bz), 3) = &
584 : rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 3) + &
585 63921 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
586 63921 : IF (az > 0) as(coset(ax, ay, az), coset(0, 1, bz), 1) = &
587 : as(coset(ax, ay, az), coset(0, 1, bz), 1) - &
588 33142 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
589 63921 : IF (ax > 0) as(coset(ax, ay, az), coset(0, 1, bz), 3) = &
590 : as(coset(ax, ay, az), coset(0, 1, bz), 3) + &
591 33142 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
592 134562 : DO by = 2, lb
593 70641 : bz = lb - by
594 70641 : f3 = f2*REAL(by - 1, dp)
595 : s(coset(ax, ay, az), coset(0, by, bz)) = &
596 : rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz)) + &
597 70641 : f3*s(coset(ax, ay, az), coset(0, by - 2, bz))
598 : as(coset(ax, ay, az), coset(0, by, bz), 1) = &
599 : rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 1) + &
600 : f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 1) + &
601 70641 : ac2(1)*s(coset(ax, ay, az), coset(0, by - 1, bz))
602 : as(coset(ax, ay, az), coset(0, by, bz), 2) = &
603 : rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 2) + &
604 : f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 2) + &
605 70641 : ac2(2)*s(coset(ax, ay, az), coset(0, by - 1, bz))
606 : as(coset(ax, ay, az), coset(0, by, bz), 3) = &
607 : rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 3) + &
608 : f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 3) + &
609 70641 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
610 70641 : IF (az > 0) as(coset(ax, ay, az), coset(0, by, bz), 1) = &
611 : as(coset(ax, ay, az), coset(0, by, bz), 1) - &
612 36502 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
613 70641 : IF (ax > 0) as(coset(ax, ay, az), coset(0, by, bz), 3) = &
614 : as(coset(ax, ay, az), coset(0, by, bz), 3) + &
615 100423 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz))
616 : END DO
617 : ELSE
618 35505 : bz = lb - 1
619 : s(coset(ax, ay, az), coset(0, 1, bz)) = &
620 : rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz)) + &
621 35505 : fy*s(coset(ax, ay - 1, az), coset(0, 0, bz))
622 : as(coset(ax, ay, az), coset(0, 1, bz), 1) = &
623 : rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 1) + &
624 : fy*as(coset(ax, ay - 1, az), coset(0, 0, bz), 1) + &
625 35505 : ac2(1)*s(coset(ax, ay, az), coset(0, 0, bz))
626 : as(coset(ax, ay, az), coset(0, 1, bz), 2) = &
627 : rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 2) + &
628 : fy*as(coset(ax, ay - 1, az), coset(0, 0, bz), 2) + &
629 35505 : ac2(2)*s(coset(ax, ay, az), coset(0, 0, bz))
630 : as(coset(ax, ay, az), coset(0, 1, bz), 3) = &
631 : rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 3) + &
632 : fy*as(coset(ax, ay - 1, az), coset(0, 0, bz), 3) + &
633 35505 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
634 35505 : IF (az > 0) as(coset(ax, ay, az), coset(0, 1, bz), 1) = &
635 : as(coset(ax, ay, az), coset(0, 1, bz), 1) - &
636 2363 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
637 35505 : IF (ax > 0) as(coset(ax, ay, az), coset(0, 1, bz), 3) = &
638 : as(coset(ax, ay, az), coset(0, 1, bz), 3) + &
639 2363 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
640 74370 : DO by = 2, lb
641 38865 : bz = lb - by
642 38865 : f3 = f2*REAL(by - 1, dp)
643 : s(coset(ax, ay, az), coset(0, by, bz)) = &
644 : rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz)) + &
645 : fy*s(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + &
646 38865 : f3*s(coset(ax, ay, az), coset(0, by - 2, bz))
647 : as(coset(ax, ay, az), coset(0, by, bz), 1) = &
648 : rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 1) + &
649 : fy*as(coset(ax, ay - 1, az), coset(0, by - 1, bz), 1) + &
650 : f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 1) + &
651 38865 : ac2(1)*s(coset(ax, ay, az), coset(0, by - 1, bz))
652 : as(coset(ax, ay, az), coset(0, by, bz), 2) = &
653 : rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 2) + &
654 : fy*as(coset(ax, ay - 1, az), coset(0, by - 1, bz), 2) + &
655 : f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 2) + &
656 38865 : ac2(2)*s(coset(ax, ay, az), coset(0, by - 1, bz))
657 : as(coset(ax, ay, az), coset(0, by, bz), 3) = &
658 : rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 3) + &
659 : fy*as(coset(ax, ay - 1, az), coset(0, by - 1, bz), 3) + &
660 : f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 3) + &
661 38865 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
662 38865 : IF (az > 0) as(coset(ax, ay, az), coset(0, by, bz), 1) = &
663 : as(coset(ax, ay, az), coset(0, by, bz), 1) - &
664 2363 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
665 38865 : IF (ax > 0) as(coset(ax, ay, az), coset(0, by, bz), 3) = &
666 : as(coset(ax, ay, az), coset(0, by, bz), 3) + &
667 37868 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz))
668 : END DO
669 : END IF
670 :
671 : ! *** Increase the angular momentum component x of function b ***
672 :
673 163347 : IF (ax == 0) THEN
674 198483 : DO by = 0, lb - 1
675 134562 : bz = lb - 1 - by
676 : s(coset(ax, ay, az), coset(1, by, bz)) = &
677 134562 : rbp(1)*s(coset(ax, ay, az), coset(0, by, bz))
678 : as(coset(ax, ay, az), coset(1, by, bz), 1) = &
679 : rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 1) + &
680 134562 : ac1(1)*s(coset(ax, ay, az), coset(0, by, bz))
681 : as(coset(ax, ay, az), coset(1, by, bz), 2) = &
682 : rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 2) + &
683 134562 : ac1(2)*s(coset(ax, ay, az), coset(0, by, bz))
684 : as(coset(ax, ay, az), coset(1, by, bz), 3) = &
685 : rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 3) + &
686 134562 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
687 134562 : IF (az > 0) as(coset(ax, ay, az), coset(1, by, bz), 2) = &
688 : as(coset(ax, ay, az), coset(1, by, bz), 2) + &
689 69644 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
690 134562 : IF (ay > 0) as(coset(ax, ay, az), coset(1, by, bz), 3) = &
691 : as(coset(ax, ay, az), coset(1, by, bz), 3) - &
692 133565 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
693 : END DO
694 134562 : DO bx = 2, lb
695 70641 : f3 = f2*REAL(bx - 1, dp)
696 213603 : DO by = 0, lb - bx
697 79041 : bz = lb - bx - by
698 : s(coset(ax, ay, az), coset(bx, by, bz)) = &
699 : rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) + &
700 79041 : f3*s(coset(ax, ay, az), coset(bx - 2, by, bz))
701 : as(coset(ax, ay, az), coset(bx, by, bz), 1) = &
702 : rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 1) + &
703 : f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 1) + &
704 79041 : ac1(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
705 : as(coset(ax, ay, az), coset(bx, by, bz), 2) = &
706 : rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 2) + &
707 : f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 2) + &
708 79041 : ac1(2)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
709 : as(coset(ax, ay, az), coset(bx, by, bz), 3) = &
710 : rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 3) + &
711 : f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 3) + &
712 79041 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
713 79041 : IF (az > 0) as(coset(ax, ay, az), coset(bx, by, bz), 2) = &
714 : as(coset(ax, ay, az), coset(bx, by, bz), 2) + &
715 40702 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
716 79041 : IF (ay > 0) as(coset(ax, ay, az), coset(bx, by, bz), 3) = &
717 : as(coset(ax, ay, az), coset(bx, by, bz), 3) - &
718 111343 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz))
719 : END DO
720 : END DO
721 : ELSE
722 109875 : DO by = 0, lb - 1
723 74370 : bz = lb - 1 - by
724 : s(coset(ax, ay, az), coset(1, by, bz)) = &
725 : rbp(1)*s(coset(ax, ay, az), coset(0, by, bz)) + &
726 74370 : fx*s(coset(ax - 1, ay, az), coset(0, by, bz))
727 : as(coset(ax, ay, az), coset(1, by, bz), 1) = &
728 : rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 1) + &
729 : fx*as(coset(ax - 1, ay, az), coset(0, by, bz), 1) + &
730 74370 : ac1(1)*s(coset(ax, ay, az), coset(0, by, bz))
731 : as(coset(ax, ay, az), coset(1, by, bz), 2) = &
732 : rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 2) + &
733 : fx*as(coset(ax - 1, ay, az), coset(0, by, bz), 2) + &
734 74370 : ac1(2)*s(coset(ax, ay, az), coset(0, by, bz))
735 : as(coset(ax, ay, az), coset(1, by, bz), 3) = &
736 : rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 3) + &
737 : fx*as(coset(ax - 1, ay, az), coset(0, by, bz), 3) + &
738 74370 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
739 74370 : IF (az > 0) as(coset(ax, ay, az), coset(1, by, bz), 2) = &
740 : as(coset(ax, ay, az), coset(1, by, bz), 2) + &
741 4726 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
742 74370 : IF (ay > 0) as(coset(ax, ay, az), coset(1, by, bz), 3) = &
743 : as(coset(ax, ay, az), coset(1, by, bz), 3) - &
744 40231 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
745 : END DO
746 74370 : DO bx = 2, lb
747 38865 : f3 = f2*REAL(bx - 1, dp)
748 117435 : DO by = 0, lb - bx
749 43065 : bz = lb - bx - by
750 : s(coset(ax, ay, az), coset(bx, by, bz)) = &
751 : rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) + &
752 : fx*s(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + &
753 43065 : f3*s(coset(ax, ay, az), coset(bx - 2, by, bz))
754 : as(coset(ax, ay, az), coset(bx, by, bz), 1) = &
755 : rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 1) + &
756 : fx*as(coset(ax - 1, ay, az), coset(bx - 1, by, bz), 1) + &
757 : f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 1) + &
758 43065 : ac1(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
759 : as(coset(ax, ay, az), coset(bx, by, bz), 2) = &
760 : rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 2) + &
761 : fx*as(coset(ax - 1, ay, az), coset(bx - 1, by, bz), 2) + &
762 : f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 2) + &
763 43065 : ac1(2)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
764 : as(coset(ax, ay, az), coset(bx, by, bz), 3) = &
765 : rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 3) + &
766 : fx*as(coset(ax - 1, ay, az), coset(bx - 1, by, bz), 3) + &
767 : f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 3) + &
768 43065 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
769 43065 : IF (az > 0) as(coset(ax, ay, az), coset(bx, by, bz), 2) = &
770 : as(coset(ax, ay, az), coset(bx, by, bz), 2) + &
771 2363 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
772 43065 : IF (ay > 0) as(coset(ax, ay, az), coset(bx, by, bz), 3) = &
773 : as(coset(ax, ay, az), coset(bx, by, bz), 3) - &
774 41228 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz))
775 : END DO
776 : END DO
777 : END IF
778 :
779 : END DO
780 : END DO
781 :
782 : END DO
783 :
784 : END IF
785 :
786 : ELSE
787 :
788 214875 : IF (lb_max > 0) THEN
789 :
790 : ! *** Vertical recurrence steps: [s|L|s] -> [s|L|b] ***
791 :
792 290092 : rbp(:) = (f1 - 1.0_dp)*rab(:)
793 :
794 : ! *** [s|p] = (Pi - Bi)*[s|s] ***
795 : ! *** [s|L|p] = (Pi - Bi)*[s|L|s] + xa/(xa+xb)*(AC x 1i)*[s|s] ***
796 :
797 72523 : s(1, 2) = rbp(1)*s(1, 1)
798 72523 : s(1, 3) = rbp(2)*s(1, 1)
799 72523 : s(1, 4) = rbp(3)*s(1, 1)
800 72523 : as(1, 2, 1) = rbp(1)*as(1, 1, 1) + ac1(1)*s(1, 1)
801 72523 : as(1, 2, 2) = rbp(1)*as(1, 1, 2) + ac1(2)*s(1, 1)
802 72523 : as(1, 2, 3) = rbp(1)*as(1, 1, 3) + ac1(3)*s(1, 1)
803 72523 : as(1, 3, 1) = rbp(2)*as(1, 1, 1) + ac2(1)*s(1, 1)
804 72523 : as(1, 3, 2) = rbp(2)*as(1, 1, 2) + ac2(2)*s(1, 1)
805 72523 : as(1, 3, 3) = rbp(2)*as(1, 1, 3) + ac2(3)*s(1, 1)
806 72523 : as(1, 4, 1) = rbp(3)*as(1, 1, 1) + ac3(1)*s(1, 1)
807 72523 : as(1, 4, 2) = rbp(3)*as(1, 1, 2) + ac3(2)*s(1, 1)
808 72523 : as(1, 4, 3) = rbp(3)*as(1, 1, 3) + ac3(3)*s(1, 1)
809 :
810 : ! *** [s|b] = (Pi - Bi)*[s|b-1i] + f2*Ni(b-1i)*[s|b-2i] ***
811 : ! *** [s|L|b] = (Pi - Bi)*[s|L|b-1i] + f2*Ni(b-1i)*[s|L|b-2i] ***
812 : ! *** + xa/(xa+xb)*(AC x 1i)*[s|s-1i] ***
813 :
814 80609 : DO lb = 2, lb_max
815 :
816 : ! *** Increase the angular momentum component z of function b ***
817 :
818 : s(1, coset(0, 0, lb)) = rbp(3)*s(1, coset(0, 0, lb - 1)) + &
819 8086 : f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2))
820 : as(1, coset(0, 0, lb), 1) = rbp(3)*as(1, coset(0, 0, lb - 1), 1) + &
821 : f2*REAL(lb - 1, dp)*as(1, coset(0, 0, lb - 2), 1) + &
822 8086 : ac3(1)*s(1, coset(0, 0, lb - 1))
823 : as(1, coset(0, 0, lb), 2) = rbp(3)*as(1, coset(0, 0, lb - 1), 2) + &
824 : f2*REAL(lb - 1, dp)*as(1, coset(0, 0, lb - 2), 2) + &
825 8086 : ac3(2)*s(1, coset(0, 0, lb - 1))
826 : as(1, coset(0, 0, lb), 3) = rbp(3)*as(1, coset(0, 0, lb - 1), 3) + &
827 : f2*REAL(lb - 1, dp)*as(1, coset(0, 0, lb - 2), 3) + &
828 8086 : ac3(3)*s(1, coset(0, 0, lb - 1))
829 :
830 : ! *** Increase the angular momentum component y of function b ***
831 :
832 8086 : bz = lb - 1
833 8086 : s(1, coset(0, 1, bz)) = rbp(2)*s(1, coset(0, 0, bz))
834 : as(1, coset(0, 1, bz), 1) = rbp(2)*as(1, coset(0, 0, bz), 1) + &
835 8086 : ac2(1)*s(1, coset(0, 0, bz))
836 : as(1, coset(0, 1, bz), 2) = rbp(2)*as(1, coset(0, 0, bz), 2) + &
837 8086 : ac2(2)*s(1, coset(0, 0, bz))
838 : as(1, coset(0, 1, bz), 3) = rbp(2)*as(1, coset(0, 0, bz), 3) + &
839 8086 : ac2(3)*s(1, coset(0, 0, bz))
840 :
841 18492 : DO by = 2, lb
842 10406 : bz = lb - by
843 : s(1, coset(0, by, bz)) = rbp(2)*s(1, coset(0, by - 1, bz)) + &
844 10406 : f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz))
845 : as(1, coset(0, by, bz), 1) = rbp(2)*as(1, coset(0, by - 1, bz), 1) + &
846 : f2*REAL(by - 1, dp)*as(1, coset(0, by - 2, bz), 1) + &
847 10406 : ac2(1)*s(1, coset(0, by - 1, bz))
848 : as(1, coset(0, by, bz), 2) = rbp(2)*as(1, coset(0, by - 1, bz), 2) + &
849 : f2*REAL(by - 1, dp)*as(1, coset(0, by - 2, bz), 2) + &
850 10406 : ac2(2)*s(1, coset(0, by - 1, bz))
851 : as(1, coset(0, by, bz), 3) = rbp(2)*as(1, coset(0, by - 1, bz), 3) + &
852 : f2*REAL(by - 1, dp)*as(1, coset(0, by - 2, bz), 3) + &
853 18492 : ac2(3)*s(1, coset(0, by - 1, bz))
854 : END DO
855 :
856 : ! *** Increase the angular momentum component x of function b ***
857 :
858 26578 : DO by = 0, lb - 1
859 18492 : bz = lb - 1 - by
860 18492 : s(1, coset(1, by, bz)) = rbp(1)*s(1, coset(0, by, bz))
861 : as(1, coset(1, by, bz), 1) = rbp(1)*as(1, coset(0, by, bz), 1) + &
862 18492 : ac1(1)*s(1, coset(0, by, bz))
863 : as(1, coset(1, by, bz), 2) = rbp(1)*as(1, coset(0, by, bz), 2) + &
864 18492 : ac1(2)*s(1, coset(0, by, bz))
865 : as(1, coset(1, by, bz), 3) = rbp(1)*as(1, coset(0, by, bz), 3) + &
866 26578 : ac1(3)*s(1, coset(0, by, bz))
867 : END DO
868 :
869 91015 : DO bx = 2, lb
870 10406 : f3 = f2*REAL(bx - 1, dp)
871 31798 : DO by = 0, lb - bx
872 13306 : bz = lb - bx - by
873 : s(1, coset(bx, by, bz)) = rbp(1)*s(1, coset(bx - 1, by, bz)) + &
874 13306 : f3*s(1, coset(bx - 2, by, bz))
875 : as(1, coset(bx, by, bz), 1) = rbp(1)*as(1, coset(bx - 1, by, bz), 1) + &
876 : f3*as(1, coset(bx - 2, by, bz), 1) + &
877 13306 : ac1(1)*s(1, coset(bx - 1, by, bz))
878 : as(1, coset(bx, by, bz), 2) = rbp(1)*as(1, coset(bx - 1, by, bz), 2) + &
879 : f3*as(1, coset(bx - 2, by, bz), 2) + &
880 13306 : ac1(2)*s(1, coset(bx - 1, by, bz))
881 : as(1, coset(bx, by, bz), 3) = rbp(1)*as(1, coset(bx - 1, by, bz), 3) + &
882 : f3*as(1, coset(bx - 2, by, bz), 3) + &
883 23712 : ac1(3)*s(1, coset(bx - 1, by, bz))
884 : END DO
885 : END DO
886 :
887 : END DO
888 :
889 : END IF
890 :
891 : END IF
892 :
893 2332695 : DO j = 1, ncoset(lb_max)
894 8749135 : DO i = 1, ncoset(la_max)
895 6416440 : angab(na + i, nb + j, 1) = as(i, j, 1)
896 6416440 : angab(na + i, nb + j, 2) = as(i, j, 2)
897 8205959 : angab(na + i, nb + j, 3) = as(i, j, 3)
898 : END DO
899 : END DO
900 :
901 1058321 : nb = nb + ncoset(lb_max)
902 :
903 : END DO
904 :
905 726666 : na = na + ncoset(la_max)
906 :
907 : END DO
908 :
909 211521 : END SUBROUTINE angmom
910 :
911 : END MODULE ai_angmom
|
