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
|
