(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
17MODULE 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
33CONTAINS
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
911END 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