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 10760 : SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, &
51 2152 : 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 4304 : DIMENSION(ncoset(la_max), ncoset(lb_max)) :: s
68 : REAL(KIND=dp), &
69 2152 : 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 8608 : rab = rbc - rac
87 2152 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
88 2152 : dab = SQRT(rab2)
89 :
90 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
91 7252276 : angab = 0.0_dp
92 147906 : s = 0.0_dp
93 445870 : as = 0.0_dp
94 :
95 2152 : na = 0
96 :
97 8758 : DO ipgf = 1, npgfa
98 :
99 6606 : nb = 0
100 :
101 16731 : DO jpgf = 1, npgfb
102 :
103 621311 : s = 0.0_dp
104 1874058 : as = 0.0_dp
105 :
106 : ! *** Screening ***
107 :
108 10125 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109 52913 : DO j = nb + 1, nb + ncoset(lb_max)
110 245669 : DO i = na + 1, na + ncoset(la_max)
111 192756 : angab(i, j, 1) = 0.0_dp
112 192756 : angab(i, j, 2) = 0.0_dp
113 243300 : angab(i, j, 3) = 0.0_dp
114 : END DO
115 : END DO
116 2369 : nb = nb + ncoset(lb_max)
117 2369 : CYCLE
118 : END IF
119 :
120 : ! *** Calculate some prefactors ***
121 :
122 7756 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
123 :
124 7756 : f0 = (pi*zetp)**1.5_dp
125 7756 : f1 = zetb(jpgf)*zetp
126 7756 : f2 = 0.5_dp*zetp
127 : !
128 7756 : bc1(1) = 0._dp
129 7756 : bc1(2) = -f1*rbc(3)
130 7756 : bc1(3) = f1*rbc(2)
131 : !
132 7756 : bc2(1) = f1*rbc(3)
133 7756 : bc2(2) = 0._dp
134 7756 : bc2(3) = -f1*rbc(1)
135 : !
136 7756 : bc3(1) = -f1*rbc(2)
137 7756 : bc3(2) = f1*rbc(1)
138 7756 : bc3(3) = 0._dp
139 : !
140 7756 : ac1(1) = 0._dp
141 7756 : ac1(2) = zeta(ipgf)*zetp*rac(3)
142 7756 : ac1(3) = -zeta(ipgf)*zetp*rac(2)
143 : !
144 7756 : ac2(1) = -zeta(ipgf)*zetp*rac(3)
145 7756 : ac2(2) = 0._dp
146 7756 : ac2(3) = zeta(ipgf)*zetp*rac(1)
147 : !
148 7756 : ac3(1) = zeta(ipgf)*zetp*rac(2)
149 7756 : ac3(2) = -zeta(ipgf)*zetp*rac(1)
150 7756 : 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 7756 : s(1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2)
156 7756 : as(1, 1, 1) = 2._dp*zeta(ipgf)*f1*(rac(2)*rbc(3) - rac(3)*rbc(2))*s(1, 1)
157 7756 : as(1, 1, 2) = 2._dp*zeta(ipgf)*f1*(rac(3)*rbc(1) - rac(1)*rbc(3))*s(1, 1)
158 7756 : 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 7756 : IF (la_max > 0) THEN
163 :
164 : ! *** Vertical recurrence steps: [s|L|s] -> [a|L|s] ***
165 :
166 21296 : 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 5324 : s(2, 1) = rap(1)*s(1, 1)
172 5324 : s(3, 1) = rap(2)*s(1, 1)
173 5324 : s(4, 1) = rap(3)*s(1, 1)
174 5324 : as(2, 1, 1) = rap(1)*as(1, 1, 1) + bc1(1)*s(1, 1)
175 5324 : as(2, 1, 2) = rap(1)*as(1, 1, 2) + bc1(2)*s(1, 1)
176 5324 : as(2, 1, 3) = rap(1)*as(1, 1, 3) + bc1(3)*s(1, 1)
177 5324 : as(3, 1, 1) = rap(2)*as(1, 1, 1) + bc2(1)*s(1, 1)
178 5324 : as(3, 1, 2) = rap(2)*as(1, 1, 2) + bc2(2)*s(1, 1)
179 5324 : as(3, 1, 3) = rap(2)*as(1, 1, 3) + bc2(3)*s(1, 1)
180 5324 : as(4, 1, 1) = rap(3)*as(1, 1, 1) + bc3(1)*s(1, 1)
181 5324 : as(4, 1, 2) = rap(3)*as(1, 1, 2) + bc3(2)*s(1, 1)
182 5324 : 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 6234 : 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 7144 : 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 5324 : IF (lb_max > 0) THEN
266 :
267 56738 : DO j = 2, ncoset(lb_max)
268 293334 : DO i = 1, ncoset(la_max)
269 236596 : s(i, j) = 0.0_dp
270 236596 : as(i, j, 1) = 0.0_dp
271 236596 : as(i, j, 2) = 0.0_dp
272 288662 : as(i, j, 3) = 0.0_dp
273 : END DO
274 : END DO
275 :
276 : ! *** Horizontal recurrence steps ***
277 :
278 18688 : 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 4672 : IF (lb_max == 1) THEN
284 1982 : la_start = la_min
285 : ELSE
286 2690 : la_start = MAX(0, la_min - 1)
287 : END IF
288 :
289 9823 : DO la = la_start, la_max - 1
290 15703 : DO ax = 0, la
291 17640 : DO ay = 0, la - ax
292 6609 : az = la - ax - ay
293 : s(coset(ax, ay, az), 2) = s(coset(ax + 1, ay, az), 1) - &
294 6609 : rab(1)*s(coset(ax, ay, az), 1)
295 : s(coset(ax, ay, az), 3) = s(coset(ax, ay + 1, az), 1) - &
296 6609 : rab(2)*s(coset(ax, ay, az), 1)
297 : s(coset(ax, ay, az), 4) = s(coset(ax, ay, az + 1), 1) - &
298 6609 : 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 6609 : 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 6609 : - 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 6609 : + 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 6609 : + 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 6609 : 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 6609 : - 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 6609 : - 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 6609 : + 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 12489 : 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 14842 : DO ax = 0, la_max
335 10170 : fx = f2*REAL(ax, dp)
336 31336 : DO ay = 0, la_max - ax
337 16494 : fy = f2*REAL(ay, dp)
338 16494 : az = la_max - ax - ay
339 16494 : fz = f2*REAL(az, dp)
340 16494 : IF (ax == 0) THEN
341 10170 : 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 10170 : 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 10170 : 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 10170 : 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 6324 : 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 6324 : 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 6324 : 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 6324 : ac1(3)*s(coset(ax, ay, az), 1)
360 : END IF
361 16494 : IF (az > 0) as(coset(ax, ay, az), 2, 2) = as(coset(ax, ay, az), 2, 2) + &
362 6324 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
363 16494 : IF (ay > 0) as(coset(ax, ay, az), 2, 3) = as(coset(ax, ay, az), 2, 3) - &
364 6324 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
365 16494 : IF (ay == 0) THEN
366 10170 : 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 10170 : 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 10170 : 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 10170 : 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 6324 : 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 6324 : 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 6324 : 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 6324 : ac2(3)*s(coset(ax, ay, az), 1)
385 : END IF
386 16494 : IF (az > 0) as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay, az), 3, 1) - &
387 6324 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
388 16494 : IF (ax > 0) as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay, az), 3, 3) + &
389 6324 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1)
390 16494 : IF (az == 0) THEN
391 10170 : 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 10170 : 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 10170 : 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 10170 : 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 6324 : 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 6324 : 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 6324 : 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 6324 : ac3(3)*s(coset(ax, ay, az), 1)
410 : END IF
411 16494 : IF (ay > 0) as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az), 4, 1) + &
412 6324 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
413 16494 : IF (ax > 0) as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az), 4, 2) - &
414 16494 : 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 9240 : 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 4568 : IF (lb == lb_max) THEN
428 2690 : la_start = la_min
429 : ELSE
430 1878 : la_start = MAX(0, la_min - 1)
431 : END IF
432 :
433 9261 : DO la = la_start, la_max - 1
434 14309 : DO ax = 0, la
435 15144 : DO ay = 0, la - ax
436 5403 : 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 5403 : 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 5403 : + 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 5403 : - 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 5403 : 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 18521 : DO by = 1, lb
460 13118 : 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 13118 : 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 13118 : - 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 13118 : 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 18521 : + 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 23569 : DO bx = 1, lb
482 42228 : DO by = 0, lb - bx
483 23707 : 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 23707 : 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 23707 : 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 23707 : + 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 36825 : - 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 18750 : DO ax = 0, la_max
517 9510 : fx = f2*REAL(ax, dp)
518 28904 : DO ay = 0, la_max - ax
519 14826 : fy = f2*REAL(ay, dp)
520 14826 : az = la_max - ax - ay
521 14826 : fz = f2*REAL(az, dp)
522 :
523 : ! *** Increase the angular momentum component z of function b ***
524 :
525 14826 : f3 = f2*REAL(lb - 1, dp)
526 :
527 14826 : 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 9510 : 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 9510 : 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 9510 : 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 9510 : 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 5316 : 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 5316 : 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 5316 : 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 5316 : ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
563 : END IF
564 14826 : 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 5316 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, 0, lb - 1))
567 14826 : 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 5316 : 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 14826 : IF (ay == 0) THEN
574 9510 : bz = lb - 1
575 : s(coset(ax, ay, az), coset(0, 1, bz)) = &
576 9510 : 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 9510 : 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 9510 : 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 9510 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
586 9510 : 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 4942 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
589 9510 : 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 4942 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
592 24028 : DO by = 2, lb
593 14518 : bz = lb - by
594 14518 : 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 14518 : 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 14518 : 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 14518 : 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 14518 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
610 14518 : 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 7446 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
613 14518 : 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 16956 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz))
616 : END DO
617 : ELSE
618 5316 : 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 5316 : 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 5316 : 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 5316 : 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 5316 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
634 5316 : 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 5316 : 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 13136 : DO by = 2, lb
641 7820 : bz = lb - by
642 7820 : 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 7820 : 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 7820 : 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 7820 : 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 7820 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
662 7820 : 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 7820 : 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 5690 : 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 24336 : IF (ax == 0) THEN
674 33538 : DO by = 0, lb - 1
675 24028 : bz = lb - 1 - by
676 : s(coset(ax, ay, az), coset(1, by, bz)) = &
677 24028 : 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 24028 : 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 24028 : 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 24028 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
687 24028 : 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 12388 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
690 24028 : 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 21898 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
693 : END DO
694 24028 : DO bx = 2, lb
695 14518 : f3 = f2*REAL(bx - 1, dp)
696 44806 : DO by = 0, lb - bx
697 20778 : 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 20778 : 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 20778 : 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 20778 : 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 20778 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
713 20778 : 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 10576 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
716 20778 : 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 25094 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz))
719 : END DO
720 : END DO
721 : ELSE
722 18452 : DO by = 0, lb - 1
723 13136 : 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 13136 : 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 13136 : 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 13136 : 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 13136 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
739 13136 : 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 13136 : 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 6064 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
745 : END DO
746 13136 : DO bx = 2, lb
747 7820 : f3 = f2*REAL(bx - 1, dp)
748 24086 : DO by = 0, lb - bx
749 10950 : 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 10950 : 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 10950 : 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 10950 : 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 10950 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
769 10950 : 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 10950 : 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 8194 : 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 2432 : IF (lb_max > 0) THEN
789 :
790 : ! *** Vertical recurrence steps: [s|L|s] -> [s|L|b] ***
791 :
792 6096 : 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 1524 : s(1, 2) = rbp(1)*s(1, 1)
798 1524 : s(1, 3) = rbp(2)*s(1, 1)
799 1524 : s(1, 4) = rbp(3)*s(1, 1)
800 1524 : as(1, 2, 1) = rbp(1)*as(1, 1, 1) + ac1(1)*s(1, 1)
801 1524 : as(1, 2, 2) = rbp(1)*as(1, 1, 2) + ac1(2)*s(1, 1)
802 1524 : as(1, 2, 3) = rbp(1)*as(1, 1, 3) + ac1(3)*s(1, 1)
803 1524 : as(1, 3, 1) = rbp(2)*as(1, 1, 1) + ac2(1)*s(1, 1)
804 1524 : as(1, 3, 2) = rbp(2)*as(1, 1, 2) + ac2(2)*s(1, 1)
805 1524 : as(1, 3, 3) = rbp(2)*as(1, 1, 3) + ac2(3)*s(1, 1)
806 1524 : as(1, 4, 1) = rbp(3)*as(1, 1, 1) + ac3(1)*s(1, 1)
807 1524 : as(1, 4, 2) = rbp(3)*as(1, 1, 2) + ac3(2)*s(1, 1)
808 1524 : 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 3424 : 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 1900 : 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 1900 : 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 1900 : 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 1900 : ac3(3)*s(1, coset(0, 0, lb - 1))
829 :
830 : ! *** Increase the angular momentum component y of function b ***
831 :
832 1900 : bz = lb - 1
833 1900 : 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 1900 : 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 1900 : 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 1900 : ac2(3)*s(1, coset(0, 0, bz))
840 :
841 5016 : DO by = 2, lb
842 3116 : bz = lb - by
843 : s(1, coset(0, by, bz)) = rbp(2)*s(1, coset(0, by - 1, bz)) + &
844 3116 : 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 3116 : 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 3116 : 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 5016 : 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 6916 : DO by = 0, lb - 1
859 5016 : bz = lb - 1 - by
860 5016 : 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 5016 : 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 5016 : 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 6916 : ac1(3)*s(1, coset(0, by, bz))
867 : END DO
868 :
869 6540 : DO bx = 2, lb
870 3116 : f3 = f2*REAL(bx - 1, dp)
871 9652 : DO by = 0, lb - bx
872 4636 : bz = lb - bx - by
873 : s(1, coset(bx, by, bz)) = rbp(1)*s(1, coset(bx - 1, by, bz)) + &
874 4636 : 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 4636 : 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 4636 : 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 7752 : 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 88718 : DO j = 1, ncoset(lb_max)
894 375642 : DO i = 1, ncoset(la_max)
895 286924 : angab(na + i, nb + j, 1) = as(i, j, 1)
896 286924 : angab(na + i, nb + j, 2) = as(i, j, 2)
897 367886 : angab(na + i, nb + j, 3) = as(i, j, 3)
898 : END DO
899 : END DO
900 :
901 14362 : nb = nb + ncoset(lb_max)
902 :
903 : END DO
904 :
905 8758 : na = na + ncoset(la_max)
906 :
907 : END DO
908 :
909 2152 : END SUBROUTINE angmom
910 :
911 : END MODULE ai_angmom
|
