Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 8990 : SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, &
51 1798 : 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 3596 : DIMENSION(ncoset(la_max), ncoset(lb_max)) :: s
68 : REAL(KIND=dp), &
69 1798 : 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 7192 : rab = rbc - rac
87 1798 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
88 1798 : dab = SQRT(rab2)
89 :
90 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
91 5912740 : angab = 0.0_dp
92 114272 : s = 0.0_dp
93 344614 : as = 0.0_dp
94 :
95 1798 : na = 0
96 :
97 6982 : DO ipgf = 1, npgfa
98 :
99 5184 : nb = 0
100 :
101 13887 : DO jpgf = 1, npgfb
102 :
103 473249 : s = 0.0_dp
104 1428450 : as = 0.0_dp
105 :
106 : ! *** Screening ***
107 :
108 8703 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109 35641 : DO j = nb + 1, nb + ncoset(lb_max)
110 161187 : DO i = na + 1, na + ncoset(la_max)
111 125546 : angab(i, j, 1) = 0.0_dp
112 125546 : angab(i, j, 2) = 0.0_dp
113 159580 : angab(i, j, 3) = 0.0_dp
114 : END DO
115 : END DO
116 1607 : nb = nb + ncoset(lb_max)
117 1607 : CYCLE
118 : END IF
119 :
120 : ! *** Calculate some prefactors ***
121 :
122 7096 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
123 :
124 7096 : f0 = (pi*zetp)**1.5_dp
125 7096 : f1 = zetb(jpgf)*zetp
126 7096 : f2 = 0.5_dp*zetp
127 : !
128 7096 : bc1(1) = 0._dp
129 7096 : bc1(2) = -f1*rbc(3)
130 7096 : bc1(3) = f1*rbc(2)
131 : !
132 7096 : bc2(1) = f1*rbc(3)
133 7096 : bc2(2) = 0._dp
134 7096 : bc2(3) = -f1*rbc(1)
135 : !
136 7096 : bc3(1) = -f1*rbc(2)
137 7096 : bc3(2) = f1*rbc(1)
138 7096 : bc3(3) = 0._dp
139 : !
140 7096 : ac1(1) = 0._dp
141 7096 : ac1(2) = zeta(ipgf)*zetp*rac(3)
142 7096 : ac1(3) = -zeta(ipgf)*zetp*rac(2)
143 : !
144 7096 : ac2(1) = -zeta(ipgf)*zetp*rac(3)
145 7096 : ac2(2) = 0._dp
146 7096 : ac2(3) = zeta(ipgf)*zetp*rac(1)
147 : !
148 7096 : ac3(1) = zeta(ipgf)*zetp*rac(2)
149 7096 : ac3(2) = -zeta(ipgf)*zetp*rac(1)
150 7096 : 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 7096 : s(1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2)
156 7096 : as(1, 1, 1) = 2._dp*zeta(ipgf)*f1*(rac(2)*rbc(3) - rac(3)*rbc(2))*s(1, 1)
157 7096 : as(1, 1, 2) = 2._dp*zeta(ipgf)*f1*(rac(3)*rbc(1) - rac(1)*rbc(3))*s(1, 1)
158 7096 : 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 7096 : IF (la_max > 0) THEN
163 :
164 : ! *** Vertical recurrence steps: [s|L|s] -> [a|L|s] ***
165 :
166 19184 : 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 4796 : s(2, 1) = rap(1)*s(1, 1)
172 4796 : s(3, 1) = rap(2)*s(1, 1)
173 4796 : s(4, 1) = rap(3)*s(1, 1)
174 4796 : as(2, 1, 1) = rap(1)*as(1, 1, 1) + bc1(1)*s(1, 1)
175 4796 : as(2, 1, 2) = rap(1)*as(1, 1, 2) + bc1(2)*s(1, 1)
176 4796 : as(2, 1, 3) = rap(1)*as(1, 1, 3) + bc1(3)*s(1, 1)
177 4796 : as(3, 1, 1) = rap(2)*as(1, 1, 1) + bc2(1)*s(1, 1)
178 4796 : as(3, 1, 2) = rap(2)*as(1, 1, 2) + bc2(2)*s(1, 1)
179 4796 : as(3, 1, 3) = rap(2)*as(1, 1, 3) + bc2(3)*s(1, 1)
180 4796 : as(4, 1, 1) = rap(3)*as(1, 1, 1) + bc3(1)*s(1, 1)
181 4796 : as(4, 1, 2) = rap(3)*as(1, 1, 2) + bc3(2)*s(1, 1)
182 4796 : 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 5706 : 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 910 : 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 910 : 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 910 : 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 910 : bc3(3)*s(coset(0, 0, la - 1), 1)
203 :
204 : ! *** Increase the angular momentum component y of function a ***
205 :
206 910 : az = la - 1
207 910 : 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 910 : 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 910 : 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 910 : bc2(3)*s(coset(0, 0, az), 1)
214 :
215 1820 : DO ay = 2, la
216 910 : az = la - ay
217 : s(coset(0, ay, az), 1) = rap(2)*s(coset(0, ay - 1, az), 1) + &
218 910 : 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 910 : 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 910 : 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 1820 : 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 2730 : DO ay = 0, la - 1
233 1820 : az = la - 1 - ay
234 1820 : 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 1820 : 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 1820 : 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 2730 : bc1(3)*s(coset(0, ay, az), 1)
241 : END DO
242 :
243 6616 : DO ax = 2, la
244 910 : f3 = f2*REAL(ax - 1, dp)
245 2730 : DO ay = 0, la - ax
246 910 : az = la - ax - ay
247 : s(coset(ax, ay, az), 1) = rap(1)*s(coset(ax - 1, ay, az), 1) + &
248 910 : 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 910 : 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 910 : 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 1820 : 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 4796 : IF (lb_max > 0) THEN
266 :
267 45298 : DO j = 2, ncoset(lb_max)
268 238246 : DO i = 1, ncoset(la_max)
269 192948 : s(i, j) = 0.0_dp
270 192948 : as(i, j, 1) = 0.0_dp
271 192948 : as(i, j, 2) = 0.0_dp
272 234102 : as(i, j, 3) = 0.0_dp
273 : END DO
274 : END DO
275 :
276 : ! *** Horizontal recurrence steps ***
277 :
278 16576 : 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 4144 : IF (lb_max == 1) THEN
284 1982 : la_start = la_min
285 : ELSE
286 2162 : la_start = MAX(0, la_min - 1)
287 : END IF
288 :
289 8767 : DO la = la_start, la_max - 1
290 14119 : DO ax = 0, la
291 16056 : DO ay = 0, la - ax
292 6081 : az = la - ax - ay
293 : s(coset(ax, ay, az), 2) = s(coset(ax + 1, ay, az), 1) - &
294 6081 : rab(1)*s(coset(ax, ay, az), 1)
295 : s(coset(ax, ay, az), 3) = s(coset(ax, ay + 1, az), 1) - &
296 6081 : rab(2)*s(coset(ax, ay, az), 1)
297 : s(coset(ax, ay, az), 4) = s(coset(ax, ay, az + 1), 1) - &
298 6081 : 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 6081 : 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 6081 : - 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 6081 : + 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 6081 : + 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 6081 : 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 6081 : - 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 6081 : - 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 6081 : + 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 11433 : 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 13258 : DO ax = 0, la_max
335 9114 : fx = f2*REAL(ax, dp)
336 28168 : DO ay = 0, la_max - ax
337 14910 : fy = f2*REAL(ay, dp)
338 14910 : az = la_max - ax - ay
339 14910 : fz = f2*REAL(az, dp)
340 14910 : IF (ax == 0) THEN
341 9114 : 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 9114 : 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 9114 : 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 9114 : 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 5796 : 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 5796 : 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 5796 : 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 5796 : ac1(3)*s(coset(ax, ay, az), 1)
360 : END IF
361 14910 : IF (az > 0) as(coset(ax, ay, az), 2, 2) = as(coset(ax, ay, az), 2, 2) + &
362 5796 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
363 14910 : IF (ay > 0) as(coset(ax, ay, az), 2, 3) = as(coset(ax, ay, az), 2, 3) - &
364 5796 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
365 14910 : IF (ay == 0) THEN
366 9114 : 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 9114 : 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 9114 : 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 9114 : 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 5796 : 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 5796 : 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 5796 : 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 5796 : ac2(3)*s(coset(ax, ay, az), 1)
385 : END IF
386 14910 : IF (az > 0) as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay, az), 3, 1) - &
387 5796 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
388 14910 : IF (ax > 0) as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay, az), 3, 3) + &
389 5796 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1)
390 14910 : IF (az == 0) THEN
391 9114 : 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 9114 : 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 9114 : 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 9114 : 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 5796 : 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 5796 : 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 5796 : 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 5796 : ac3(3)*s(coset(ax, ay, az), 1)
410 : END IF
411 14910 : IF (ay > 0) as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az), 4, 1) + &
412 5796 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
413 14910 : IF (ax > 0) as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az), 4, 2) - &
414 14910 : 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 7656 : 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 3512 : IF (lb == lb_max) THEN
428 2162 : la_start = la_min
429 : ELSE
430 1350 : la_start = MAX(0, la_min - 1)
431 : END IF
432 :
433 7167 : DO la = la_start, la_max - 1
434 11177 : DO ax = 0, la
435 12030 : DO ay = 0, la - ax
436 4365 : 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 4365 : 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 4365 : + 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 4365 : - 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 4365 : 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 14721 : DO by = 1, lb
460 10356 : 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 10356 : 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 10356 : - 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 10356 : 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 14721 : + 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 18731 : DO bx = 1, lb
482 33086 : DO by = 0, lb - bx
483 18365 : 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 18365 : 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 18365 : 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 18365 : + 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 28721 : - 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 15054 : DO ax = 0, la_max
517 7398 : fx = f2*REAL(ax, dp)
518 22568 : DO ay = 0, la_max - ax
519 11658 : fy = f2*REAL(ay, dp)
520 11658 : az = la_max - ax - ay
521 11658 : fz = f2*REAL(az, dp)
522 :
523 : ! *** Increase the angular momentum component z of function b ***
524 :
525 11658 : f3 = f2*REAL(lb - 1, dp)
526 :
527 11658 : 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 7398 : 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 7398 : 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 7398 : 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 7398 : 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 4260 : 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 4260 : 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 4260 : 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 4260 : ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
563 : END IF
564 11658 : 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 4260 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, 0, lb - 1))
567 11658 : 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 4260 : 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 11658 : IF (ay == 0) THEN
574 7398 : bz = lb - 1
575 : s(coset(ax, ay, az), coset(0, 1, bz)) = &
576 7398 : 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 7398 : 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 7398 : 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 7398 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
586 7398 : 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 3886 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
589 7398 : 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 3886 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
592 18396 : DO by = 2, lb
593 10998 : bz = lb - by
594 10998 : 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 10998 : 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 10998 : 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 10998 : 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 10998 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
610 10998 : 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 5686 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
613 10998 : 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 13084 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz))
616 : END DO
617 : ELSE
618 4260 : 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 4260 : 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 4260 : 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 4260 : 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 4260 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
634 4260 : 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 374 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
637 4260 : 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 374 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
640 10320 : DO by = 2, lb
641 6060 : bz = lb - by
642 6060 : 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 6060 : 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 6060 : 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 6060 : 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 6060 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
662 6060 : 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 374 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
665 6060 : 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 4634 : 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 19056 : IF (ax == 0) THEN
674 25794 : DO by = 0, lb - 1
675 18396 : bz = lb - 1 - by
676 : s(coset(ax, ay, az), coset(1, by, bz)) = &
677 18396 : 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 18396 : 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 18396 : 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 18396 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
687 18396 : 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 9572 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
690 18396 : 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 16970 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
693 : END DO
694 18396 : DO bx = 2, lb
695 10998 : f3 = f2*REAL(bx - 1, dp)
696 33894 : DO by = 0, lb - bx
697 15498 : 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 15498 : 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 15498 : 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 15498 : 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 15498 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
713 15498 : 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 7936 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
716 15498 : 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 18934 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz))
719 : END DO
720 : END DO
721 : ELSE
722 14580 : DO by = 0, lb - 1
723 10320 : 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 10320 : 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 10320 : 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 10320 : 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 10320 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
739 10320 : 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 748 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
742 10320 : 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 5008 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
745 : END DO
746 10320 : DO bx = 2, lb
747 6060 : f3 = f2*REAL(bx - 1, dp)
748 18630 : DO by = 0, lb - bx
749 8310 : 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 8310 : 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 8310 : 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 8310 : 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 8310 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
769 8310 : 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 374 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
772 8310 : 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 6434 : 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 2300 : IF (lb_max > 0) THEN
789 :
790 : ! *** Vertical recurrence steps: [s|L|s] -> [s|L|b] ***
791 :
792 5568 : 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 1392 : s(1, 2) = rbp(1)*s(1, 1)
798 1392 : s(1, 3) = rbp(2)*s(1, 1)
799 1392 : s(1, 4) = rbp(3)*s(1, 1)
800 1392 : as(1, 2, 1) = rbp(1)*as(1, 1, 1) + ac1(1)*s(1, 1)
801 1392 : as(1, 2, 2) = rbp(1)*as(1, 1, 2) + ac1(2)*s(1, 1)
802 1392 : as(1, 2, 3) = rbp(1)*as(1, 1, 3) + ac1(3)*s(1, 1)
803 1392 : as(1, 3, 1) = rbp(2)*as(1, 1, 1) + ac2(1)*s(1, 1)
804 1392 : as(1, 3, 2) = rbp(2)*as(1, 1, 2) + ac2(2)*s(1, 1)
805 1392 : as(1, 3, 3) = rbp(2)*as(1, 1, 3) + ac2(3)*s(1, 1)
806 1392 : as(1, 4, 1) = rbp(3)*as(1, 1, 1) + ac3(1)*s(1, 1)
807 1392 : as(1, 4, 2) = rbp(3)*as(1, 1, 2) + ac3(2)*s(1, 1)
808 1392 : 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 3028 : 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 1636 : 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 1636 : 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 1636 : 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 1636 : ac3(3)*s(1, coset(0, 0, lb - 1))
829 :
830 : ! *** Increase the angular momentum component y of function b ***
831 :
832 1636 : bz = lb - 1
833 1636 : 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 1636 : 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 1636 : 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 1636 : ac2(3)*s(1, coset(0, 0, bz))
840 :
841 4312 : DO by = 2, lb
842 2676 : bz = lb - by
843 : s(1, coset(0, by, bz)) = rbp(2)*s(1, coset(0, by - 1, bz)) + &
844 2676 : 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 2676 : 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 2676 : 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 4312 : 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 5948 : DO by = 0, lb - 1
859 4312 : bz = lb - 1 - by
860 4312 : 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 4312 : 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 4312 : 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 5948 : ac1(3)*s(1, coset(0, by, bz))
867 : END DO
868 :
869 5704 : DO bx = 2, lb
870 2676 : f3 = f2*REAL(bx - 1, dp)
871 8288 : DO by = 0, lb - bx
872 3976 : bz = lb - bx - by
873 : s(1, coset(bx, by, bz)) = rbp(1)*s(1, coset(bx - 1, by, bz)) + &
874 3976 : 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 3976 : 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 3976 : 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 6652 : 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 73758 : DO j = 1, ncoset(lb_max)
894 312062 : DO i = 1, ncoset(la_max)
895 238304 : angab(na + i, nb + j, 1) = as(i, j, 1)
896 238304 : angab(na + i, nb + j, 2) = as(i, j, 2)
897 304966 : angab(na + i, nb + j, 3) = as(i, j, 3)
898 : END DO
899 : END DO
900 :
901 12280 : nb = nb + ncoset(lb_max)
902 :
903 : END DO
904 :
905 6982 : na = na + ncoset(la_max)
906 :
907 : END DO
908 :
909 1798 : END SUBROUTINE angmom
910 :
911 : END MODULE ai_angmom
|
