LCOV - code coverage report
Current view: top level - xc-pot - gaunt.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 51 53 96.2 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             : MODULE m_gaunt
       2             : !*********************************************************************
       3             : !     Modified module to include old gaunt2 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
       9             :    REAL,SAVE,ALLOCATABLE::w(:),yr(:,:)
      10             :    PUBLIC gaunt1,gaunt2
      11             : CONTAINS
      12     9666420 :    REAL FUNCTION gaunt1(lp,l,ls,mp,m,ms,lmaxd)
      13             : !*********************************************************************
      14             : !     gaunt computes the integral of conjg(y(lp,mp))*y(l,m)*y(ls,ms)
      15             : !     for lp+l+ls .lt. 2*ngntd
      16             : !     using gaussian quadrature as given by
      17             : !     m. abramowitz and i.a. stegun, handbook of mathematical functions,
      18             : !     nbs applied mathematics series 55 (1968), pages 887 and 916
      19             : !     m. weinert and e. wimmer
      20             : !     northwestern university march 1980
      21             : !     modified to use calculated points and weights
      22             : !     to make it dynamic.   (m.w.  jan. 1982)
      23             : !*********************************************************************
      24             :       IMPLICIT NONE
      25             : !     ..
      26             : !     .. Scalar Arguments ..
      27             :       INTEGER,INTENT(IN):: l,lp,ls,m,mp,ms,lmaxd
      28             : !     ..
      29             : !     .. Local Scalars ..
      30             :       REAL :: zero
      31             :       INTEGER :: i,il,ilp,ils,n
      32             : !     ..
      33             : !     .. Intrinsic Functions ..
      34             :       INTRINSIC mod
      35             : !     ..
      36             : !     .. Data statements ..
      37             :       DATA zero/0.0e0/
      38             : !     ..
      39     9666420 :       n= (3*lmaxd)/4+1
      40             : 
      41             : ! heck if this is first call to subroutine
      42     9666420 :       IF ( .NOT. ALLOCATED(YR)) CALL gaunt2(lmaxd)
      43             : ! heck if the previous call of the subroutine was with the same lmaxd
      44     9666420 :       IF( lmaxd /= lmaxdp ) THEN
      45           0 :          DEALLOCATE(yr,w)
      46           0 :          CALL gaunt2(lmaxd)
      47             :       END IF
      48             : 
      49     9666420 :       gaunt1 = zero
      50     9666420 :       IF (mp /= (m+ms)) RETURN
      51     9440444 :       IF (MOD((l+lp+ls),2) == 1) RETURN
      52     9432436 :       IF ((l+lp-ls) < 0) RETURN
      53     9430772 :       IF ((l-lp+ls) < 0) RETURN
      54     9430772 :       IF ((lp-l+ls) < 0) RETURN
      55     9429108 :       il = l* (l+1) + m + 1
      56     9429108 :       ilp = lp* (lp+1) + mp + 1
      57     9429108 :       ils = ls* (ls+1) + ms + 1
      58    80536516 :       DO i = 1,n
      59    80536516 :          gaunt1 = gaunt1 + w(i)*yr(i,ilp)*yr(i,il)*yr(i,ils)
      60             :       END DO
      61             :       RETURN
      62             :    END FUNCTION
      63             : 
      64             : !     private subroutine for initializing the private arrays!
      65         608 :    SUBROUTINE gaunt2( &
      66             :       lmaxd)
      67             : !**********************************************************************
      68             : !     sets up values needed for gaunt1
      69             : !        m. weinert  january 1982
      70             : !**********************************************************************
      71             :       USE m_constants, ONLY : pimach
      72             :       USE m_grule
      73             :       USE m_juDFT_stop
      74             : !$    USE omp_lib
      75             :       IMPLICIT NONE
      76             : 
      77             :       INTEGER, INTENT (IN)  :: lmaxd
      78             : !     ..
      79             : !     .. Local Scalars ..
      80             :       REAL :: a,cd,cth,fac,fpi,rf,sgm,sth,t
      81             :       INTEGER :: k,l,lm,lomax,m,nn
      82             :       INTEGER :: n,lmax1d
      83             : !     ..
      84             : !     .. Local Arrays ..
      85        1216 :       REAL :: p(0:lmaxd+1,0:lmaxd+1),x((3*lmaxd)/4+1)
      86             : !     ..
      87             : !     .. Intrinsic Functions ..
      88             :       INTRINSIC sqrt
      89             : !     ..
      90         608 :       if (allocated(w)) return
      91          66 : !$    if (omp_in_parallel() .and. omp_get_num_threads() > 1) call juDFT_error("BUG IN GAUNT!!")
      92          66 :       ALLOCATE(w((3*lmaxd)/4+1),yr((3*lmaxd)/4+1,(lmaxd+1)**2))
      93             : 
      94          66 :       n = (3*lmaxd)/4+1
      95          66 :       lmaxdp = lmaxd
      96          66 :       lmax1d = lmaxd+1
      97             : 
      98          66 :       fpi = 4.0 * pimach()
      99          66 :       rf = fpi** (1./3.)
     100          66 :       lomax = lmax1d - 1
     101             : !--->    obtain gauss-legendre points and weights
     102          66 :       nn = 2*n
     103          66 :       CALL grule(nn,x,w)
     104             : !--->    generate associated legendre functions for m.ge.0
     105         612 :       DO  k = 1,n
     106         546 :          cth = x(k)
     107         546 :          sth = sqrt(1.0-cth*cth)
     108         546 :          fac = 1.0
     109             :          !--->    loop over m values
     110        6596 :          DO  m = 0,lomax
     111        6050 :             fac = - (2*m-1)*fac
     112        6050 :             p(m,m) = fac
     113        6050 :             p(m+1,m) = (2*m+1)*cth*fac
     114             :             !--->    recurse upward in l
     115       32034 :             DO  l = m + 2,lomax
     116       32034 :                p(l,m) = ((2*l-1)*cth*p(l-1,m)- (l+m-1)*p(l-2,m))/ (l-m)
     117             :             ENDDO
     118        6596 :             fac = fac*sth
     119             :          ENDDO
     120             :          !--->    multiply in the normalization factors
     121        6662 :          DO  l = 0,lomax
     122        6050 :             a = rf*sqrt((2*l+1)/fpi)
     123        6050 :             cd = 1
     124        6050 :             lm = l* (l+1) + 1
     125        6050 :             yr(k,lm) = a*p(l,0)
     126        6050 :             sgm = -1.
     127       37538 :             DO  m = 1,l
     128       31488 :                t = (l+1-m)* (l+m)
     129       31488 :                cd = cd/t
     130       31488 :                yr(k,lm+m) = a*sqrt(cd)*p(l,m)
     131       31488 :                yr(k,lm-m) = sgm*a*sqrt(cd)*p(l,m)
     132       37538 :                sgm = -sgm
     133             :             ENDDO
     134             :          ENDDO
     135             :       ENDDO
     136             :       RETURN
     137             :    END SUBROUTINE
     138             : 
     139             : END MODULE

Generated by: LCOV version 1.13