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 12110 : SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, &
51 2422 : 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 4844 : DIMENSION(ncoset(la_max), ncoset(lb_max)) :: s
68 : REAL(KIND=dp), &
69 2422 : 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 9688 : rab = rbc - rac
87 2422 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
88 2422 : dab = SQRT(rab2)
89 :
90 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
91 8273956 : angab = 0.0_dp
92 167286 : s = 0.0_dp
93 504280 : as = 0.0_dp
94 :
95 2422 : na = 0
96 :
97 9658 : DO ipgf = 1, npgfa
98 :
99 7236 : nb = 0
100 :
101 17991 : DO jpgf = 1, npgfb
102 :
103 673811 : s = 0.0_dp
104 2032188 : as = 0.0_dp
105 :
106 : ! *** Screening ***
107 :
108 10755 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109 57537 : DO j = nb + 1, nb + ncoset(lb_max)
110 264073 : DO i = na + 1, na + ncoset(la_max)
111 206536 : angab(i, j, 1) = 0.0_dp
112 206536 : angab(i, j, 2) = 0.0_dp
113 261500 : angab(i, j, 3) = 0.0_dp
114 : END DO
115 : END DO
116 2573 : nb = nb + ncoset(lb_max)
117 2573 : CYCLE
118 : END IF
119 :
120 : ! *** Calculate some prefactors ***
121 :
122 8182 : zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
123 :
124 8182 : f0 = (pi*zetp)**1.5_dp
125 8182 : f1 = zetb(jpgf)*zetp
126 8182 : f2 = 0.5_dp*zetp
127 : !
128 8182 : bc1(1) = 0._dp
129 8182 : bc1(2) = -f1*rbc(3)
130 8182 : bc1(3) = f1*rbc(2)
131 : !
132 8182 : bc2(1) = f1*rbc(3)
133 8182 : bc2(2) = 0._dp
134 8182 : bc2(3) = -f1*rbc(1)
135 : !
136 8182 : bc3(1) = -f1*rbc(2)
137 8182 : bc3(2) = f1*rbc(1)
138 8182 : bc3(3) = 0._dp
139 : !
140 8182 : ac1(1) = 0._dp
141 8182 : ac1(2) = zeta(ipgf)*zetp*rac(3)
142 8182 : ac1(3) = -zeta(ipgf)*zetp*rac(2)
143 : !
144 8182 : ac2(1) = -zeta(ipgf)*zetp*rac(3)
145 8182 : ac2(2) = 0._dp
146 8182 : ac2(3) = zeta(ipgf)*zetp*rac(1)
147 : !
148 8182 : ac3(1) = zeta(ipgf)*zetp*rac(2)
149 8182 : ac3(2) = -zeta(ipgf)*zetp*rac(1)
150 8182 : 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 8182 : s(1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2)
156 8182 : as(1, 1, 1) = 2._dp*zeta(ipgf)*f1*(rac(2)*rbc(3) - rac(3)*rbc(2))*s(1, 1)
157 8182 : as(1, 1, 2) = 2._dp*zeta(ipgf)*f1*(rac(3)*rbc(1) - rac(1)*rbc(3))*s(1, 1)
158 8182 : 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 8182 : IF (la_max > 0) THEN
163 :
164 : ! *** Vertical recurrence steps: [s|L|s] -> [a|L|s] ***
165 :
166 22232 : 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 5558 : s(2, 1) = rap(1)*s(1, 1)
172 5558 : s(3, 1) = rap(2)*s(1, 1)
173 5558 : s(4, 1) = rap(3)*s(1, 1)
174 5558 : as(2, 1, 1) = rap(1)*as(1, 1, 1) + bc1(1)*s(1, 1)
175 5558 : as(2, 1, 2) = rap(1)*as(1, 1, 2) + bc1(2)*s(1, 1)
176 5558 : as(2, 1, 3) = rap(1)*as(1, 1, 3) + bc1(3)*s(1, 1)
177 5558 : as(3, 1, 1) = rap(2)*as(1, 1, 1) + bc2(1)*s(1, 1)
178 5558 : as(3, 1, 2) = rap(2)*as(1, 1, 2) + bc2(2)*s(1, 1)
179 5558 : as(3, 1, 3) = rap(2)*as(1, 1, 3) + bc2(3)*s(1, 1)
180 5558 : as(4, 1, 1) = rap(3)*as(1, 1, 1) + bc3(1)*s(1, 1)
181 5558 : as(4, 1, 2) = rap(3)*as(1, 1, 2) + bc3(2)*s(1, 1)
182 5558 : 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 6468 : 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 7378 : 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 5558 : IF (lb_max > 0) THEN
266 :
267 61808 : DO j = 2, ncoset(lb_max)
268 317748 : DO i = 1, ncoset(la_max)
269 255940 : s(i, j) = 0.0_dp
270 255940 : as(i, j, 1) = 0.0_dp
271 255940 : as(i, j, 2) = 0.0_dp
272 312842 : as(i, j, 3) = 0.0_dp
273 : END DO
274 : END DO
275 :
276 : ! *** Horizontal recurrence steps ***
277 :
278 19624 : 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 4906 : IF (lb_max == 1) THEN
284 1982 : la_start = la_min
285 : ELSE
286 2924 : la_start = MAX(0, la_min - 1)
287 : END IF
288 :
289 10291 : DO la = la_start, la_max - 1
290 16405 : DO ax = 0, la
291 18342 : DO ay = 0, la - ax
292 6843 : az = la - ax - ay
293 : s(coset(ax, ay, az), 2) = s(coset(ax + 1, ay, az), 1) - &
294 6843 : rab(1)*s(coset(ax, ay, az), 1)
295 : s(coset(ax, ay, az), 3) = s(coset(ax, ay + 1, az), 1) - &
296 6843 : rab(2)*s(coset(ax, ay, az), 1)
297 : s(coset(ax, ay, az), 4) = s(coset(ax, ay, az + 1), 1) - &
298 6843 : 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 6843 : 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 6843 : - 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 6843 : + 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 6843 : + 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 6843 : 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 6843 : - 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 6843 : - 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 6843 : + 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 12957 : 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 15544 : DO ax = 0, la_max
335 10638 : fx = f2*REAL(ax, dp)
336 32740 : DO ay = 0, la_max - ax
337 17196 : fy = f2*REAL(ay, dp)
338 17196 : az = la_max - ax - ay
339 17196 : fz = f2*REAL(az, dp)
340 17196 : IF (ax == 0) THEN
341 10638 : 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 10638 : 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 10638 : 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 10638 : 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 6558 : 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 6558 : 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 6558 : 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 6558 : ac1(3)*s(coset(ax, ay, az), 1)
360 : END IF
361 17196 : IF (az > 0) as(coset(ax, ay, az), 2, 2) = as(coset(ax, ay, az), 2, 2) + &
362 6558 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
363 17196 : IF (ay > 0) as(coset(ax, ay, az), 2, 3) = as(coset(ax, ay, az), 2, 3) - &
364 6558 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
365 17196 : IF (ay == 0) THEN
366 10638 : 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 10638 : 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 10638 : 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 10638 : 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 6558 : 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 6558 : 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 6558 : 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 6558 : ac2(3)*s(coset(ax, ay, az), 1)
385 : END IF
386 17196 : IF (az > 0) as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay, az), 3, 1) - &
387 6558 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1)
388 17196 : IF (ax > 0) as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay, az), 3, 3) + &
389 6558 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1)
390 17196 : IF (az == 0) THEN
391 10638 : 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 10638 : 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 10638 : 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 10638 : 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 6558 : 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 6558 : 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 6558 : 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 6558 : ac3(3)*s(coset(ax, ay, az), 1)
410 : END IF
411 17196 : IF (ay > 0) as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az), 4, 1) + &
412 6558 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1)
413 17196 : IF (ax > 0) as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az), 4, 2) - &
414 17196 : 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 9942 : 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 5036 : IF (lb == lb_max) THEN
428 2924 : la_start = la_min
429 : ELSE
430 2112 : la_start = MAX(0, la_min - 1)
431 : END IF
432 :
433 10149 : DO la = la_start, la_max - 1
434 15617 : DO ax = 0, la
435 16404 : DO ay = 0, la - ax
436 5823 : 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 5823 : 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 5823 : + 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 5823 : - 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 5823 : 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 20045 : DO by = 1, lb
460 14222 : 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 14222 : 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 14222 : - 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 14222 : 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 20045 : + 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 25513 : DO bx = 1, lb
482 45866 : DO by = 0, lb - bx
483 25821 : 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 25821 : 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 25821 : 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 25821 : + 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 40043 : - 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 20388 : DO ax = 0, la_max
517 10446 : fx = f2*REAL(ax, dp)
518 31712 : DO ay = 0, la_max - ax
519 16230 : fy = f2*REAL(ay, dp)
520 16230 : az = la_max - ax - ay
521 16230 : fz = f2*REAL(az, dp)
522 :
523 : ! *** Increase the angular momentum component z of function b ***
524 :
525 16230 : f3 = f2*REAL(lb - 1, dp)
526 :
527 16230 : 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 10446 : 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 10446 : 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 10446 : 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 10446 : 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 5784 : 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 5784 : 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 5784 : 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 5784 : ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
563 : END IF
564 16230 : 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 5784 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, 0, lb - 1))
567 16230 : 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 5784 : 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 16230 : IF (ay == 0) THEN
574 10446 : bz = lb - 1
575 : s(coset(ax, ay, az), coset(0, 1, bz)) = &
576 10446 : 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 10446 : 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 10446 : 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 10446 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
586 10446 : 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 5410 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
589 10446 : 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 5410 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
592 26524 : DO by = 2, lb
593 16078 : bz = lb - by
594 16078 : 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 16078 : 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 16078 : 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 16078 : 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 16078 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
610 16078 : 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 8226 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
613 16078 : 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 18672 : f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz))
616 : END DO
617 : ELSE
618 5784 : 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 5784 : 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 5784 : 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 5784 : 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 5784 : ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
634 5784 : 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 5784 : 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 14384 : DO by = 2, lb
641 8600 : bz = lb - by
642 8600 : 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 8600 : 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 8600 : 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 8600 : 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 8600 : ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
662 8600 : 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 8600 : 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 6158 : 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 26676 : IF (ax == 0) THEN
674 36970 : DO by = 0, lb - 1
675 26524 : bz = lb - 1 - by
676 : s(coset(ax, ay, az), coset(1, by, bz)) = &
677 26524 : 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 26524 : 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 26524 : 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 26524 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
687 26524 : 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 13636 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
690 26524 : 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 24082 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
693 : END DO
694 26524 : DO bx = 2, lb
695 16078 : f3 = f2*REAL(bx - 1, dp)
696 49642 : DO by = 0, lb - bx
697 23118 : 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 23118 : 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 23118 : 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 23118 : 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 23118 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
713 23118 : 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 11746 : f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
716 23118 : 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 27824 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz))
719 : END DO
720 : END DO
721 : ELSE
722 20168 : DO by = 0, lb - 1
723 14384 : 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 14384 : 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 14384 : 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 14384 : 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 14384 : ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
739 14384 : 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 14384 : 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 6532 : f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
745 : END DO
746 14384 : DO bx = 2, lb
747 8600 : f3 = f2*REAL(bx - 1, dp)
748 26504 : DO by = 0, lb - bx
749 12120 : 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 12120 : 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 12120 : 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 12120 : 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 12120 : ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
769 12120 : 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 12120 : 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 8974 : 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 2624 : IF (lb_max > 0) THEN
789 :
790 : ! *** Vertical recurrence steps: [s|L|s] -> [s|L|b] ***
791 :
792 6864 : 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 1716 : s(1, 2) = rbp(1)*s(1, 1)
798 1716 : s(1, 3) = rbp(2)*s(1, 1)
799 1716 : s(1, 4) = rbp(3)*s(1, 1)
800 1716 : as(1, 2, 1) = rbp(1)*as(1, 1, 1) + ac1(1)*s(1, 1)
801 1716 : as(1, 2, 2) = rbp(1)*as(1, 1, 2) + ac1(2)*s(1, 1)
802 1716 : as(1, 2, 3) = rbp(1)*as(1, 1, 3) + ac1(3)*s(1, 1)
803 1716 : as(1, 3, 1) = rbp(2)*as(1, 1, 1) + ac2(1)*s(1, 1)
804 1716 : as(1, 3, 2) = rbp(2)*as(1, 1, 2) + ac2(2)*s(1, 1)
805 1716 : as(1, 3, 3) = rbp(2)*as(1, 1, 3) + ac2(3)*s(1, 1)
806 1716 : as(1, 4, 1) = rbp(3)*as(1, 1, 1) + ac3(1)*s(1, 1)
807 1716 : as(1, 4, 2) = rbp(3)*as(1, 1, 2) + ac3(2)*s(1, 1)
808 1716 : 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 4000 : 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 2284 : 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 2284 : 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 2284 : 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 2284 : ac3(3)*s(1, coset(0, 0, lb - 1))
829 :
830 : ! *** Increase the angular momentum component y of function b ***
831 :
832 2284 : bz = lb - 1
833 2284 : 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 2284 : 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 2284 : 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 2284 : ac2(3)*s(1, coset(0, 0, bz))
840 :
841 6040 : DO by = 2, lb
842 3756 : bz = lb - by
843 : s(1, coset(0, by, bz)) = rbp(2)*s(1, coset(0, by - 1, bz)) + &
844 3756 : 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 3756 : 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 3756 : 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 6040 : 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 8324 : DO by = 0, lb - 1
859 6040 : bz = lb - 1 - by
860 6040 : 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 6040 : 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 6040 : 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 8324 : ac1(3)*s(1, coset(0, by, bz))
867 : END DO
868 :
869 7756 : DO bx = 2, lb
870 3756 : f3 = f2*REAL(bx - 1, dp)
871 11636 : DO by = 0, lb - bx
872 5596 : bz = lb - bx - by
873 : s(1, coset(bx, by, bz)) = rbp(1)*s(1, coset(bx - 1, by, bz)) + &
874 5596 : 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 5596 : 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 5596 : 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 9352 : 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 98374 : DO j = 1, ncoset(lb_max)
894 409738 : DO i = 1, ncoset(la_max)
895 311364 : angab(na + i, nb + j, 1) = as(i, j, 1)
896 311364 : angab(na + i, nb + j, 2) = as(i, j, 2)
897 401556 : angab(na + i, nb + j, 3) = as(i, j, 3)
898 : END DO
899 : END DO
900 :
901 15418 : nb = nb + ncoset(lb_max)
902 :
903 : END DO
904 :
905 9658 : na = na + ncoset(la_max)
906 :
907 : END DO
908 :
909 2422 : END SUBROUTINE angmom
910 :
911 : END MODULE ai_angmom
|
