gaunt.f90 3.77 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1 2
MODULE m_gaunt
!*********************************************************************
Matthias Redies's avatar
Matthias Redies committed
3
!     Modified module to include old gaunt_init subroutine
Matthias Redies's avatar
Matthias Redies committed
4 5 6 7 8 9
!     the private arrays are allocated and computed in the first call to gaunt1
!                                            Daniel Wortmann
!*********************************************************************
   PRIVATE
   INTEGER,SAVE         :: lmaxdp
   REAL,SAVE,ALLOCATABLE::w(:),yr(:,:)
Matthias Redies's avatar
Matthias Redies committed
10
   PUBLIC gaunt1, gaunt_init
Matthias Redies's avatar
Matthias Redies committed
11
CONTAINS
Matthias Redies's avatar
Matthias Redies committed
12
   FUNCTION gaunt1(lp,l,ls,mp,m,ms,lmaxd)
Matthias Redies's avatar
Matthias Redies committed
13 14 15 16 17 18 19 20 21 22 23
!*********************************************************************
!     gaunt computes the integral of conjg(y(lp,mp))*y(l,m)*y(ls,ms)
!     for lp+l+ls .lt. 2*ngntd
!     using gaussian quadrature as given by
!     m. abramowitz and i.a. stegun, handbook of mathematical functions,
!     nbs applied mathematics series 55 (1968), pages 887 and 916
!     m. weinert and e. wimmer
!     northwestern university march 1980
!     modified to use calculated points and weights
!     to make it dynamic.   (m.w.  jan. 1982)
!*********************************************************************
Matthias Redies's avatar
Matthias Redies committed
24
      USE m_judft
Matthias Redies's avatar
Matthias Redies committed
25
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
26 27
      INTEGER,INTENT(IN) :: l,lp,ls,m,mp,ms,lmaxd
      REAL               :: gaunt1
Matthias Redies's avatar
Matthias Redies committed
28 29
      INTEGER :: i,il,ilp,ils,n

30 31

      n= (3*lmaxd)/4+1
Matthias Redies's avatar
Matthias Redies committed
32
! heck if this is first call to subroutine
Matthias Redies's avatar
Matthias Redies committed
33
      IF(.NOT. ALLOCATED(YR)) CALL gaunt_init(lmaxd)
Matthias Redies's avatar
Matthias Redies committed
34
! heck if the previous call of the subroutine was with the same lmaxd
Matthias Redies's avatar
Matthias Redies committed
35
      IF(lmaxd > lmaxdp) call juDFT_error("Can't calc gaunt. lmaxd too high")
Matthias Redies's avatar
Matthias Redies committed
36

Matthias Redies's avatar
Matthias Redies committed
37
      gaunt1 = 0.0
Matthias Redies's avatar
Matthias Redies committed
38 39 40 41 42
      IF (mp /= (m+ms)) RETURN
      IF (MOD((l+lp+ls),2) == 1) RETURN
      IF ((l+lp-ls) < 0) RETURN
      IF ((l-lp+ls) < 0) RETURN
      IF ((lp-l+ls) < 0) RETURN
43 44 45 46 47 48

      il  = l  * (l  + 1) + m  + 1
      ilp = lp * (lp + 1) + mp + 1
      ils = ls * (ls + 1) + ms + 1

      gaunt1 = dot_product(w, yr(:,ilp)*yr(:,il)*yr(:,ils))
Matthias Redies's avatar
Matthias Redies committed
49 50 51
   END FUNCTION

!     private subroutine for initializing the private arrays!
Matthias Redies's avatar
Matthias Redies committed
52
   SUBROUTINE gaunt_init(lmaxd)
Matthias Redies's avatar
Matthias Redies committed
53 54 55 56 57 58 59 60 61 62 63
!**********************************************************************
!     sets up values needed for gaunt1
!        m. weinert  january 1982
!**********************************************************************
      USE m_constants, ONLY : pimach
      USE m_grule
      USE m_juDFT_stop
      IMPLICIT NONE

      INTEGER, INTENT (IN)  :: lmaxd
      REAL :: a,cd,cth,fac,fpi,rf,sgm,sth,t
Matthias Redies's avatar
Matthias Redies committed
64
      INTEGER :: k,l,lm,lomax,m
Matthias Redies's avatar
Matthias Redies committed
65 66 67
      INTEGER :: n,lmax1d
      REAL :: p(0:lmaxd+1,0:lmaxd+1),x((3*lmaxd)/4+1)

68
      if (allocated(w)) return
Matthias Redies's avatar
Matthias Redies committed
69
      n = (3*lmaxd)/4+1
70 71
      ALLOCATE(w(n),  source=0.0)
      ALLOCATE(yr(n,(lmaxd+1)**2), source=0.0)
Matthias Redies's avatar
Matthias Redies committed
72 73 74 75 76 77 78
      lmaxdp = lmaxd
      lmax1d = lmaxd+1

      fpi = 4.0 * pimach()
      rf = fpi** (1./3.)
      lomax = lmax1d - 1
!--->    obtain gauss-legendre points and weights
Matthias Redies's avatar
Matthias Redies committed
79
      CALL grule(2*n,x,w)
Matthias Redies's avatar
Matthias Redies committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
!--->    generate associated legendre functions for m.ge.0
      DO  k = 1,n
         cth = x(k)
         sth = sqrt(1.0-cth*cth)
         fac = 1.0
         !--->    loop over m values
         DO  m = 0,lomax
            fac = - (2*m-1)*fac
            p(m,m) = fac
            p(m+1,m) = (2*m+1)*cth*fac
            !--->    recurse upward in l
            DO  l = m + 2,lomax
               p(l,m) = ((2*l-1)*cth*p(l-1,m)- (l+m-1)*p(l-2,m))/ (l-m)
            ENDDO
            fac = fac*sth
         ENDDO
         !--->    multiply in the normalization factors
         DO  l = 0,lomax
            a = rf*sqrt((2*l+1)/fpi)
            cd = 1
            lm = l* (l+1) + 1
            yr(k,lm) = a*p(l,0)
            sgm = -1.
            DO  m = 1,l
               t = (l+1-m)* (l+m)
               cd = cd/t
               yr(k,lm+m) = a*sqrt(cd)*p(l,m)
               yr(k,lm-m) = sgm*a*sqrt(cd)*p(l,m)
               sgm = -sgm
            ENDDO
         ENDDO
      ENDDO
   END SUBROUTINE
END MODULE