(git:b279b6b)
ai_angmom.F
Go to the documentation of this file.
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  SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, &
51  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  DIMENSION(ncoset(la_max), ncoset(lb_max)) :: s
68  REAL(kind=dp), &
69  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  rab = rbc - rac
87  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
88  dab = sqrt(rab2)
89 
90 ! *** Loop over all pairs of primitive Gaussian-type functions ***
91  angab = 0.0_dp
92  s = 0.0_dp
93  as = 0.0_dp
94 
95  na = 0
96 
97  DO ipgf = 1, npgfa
98 
99  nb = 0
100 
101  DO jpgf = 1, npgfb
102 
103  s = 0.0_dp
104  as = 0.0_dp
105 
106 ! *** Screening ***
107 
108  IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
109  DO j = nb + 1, nb + ncoset(lb_max)
110  DO i = na + 1, na + ncoset(la_max)
111  angab(i, j, 1) = 0.0_dp
112  angab(i, j, 2) = 0.0_dp
113  angab(i, j, 3) = 0.0_dp
114  END DO
115  END DO
116  nb = nb + ncoset(lb_max)
117  cycle
118  END IF
119 
120 ! *** Calculate some prefactors ***
121 
122  zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
123 
124  f0 = (pi*zetp)**1.5_dp
125  f1 = zetb(jpgf)*zetp
126  f2 = 0.5_dp*zetp
127  !
128  bc1(1) = 0._dp
129  bc1(2) = -f1*rbc(3)
130  bc1(3) = f1*rbc(2)
131  !
132  bc2(1) = f1*rbc(3)
133  bc2(2) = 0._dp
134  bc2(3) = -f1*rbc(1)
135  !
136  bc3(1) = -f1*rbc(2)
137  bc3(2) = f1*rbc(1)
138  bc3(3) = 0._dp
139  !
140  ac1(1) = 0._dp
141  ac1(2) = zeta(ipgf)*zetp*rac(3)
142  ac1(3) = -zeta(ipgf)*zetp*rac(2)
143  !
144  ac2(1) = -zeta(ipgf)*zetp*rac(3)
145  ac2(2) = 0._dp
146  ac2(3) = zeta(ipgf)*zetp*rac(1)
147  !
148  ac3(1) = zeta(ipgf)*zetp*rac(2)
149  ac3(2) = -zeta(ipgf)*zetp*rac(1)
150  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  s(1, 1) = f0*exp(-zeta(ipgf)*f1*rab2)
156  as(1, 1, 1) = 2._dp*zeta(ipgf)*f1*(rac(2)*rbc(3) - rac(3)*rbc(2))*s(1, 1)
157  as(1, 1, 2) = 2._dp*zeta(ipgf)*f1*(rac(3)*rbc(1) - rac(1)*rbc(3))*s(1, 1)
158  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  IF (la_max > 0) THEN
163 
164 ! *** Vertical recurrence steps: [s|L|s] -> [a|L|s] ***
165 
166  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  s(2, 1) = rap(1)*s(1, 1)
172  s(3, 1) = rap(2)*s(1, 1)
173  s(4, 1) = rap(3)*s(1, 1)
174  as(2, 1, 1) = rap(1)*as(1, 1, 1) + bc1(1)*s(1, 1)
175  as(2, 1, 2) = rap(1)*as(1, 1, 2) + bc1(2)*s(1, 1)
176  as(2, 1, 3) = rap(1)*as(1, 1, 3) + bc1(3)*s(1, 1)
177  as(3, 1, 1) = rap(2)*as(1, 1, 1) + bc2(1)*s(1, 1)
178  as(3, 1, 2) = rap(2)*as(1, 1, 2) + bc2(2)*s(1, 1)
179  as(3, 1, 3) = rap(2)*as(1, 1, 3) + bc2(3)*s(1, 1)
180  as(4, 1, 1) = rap(3)*as(1, 1, 1) + bc3(1)*s(1, 1)
181  as(4, 1, 2) = rap(3)*as(1, 1, 2) + bc3(2)*s(1, 1)
182  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  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  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  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  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  bc3(3)*s(coset(0, 0, la - 1), 1)
203 
204 ! *** Increase the angular momentum component y of function a ***
205 
206  az = la - 1
207  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  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  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  bc2(3)*s(coset(0, 0, az), 1)
214 
215  DO ay = 2, la
216  az = la - ay
217  s(coset(0, ay, az), 1) = rap(2)*s(coset(0, ay - 1, az), 1) + &
218  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  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  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  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  DO ay = 0, la - 1
233  az = la - 1 - ay
234  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  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  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  bc1(3)*s(coset(0, ay, az), 1)
241  END DO
242 
243  DO ax = 2, la
244  f3 = f2*real(ax - 1, dp)
245  DO ay = 0, la - ax
246  az = la - ax - ay
247  s(coset(ax, ay, az), 1) = rap(1)*s(coset(ax - 1, ay, az), 1) + &
248  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  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  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  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  IF (lb_max > 0) THEN
266 
267  DO j = 2, ncoset(lb_max)
268  DO i = 1, ncoset(la_max)
269  s(i, j) = 0.0_dp
270  as(i, j, 1) = 0.0_dp
271  as(i, j, 2) = 0.0_dp
272  as(i, j, 3) = 0.0_dp
273  END DO
274  END DO
275 
276 ! *** Horizontal recurrence steps ***
277 
278  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  IF (lb_max == 1) THEN
284  la_start = la_min
285  ELSE
286  la_start = max(0, la_min - 1)
287  END IF
288 
289  DO la = la_start, la_max - 1
290  DO ax = 0, la
291  DO ay = 0, la - ax
292  az = la - ax - ay
293  s(coset(ax, ay, az), 2) = s(coset(ax + 1, ay, az), 1) - &
294  rab(1)*s(coset(ax, ay, az), 1)
295  s(coset(ax, ay, az), 3) = s(coset(ax, ay + 1, az), 1) - &
296  rab(2)*s(coset(ax, ay, az), 1)
297  s(coset(ax, ay, az), 4) = s(coset(ax, ay, az + 1), 1) - &
298  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  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  - 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  + 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  + 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  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  - 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  - 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  + 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  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  DO ax = 0, la_max
335  fx = f2*real(ax, dp)
336  DO ay = 0, la_max - ax
337  fy = f2*real(ay, dp)
338  az = la_max - ax - ay
339  fz = f2*real(az, dp)
340  IF (ax == 0) THEN
341  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  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  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  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  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  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  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  ac1(3)*s(coset(ax, ay, az), 1)
360  END IF
361  IF (az > 0) as(coset(ax, ay, az), 2, 2) = as(coset(ax, ay, az), 2, 2) + &
362  f2*real(az, dp)*s(coset(ax, ay, az - 1), 1)
363  IF (ay > 0) as(coset(ax, ay, az), 2, 3) = as(coset(ax, ay, az), 2, 3) - &
364  f2*real(ay, dp)*s(coset(ax, ay - 1, az), 1)
365  IF (ay == 0) THEN
366  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  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  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  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  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  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  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  ac2(3)*s(coset(ax, ay, az), 1)
385  END IF
386  IF (az > 0) as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay, az), 3, 1) - &
387  f2*real(az, dp)*s(coset(ax, ay, az - 1), 1)
388  IF (ax > 0) as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay, az), 3, 3) + &
389  f2*real(ax, dp)*s(coset(ax - 1, ay, az), 1)
390  IF (az == 0) THEN
391  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  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  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  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  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  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  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  ac3(3)*s(coset(ax, ay, az), 1)
410  END IF
411  IF (ay > 0) as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az), 4, 1) + &
412  f2*real(ay, dp)*s(coset(ax, ay - 1, az), 1)
413  IF (ax > 0) as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az), 4, 2) - &
414  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  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  IF (lb == lb_max) THEN
428  la_start = la_min
429  ELSE
430  la_start = max(0, la_min - 1)
431  END IF
432 
433  DO la = la_start, la_max - 1
434  DO ax = 0, la
435  DO ay = 0, la - ax
436  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  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  + 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  - 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  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  DO by = 1, lb
460  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  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  - 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  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  + 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  DO bx = 1, lb
482  DO by = 0, lb - bx
483  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  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  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  + 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  - 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  DO ax = 0, la_max
517  fx = f2*real(ax, dp)
518  DO ay = 0, la_max - ax
519  fy = f2*real(ay, dp)
520  az = la_max - ax - ay
521  fz = f2*real(az, dp)
522 
523 ! *** Increase the angular momentum component z of function b ***
524 
525  f3 = f2*real(lb - 1, dp)
526 
527  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  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  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  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  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  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  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  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  ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1))
563  END IF
564  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  f2*real(ay, dp)*s(coset(ax, ay - 1, az), coset(0, 0, lb - 1))
567  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  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  IF (ay == 0) THEN
574  bz = lb - 1
575  s(coset(ax, ay, az), coset(0, 1, bz)) = &
576  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  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  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  ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
586  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
589  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  f2*real(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
592  DO by = 2, lb
593  bz = lb - by
594  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  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  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  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  ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
610  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
613  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  f2*real(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz))
616  END DO
617  ELSE
618  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  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  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  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  ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz))
634  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz))
637  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  f2*real(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz))
640  DO by = 2, lb
641  bz = lb - by
642  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  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  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  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  ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz))
662  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz))
665  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  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  IF (ax == 0) THEN
674  DO by = 0, lb - 1
675  bz = lb - 1 - by
676  s(coset(ax, ay, az), coset(1, by, bz)) = &
677  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  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  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  ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
687  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
690  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  f2*real(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
693  END DO
694  DO bx = 2, lb
695  f3 = f2*real(bx - 1, dp)
696  DO by = 0, lb - bx
697  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  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  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  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  ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
713  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
716  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  f2*real(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz))
719  END DO
720  END DO
721  ELSE
722  DO by = 0, lb - 1
723  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  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  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  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  ac1(3)*s(coset(ax, ay, az), coset(0, by, bz))
739  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz))
742  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  f2*real(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz))
745  END DO
746  DO bx = 2, lb
747  f3 = f2*real(bx - 1, dp)
748  DO by = 0, lb - bx
749  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  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  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  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  ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz))
769  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  f2*real(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz))
772  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  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  IF (lb_max > 0) THEN
789 
790 ! *** Vertical recurrence steps: [s|L|s] -> [s|L|b] ***
791 
792  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  s(1, 2) = rbp(1)*s(1, 1)
798  s(1, 3) = rbp(2)*s(1, 1)
799  s(1, 4) = rbp(3)*s(1, 1)
800  as(1, 2, 1) = rbp(1)*as(1, 1, 1) + ac1(1)*s(1, 1)
801  as(1, 2, 2) = rbp(1)*as(1, 1, 2) + ac1(2)*s(1, 1)
802  as(1, 2, 3) = rbp(1)*as(1, 1, 3) + ac1(3)*s(1, 1)
803  as(1, 3, 1) = rbp(2)*as(1, 1, 1) + ac2(1)*s(1, 1)
804  as(1, 3, 2) = rbp(2)*as(1, 1, 2) + ac2(2)*s(1, 1)
805  as(1, 3, 3) = rbp(2)*as(1, 1, 3) + ac2(3)*s(1, 1)
806  as(1, 4, 1) = rbp(3)*as(1, 1, 1) + ac3(1)*s(1, 1)
807  as(1, 4, 2) = rbp(3)*as(1, 1, 2) + ac3(2)*s(1, 1)
808  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  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  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  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  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  ac3(3)*s(1, coset(0, 0, lb - 1))
829 
830 ! *** Increase the angular momentum component y of function b ***
831 
832  bz = lb - 1
833  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  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  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  ac2(3)*s(1, coset(0, 0, bz))
840 
841  DO by = 2, lb
842  bz = lb - by
843  s(1, coset(0, by, bz)) = rbp(2)*s(1, coset(0, by - 1, bz)) + &
844  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  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  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  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  DO by = 0, lb - 1
859  bz = lb - 1 - by
860  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  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  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  ac1(3)*s(1, coset(0, by, bz))
867  END DO
868 
869  DO bx = 2, lb
870  f3 = f2*real(bx - 1, dp)
871  DO by = 0, lb - bx
872  bz = lb - bx - by
873  s(1, coset(bx, by, bz)) = rbp(1)*s(1, coset(bx - 1, by, bz)) + &
874  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  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  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  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  DO j = 1, ncoset(lb_max)
894  DO i = 1, ncoset(la_max)
895  angab(na + i, nb + j, 1) = as(i, j, 1)
896  angab(na + i, nb + j, 2) = as(i, j, 2)
897  angab(na + i, nb + j, 3) = as(i, j, 3)
898  END DO
899  END DO
900 
901  nb = nb + ncoset(lb_max)
902 
903  END DO
904 
905  na = na + ncoset(la_max)
906 
907  END DO
908 
909  END SUBROUTINE angmom
910 
911 END MODULE ai_angmom
ai_angmom
Calculation of the angular momentum integrals over Cartesian Gaussian-type functions.
Definition: ai_angmom.F:17
ai_angmom::angmom
subroutine, public angmom(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
...
Definition: ai_angmom.F:52
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition: orbital_pointers.F:45