LCOV - code coverage report
Current view: top level - src/aobasis - ai_angmom.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 349 349 100.0 %
Date: 2025-02-18 08:24:35 Functions: 1 1 100.0 %

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

Generated by: LCOV version 1.15