Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : module m_wann_ujugaunt
8 : use m_juDFT
9 : contains
10 2 : SUBROUTINE wann_ujugaunt(
11 : c***********************************************************************
12 : c wann_ujugaunt calculates integrals of radial wave functions with
13 : c bessel functions and multiplies them with an angular factor.
14 : c Calculating them only once gives some speed-up of wann_mmkb_sph.
15 : c Frank Freimuth, October 2006
16 : c***********************************************************************
17 8 : > llod,nntot,kdiff,lmax,
18 2 : > ntype,ntypd,bbmat,bmat,
19 4 : > nlod,nlo,llo,flo,flo_b,f,f_b,g,g_b,
20 4 : > jri,rmsh,dx,jmtd,
21 : > lmaxd,lmd,
22 2 : < ujug,ujdg,djug,djdg,ujulog,djulog,
23 2 : < ulojug,ulojdg,ulojulog,l_q,sign_q)
24 : use m_constants, only : pimach
25 : use m_sphbes
26 : use m_ylm
27 : use m_intgr, only : intgr3
28 : use m_gaunt, only: gaunt1
29 :
30 : IMPLICIT NONE
31 : integer, intent (in) :: llod
32 : INTEGER, INTENT (IN) :: nntot,ntype,ntypd
33 : INTEGER, INTENT (IN) :: lmaxd,jmtd,lmd
34 : real, intent (in) :: kdiff(:,:) !(3,nntot)
35 : real, intent (in) :: bbmat(:,:) !(3,3)
36 : real, intent (in) :: bmat(:,:) !(3,3)
37 : integer, intent (in) :: lmax(:) !(ntypd)
38 : integer, intent (in) :: nlod
39 : integer, intent (in) :: jri(:) !(ntypd)
40 : integer, intent (in) :: nlo(:) !(ntypd)
41 : integer, intent (in) :: llo(:,:) !(nlod,ntypd)
42 : real, intent (in) :: f(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
43 : real, intent (in) :: f_b(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
44 : real, intent (in) :: g(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
45 : real, intent (in) :: g_b(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
46 : real, intent (in) :: flo(:,:,:,:) !(ntypd,jmtd,2,nlod)
47 : real, intent (in) :: flo_b(:,:,:,:) !(ntypd,jmtd,2,nlod)
48 : real, intent (in) :: rmsh(:,:) !(jmtd,ntypd)
49 : real, intent (in) :: dx(:) !(ntypd)
50 :
51 : logical, intent (in) :: l_q ! if true, we deal with q points
52 : integer, intent (in) :: sign_q ! if we deal with q points, we might pick up an additional sign
53 :
54 : complex, intent (out) :: ujug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
55 : complex, intent (out) :: ujdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
56 : complex, intent (out) :: djug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
57 : complex, intent (out) :: djdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
58 : complex, intent (out) :: ujulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
59 : complex, intent (out) :: djulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
60 : complex, intent (out) :: ulojug(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
61 : complex, intent (out) :: ulojdg(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
62 : complex, intent (out) :: ulojulog(:,-llod:,:,-llod:,:,:) !(1:nlod,-llod:llod,1:nlod,-llod:llod,1:ntype,1:nntot)
63 :
64 2 : real, allocatable :: djd(:,:,:),ujd(:,:,:),uju(:,:,:)
65 2 : real, allocatable :: dju(:,:,:)
66 2 : real, allocatable :: ujulo(:,:,:),djulo(:,:,:),ulojulo(:,:,:)
67 2 : real, allocatable :: uloju(:,:,:),ulojd(:,:,:)
68 : integer :: ikpt_b,i,lwn,n,lpp,lop,lo,l,lp
69 : integer :: lmini,lmaxi,m,mp,llpp,mpp
70 : integer :: lmpp,lminp,lmaxp,lm,lpmp
71 2 : real :: rk,bpt(3),gs,jlpp(0:lmaxd)
72 2 : real :: jj(0:lmaxd,jmtd),x(jmtd)
73 : real :: bkrot(3)
74 2 : complex :: ylmpp((lmaxd+1)**2),factor,ic
75 : complex :: factor2
76 :
77 2 : call timestart("wann_ujugaunt")
78 2 : ic = cmplx(0.,1.)
79 :
80 10 : allocate( djd(0:lmaxd,0:lmaxd,0:lmaxd) )
81 8 : allocate( ujd(0:lmaxd,0:lmaxd,0:lmaxd) )
82 8 : allocate( dju(0:lmaxd,0:lmaxd,0:lmaxd) )
83 8 : allocate( uju(0:lmaxd,0:lmaxd,0:lmaxd) )
84 :
85 10 : allocate( ujulo(nlod,0:lmaxd,0:lmaxd) )
86 8 : allocate( djulo(nlod,0:lmaxd,0:lmaxd) )
87 8 : allocate( uloju(nlod,0:lmaxd,0:lmaxd) )
88 8 : allocate( ulojd(nlod,0:lmaxd,0:lmaxd) )
89 :
90 10 : allocate( ulojulo(nlod,nlod,0:lmaxd) )
91 :
92 39234 : ujug = cmplx(0.0,0.0)
93 39234 : ujdg = cmplx(0.0,0.0)
94 39234 : djug = cmplx(0.0,0.0)
95 39234 : djdg = cmplx(0.0,0.0)
96 :
97 50 : ujulog = cmplx(0.0,0.0)
98 50 : djulog = cmplx(0.0,0.0)
99 50 : ulojug = cmplx(0.0,0.0)
100 50 : ulojdg = cmplx(0.0,0.0)
101 :
102 50 : ulojulog = cmplx(0.0,0.0)
103 :
104 18 : do ikpt_b=1,nntot
105 64 : bpt(:)=kdiff(:,ikpt_b)
106 336 : rk = sqrt(dot_product(bpt,matmul(bbmat,bpt)))
107 : !write(*,*)'ujugaunt rk',rk
108 :
109 34 : do n=1,ntype
110 16 : lwn = lmax(n)
111 : c...generate the j_lpp(br) on the radial grid
112 7552 : do i = 1,jri(n)
113 7536 : gs = rk*rmsh(i,n)
114 7536 : call sphbes(lwn,gs,jlpp)
115 60304 : jj(:,i) = jlpp(:)
116 : enddo
117 128 : do lpp = 0,lwn ! lpp is the ang. momentum of the bessel function
118 : c***************************************************************************
119 : c...the local orbitals overlaps
120 : c***************************************************************************
121 112 : if (nlo(n).GE.1) then
122 0 : do lop = 1,nlo(n)
123 0 : do lo = 1,nlo(n)
124 0 : l = llo(lo,n)
125 0 : lp = llo(lop,n)
126 0 : lmini = abs(lp - l)
127 0 : lmaxi = lp + l
128 : c..the gaunt conditions
129 0 : if ((mod(l+lp+lpp,2).eq.1) .or. (lpp.LT.lmini) .or.
130 0 : + (lpp.gt.lmaxi)) then
131 0 : ulojulo(lo,lop,lpp) = 0.
132 : else
133 0 : do i = 1,jri(n)
134 : x(i) = ( flo(n,i,1,lo)*flo_b(n,i,1,lop)+
135 0 : + flo(n,i,2,lo)*flo_b(n,i,2,lop) )*jj(lpp,i)
136 : enddo
137 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),ulojulo(lo,lop,lpp))
138 : endif
139 : enddo
140 : enddo
141 : endif ! local orbitals
142 : c**************************************************************************
143 : c...overlaps of the apws only
144 : c**************************************************************************
145 912 : do lp = 0,lwn
146 6272 : do l = 0,lwn
147 5488 : lmini = abs(lp-l)
148 5488 : lmaxi = lp + l
149 : c..gaunt conditions
150 5488 : if ((mod(l+lp+lpp,2).eq.1) .or. (lpp.LT.lmini) .or.
151 784 : + (lpp.gt.lmaxi)) then
152 3792 : uju(l,lp,lpp) = 0.
153 3792 : ujd(l,lp,lpp) = 0.
154 3792 : dju(l,lp,lpp) = 0.
155 3792 : djd(l,lp,lpp) = 0.
156 : else
157 800512 : do i = 1,jri(n)
158 : x(i) = ( f(n,i,1,l)*f_b(n,i,1,lp)+
159 800512 : + f(n,i,2,l)*f_b(n,i,2,lp) )*jj(lpp,i)
160 : enddo
161 1696 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),uju(l,lp,lpp))
162 :
163 800512 : do i = 1,jri(n)
164 : x(i) = ( f(n,i,1,l)*g_b(n,i,1,lp)+
165 800512 : + f(n,i,2,l)*g_b(n,i,2,lp) )*jj(lpp,i)
166 : enddo
167 1696 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujd(l,lp,lpp))
168 :
169 800512 : do i = 1,jri(n)
170 : x(i) = ( g(n,i,1,l)*f_b(n,i,1,lp)+
171 800512 : + g(n,i,2,l)*f_b(n,i,2,lp) )*jj(lpp,i)
172 : enddo
173 1696 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),dju(l,lp,lpp))
174 :
175 800512 : do i = 1,jri(n)
176 : x(i) = ( g(n,i,1,l)*g_b(n,i,1,lp)+
177 800512 : + g(n,i,2,l)*g_b(n,i,2,lp) )*jj(lpp,i)
178 : enddo
179 1696 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),djd(l,lp,lpp))
180 : endif
181 : enddo ! l
182 :
183 : c********************************************************************
184 : c...overlaps of the lo's with the apws
185 : c********************************************************************
186 896 : if (nlo(n).GE.1) then
187 0 : do lo = 1,nlo(n)
188 0 : l = llo(lo,n)
189 0 : lmini = abs(lp-l)
190 0 : lmaxi = lp + l
191 : c..gaunt conditions
192 0 : if ((mod(l+lp+lpp,2).eq.1) .OR. (lpp.lt.lmini) .or.
193 0 : + (lpp.gt.lmaxi)) then
194 0 : ujulo(lo,lp,lpp) = 0.
195 0 : djulo(lo,lp,lpp) = 0.
196 0 : uloju(lo,lp,lpp) = 0.
197 0 : ulojd(lo,lp,lpp) = 0.
198 : else
199 0 : do i = 1,jri(n)
200 : x(i) = ( flo(n,i,1,lo)*f_b(n,i,1,lp)+
201 0 : + flo(n,i,2,lo)*f_b(n,i,2,lp) )*jj(lpp,i)
202 : enddo
203 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujulo(lo,lp,lpp))
204 0 : do i = 1,jri(n)
205 : x(i) = ( flo(n,i,1,lo)*g_b(n,i,1,lp)+
206 0 : + flo(n,i,2,lo)*g_b(n,i,2,lp) )*jj(lpp,i)
207 : enddo
208 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),djulo(lo,lp,lpp))
209 0 : do i = 1,jri(n)
210 : x(i) = ( flo_b(n,i,1,lo)*f(n,i,1,lp)+
211 0 : + flo_b(n,i,2,lo)*f(n,i,2,lp) )*jj(lpp,i)
212 : enddo
213 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),uloju(lo,lp,lpp))
214 0 : do i = 1,jri(n)
215 : x(i) = ( flo_b(n,i,1,lo)*g(n,i,1,lp)+
216 0 : + flo_b(n,i,2,lo)*g(n,i,2,lp) )*jj(lpp,i)
217 : enddo
218 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),ulojd(lo,lp,lpp))
219 : endif
220 : enddo !lo
221 : endif ! local orbitals
222 : enddo !lp
223 : enddo !lpp
224 : c********************************************************************
225 : c multiply with gaunt-coefficient (apw-apw)
226 : c********************************************************************
227 208 : bkrot=matmul(bpt,bmat)
228 16 : call ylm4(lwn,bkrot,ylmpp)
229 128 : do l = 0,lwn
230 912 : do m = -l,l
231 784 : lm=l*(l+1)+m
232 6384 : do lp = 0,lwn
233 44688 : do mp = -lp,lp
234 38416 : lpmp=lp*(lp+1)+mp
235 312816 : do lpp = 0,lwn
236 268912 : llpp = lpp*(lpp+1)
237 :
238 268912 : mpp = mp - m
239 :
240 268912 : lmpp = llpp + mpp
241 268912 : lmini = abs(l-lpp)
242 268912 : lmaxi = l+lpp
243 : if ((lmini.le.lp).and.(lp.le.lmaxi).and.
244 307328 : & (mod(l+lp+lpp,2).eq.0).and.(abs(mpp).LE.lpp))then
245 :
246 : factor=conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
247 66208 : * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
248 :
249 : c if(factor.ne.0 .and. uju(l,lp,lpp).ne.0)
250 : c > write(*,*)lpp,lp,l,mp,m
251 :
252 66208 : if(l_q) then
253 0 : factor=(sign_q**lpp)*factor ! additional sign for q points
254 : endif
255 :
256 : ujug(lpmp,lm,n,ikpt_b)=ujug(lpmp,lm,n,ikpt_b)+
257 66208 : + factor*uju(l,lp,lpp)
258 : ujdg(lpmp,lm,n,ikpt_b)=ujdg(lpmp,lm,n,ikpt_b)+
259 66208 : + factor*ujd(l,lp,lpp)
260 : djug(lpmp,lm,n,ikpt_b)=djug(lpmp,lm,n,ikpt_b)+
261 66208 : + factor*dju(l,lp,lpp)
262 : djdg(lpmp,lm,n,ikpt_b)=djdg(lpmp,lm,n,ikpt_b)+
263 66208 : + factor*djd(l,lp,lpp)
264 :
265 : endif
266 : enddo ! lpp
267 : enddo ! mp
268 : enddo ! lp
269 : enddo ! m
270 : enddo ! l
271 : c******************************************************************
272 : c multiply with the gaunt-coefficient (apw-lo)
273 : c******************************************************************
274 32 : if (nlo(n).ge.1) then
275 0 : do lo = 1,nlo(n)
276 0 : l = llo(lo,n)
277 0 : do m = -l,l
278 0 : lm=l*(l+1)+m
279 0 : do lp = 0,lwn
280 0 : do mp = -lp,lp
281 0 : lpmp=lp*(lp+1)+mp
282 0 : do lpp = 0,lwn
283 0 : llpp = lpp*(lpp+1)
284 0 : lmini = abs(l-lpp)
285 0 : lmaxi = l+lpp
286 0 : lminp = abs(lp-lpp)
287 0 : lmaxp = lp+lpp
288 : if ((lmini.le.lp).and.(lp.le.lmaxi).and.
289 0 : & (mod(l+lp+lpp,2).eq.0).and.(abs(mp-m).LE.lpp)) then
290 0 : mpp = mp - m
291 0 : lmpp = llpp + mpp
292 : factor=conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
293 0 : * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
294 0 : if(l_q) then
295 0 : factor=(sign_q**lpp)*factor ! additional sign for q points
296 : endif
297 :
298 : ujulog(lpmp,lo,m,n,ikpt_b)=ujulog(lpmp,lo,m,n,ikpt_b)+
299 0 : + factor*ujulo(lo,lp,lpp)
300 : djulog(lpmp,lo,m,n,ikpt_b)=djulog(lpmp,lo,m,n,ikpt_b)+
301 0 : + factor*djulo(lo,lp,lpp)
302 : endif
303 :
304 : if ((lminp.le.l).and.(l.le.lmaxp).and.
305 0 : & (mod(l+lp+lpp,2).eq.0).and.(abs(m-mp).LE.lpp)) then
306 0 : mpp = m - mp
307 0 : lmpp = llpp + mpp
308 : factor=conjg(ylmpp(lmpp+1))*(ic**(lp+lpp-l))*
309 0 : * gaunt1(l,lpp,lp,m,mpp,mp,lmaxd)
310 0 : if(l_q) then
311 0 : factor=(sign_q**lpp)*factor
312 : endif
313 :
314 : ulojug(lpmp,lo,m,n,ikpt_b)=ulojug(lpmp,lo,m,n,ikpt_b)+
315 0 : + factor*uloju(lo,lp,lpp)
316 : ulojdg(lpmp,lo,m,n,ikpt_b)=ulojdg(lpmp,lo,m,n,ikpt_b)+
317 0 : + factor*ulojd(lo,lp,lpp)
318 : endif
319 : enddo ! lpp
320 : enddo ! mp
321 : enddo ! lp
322 : enddo ! m lo
323 : enddo ! lo
324 : c*************************************************************
325 : c multiply with the gaunt-coefficient (lo-lo)
326 : c*************************************************************
327 0 : do lo = 1,nlo(n)
328 0 : l = llo(lo,n)
329 0 : do m = -l,l
330 0 : lm=l*(l+1)+m
331 0 : do lop = 1,nlo(n)
332 0 : lp = llo(lop,n)
333 0 : do mp = -lp,lp
334 0 : lpmp=lp*(lp+1)+mp
335 0 : do lpp = 0,lwn
336 0 : llpp = lpp*(lpp+1)
337 0 : mpp = mp - m
338 0 : lmpp = llpp + mpp
339 0 : lmini = abs(l-lpp)
340 0 : lmaxi = l+lpp
341 : if ((lmini.le.lp).and.(lp.le.lmaxi).and.
342 0 : & (mod(l+lp+lpp,2).eq.0).and.(abs(mpp).LE.lpp))then
343 :
344 : factor= conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
345 0 : * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
346 0 : if(l_q) then
347 0 : factor=(sign_q**lpp)*factor ! additional sign for q points
348 : endif
349 :
350 : ulojulog(lop,mp,lo,m,n,ikpt_b)=
351 : = ulojulog(lop,mp,lo,m,n,ikpt_b)+
352 0 : + ulojulo(lo,lop,lpp)*factor
353 : ! + conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
354 : ! * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)*
355 : ! * ulojulo(lo,lop,lpp)
356 : endif
357 : enddo ! lpp
358 : enddo ! mp lop
359 : enddo ! lop
360 : enddo ! m lo
361 : enddo ! lo
362 : endif ! local orbitals on this atom
363 : enddo !ntype
364 : enddo !ikpt_b
365 2 : deallocate(djd,ujd,uju,ujulo,djulo,ulojulo)
366 2 : deallocate(dju,uloju,ulojd)
367 :
368 2 : call timestop("wann_ujugaunt")
369 2 : end subroutine wann_ujugaunt
370 16 : end module m_wann_ujugaunt
|
