Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of 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 9410 : SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, &
51 1882 : 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 3764 : DIMENSION(ncoset(la_max), ncoset(lb_max)) :: s
68 : REAL(KIND=dp), &
69 1882 : 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 7528 : rab = rbc - rac
87 1882 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
88 1882 : dab = SQRT(rab2)
89 :
90 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
91 6230596 : angab = 0.0_dp
92 128526 : s = 0.0_dp
93 387460 : as = 0.0_dp
94 :
95 1882 : na = 0
96 :
97 7858 : DO ipgf = 1, npgfa
98 :
99 5976 : nb = 0
100 :
101 15471 : DO jpgf = 1, npgfb
102 :
103 568811 : s = 0.0_dp
104 1715928 : as = 0.0_dp
105 :
106 : ! *** Screening ***
107 :
108 9495 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109 48289 : DO j = nb + 1, nb + ncoset(lb_max)
110 227265 : DO i = na + 1, na + ncoset(la_max)
111 178976 : angab(i, j, 1) = 0.0_dp
112 178976 : angab(i, j, 2) = 0.0_dp
113 225100 : angab(i, j, 3) = 0.0_dp
114 : END DO
115 : END DO
116 2165 : nb = nb + ncoset(lb_max)
117 2165 : CYCLE
118 : END IF
119 :
120 : ! *** Calculate some prefactors ***
121 :
122 7330 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
123 :
124 7330 : f0 = (pi*zetp)**1.5_dp
125 7330 : f1 = zetb(jpgf)*zetp
126 7330 : f2 = 0.5_dp*zetp
127 : !
128 7330 : bc1(1) = 0._dp
129 7330 : bc1(2) = -f1*rbc(3)
130 7330 : bc1(3) = f1*rbc(2)
131 : !
132 7330 : bc2(1) = f1*rbc(3)
133 7330 : bc2(2) = 0._dp
134 7330 : bc2(3) = -f1*rbc(1)
135 : !
136 7330 : bc3(1) = -f1*rbc(2)
137 7330 : bc3(2) = f1*rbc(1)
138 7330 : bc3(3) = 0._dp
139 : !
140 7330 : ac1(1) = 0._dp
141 7330 : ac1(2) = zeta(ipgf)*zetp*rac(3)
142 7330 : ac1(3) = -zeta(ipgf)*zetp*rac(2)
143 : !
144 7330 : ac2(1) = -zeta(ipgf)*zetp*rac(3)
145 7330 : ac2(2) = 0._dp
146 7330 : ac2(3) = zeta(ipgf)*zetp*rac(1)
147 : !
148 7330 : ac3(1) = zeta(ipgf)*zetp*rac(2)
149 7330 : ac3(2) = -zeta(ipgf)*zetp*rac(1)
150 7330 : 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 7330 : s(1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2)
156 7330 : as(1, 1, 1) = 2._dp*zeta(ipgf)*f1*(rac(2)*rbc(3) - rac(3)*rbc(2))*s(1, 1)
157 7330 : as(1, 1, 2) = 2._dp*zeta(ipgf)*f1*(rac(3)*rbc(1) - rac(1)*rbc(3))*s(1, 1)
158 7330 : 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 7330 : IF (la_max > 0) THEN
163 :
164 : ! *** Vertical recurrence steps: [s|L|s] -> [a|L|s] ***
165 :
166 20360 : 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 5090 : s(2, 1) = rap(1)*s(1, 1)
172 5090 : s(3, 1) = rap(2)*s(1, 1)
173 5090 : s(4, 1) = rap(3)*s(1, 1)
174 5090 : as(2, 1, 1) = rap(1)*as(1, 1, 1) + bc1(1)*s(1, 1)
175 5090 : as(2, 1, 2) = rap(1)*as(1, 1, 2) + bc1(2)*s(1, 1)
176 5090 : as(2, 1, 3) = rap(1)*as(1, 1, 3) + bc1(3)*s(1, 1)
177 5090 : as(3, 1, 1) = rap(2)*as(1, 1, 1) + bc2(1)*s(1, 1)
178 5090 : as(3, 1, 2) = rap(2)*as(1, 1, 2) + bc2(2)*s(1, 1)
179 5090 : as(3, 1, 3) = rap(2)*as(1, 1, 3) + bc2(3)*s(1, 1)
180 5090 : as(4, 1, 1) = rap(3)*as(1, 1, 1) + bc3(1)*s(1, 1)
181 5090 : as(4, 1, 2) = rap(3)*as(1, 1, 2) + bc3(2)*s(1, 1)
182 5090 : 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 6000 : 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 6910 : 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 5090 : IF (lb_max > 0) THEN
266 :
267 51668 : DO j = 2, ncoset(lb_max)
268 268920 : DO i = 1, ncoset(la_max)
269 217252 : s(i, j) = 0.0_dp
270 217252 : as(i, j, 1) = 0.0_dp
271 217252 : as(i, j, 2) = 0.0_dp
272 264482 : as(i, j, 3) = 0.0_dp
273 : END DO
274 : END DO
275 :
276 : ! *** Horizontal recurrence steps ***
277 :
278 17752 : 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 4438 : IF (lb_max == 1) THEN
284 1982 : la_start = la_min
285 : ELSE
286 2456 : la_start = MAX(0, la_min - 1)
287 : END IF
288 :
289 9355 : DO la = la_start, la_max - 1
290 15001 : DO ax = 0, la
291 16938 : DO ay = 0, la - ax
292 6375 : az = la - ax - ay
293 : s(coset(ax, ay, az), 2) = s(coset(ax + 1, ay, az), 1) - &
294 6375 : rab(1)*s(coset(ax, ay, az), 1)
295 : s(coset(ax, ay, az), 3) = s(coset(ax, ay + 1, az), 1) - &
296 6375 : rab(2)*s(coset(ax, ay, az), 1)
297 : s(coset(ax, ay, az), 4) = s(coset(ax, ay, az + 1), 1) - &
298 6375 : 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 6375 : 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 6375 : - 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 6375 : + 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 6375 : + 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 6375 : 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 6375 : - 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 6375 : - 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 6375 : + 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 12021 : 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 14140 : DO ax = 0, la_max
335 9702 : fx = f2*REAL(ax, dp)
336 29932 : DO ay = 0, la_max - ax
337 15792 : fy = f2*REAL(ay, dp)
338 15792 : az = la_max - ax - ay
339 15792 : fz = f2*REAL(az, dp)
340 15792 : IF (ax == 0) THEN
341 9702 : 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 9702 : 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 9702 : 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 9702 : 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 6090 : 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 6090 : 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 6090 : 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 6090 : ac1(3)*s(coset(ax, ay, az), 1)
360 : END IF
361 15792 : IF (az > 0) as(coset(ax, ay, az), 2, 2) = as(coset(ax, ay, az), 2, 2) + &
362 6090 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
363 15792 : IF (ay > 0) as(coset(ax, ay, az), 2, 3) = as(coset(ax, ay, az), 2, 3) - &
364 6090 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
365 15792 : IF (ay == 0) THEN
366 9702 : 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 9702 : 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 9702 : 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 9702 : 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 6090 : 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 6090 : 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 6090 : 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 6090 : ac2(3)*s(coset(ax, ay, az), 1)
385 : END IF
386 15792 : IF (az > 0) as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay, az), 3, 1) - &
387 6090 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
388 15792 : IF (ax > 0) as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay, az), 3, 3) + &
389 6090 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1)
390 15792 : IF (az == 0) THEN
391 9702 : 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 9702 : 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 9702 : 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 9702 : 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 6090 : 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 6090 : 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 6090 : 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 6090 : ac3(3)*s(coset(ax, ay, az), 1)
410 : END IF
411 15792 : IF (ay > 0) as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az), 4, 1) + &
412 6090 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
413 15792 : IF (ax > 0) as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az), 4, 2) - &
414 15792 : 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 8538 : 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 4100 : IF (lb == lb_max) THEN
428 2456 : la_start = la_min
429 : ELSE
430 1644 : la_start = MAX(0, la_min - 1)
431 : END IF
432 :
433 8373 : DO la = la_start, la_max - 1
434 13001 : DO ax = 0, la
435 13884 : DO ay = 0, la - ax
436 4983 : 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 4983 : 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 4983 : + 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 4983 : - 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 4983 : 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 16997 : DO by = 1, lb
460 12014 : 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 12014 : 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 12014 : - 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 12014 : 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 16997 : + 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 21625 : DO bx = 1, lb
482 38590 : DO by = 0, lb - bx
483 21593 : 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 21593 : 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 21593 : 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 21593 : + 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 33607 : - 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 17112 : DO ax = 0, la_max
517 8574 : fx = f2*REAL(ax, dp)
518 26096 : DO ay = 0, la_max - ax
519 13422 : fy = f2*REAL(ay, dp)
520 13422 : az = la_max - ax - ay
521 13422 : fz = f2*REAL(az, dp)
522 :
523 : ! *** Increase the angular momentum component z of function b ***
524 :
525 13422 : f3 = f2*REAL(lb - 1, dp)
526 :
527 13422 : 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 8574 : 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 8574 : 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 8574 : 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 8574 : 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 4848 : 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 4848 : 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 4848 : 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 4848 : ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
563 : END IF
564 13422 : 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 4848 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, 0, lb - 1))
567 13422 : 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 4848 : 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 13422 : IF (ay == 0) THEN
574 8574 : bz = lb - 1
575 : s(coset(ax, ay, az), coset(0, 1, bz)) = &
576 8574 : 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 8574 : 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 8574 : 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 8574 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
586 8574 : 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 4474 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
589 8574 : 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 4474 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
592 21532 : DO by = 2, lb
593 12958 : bz = lb - by
594 12958 : 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 12958 : 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 12958 : 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 12958 : 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 12958 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
610 12958 : 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 6666 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
613 12958 : 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 15240 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz))
616 : END DO
617 : ELSE
618 4848 : 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 4848 : 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 4848 : 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 4848 : 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 4848 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
634 4848 : 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 4848 : 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 11888 : DO by = 2, lb
641 7040 : bz = lb - by
642 7040 : 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 7040 : 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 7040 : 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 7040 : 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 7040 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
662 7040 : 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 7040 : 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 5222 : 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 21996 : IF (ax == 0) THEN
674 30106 : DO by = 0, lb - 1
675 21532 : bz = lb - 1 - by
676 : s(coset(ax, ay, az), coset(1, by, bz)) = &
677 21532 : 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 21532 : 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 21532 : 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 21532 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
687 21532 : 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 11140 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
690 21532 : 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 19714 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
693 : END DO
694 21532 : DO bx = 2, lb
695 12958 : f3 = f2*REAL(bx - 1, dp)
696 39970 : DO by = 0, lb - bx
697 18438 : 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 18438 : 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 18438 : 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 18438 : 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 18438 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
713 18438 : 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 9406 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
716 18438 : 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 22364 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz))
719 : END DO
720 : END DO
721 : ELSE
722 16736 : DO by = 0, lb - 1
723 11888 : 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 11888 : 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 11888 : 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 11888 : 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 11888 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
739 11888 : 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 11888 : 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 5596 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
745 : END DO
746 11888 : DO bx = 2, lb
747 7040 : f3 = f2*REAL(bx - 1, dp)
748 21668 : DO by = 0, lb - bx
749 9780 : 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 9780 : 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 9780 : 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 9780 : 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 9780 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
769 9780 : 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 9780 : 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 7414 : 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 2240 : IF (lb_max > 0) THEN
789 :
790 : ! *** Vertical recurrence steps: [s|L|s] -> [s|L|b] ***
791 :
792 5328 : 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 1332 : s(1, 2) = rbp(1)*s(1, 1)
798 1332 : s(1, 3) = rbp(2)*s(1, 1)
799 1332 : s(1, 4) = rbp(3)*s(1, 1)
800 1332 : as(1, 2, 1) = rbp(1)*as(1, 1, 1) + ac1(1)*s(1, 1)
801 1332 : as(1, 2, 2) = rbp(1)*as(1, 1, 2) + ac1(2)*s(1, 1)
802 1332 : as(1, 2, 3) = rbp(1)*as(1, 1, 3) + ac1(3)*s(1, 1)
803 1332 : as(1, 3, 1) = rbp(2)*as(1, 1, 1) + ac2(1)*s(1, 1)
804 1332 : as(1, 3, 2) = rbp(2)*as(1, 1, 2) + ac2(2)*s(1, 1)
805 1332 : as(1, 3, 3) = rbp(2)*as(1, 1, 3) + ac2(3)*s(1, 1)
806 1332 : as(1, 4, 1) = rbp(3)*as(1, 1, 1) + ac3(1)*s(1, 1)
807 1332 : as(1, 4, 2) = rbp(3)*as(1, 1, 2) + ac3(2)*s(1, 1)
808 1332 : 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 2848 : 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 1516 : 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 1516 : 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 1516 : 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 1516 : ac3(3)*s(1, coset(0, 0, lb - 1))
829 :
830 : ! *** Increase the angular momentum component y of function b ***
831 :
832 1516 : bz = lb - 1
833 1516 : 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 1516 : 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 1516 : 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 1516 : ac2(3)*s(1, coset(0, 0, bz))
840 :
841 3992 : DO by = 2, lb
842 2476 : bz = lb - by
843 : s(1, coset(0, by, bz)) = rbp(2)*s(1, coset(0, by - 1, bz)) + &
844 2476 : 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 2476 : 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 2476 : 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 3992 : 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 5508 : DO by = 0, lb - 1
859 3992 : bz = lb - 1 - by
860 3992 : 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 3992 : 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 3992 : 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 5508 : ac1(3)*s(1, coset(0, by, bz))
867 : END DO
868 :
869 5324 : DO bx = 2, lb
870 2476 : f3 = f2*REAL(bx - 1, dp)
871 7668 : DO by = 0, lb - bx
872 3676 : bz = lb - bx - by
873 : s(1, coset(bx, by, bz)) = rbp(1)*s(1, coset(bx - 1, by, bz)) + &
874 3676 : 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 3676 : 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 3676 : 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 6152 : 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 79062 : DO j = 1, ncoset(lb_max)
894 341546 : DO i = 1, ncoset(la_max)
895 262484 : angab(na + i, nb + j, 1) = as(i, j, 1)
896 262484 : angab(na + i, nb + j, 2) = as(i, j, 2)
897 334216 : angab(na + i, nb + j, 3) = as(i, j, 3)
898 : END DO
899 : END DO
900 :
901 13306 : nb = nb + ncoset(lb_max)
902 :
903 : END DO
904 :
905 7858 : na = na + ncoset(la_max)
906 :
907 : END DO
908 :
909 1882 : END SUBROUTINE angmom
910 :
911 : END MODULE ai_angmom
|
