LCOV - code coverage report
Current view: top level - xc-pot - gaunt.f90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 100.0 % 98 98
Test Date: 2025-07-19 04:34:44 Functions: 100.0 % 4 4

            Line data    Source code
       1              : MODULE m_gaunt
       2              : !*********************************************************************
       3              : !     Modified module to include old gaunt_init subroutine
       4              : !     the private arrays are allocated and computed in the first call to gaunt1
       5              : !                                            Daniel Wortmann
       6              : !*********************************************************************
       7              :    PRIVATE
       8              :    INTEGER,SAVE         :: lmaxdp,lmaxdp2
       9              :    REAL,SAVE,ALLOCATABLE::w(:),yr(:,:)
      10              :    REAL,SAVE,ALLOCATABLE::w2(:),yr2(:,:)
      11              :    PUBLIC gaunt1, gaunt2, gaunt_init, gaunt_init2
      12              : CONTAINS
      13     90703136 :    FUNCTION gaunt1(lp,l,ls,mp,m,ms,lmaxd)
      14              :       !! This is a test for the new FORD documentation.
      15              :       !! In this function we calculate the Gaunt coefficients
      16              :       !! $$ G_{l',l,l''}^{m',m,m''}\equiv\int d\Omega ~ Y_{l'}^{m'*}Y_{l}^{m}Y_{l''}^{m''} $$
      17              : 
      18              : !*********************************************************************
      19              : !     gaunt computes the integral of conjg(y(lp,mp))*y(l,m)*y(ls,ms)
      20              : !     for lp+l+ls .lt. 2*ngntd
      21              : !     using gaussian quadrature as given by
      22              : !     m. abramowitz and i.a. stegun, handbook of mathematical functions,
      23              : !     nbs applied mathematics series 55 (1968), pages 887 and 916
      24              : !     m. weinert and e. wimmer
      25              : !     northwestern university march 1980
      26              : !     modified to use calculated points and weights
      27              : !     to make it dynamic.   (m.w.  jan. 1982)
      28              : !*********************************************************************
      29              :       USE m_judft
      30              :       IMPLICIT NONE
      31              :       INTEGER,INTENT(IN) :: l,lp,ls,m,mp,ms,lmaxd
      32              :       REAL               :: gaunt1
      33              :       INTEGER :: i,il,ilp,ils,n
      34              : 
      35              : 
      36     90703136 :       n= (3*lmaxd)/4+1
      37              : ! heck if this is first call to subroutine
      38     90703136 :       IF(.NOT. ALLOCATED(YR)) CALL gaunt_init(lmaxd)
      39              : ! heck if the previous call of the subroutine was with the same lmaxd
      40     90703136 :       IF(lmaxd > lmaxdp) call juDFT_error("Can't calc gaunt. lmaxd too high")
      41              : 
      42     90703136 :       gaunt1 = 0.0
      43     90703136 :       IF (mp /= (m+ms)) RETURN
      44     89585508 :       IF (MOD((l+lp+ls),2) == 1) RETURN
      45     89548762 :       IF ((l+lp-ls) < 0) RETURN
      46     89537958 :       IF ((l-lp+ls) < 0) RETURN
      47     89537510 :       IF ((lp-l+ls) < 0) RETURN
      48              : 
      49     89527154 :       il  = l  * (l  + 1) + m  + 1
      50     89527154 :       ilp = lp * (lp + 1) + mp + 1
      51     89527154 :       ils = ls * (ls + 1) + ms + 1
      52              : 
      53    756458918 :       gaunt1 = dot_product(w, yr(:,ilp)*yr(:,il)*yr(:,ils))
      54              :    END FUNCTION
      55              : 
      56        38961 :    FUNCTION gaunt2(lp,l,ls,mp,m,ms,lmaxd)
      57              :       USE m_judft
      58              :       IMPLICIT NONE
      59              :       INTEGER,INTENT(IN) :: l,lp,ls,m,mp,ms,lmaxd
      60              :       REAL               :: gaunt2
      61              :       INTEGER :: i,il,ilp,ils,n
      62              : 
      63              : 
      64        38961 :       n= (3*lmaxd)/4+1
      65              : ! heck if this is first call to subroutine
      66        38961 :       IF(.NOT. ALLOCATED(YR2)) CALL gaunt_init2(lmaxd)
      67              : ! heck if the previous call of the subroutine was with the same lmaxd
      68        38961 :       IF(lmaxd > lmaxdp2) call juDFT_error("Can't calc gaunt. lmaxd too high")
      69              : 
      70        38961 :       gaunt2 = 0.0
      71        38961 :       IF (mp /= (m+ms)) RETURN
      72        38961 :       IF (MOD((l+lp+ls),2) == 1) RETURN
      73        38961 :       IF ((l+lp-ls) < 0) RETURN
      74        38961 :       IF ((l-lp+ls) < 0) RETURN
      75        38961 :       IF ((lp-l+ls) < 0) RETURN
      76              : 
      77        38961 :       il  = l  * (l  + 1) + m  + 1
      78        38961 :       ilp = lp * (lp + 1) + mp + 1
      79        38961 :       ils = ls * (ls + 1) + ms + 1
      80              : 
      81       545454 :       gaunt2 = dot_product(w2, yr2(:,ilp)*yr2(:,il)*yr2(:,ils))
      82              :    END FUNCTION
      83              : 
      84              : !     private subroutine for initializing the private arrays!
      85          162 :    SUBROUTINE gaunt_init(lmaxd)
      86              : !**********************************************************************
      87              : !     sets up values needed for gaunt1
      88              : !        m. weinert  january 1982
      89              : !**********************************************************************
      90              :       USE m_constants, ONLY : pimach
      91              :       USE m_grule
      92              :       USE m_juDFT_stop
      93              :       IMPLICIT NONE
      94              : 
      95              :       INTEGER, INTENT (IN)  :: lmaxd
      96              :       REAL :: a,cd,cth,fac,fpi,rf,sgm,sth,t
      97              :       INTEGER :: k,l,lm,lomax,m
      98              :       INTEGER :: n,lmax1d
      99          162 :       REAL :: p(0:lmaxd+1,0:lmaxd+1),x((3*lmaxd)/4+1)
     100              : 
     101          162 :       if (allocated(w)) return
     102          162 :       n = (3*lmaxd)/4+1
     103         1408 :       ALLOCATE(w(n),  source=0.0)
     104       172410 :       ALLOCATE(yr(n,(lmaxd+1)**2), source=0.0)
     105          162 :       lmaxdp = lmaxd
     106          162 :       lmax1d = lmaxd+1
     107              : 
     108          162 :       fpi = 4.0 * pimach()
     109          162 :       rf = fpi** (1./3.)
     110          162 :       lomax = lmax1d - 1
     111              : !--->    obtain gauss-legendre points and weights
     112          162 :       CALL grule(2*n,x,w)
     113              : !--->    generate associated legendre functions for m.ge.0
     114         1408 :       DO  k = 1,n
     115         1246 :          cth = x(k)
     116         1246 :          sth = sqrt(1.0-cth*cth)
     117         1246 :          fac = 1.0
     118              :          !--->    loop over m values
     119        14918 :          DO  m = 0,lomax
     120        13672 :             fac = - (2*m-1)*fac
     121        13672 :             p(m,m) = fac
     122        13672 :             p(m+1,m) = (2*m+1)*cth*fac
     123              :             !--->    recurse upward in l
     124        71002 :             DO  l = m + 2,lomax
     125        71002 :                p(l,m) = ((2*l-1)*cth*p(l-1,m)- (l+m-1)*p(l-2,m))/ (l-m)
     126              :             ENDDO
     127        14918 :             fac = fac*sth
     128              :          ENDDO
     129              :          !--->    multiply in the normalization factors
     130        15080 :          DO  l = 0,lomax
     131        13672 :             a = rf*sqrt((2*l+1)/fpi)
     132        13672 :             cd = 1
     133        13672 :             lm = l* (l+1) + 1
     134        13672 :             yr(k,lm) = a*p(l,0)
     135        13672 :             sgm = -1.
     136        84674 :             DO  m = 1,l
     137        69756 :                t = (l+1-m)* (l+m)
     138        69756 :                cd = cd/t
     139        69756 :                yr(k,lm+m) = a*sqrt(cd)*p(l,m)
     140        69756 :                yr(k,lm-m) = sgm*a*sqrt(cd)*p(l,m)
     141        83428 :                sgm = -sgm
     142              :             ENDDO
     143              :          ENDDO
     144              :       ENDDO
     145              :    END SUBROUTINE
     146              : 
     147            2 :    SUBROUTINE gaunt_init2(lmaxd)
     148              :       USE m_constants, ONLY : pimach
     149              :       USE m_grule
     150              :       USE m_juDFT_stop
     151              :       IMPLICIT NONE
     152              : 
     153              :       INTEGER, INTENT (IN)  :: lmaxd
     154              :       REAL :: a,cd,cth,fac,fpi,rf,sgm,sth,t
     155              :       INTEGER :: k,l,lm,lomax,m
     156              :       INTEGER :: n,lmax1d
     157            2 :       REAL :: p(0:lmaxd+1,0:lmaxd+1),x((3*lmaxd)/4+1)
     158              : 
     159            2 :       if (allocated(w2)) return
     160            2 :       n = (3*lmaxd)/4+1
     161           28 :       ALLOCATE(w2(n),  source=0.0)
     162         8094 :       ALLOCATE(yr2(n,(lmaxd+1)**2), source=0.0)
     163            2 :       lmaxdp2 = lmaxd
     164            2 :       lmax1d = lmaxd+1
     165              : 
     166            2 :       fpi = 4.0 * pimach()
     167            2 :       rf = fpi** (1./3.)
     168            2 :       lomax = lmax1d - 1
     169              : !--->    obtain gauss-legendre points and weights
     170            2 :       CALL grule(2*n,x,w2)
     171              : !--->    generate associated legendre functions for m.ge.0
     172           28 :       DO  k = 1,n
     173           26 :          cth = x(k)
     174           26 :          sth = sqrt(1.0-cth*cth)
     175           26 :          fac = 1.0
     176              :          !--->    loop over m values
     177          468 :          DO  m = 0,lomax
     178          442 :             fac = - (2*m-1)*fac
     179          442 :             p(m,m) = fac
     180          442 :             p(m+1,m) = (2*m+1)*cth*fac
     181              :             !--->    recurse upward in l
     182         3562 :             DO  l = m + 2,lomax
     183         3562 :                p(l,m) = ((2*l-1)*cth*p(l-1,m)- (l+m-1)*p(l-2,m))/ (l-m)
     184              :             ENDDO
     185          468 :             fac = fac*sth
     186              :          ENDDO
     187              :          !--->    multiply in the normalization factors
     188          470 :          DO  l = 0,lomax
     189          442 :             a = rf*sqrt((2*l+1)/fpi)
     190          442 :             cd = 1
     191          442 :             lm = l* (l+1) + 1
     192          442 :             yr2(k,lm) = a*p(l,0)
     193          442 :             sgm = -1.
     194         4004 :             DO  m = 1,l
     195         3536 :                t = (l+1-m)* (l+m)
     196         3536 :                cd = cd/t
     197         3536 :                yr2(k,lm+m) = a*sqrt(cd)*p(l,m)
     198         3536 :                yr2(k,lm-m) = sgm*a*sqrt(cd)*p(l,m)
     199         3978 :                sgm = -sgm
     200              :             ENDDO
     201              :          ENDDO
     202              :       ENDDO
     203              :    END SUBROUTINE
     204              : END MODULE
        

Generated by: LCOV version 2.0-1