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 90720185 : 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 90720185 : n= (3*lmaxd)/4+1
37 : ! heck if this is first call to subroutine
38 90720185 : IF(.NOT. ALLOCATED(YR)) CALL gaunt_init(lmaxd)
39 : ! heck if the previous call of the subroutine was with the same lmaxd
40 90720185 : IF(lmaxd > lmaxdp) call juDFT_error("Can't calc gaunt. lmaxd too high")
41 :
42 90720185 : gaunt1 = 0.0
43 90720185 : IF (mp /= (m+ms)) RETURN
44 89635173 : IF (MOD((l+lp+ls),2) == 1) RETURN
45 89600107 : IF ((l+lp-ls) < 0) RETURN
46 89589539 : IF ((l-lp+ls) < 0) RETURN
47 89589539 : IF ((lp-l+ls) < 0) RETURN
48 :
49 89578971 : il = l * (l + 1) + m + 1
50 89578971 : ilp = lp * (lp + 1) + mp + 1
51 89578971 : ils = ls * (ls + 1) + ms + 1
52 :
53 756929822 : 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 160 : 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 160 : REAL :: p(0:lmaxd+1,0:lmaxd+1),x((3*lmaxd)/4+1)
100 :
101 160 : if (allocated(w)) return
102 160 : n = (3*lmaxd)/4+1
103 1714 : ALLOCATE(w(n), source=0.0)
104 171992 : ALLOCATE(yr(n,(lmaxd+1)**2), source=0.0)
105 160 : lmaxdp = lmaxd
106 160 : lmax1d = lmaxd+1
107 :
108 160 : fpi = 4.0 * pimach()
109 160 : rf = fpi** (1./3.)
110 160 : lomax = lmax1d - 1
111 : !---> obtain gauss-legendre points and weights
112 160 : CALL grule(2*n,x,w)
113 : !---> generate associated legendre functions for m.ge.0
114 1394 : DO k = 1,n
115 1234 : cth = x(k)
116 1234 : sth = sqrt(1.0-cth*cth)
117 1234 : fac = 1.0
118 : !---> loop over m values
119 14810 : DO m = 0,lomax
120 13576 : fac = - (2*m-1)*fac
121 13576 : p(m,m) = fac
122 13576 : p(m+1,m) = (2*m+1)*cth*fac
123 : !---> recurse upward in l
124 70654 : DO l = m + 2,lomax
125 70654 : p(l,m) = ((2*l-1)*cth*p(l-1,m)- (l+m-1)*p(l-2,m))/ (l-m)
126 : ENDDO
127 14810 : fac = fac*sth
128 : ENDDO
129 : !---> multiply in the normalization factors
130 14970 : DO l = 0,lomax
131 13576 : a = rf*sqrt((2*l+1)/fpi)
132 13576 : cd = 1
133 13576 : lm = l* (l+1) + 1
134 13576 : yr(k,lm) = a*p(l,0)
135 13576 : sgm = -1.
136 84230 : DO m = 1,l
137 69420 : t = (l+1-m)* (l+m)
138 69420 : cd = cd/t
139 69420 : yr(k,lm+m) = a*sqrt(cd)*p(l,m)
140 69420 : yr(k,lm-m) = sgm*a*sqrt(cd)*p(l,m)
141 82996 : 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 32 : ALLOCATE(w2(n), source=0.0)
162 8100 : 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
|
