1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of the angular momentum integrals over 8!> Cartesian Gaussian-type functions. 9!> \par Literature 10!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986) 11!> \par History 12!> none 13!> \author JGH (20.02.2005) 14! ************************************************************************************************** 15MODULE ai_angmom 16 17 USE kinds, ONLY: dp 18 USE mathconstants, ONLY: pi 19 USE orbital_pointers, ONLY: coset,& 20 ncoset 21#include "../base/base_uses.f90" 22 23 IMPLICIT NONE 24 25 PRIVATE 26 27! *** Public subroutines *** 28 29 PUBLIC :: angmom 30 31CONTAINS 32 33! ************************************************************************************************** 34!> \brief ... 35!> \param la_max ... 36!> \param npgfa ... 37!> \param zeta ... 38!> \param rpgfa ... 39!> \param la_min ... 40!> \param lb_max ... 41!> \param npgfb ... 42!> \param zetb ... 43!> \param rpgfb ... 44!> \param rac ... 45!> \param rbc ... 46!> \param angab ... 47! ************************************************************************************************** 48 SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, & 49 lb_max, npgfb, zetb, rpgfb, rac, rbc, angab) 50 51 INTEGER, INTENT(IN) :: la_max, npgfa 52 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa 53 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb 54 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb 55 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc 56 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: angab 57 58 INTEGER :: ax, ay, az, bx, by, bz, i, ipgf, j, & 59 jpgf, la, la_start, lb, na, nb 60 REAL(KIND=dp) :: dab, f0, f1, f2, f3, fx, fy, fz, rab2, & 61 zetp 62 REAL(KIND=dp), DIMENSION(3) :: ac1, ac2, ac3, bc1, bc2, bc3, rab, rap, & 63 rbp 64 REAL(KIND=dp), & 65 DIMENSION(ncoset(la_max), ncoset(lb_max)) :: s 66 REAL(KIND=dp), & 67 DIMENSION(ncoset(la_max), ncoset(lb_max), 3) :: as 68 69! ax,ay,az : Angular momentum index numbers of orbital a. 70! bx,by,bz : Angular momentum index numbers of orbital b. 71! coset : Cartesian orbital set pointer. 72! dab : Distance between the atomic centers a and b. 73! l{a,b} : Angular momentum quantum number of shell a or b. 74! l{a,b}_max: Maximum angular momentum quantum number of shell a or b. 75! l{a,b}_min: Minimum angular momentum quantum number of shell a or b. 76! rac : Distance vector between the atomic center a and C. 77! rbc : Distance vector between the atomic center b and C. 78! rpgf{a,b} : Radius of the primitive Gaussian-type function a or b. 79! angab : Shell set of angular momentum integrals. 80! zet{a,b} : Exponents of the Gaussian-type functions a or b. 81! zetp : Reciprocal of the sum of the exponents of orbital a and b. 82! *** Calculate the distance between the centers a and b *** 83 84 rab = rbc - rac 85 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 86 dab = SQRT(rab2) 87 88! *** Loop over all pairs of primitive Gaussian-type functions *** 89 angab = 0.0_dp 90 s = 0.0_dp 91 as = 0.0_dp 92 93 na = 0 94 95 DO ipgf = 1, npgfa 96 97 nb = 0 98 99 DO jpgf = 1, npgfb 100 101 s = 0.0_dp 102 as = 0.0_dp 103 104! *** Screening *** 105 106 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN 107 DO j = nb + 1, nb + ncoset(lb_max) 108 DO i = na + 1, na + ncoset(la_max) 109 angab(i, j, 1) = 0.0_dp 110 angab(i, j, 2) = 0.0_dp 111 angab(i, j, 3) = 0.0_dp 112 END DO 113 END DO 114 nb = nb + ncoset(lb_max) 115 CYCLE 116 END IF 117 118! *** Calculate some prefactors *** 119 120 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf)) 121 122 f0 = (pi*zetp)**1.5_dp 123 f1 = zetb(jpgf)*zetp 124 f2 = 0.5_dp*zetp 125 ! 126 bc1(1) = 0._dp 127 bc1(2) = -f1*rbc(3) 128 bc1(3) = f1*rbc(2) 129 ! 130 bc2(1) = f1*rbc(3) 131 bc2(2) = 0._dp 132 bc2(3) = -f1*rbc(1) 133 ! 134 bc3(1) = -f1*rbc(2) 135 bc3(2) = f1*rbc(1) 136 bc3(3) = 0._dp 137 ! 138 ac1(1) = 0._dp 139 ac1(2) = zeta(ipgf)*zetp*rac(3) 140 ac1(3) = -zeta(ipgf)*zetp*rac(2) 141 ! 142 ac2(1) = -zeta(ipgf)*zetp*rac(3) 143 ac2(2) = 0._dp 144 ac2(3) = zeta(ipgf)*zetp*rac(1) 145 ! 146 ac3(1) = zeta(ipgf)*zetp*rac(2) 147 ac3(2) = -zeta(ipgf)*zetp*rac(1) 148 ac3(3) = 0._dp 149 ! 150! *** Calculate the basic two-center overlap integral [s|s] *** 151! *** Calculate the basic two-center angular momentum integral [s|L|s] *** 152 153 s(1, 1) = f0*EXP(-zeta(ipgf)*f1*rab2) 154 as(1, 1, 1) = 2._dp*zeta(ipgf)*f1*(rac(2)*rbc(3) - rac(3)*rbc(2))*s(1, 1) 155 as(1, 1, 2) = 2._dp*zeta(ipgf)*f1*(rac(3)*rbc(1) - rac(1)*rbc(3))*s(1, 1) 156 as(1, 1, 3) = 2._dp*zeta(ipgf)*f1*(rac(1)*rbc(2) - rac(2)*rbc(1))*s(1, 1) 157 158! *** Recurrence steps: [s|L|s] -> [a|L|b] *** 159 160 IF (la_max > 0) THEN 161 162! *** Vertical recurrence steps: [s|L|s] -> [a|L|s] *** 163 164 rap(:) = f1*rab(:) 165 166! *** [p|s] = (Pi - Ai)*[s|s] (i = x,y,z) *** 167! *** [p|Ln|s] = (Pi - Ai)*[s|Ln|s]+xb/(xa+xb)*(1i x BC)*[s|s] (i = x,y,z) *** 168 169 s(2, 1) = rap(1)*s(1, 1) 170 s(3, 1) = rap(2)*s(1, 1) 171 s(4, 1) = rap(3)*s(1, 1) 172 as(2, 1, 1) = rap(1)*as(1, 1, 1) + bc1(1)*s(1, 1) 173 as(2, 1, 2) = rap(1)*as(1, 1, 2) + bc1(2)*s(1, 1) 174 as(2, 1, 3) = rap(1)*as(1, 1, 3) + bc1(3)*s(1, 1) 175 as(3, 1, 1) = rap(2)*as(1, 1, 1) + bc2(1)*s(1, 1) 176 as(3, 1, 2) = rap(2)*as(1, 1, 2) + bc2(2)*s(1, 1) 177 as(3, 1, 3) = rap(2)*as(1, 1, 3) + bc2(3)*s(1, 1) 178 as(4, 1, 1) = rap(3)*as(1, 1, 1) + bc3(1)*s(1, 1) 179 as(4, 1, 2) = rap(3)*as(1, 1, 2) + bc3(2)*s(1, 1) 180 as(4, 1, 3) = rap(3)*as(1, 1, 3) + bc3(3)*s(1, 1) 181 182! *** [a|s] = (Pi - Ai)*[a-1i|s] + f2*Ni(a-1i)*[a-2i|s] *** 183! *** [a|Ln|s] = (Pi - Ai)*[a-1i|Ln|s] + f2*Ni(a-1i)*[a-2i|Ln|s] *** 184! *** + xb/(xa+xb)*{(1i x BC)}_n*[a-1i|s] *** 185 186 DO la = 2, la_max 187 188! *** Increase the angular momentum component z of function a *** 189 190 s(coset(0, 0, la), 1) = rap(3)*s(coset(0, 0, la - 1), 1) + & 191 f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1) 192 as(coset(0, 0, la), 1, 1) = rap(3)*as(coset(0, 0, la - 1), 1, 1) + & 193 f2*REAL(la - 1, dp)*as(coset(0, 0, la - 2), 1, 1) + & 194 bc3(1)*s(coset(0, 0, la - 1), 1) 195 as(coset(0, 0, la), 1, 2) = rap(3)*as(coset(0, 0, la - 1), 1, 2) + & 196 f2*REAL(la - 1, dp)*as(coset(0, 0, la - 2), 1, 2) + & 197 bc3(2)*s(coset(0, 0, la - 1), 1) 198 as(coset(0, 0, la), 1, 3) = rap(3)*as(coset(0, 0, la - 1), 1, 3) + & 199 f2*REAL(la - 1, dp)*as(coset(0, 0, la - 2), 1, 3) + & 200 bc3(3)*s(coset(0, 0, la - 1), 1) 201 202! *** Increase the angular momentum component y of function a *** 203 204 az = la - 1 205 s(coset(0, 1, az), 1) = rap(2)*s(coset(0, 0, az), 1) 206 as(coset(0, 1, az), 1, 1) = rap(2)*as(coset(0, 0, az), 1, 1) + & 207 bc2(1)*s(coset(0, 0, az), 1) 208 as(coset(0, 1, az), 1, 2) = rap(2)*as(coset(0, 0, az), 1, 2) + & 209 bc2(2)*s(coset(0, 0, az), 1) 210 as(coset(0, 1, az), 1, 3) = rap(2)*as(coset(0, 0, az), 1, 3) + & 211 bc2(3)*s(coset(0, 0, az), 1) 212 213 DO ay = 2, la 214 az = la - ay 215 s(coset(0, ay, az), 1) = rap(2)*s(coset(0, ay - 1, az), 1) + & 216 f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1) 217 as(coset(0, ay, az), 1, 1) = rap(2)*as(coset(0, ay - 1, az), 1, 1) + & 218 f2*REAL(ay - 1, dp)*as(coset(0, ay - 2, az), 1, 1) + & 219 bc2(1)*s(coset(0, ay - 1, az), 1) 220 as(coset(0, ay, az), 1, 2) = rap(2)*as(coset(0, ay - 1, az), 1, 2) + & 221 f2*REAL(ay - 1, dp)*as(coset(0, ay - 2, az), 1, 2) + & 222 bc2(2)*s(coset(0, ay - 1, az), 1) 223 as(coset(0, ay, az), 1, 3) = rap(2)*as(coset(0, ay - 1, az), 1, 3) + & 224 f2*REAL(ay - 1, dp)*as(coset(0, ay - 2, az), 1, 3) + & 225 bc2(3)*s(coset(0, ay - 1, az), 1) 226 END DO 227 228! *** Increase the angular momentum component x of function a *** 229 230 DO ay = 0, la - 1 231 az = la - 1 - ay 232 s(coset(1, ay, az), 1) = rap(1)*s(coset(0, ay, az), 1) 233 as(coset(1, ay, az), 1, 1) = rap(1)*as(coset(0, ay, az), 1, 1) + & 234 bc1(1)*s(coset(0, ay, az), 1) 235 as(coset(1, ay, az), 1, 2) = rap(1)*as(coset(0, ay, az), 1, 2) + & 236 bc1(2)*s(coset(0, ay, az), 1) 237 as(coset(1, ay, az), 1, 3) = rap(1)*as(coset(0, ay, az), 1, 3) + & 238 bc1(3)*s(coset(0, ay, az), 1) 239 END DO 240 241 DO ax = 2, la 242 f3 = f2*REAL(ax - 1, dp) 243 DO ay = 0, la - ax 244 az = la - ax - ay 245 s(coset(ax, ay, az), 1) = rap(1)*s(coset(ax - 1, ay, az), 1) + & 246 f3*s(coset(ax - 2, ay, az), 1) 247 as(coset(ax, ay, az), 1, 1) = rap(1)*as(coset(ax - 1, ay, az), 1, 1) + & 248 f3*as(coset(ax - 2, ay, az), 1, 1) + & 249 bc1(1)*s(coset(ax - 1, ay, az), 1) 250 as(coset(ax, ay, az), 1, 2) = rap(1)*as(coset(ax - 1, ay, az), 1, 2) + & 251 f3*as(coset(ax - 2, ay, az), 1, 2) + & 252 bc1(2)*s(coset(ax - 1, ay, az), 1) 253 as(coset(ax, ay, az), 1, 3) = rap(1)*as(coset(ax - 1, ay, az), 1, 3) + & 254 f3*as(coset(ax - 2, ay, az), 1, 3) + & 255 bc1(3)*s(coset(ax - 1, ay, az), 1) 256 END DO 257 END DO 258 259 END DO 260 261! *** Recurrence steps: [a|L|s] -> [a|L|b] *** 262 263 IF (lb_max > 0) THEN 264 265 DO j = 2, ncoset(lb_max) 266 DO i = 1, ncoset(la_max) 267 s(i, j) = 0.0_dp 268 as(i, j, 1) = 0.0_dp 269 as(i, j, 2) = 0.0_dp 270 as(i, j, 3) = 0.0_dp 271 END DO 272 END DO 273 274! *** Horizontal recurrence steps *** 275 276 rbp(:) = rap(:) - rab(:) 277 278! *** [a|L|p] = [a+1i|Lm|s] - (Bi - Ai)*[a|Lm|s] *** 279! *** + [a+1k|s] + (Ak - Ck)*[a|s] eps(i,m,k) 280 281 IF (lb_max == 1) THEN 282 la_start = la_min 283 ELSE 284 la_start = MAX(0, la_min - 1) 285 END IF 286 287 DO la = la_start, la_max - 1 288 DO ax = 0, la 289 DO ay = 0, la - ax 290 az = la - ax - ay 291 s(coset(ax, ay, az), 2) = s(coset(ax + 1, ay, az), 1) - & 292 rab(1)*s(coset(ax, ay, az), 1) 293 s(coset(ax, ay, az), 3) = s(coset(ax, ay + 1, az), 1) - & 294 rab(2)*s(coset(ax, ay, az), 1) 295 s(coset(ax, ay, az), 4) = s(coset(ax, ay, az + 1), 1) - & 296 rab(3)*s(coset(ax, ay, az), 1) 297 as(coset(ax, ay, az), 2, 1) = as(coset(ax + 1, ay, az), 1, 1) - & 298 rab(1)*as(coset(ax, ay, az), 1, 1) 299 as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay + 1, az), 1, 1) - & 300 rab(2)*as(coset(ax, ay, az), 1, 1) & 301 - s(coset(ax, ay, az + 1), 1) - rac(3)*s(coset(ax, ay, az), 1) 302 as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az + 1), 1, 1) - & 303 rab(3)*as(coset(ax, ay, az), 1, 1) & 304 + s(coset(ax, ay + 1, az), 1) + rac(2)*s(coset(ax, ay, az), 1) 305 as(coset(ax, ay, az), 2, 2) = as(coset(ax + 1, ay, az), 1, 2) - & 306 rab(1)*as(coset(ax, ay, az), 1, 2) & 307 + s(coset(ax, ay, az + 1), 1) + rac(3)*s(coset(ax, ay, az), 1) 308 as(coset(ax, ay, az), 3, 2) = as(coset(ax, ay + 1, az), 1, 2) - & 309 rab(2)*as(coset(ax, ay, az), 1, 2) 310 as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az + 1), 1, 2) - & 311 rab(3)*as(coset(ax, ay, az), 1, 2) & 312 - s(coset(ax + 1, ay, az), 1) - rac(1)*s(coset(ax, ay, az), 1) 313 as(coset(ax, ay, az), 2, 3) = as(coset(ax + 1, ay, az), 1, 3) - & 314 rab(1)*as(coset(ax, ay, az), 1, 3) & 315 - s(coset(ax, ay + 1, az), 1) - rac(2)*s(coset(ax, ay, az), 1) 316 as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay + 1, az), 1, 3) - & 317 rab(2)*as(coset(ax, ay, az), 1, 3) & 318 + s(coset(ax + 1, ay, az), 1) + rac(1)*s(coset(ax, ay, az), 1) 319 as(coset(ax, ay, az), 4, 3) = as(coset(ax, ay, az + 1), 1, 3) - & 320 rab(3)*as(coset(ax, ay, az), 1, 3) 321 END DO 322 END DO 323 END DO 324 325! *** Vertical recurrence step *** 326 327! *** [a|p] = (Pi - Bi)*[a|s] + f2*Ni(a)*[a-1i|s] *** 328! *** [a|L|p] = (Pi - Bi)*[a|L|s] + f2*Ni(a)*[a-1i|L|s] *** 329! *** + xa/(xa+xb)*(AC x 1i)*[a|s] *** 330! *** + 0.5*zetp*SUM_k Nk(a)*(1k x 1i)*[a-1k|s] *** 331 332 DO ax = 0, la_max 333 fx = f2*REAL(ax, dp) 334 DO ay = 0, la_max - ax 335 fy = f2*REAL(ay, dp) 336 az = la_max - ax - ay 337 fz = f2*REAL(az, dp) 338 IF (ax == 0) THEN 339 s(coset(ax, ay, az), 2) = rbp(1)*s(coset(ax, ay, az), 1) 340 as(coset(ax, ay, az), 2, 1) = rbp(1)*as(coset(ax, ay, az), 1, 1) + & 341 ac1(1)*s(coset(ax, ay, az), 1) 342 as(coset(ax, ay, az), 2, 2) = rbp(1)*as(coset(ax, ay, az), 1, 2) + & 343 ac1(2)*s(coset(ax, ay, az), 1) 344 as(coset(ax, ay, az), 2, 3) = rbp(1)*as(coset(ax, ay, az), 1, 3) + & 345 ac1(3)*s(coset(ax, ay, az), 1) 346 ELSE 347 s(coset(ax, ay, az), 2) = rbp(1)*s(coset(ax, ay, az), 1) + & 348 fx*s(coset(ax - 1, ay, az), 1) 349 as(coset(ax, ay, az), 2, 1) = rbp(1)*as(coset(ax, ay, az), 1, 1) + & 350 fx*as(coset(ax - 1, ay, az), 1, 1) + & 351 ac1(1)*s(coset(ax, ay, az), 1) 352 as(coset(ax, ay, az), 2, 2) = rbp(1)*as(coset(ax, ay, az), 1, 2) + & 353 fx*as(coset(ax - 1, ay, az), 1, 2) + & 354 ac1(2)*s(coset(ax, ay, az), 1) 355 as(coset(ax, ay, az), 2, 3) = rbp(1)*as(coset(ax, ay, az), 1, 3) + & 356 fx*as(coset(ax - 1, ay, az), 1, 3) + & 357 ac1(3)*s(coset(ax, ay, az), 1) 358 END IF 359 IF (az > 0) as(coset(ax, ay, az), 2, 2) = as(coset(ax, ay, az), 2, 2) + & 360 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1) 361 IF (ay > 0) as(coset(ax, ay, az), 2, 3) = as(coset(ax, ay, az), 2, 3) - & 362 f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1) 363 IF (ay == 0) THEN 364 s(coset(ax, ay, az), 3) = rbp(2)*s(coset(ax, ay, az), 1) 365 as(coset(ax, ay, az), 3, 1) = rbp(2)*as(coset(ax, ay, az), 1, 1) + & 366 ac2(1)*s(coset(ax, ay, az), 1) 367 as(coset(ax, ay, az), 3, 2) = rbp(2)*as(coset(ax, ay, az), 1, 2) + & 368 ac2(2)*s(coset(ax, ay, az), 1) 369 as(coset(ax, ay, az), 3, 3) = rbp(2)*as(coset(ax, ay, az), 1, 3) + & 370 ac2(3)*s(coset(ax, ay, az), 1) 371 ELSE 372 s(coset(ax, ay, az), 3) = rbp(2)*s(coset(ax, ay, az), 1) + & 373 fy*s(coset(ax, ay - 1, az), 1) 374 as(coset(ax, ay, az), 3, 1) = rbp(2)*as(coset(ax, ay, az), 1, 1) + & 375 fy*as(coset(ax, ay - 1, az), 1, 1) + & 376 ac2(1)*s(coset(ax, ay, az), 1) 377 as(coset(ax, ay, az), 3, 2) = rbp(2)*as(coset(ax, ay, az), 1, 2) + & 378 fy*as(coset(ax, ay - 1, az), 1, 2) + & 379 ac2(2)*s(coset(ax, ay, az), 1) 380 as(coset(ax, ay, az), 3, 3) = rbp(2)*as(coset(ax, ay, az), 1, 3) + & 381 fy*as(coset(ax, ay - 1, az), 1, 3) + & 382 ac2(3)*s(coset(ax, ay, az), 1) 383 END IF 384 IF (az > 0) as(coset(ax, ay, az), 3, 1) = as(coset(ax, ay, az), 3, 1) - & 385 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), 1) 386 IF (ax > 0) as(coset(ax, ay, az), 3, 3) = as(coset(ax, ay, az), 3, 3) + & 387 f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1) 388 IF (az == 0) THEN 389 s(coset(ax, ay, az), 4) = rbp(3)*s(coset(ax, ay, az), 1) 390 as(coset(ax, ay, az), 4, 1) = rbp(3)*as(coset(ax, ay, az), 1, 1) + & 391 ac3(1)*s(coset(ax, ay, az), 1) 392 as(coset(ax, ay, az), 4, 2) = rbp(3)*as(coset(ax, ay, az), 1, 2) + & 393 ac3(2)*s(coset(ax, ay, az), 1) 394 as(coset(ax, ay, az), 4, 3) = rbp(3)*as(coset(ax, ay, az), 1, 3) + & 395 ac3(3)*s(coset(ax, ay, az), 1) 396 ELSE 397 s(coset(ax, ay, az), 4) = rbp(3)*s(coset(ax, ay, az), 1) + & 398 fz*s(coset(ax, ay, az - 1), 1) 399 as(coset(ax, ay, az), 4, 1) = rbp(3)*as(coset(ax, ay, az), 1, 1) + & 400 fz*as(coset(ax, ay, az - 1), 1, 1) + & 401 ac3(1)*s(coset(ax, ay, az), 1) 402 as(coset(ax, ay, az), 4, 2) = rbp(3)*as(coset(ax, ay, az), 1, 2) + & 403 fz*as(coset(ax, ay, az - 1), 1, 2) + & 404 ac3(2)*s(coset(ax, ay, az), 1) 405 as(coset(ax, ay, az), 4, 3) = rbp(3)*as(coset(ax, ay, az), 1, 3) + & 406 fz*as(coset(ax, ay, az - 1), 1, 3) + & 407 ac3(3)*s(coset(ax, ay, az), 1) 408 END IF 409 IF (ay > 0) as(coset(ax, ay, az), 4, 1) = as(coset(ax, ay, az), 4, 1) + & 410 f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), 1) 411 IF (ax > 0) as(coset(ax, ay, az), 4, 2) = as(coset(ax, ay, az), 4, 2) - & 412 f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), 1) 413 END DO 414 END DO 415 416! *** Recurrence steps: [a|L|p] -> [a|L|b] *** 417 418 DO lb = 2, lb_max 419 420! *** Horizontal recurrence steps *** 421 422! *** [a|Lm|b] = [a+1i|Lm|b-1i] - (Bi - Ai)*[a|Lm|b-1i] *** 423! *** + [a+1k|b-1i] + (Ak - Ck)*[a|b-1i] eps(i,m,k) 424 425 IF (lb == lb_max) THEN 426 la_start = la_min 427 ELSE 428 la_start = MAX(0, la_min - 1) 429 END IF 430 431 DO la = la_start, la_max - 1 432 DO ax = 0, la 433 DO ay = 0, la - ax 434 az = la - ax - ay 435 436! *** Shift of angular momentum component z from a to b *** 437 438 s(coset(ax, ay, az), coset(0, 0, lb)) = & 439 s(coset(ax, ay, az + 1), coset(0, 0, lb - 1)) - & 440 rab(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 441 as(coset(ax, ay, az), coset(0, 0, lb), 1) = & 442 as(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 1) - & 443 rab(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 1) & 444 + s(coset(ax, ay + 1, az), coset(0, 0, lb - 1)) & 445 + rac(2)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 446 as(coset(ax, ay, az), coset(0, 0, lb), 2) = & 447 as(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 2) - & 448 rab(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 2) & 449 - s(coset(ax + 1, ay, az), coset(0, 0, lb - 1)) & 450 - rac(1)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 451 as(coset(ax, ay, az), coset(0, 0, lb), 3) = & 452 as(coset(ax, ay, az + 1), coset(0, 0, lb - 1), 3) - & 453 rab(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 3) 454 455! *** Shift of angular momentum component y from a to b *** 456 457 DO by = 1, lb 458 bz = lb - by 459 s(coset(ax, ay, az), coset(0, by, bz)) = & 460 s(coset(ax, ay + 1, az), coset(0, by - 1, bz)) - & 461 rab(2)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 462 as(coset(ax, ay, az), coset(0, by, bz), 1) = & 463 as(coset(ax, ay + 1, az), coset(0, by - 1, bz), 1) - & 464 rab(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 1) & 465 - s(coset(ax, ay, az + 1), coset(0, by - 1, bz)) & 466 - rac(3)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 467 as(coset(ax, ay, az), coset(0, by, bz), 2) = & 468 as(coset(ax, ay + 1, az), coset(0, by - 1, bz), 2) - & 469 rab(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 2) 470 as(coset(ax, ay, az), coset(0, by, bz), 3) = & 471 as(coset(ax, ay + 1, az), coset(0, by - 1, bz), 3) - & 472 rab(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 3) & 473 + s(coset(ax + 1, ay, az), coset(0, by - 1, bz)) & 474 + rac(1)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 475 END DO 476 477! *** Shift of angular momentum component x from a to b *** 478 479 DO bx = 1, lb 480 DO by = 0, lb - bx 481 bz = lb - bx - by 482 s(coset(ax, ay, az), coset(bx, by, bz)) = & 483 s(coset(ax + 1, ay, az), coset(bx - 1, by, bz)) - & 484 rab(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 485 as(coset(ax, ay, az), coset(bx, by, bz), 1) = & 486 as(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 1) - & 487 rab(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 1) 488 as(coset(ax, ay, az), coset(bx, by, bz), 2) = & 489 as(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 2) - & 490 rab(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 2) & 491 + s(coset(ax, ay, az + 1), coset(bx - 1, by, bz)) & 492 + rac(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 493 as(coset(ax, ay, az), coset(bx, by, bz), 3) = & 494 as(coset(ax + 1, ay, az), coset(bx - 1, by, bz), 3) - & 495 rab(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 3) & 496 - s(coset(ax, ay + 1, az), coset(bx - 1, by, bz)) & 497 - rac(2)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 498 END DO 499 END DO 500 501 END DO 502 END DO 503 END DO 504 505! *** Vertical recurrence step *** 506 507! *** [a|b] = (Pi - Bi)*[a|b-1i] + f2*Ni(a)*[a-1i|b-1i] + *** 508! *** f2*Ni(b-1i)*[a|b-2i] *** 509! *** [a|L|b] = (Pi - Bi)*[a|L|b-1i] + f2*Ni(a)*[a-1i|L|b-1i] + *** 510! *** f2*Ni(b-1i)*[a|L|b-2i] *** 511! *** + xa/(xa+xb)*(AC x 1i)*[a|b-1i] *** 512! *** + 0.5*zetp*SUM_k Nk(a)*(1k x 1i)*[a-1k|b-1i] *** 513 514 DO ax = 0, la_max 515 fx = f2*REAL(ax, dp) 516 DO ay = 0, la_max - ax 517 fy = f2*REAL(ay, dp) 518 az = la_max - ax - ay 519 fz = f2*REAL(az, dp) 520 521! *** Increase the angular momentum component z of function b *** 522 523 f3 = f2*REAL(lb - 1, dp) 524 525 IF (az == 0) THEN 526 s(coset(ax, ay, az), coset(0, 0, lb)) = & 527 rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) + & 528 f3*s(coset(ax, ay, az), coset(0, 0, lb - 2)) 529 as(coset(ax, ay, az), coset(0, 0, lb), 1) = & 530 rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 1) + & 531 f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 1) + & 532 ac3(1)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 533 as(coset(ax, ay, az), coset(0, 0, lb), 2) = & 534 rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 2) + & 535 f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 2) + & 536 ac3(2)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 537 as(coset(ax, ay, az), coset(0, 0, lb), 3) = & 538 rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 3) + & 539 f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 3) + & 540 ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 541 ELSE 542 s(coset(ax, ay, az), coset(0, 0, lb)) = & 543 rbp(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) + & 544 fz*s(coset(ax, ay, az - 1), coset(0, 0, lb - 1)) + & 545 f3*s(coset(ax, ay, az), coset(0, 0, lb - 2)) 546 as(coset(ax, ay, az), coset(0, 0, lb), 1) = & 547 rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 1) + & 548 fz*as(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 1) + & 549 f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 1) + & 550 ac3(1)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 551 as(coset(ax, ay, az), coset(0, 0, lb), 2) = & 552 rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 2) + & 553 fz*as(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 2) + & 554 f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 2) + & 555 ac3(2)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 556 as(coset(ax, ay, az), coset(0, 0, lb), 3) = & 557 rbp(3)*as(coset(ax, ay, az), coset(0, 0, lb - 1), 3) + & 558 fz*as(coset(ax, ay, az - 1), coset(0, 0, lb - 1), 3) + & 559 f3*as(coset(ax, ay, az), coset(0, 0, lb - 2), 3) + & 560 ac3(3)*s(coset(ax, ay, az), coset(0, 0, lb - 1)) 561 END IF 562 IF (ay > 0) as(coset(ax, ay, az), coset(0, 0, lb), 1) = & 563 as(coset(ax, ay, az), coset(0, 0, lb), 1) + & 564 f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, 0, lb - 1)) 565 IF (ax > 0) as(coset(ax, ay, az), coset(0, 0, lb), 2) = & 566 as(coset(ax, ay, az), coset(0, 0, lb), 2) - & 567 f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, lb - 1)) 568 569! *** Increase the angular momentum component y of function b *** 570 571 IF (ay == 0) THEN 572 bz = lb - 1 573 s(coset(ax, ay, az), coset(0, 1, bz)) = & 574 rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz)) 575 as(coset(ax, ay, az), coset(0, 1, bz), 1) = & 576 rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 1) + & 577 ac2(1)*s(coset(ax, ay, az), coset(0, 0, bz)) 578 as(coset(ax, ay, az), coset(0, 1, bz), 2) = & 579 rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 2) + & 580 ac2(2)*s(coset(ax, ay, az), coset(0, 0, bz)) 581 as(coset(ax, ay, az), coset(0, 1, bz), 3) = & 582 rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 3) + & 583 ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz)) 584 IF (az > 0) as(coset(ax, ay, az), coset(0, 1, bz), 1) = & 585 as(coset(ax, ay, az), coset(0, 1, bz), 1) - & 586 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz)) 587 IF (ax > 0) as(coset(ax, ay, az), coset(0, 1, bz), 3) = & 588 as(coset(ax, ay, az), coset(0, 1, bz), 3) + & 589 f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz)) 590 DO by = 2, lb 591 bz = lb - by 592 f3 = f2*REAL(by - 1, dp) 593 s(coset(ax, ay, az), coset(0, by, bz)) = & 594 rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz)) + & 595 f3*s(coset(ax, ay, az), coset(0, by - 2, bz)) 596 as(coset(ax, ay, az), coset(0, by, bz), 1) = & 597 rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 1) + & 598 f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 1) + & 599 ac2(1)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 600 as(coset(ax, ay, az), coset(0, by, bz), 2) = & 601 rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 2) + & 602 f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 2) + & 603 ac2(2)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 604 as(coset(ax, ay, az), coset(0, by, bz), 3) = & 605 rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 3) + & 606 f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 3) + & 607 ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 608 IF (az > 0) as(coset(ax, ay, az), coset(0, by, bz), 1) = & 609 as(coset(ax, ay, az), coset(0, by, bz), 1) - & 610 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz)) 611 IF (ax > 0) as(coset(ax, ay, az), coset(0, by, bz), 3) = & 612 as(coset(ax, ay, az), coset(0, by, bz), 3) + & 613 f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz)) 614 END DO 615 ELSE 616 bz = lb - 1 617 s(coset(ax, ay, az), coset(0, 1, bz)) = & 618 rbp(2)*s(coset(ax, ay, az), coset(0, 0, bz)) + & 619 fy*s(coset(ax, ay - 1, az), coset(0, 0, bz)) 620 as(coset(ax, ay, az), coset(0, 1, bz), 1) = & 621 rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 1) + & 622 fy*as(coset(ax, ay - 1, az), coset(0, 0, bz), 1) + & 623 ac2(1)*s(coset(ax, ay, az), coset(0, 0, bz)) 624 as(coset(ax, ay, az), coset(0, 1, bz), 2) = & 625 rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 2) + & 626 fy*as(coset(ax, ay - 1, az), coset(0, 0, bz), 2) + & 627 ac2(2)*s(coset(ax, ay, az), coset(0, 0, bz)) 628 as(coset(ax, ay, az), coset(0, 1, bz), 3) = & 629 rbp(2)*as(coset(ax, ay, az), coset(0, 0, bz), 3) + & 630 fy*as(coset(ax, ay - 1, az), coset(0, 0, bz), 3) + & 631 ac2(3)*s(coset(ax, ay, az), coset(0, 0, bz)) 632 IF (az > 0) as(coset(ax, ay, az), coset(0, 1, bz), 1) = & 633 as(coset(ax, ay, az), coset(0, 1, bz), 1) - & 634 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, 0, bz)) 635 IF (ax > 0) as(coset(ax, ay, az), coset(0, 1, bz), 3) = & 636 as(coset(ax, ay, az), coset(0, 1, bz), 3) + & 637 f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, 0, bz)) 638 DO by = 2, lb 639 bz = lb - by 640 f3 = f2*REAL(by - 1, dp) 641 s(coset(ax, ay, az), coset(0, by, bz)) = & 642 rbp(2)*s(coset(ax, ay, az), coset(0, by - 1, bz)) + & 643 fy*s(coset(ax, ay - 1, az), coset(0, by - 1, bz)) + & 644 f3*s(coset(ax, ay, az), coset(0, by - 2, bz)) 645 as(coset(ax, ay, az), coset(0, by, bz), 1) = & 646 rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 1) + & 647 fy*as(coset(ax, ay - 1, az), coset(0, by - 1, bz), 1) + & 648 f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 1) + & 649 ac2(1)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 650 as(coset(ax, ay, az), coset(0, by, bz), 2) = & 651 rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 2) + & 652 fy*as(coset(ax, ay - 1, az), coset(0, by - 1, bz), 2) + & 653 f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 2) + & 654 ac2(2)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 655 as(coset(ax, ay, az), coset(0, by, bz), 3) = & 656 rbp(2)*as(coset(ax, ay, az), coset(0, by - 1, bz), 3) + & 657 fy*as(coset(ax, ay - 1, az), coset(0, by - 1, bz), 3) + & 658 f3*as(coset(ax, ay, az), coset(0, by - 2, bz), 3) + & 659 ac2(3)*s(coset(ax, ay, az), coset(0, by - 1, bz)) 660 IF (az > 0) as(coset(ax, ay, az), coset(0, by, bz), 1) = & 661 as(coset(ax, ay, az), coset(0, by, bz), 1) - & 662 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by - 1, bz)) 663 IF (ax > 0) as(coset(ax, ay, az), coset(0, by, bz), 3) = & 664 as(coset(ax, ay, az), coset(0, by, bz), 3) + & 665 f2*REAL(ax, dp)*s(coset(ax - 1, ay, az), coset(0, by - 1, bz)) 666 END DO 667 END IF 668 669! *** Increase the angular momentum component x of function b *** 670 671 IF (ax == 0) THEN 672 DO by = 0, lb - 1 673 bz = lb - 1 - by 674 s(coset(ax, ay, az), coset(1, by, bz)) = & 675 rbp(1)*s(coset(ax, ay, az), coset(0, by, bz)) 676 as(coset(ax, ay, az), coset(1, by, bz), 1) = & 677 rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 1) + & 678 ac1(1)*s(coset(ax, ay, az), coset(0, by, bz)) 679 as(coset(ax, ay, az), coset(1, by, bz), 2) = & 680 rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 2) + & 681 ac1(2)*s(coset(ax, ay, az), coset(0, by, bz)) 682 as(coset(ax, ay, az), coset(1, by, bz), 3) = & 683 rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 3) + & 684 ac1(3)*s(coset(ax, ay, az), coset(0, by, bz)) 685 IF (az > 0) as(coset(ax, ay, az), coset(1, by, bz), 2) = & 686 as(coset(ax, ay, az), coset(1, by, bz), 2) + & 687 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz)) 688 IF (ay > 0) as(coset(ax, ay, az), coset(1, by, bz), 3) = & 689 as(coset(ax, ay, az), coset(1, by, bz), 3) - & 690 f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz)) 691 END DO 692 DO bx = 2, lb 693 f3 = f2*REAL(bx - 1, dp) 694 DO by = 0, lb - bx 695 bz = lb - bx - by 696 s(coset(ax, ay, az), coset(bx, by, bz)) = & 697 rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) + & 698 f3*s(coset(ax, ay, az), coset(bx - 2, by, bz)) 699 as(coset(ax, ay, az), coset(bx, by, bz), 1) = & 700 rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 1) + & 701 f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 1) + & 702 ac1(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 703 as(coset(ax, ay, az), coset(bx, by, bz), 2) = & 704 rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 2) + & 705 f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 2) + & 706 ac1(2)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 707 as(coset(ax, ay, az), coset(bx, by, bz), 3) = & 708 rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 3) + & 709 f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 3) + & 710 ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 711 IF (az > 0) as(coset(ax, ay, az), coset(bx, by, bz), 2) = & 712 as(coset(ax, ay, az), coset(bx, by, bz), 2) + & 713 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz)) 714 IF (ay > 0) as(coset(ax, ay, az), coset(bx, by, bz), 3) = & 715 as(coset(ax, ay, az), coset(bx, by, bz), 3) - & 716 f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz)) 717 END DO 718 END DO 719 ELSE 720 DO by = 0, lb - 1 721 bz = lb - 1 - by 722 s(coset(ax, ay, az), coset(1, by, bz)) = & 723 rbp(1)*s(coset(ax, ay, az), coset(0, by, bz)) + & 724 fx*s(coset(ax - 1, ay, az), coset(0, by, bz)) 725 as(coset(ax, ay, az), coset(1, by, bz), 1) = & 726 rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 1) + & 727 fx*as(coset(ax - 1, ay, az), coset(0, by, bz), 1) + & 728 ac1(1)*s(coset(ax, ay, az), coset(0, by, bz)) 729 as(coset(ax, ay, az), coset(1, by, bz), 2) = & 730 rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 2) + & 731 fx*as(coset(ax - 1, ay, az), coset(0, by, bz), 2) + & 732 ac1(2)*s(coset(ax, ay, az), coset(0, by, bz)) 733 as(coset(ax, ay, az), coset(1, by, bz), 3) = & 734 rbp(1)*as(coset(ax, ay, az), coset(0, by, bz), 3) + & 735 fx*as(coset(ax - 1, ay, az), coset(0, by, bz), 3) + & 736 ac1(3)*s(coset(ax, ay, az), coset(0, by, bz)) 737 IF (az > 0) as(coset(ax, ay, az), coset(1, by, bz), 2) = & 738 as(coset(ax, ay, az), coset(1, by, bz), 2) + & 739 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(0, by, bz)) 740 IF (ay > 0) as(coset(ax, ay, az), coset(1, by, bz), 3) = & 741 as(coset(ax, ay, az), coset(1, by, bz), 3) - & 742 f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(0, by, bz)) 743 END DO 744 DO bx = 2, lb 745 f3 = f2*REAL(bx - 1, dp) 746 DO by = 0, lb - bx 747 bz = lb - bx - by 748 s(coset(ax, ay, az), coset(bx, by, bz)) = & 749 rbp(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) + & 750 fx*s(coset(ax - 1, ay, az), coset(bx - 1, by, bz)) + & 751 f3*s(coset(ax, ay, az), coset(bx - 2, by, bz)) 752 as(coset(ax, ay, az), coset(bx, by, bz), 1) = & 753 rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 1) + & 754 fx*as(coset(ax - 1, ay, az), coset(bx - 1, by, bz), 1) + & 755 f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 1) + & 756 ac1(1)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 757 as(coset(ax, ay, az), coset(bx, by, bz), 2) = & 758 rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 2) + & 759 fx*as(coset(ax - 1, ay, az), coset(bx - 1, by, bz), 2) + & 760 f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 2) + & 761 ac1(2)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 762 as(coset(ax, ay, az), coset(bx, by, bz), 3) = & 763 rbp(1)*as(coset(ax, ay, az), coset(bx - 1, by, bz), 3) + & 764 fx*as(coset(ax - 1, ay, az), coset(bx - 1, by, bz), 3) + & 765 f3*as(coset(ax, ay, az), coset(bx - 2, by, bz), 3) + & 766 ac1(3)*s(coset(ax, ay, az), coset(bx - 1, by, bz)) 767 IF (az > 0) as(coset(ax, ay, az), coset(bx, by, bz), 2) = & 768 as(coset(ax, ay, az), coset(bx, by, bz), 2) + & 769 f2*REAL(az, dp)*s(coset(ax, ay, az - 1), coset(bx - 1, by, bz)) 770 IF (ay > 0) as(coset(ax, ay, az), coset(bx, by, bz), 3) = & 771 as(coset(ax, ay, az), coset(bx, by, bz), 3) - & 772 f2*REAL(ay, dp)*s(coset(ax, ay - 1, az), coset(bx - 1, by, bz)) 773 END DO 774 END DO 775 END IF 776 777 END DO 778 END DO 779 780 END DO 781 782 END IF 783 784 ELSE 785 786 IF (lb_max > 0) THEN 787 788! *** Vertical recurrence steps: [s|L|s] -> [s|L|b] *** 789 790 rbp(:) = (f1 - 1.0_dp)*rab(:) 791 792! *** [s|p] = (Pi - Bi)*[s|s] *** 793! *** [s|L|p] = (Pi - Bi)*[s|L|s] + xa/(xa+xb)*(AC x 1i)*[s|s] *** 794 795 s(1, 2) = rbp(1)*s(1, 1) 796 s(1, 3) = rbp(2)*s(1, 1) 797 s(1, 4) = rbp(3)*s(1, 1) 798 as(1, 2, 1) = rbp(1)*as(1, 1, 1) + ac1(1)*s(1, 1) 799 as(1, 2, 2) = rbp(1)*as(1, 1, 2) + ac1(2)*s(1, 1) 800 as(1, 2, 3) = rbp(1)*as(1, 1, 3) + ac1(3)*s(1, 1) 801 as(1, 3, 1) = rbp(2)*as(1, 1, 1) + ac2(1)*s(1, 1) 802 as(1, 3, 2) = rbp(2)*as(1, 1, 2) + ac2(2)*s(1, 1) 803 as(1, 3, 3) = rbp(2)*as(1, 1, 3) + ac2(3)*s(1, 1) 804 as(1, 4, 1) = rbp(3)*as(1, 1, 1) + ac3(1)*s(1, 1) 805 as(1, 4, 2) = rbp(3)*as(1, 1, 2) + ac3(2)*s(1, 1) 806 as(1, 4, 3) = rbp(3)*as(1, 1, 3) + ac3(3)*s(1, 1) 807 808! *** [s|b] = (Pi - Bi)*[s|b-1i] + f2*Ni(b-1i)*[s|b-2i] *** 809! *** [s|L|b] = (Pi - Bi)*[s|L|b-1i] + f2*Ni(b-1i)*[s|L|b-2i] *** 810! *** + xa/(xa+xb)*(AC x 1i)*[s|s-1i] *** 811 812 DO lb = 2, lb_max 813 814! *** Increase the angular momentum component z of function b *** 815 816 s(1, coset(0, 0, lb)) = rbp(3)*s(1, coset(0, 0, lb - 1)) + & 817 f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2)) 818 as(1, coset(0, 0, lb), 1) = rbp(3)*as(1, coset(0, 0, lb - 1), 1) + & 819 f2*REAL(lb - 1, dp)*as(1, coset(0, 0, lb - 2), 1) + & 820 ac3(1)*s(1, coset(0, 0, lb - 1)) 821 as(1, coset(0, 0, lb), 2) = rbp(3)*as(1, coset(0, 0, lb - 1), 2) + & 822 f2*REAL(lb - 1, dp)*as(1, coset(0, 0, lb - 2), 2) + & 823 ac3(2)*s(1, coset(0, 0, lb - 1)) 824 as(1, coset(0, 0, lb), 3) = rbp(3)*as(1, coset(0, 0, lb - 1), 3) + & 825 f2*REAL(lb - 1, dp)*as(1, coset(0, 0, lb - 2), 3) + & 826 ac3(3)*s(1, coset(0, 0, lb - 1)) 827 828! *** Increase the angular momentum component y of function b *** 829 830 bz = lb - 1 831 s(1, coset(0, 1, bz)) = rbp(2)*s(1, coset(0, 0, bz)) 832 as(1, coset(0, 1, bz), 1) = rbp(2)*as(1, coset(0, 0, bz), 1) + & 833 ac2(1)*s(1, coset(0, 0, bz)) 834 as(1, coset(0, 1, bz), 2) = rbp(2)*as(1, coset(0, 0, bz), 2) + & 835 ac2(2)*s(1, coset(0, 0, bz)) 836 as(1, coset(0, 1, bz), 3) = rbp(2)*as(1, coset(0, 0, bz), 3) + & 837 ac2(3)*s(1, coset(0, 0, bz)) 838 839 DO by = 2, lb 840 bz = lb - by 841 s(1, coset(0, by, bz)) = rbp(2)*s(1, coset(0, by - 1, bz)) + & 842 f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz)) 843 as(1, coset(0, by, bz), 1) = rbp(2)*as(1, coset(0, by - 1, bz), 1) + & 844 f2*REAL(by - 1, dp)*as(1, coset(0, by - 2, bz), 1) + & 845 ac2(1)*s(1, coset(0, by - 1, bz)) 846 as(1, coset(0, by, bz), 2) = rbp(2)*as(1, coset(0, by - 1, bz), 2) + & 847 f2*REAL(by - 1, dp)*as(1, coset(0, by - 2, bz), 2) + & 848 ac2(2)*s(1, coset(0, by - 1, bz)) 849 as(1, coset(0, by, bz), 3) = rbp(2)*as(1, coset(0, by - 1, bz), 3) + & 850 f2*REAL(by - 1, dp)*as(1, coset(0, by - 2, bz), 3) + & 851 ac2(3)*s(1, coset(0, by - 1, bz)) 852 END DO 853 854! *** Increase the angular momentum component x of function b *** 855 856 DO by = 0, lb - 1 857 bz = lb - 1 - by 858 s(1, coset(1, by, bz)) = rbp(1)*s(1, coset(0, by, bz)) 859 as(1, coset(1, by, bz), 1) = rbp(1)*as(1, coset(0, by, bz), 1) + & 860 ac1(1)*s(1, coset(0, by, bz)) 861 as(1, coset(1, by, bz), 2) = rbp(1)*as(1, coset(0, by, bz), 2) + & 862 ac1(2)*s(1, coset(0, by, bz)) 863 as(1, coset(1, by, bz), 3) = rbp(1)*as(1, coset(0, by, bz), 3) + & 864 ac1(3)*s(1, coset(0, by, bz)) 865 END DO 866 867 DO bx = 2, lb 868 f3 = f2*REAL(bx - 1, dp) 869 DO by = 0, lb - bx 870 bz = lb - bx - by 871 s(1, coset(bx, by, bz)) = rbp(1)*s(1, coset(bx - 1, by, bz)) + & 872 f3*s(1, coset(bx - 2, by, bz)) 873 as(1, coset(bx, by, bz), 1) = rbp(1)*as(1, coset(bx - 1, by, bz), 1) + & 874 f3*as(1, coset(bx - 2, by, bz), 1) + & 875 ac1(1)*s(1, coset(bx - 1, by, bz)) 876 as(1, coset(bx, by, bz), 2) = rbp(1)*as(1, coset(bx - 1, by, bz), 2) + & 877 f3*as(1, coset(bx - 2, by, bz), 2) + & 878 ac1(2)*s(1, coset(bx - 1, by, bz)) 879 as(1, coset(bx, by, bz), 3) = rbp(1)*as(1, coset(bx - 1, by, bz), 3) + & 880 f3*as(1, coset(bx - 2, by, bz), 3) + & 881 ac1(3)*s(1, coset(bx - 1, by, bz)) 882 END DO 883 END DO 884 885 END DO 886 887 END IF 888 889 END IF 890 891 DO j = 1, ncoset(lb_max) 892 DO i = 1, ncoset(la_max) 893 angab(na + i, nb + j, 1) = as(i, j, 1) 894 angab(na + i, nb + j, 2) = as(i, j, 2) 895 angab(na + i, nb + j, 3) = as(i, j, 3) 896 END DO 897 END DO 898 899 nb = nb + ncoset(lb_max) 900 901 END DO 902 903 na = na + ncoset(la_max) 904 905 END DO 906 907 END SUBROUTINE angmom 908 909END MODULE ai_angmom 910