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_gabeffgagaunt
8 : contains
9 0 : SUBROUTINE wann_gabeffgagaunt(
10 : > vTot,
11 : > l_perpmagatlres,perpmagl,
12 0 : > neq,mlh,nlhd,nlh,ntypsy,llh,lmax,
13 0 : > nmem,ntype,ntypd,bbmat,bmat,
14 0 : > nlod,llod,nlo,llo,flo,f,g,jri,rmsh,dx,jmtd,
15 0 : > lmaxd,lmd,clnu,
16 0 : < ujug,ujdg,djug,djdg,ujulog,djulog,
17 0 : < ulojug,ulojdg,ulojulog)
18 : c*************************************************************************
19 : c wann_gabeffgagaunt calculates integrals of radial wave functions with
20 : c the exchange field and multiplies them with an angular factor.
21 : c
22 : c Frank Freimuth, 2010
23 : c*************************************************************************
24 : USE m_types
25 : use m_constants
26 : use m_sphbes
27 : use m_ylm
28 : use m_intgr, only : intgr3
29 : use m_gaunt, only: gaunt1
30 :
31 : IMPLICIT NONE
32 : TYPE(t_potden), INTENT(IN) :: vTot
33 : logical, intent (in) :: l_perpmagatlres
34 : integer, intent (in) :: perpmagl
35 : integer, intent (in) :: neq(:)
36 : INTEGER, INTENT (IN) :: mlh(:,0:,:)!(memd,0:nlhd,ntypsd)
37 : INTEGER, INTENT (IN) :: nlhd
38 : integer, intent (in) :: nlh(:)
39 : integer, intent (in) :: ntypsy(:)
40 : integer, intent (in) :: nmem(0:,:)
41 : integer, intent (in) :: ntype,ntypd
42 : integer, intent (in) :: llh(0:,:)
43 : INTEGER, INTENT (IN) :: lmaxd,jmtd,lmd
44 : real, intent (in) :: bbmat(:,:)!(3,3)
45 : real, intent (in) :: bmat(:,:)!(3,3)
46 : integer, intent (in) :: lmax(:)!(ntypd)
47 : integer, intent (in) :: nlod,llod
48 : integer, intent (in) :: jri(:)!(ntypd)
49 : integer, intent (in) :: nlo(:)!(ntypd)
50 : integer, intent (in) :: llo(:,:)!(nlod,ntypd)
51 : real, intent (in) :: f(:,:,:,0:,:)!(ntypd,jmtd,2,0:lmaxd,2)
52 : real, intent (in) :: g(:,:,:,0:,:)!(ntypd,jmtd,2,0:lmaxd,2)
53 : real, intent (in) :: flo(:,:,:,:,:)!(ntypd,jmtd,2,nlod,2)
54 : real, intent (in) :: rmsh(:,:)!(jmtd,ntypd)
55 : real, intent (in) :: dx(:)!(ntypd)
56 : complex, intent (in) :: clnu(:,0:,:)!(memd,0:nlhd,ntypsd)
57 :
58 : complex, intent (out) :: ujug(0:,0:,1:)!(0:lmd,0:lmd,1:ntype)
59 : complex, intent (out) :: ujdg(0:,0:,1:)!(0:lmd,0:lmd,1:ntype)
60 : complex, intent (out) :: djug(0:,0:,1:)!(0:lmd,0:lmd,1:ntype)
61 : complex, intent (out) :: djdg(0:,0:,1:)!(0:lmd,0:lmd,1:ntype)
62 :
63 : complex, intent (out) :: ujulog(0:,:,-llod:,:) !(0:lmd,nlod,-llod:llod,1:ntype)
64 : complex, intent (out) :: djulog(0:,:,-llod:,:) !(0:lmd,nlod,-llod:llod,1:ntype)
65 : complex, intent (out) :: ulojug(0:,:,-llod:,:) !(0:lmd,nlod,-llod:llod,1:ntype)
66 : complex, intent (out) :: ulojdg(0:,:,-llod:,:) !(0:lmd,nlod,-llod:llod,1:ntype)
67 : complex, intent (out) :: ulojulog(:,-llod:,:,-llod:,:) !(1:nlod,-llod:llod,1:nlod,-llod:llod,1:ntype)
68 :
69 0 : real, allocatable :: djd(:,:,:),ujd(:,:,:),uju(:,:,:)
70 0 : real, allocatable :: dju(:,:,:),loju(:,:,:),lojd(:,:,:)
71 0 : real, allocatable :: ujulo(:,:,:),djulo(:,:,:),ulojulo(:,:,:)
72 : integer :: i,lwn,n,lpp,lop,lo,l,lp,lmini,lmaxi,m,mp,llpp,mpp
73 : integer :: lmpp,lminp,lmaxp,lm,lpmp,na,ns,lh,j,jmem
74 0 : real :: rk,gs,jj(0:lmaxd,jmtd),x(jmtd)
75 : complex :: factor,ic
76 : logical :: l_doit
77 0 : real,allocatable :: beff(:,:,:)
78 :
79 0 : call timestart("wann_gabeffgagaunt")
80 0 : ic = cmplx(0.,1.)
81 :
82 0 : allocate( beff(jmtd,0:nlhd,ntype) )
83 :
84 : if(.false.)then !from file
85 : open(777,file='beff',form='unformatted')
86 : read(777)beff
87 : close(777)
88 : else !from vTot
89 0 : beff=(vTot%mt(:,:,:,1)-vTot%mt(:,:,:,2))/2.0
90 : ! In vgen/vgen_finalize.F90 there is the line
91 : ! vTot%mt(:atoms%jri(n),0,n,js) = atoms%rmsh(:atoms%jri(n),n)*vTot%mt(:atoms%jri(n),0,n,js)/sfp_const
92 : ! The following line removes this factor:
93 0 : do n=1,ntype
94 0 : beff(:jri(n),0,n)=beff(:jri(n),0,n)*sfp_const/rmsh(:jri(n),n)
95 : enddo
96 : endif
97 :
98 0 : if(l_perpmagatlres)then
99 : na=1
100 0 : do n=1,ntype
101 0 : ns=ntypsy(na)
102 0 : do lh = 0,nlh(ns)
103 0 : lpp=llh(lh,ns)
104 0 : if(lpp.ne.perpmagl)then
105 0 : beff(:,lh,n)=0.0
106 : endif
107 : enddo !lh
108 0 : na=na+neq(n)
109 : enddo !n
110 : endif !l_perpmagatlres
111 :
112 : allocate ( djd(0:lmaxd,0:lmaxd,0:nlhd),
113 : + dju(0:lmaxd,0:lmaxd,0:nlhd),
114 : + ujd(0:lmaxd,0:lmaxd,0:nlhd),
115 : + uju(0:lmaxd,0:lmaxd,0:nlhd),
116 :
117 : + ujulo(nlod,0:lmaxd,0:nlhd),
118 : + djulo(nlod,0:lmaxd,0:nlhd),
119 : + loju(nlod,0:lmaxd,0:nlhd),
120 : + lojd(nlod,0:lmaxd,0:nlhd),
121 0 : + ulojulo(nlod,nlod,0:nlhd) )
122 :
123 :
124 0 : na=1
125 0 : do n=1,ntype
126 0 : ns=ntypsy(na)
127 0 : lwn = lmax(n)
128 0 : do lh = 0,nlh(ns)
129 0 : lpp=llh(lh,ns)
130 : c***************************************************************************
131 : c...the local orbitals overlaps
132 : c***************************************************************************
133 0 : if (nlo(n).GE.1) then
134 0 : do lop = 1,nlo(n)
135 0 : do lo = 1,nlo(n)
136 0 : l = llo(lo,n)
137 0 : lp = llo(lop,n)
138 0 : lmini = abs(lp - l)
139 0 : lmaxi = lp + l
140 : c..the gaunt conditions
141 0 : if ((mod(l+lp+lpp,2).eq.1) .or. (lpp.LT.lmini) .or.
142 0 : + (lpp.gt.lmaxi)) then
143 0 : ulojulo(lo,lop,lh) = 0.
144 : else
145 0 : do i = 1,jri(n)
146 : x(i) = ( flo(n,i,1,lo,1)*flo(n,i,1,lop,2)+
147 0 : + flo(n,i,2,lo,1)*flo(n,i,2,lop,2) )*beff(i,lh,n)
148 : enddo
149 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),ulojulo(lo,lop,lh))
150 : endif
151 : enddo
152 : enddo
153 : endif ! local orbitals
154 : c**************************************************************************
155 : c...overlaps of the apws only
156 : c**************************************************************************
157 0 : do lp = 0,lwn
158 0 : do l = 0,lwn
159 0 : lmini = abs(lp-l)
160 0 : lmaxi = lp + l
161 : c..gaunt conditions
162 0 : if ((mod(l+lp+lpp,2).eq.1) .or. (lpp.LT.lmini) .or.
163 0 : + (lpp.gt.lmaxi)) then
164 0 : uju(l,lp,lh) = 0.
165 0 : dju(l,lp,lh) = 0.
166 0 : ujd(l,lp,lh) = 0.
167 0 : djd(l,lp,lh) = 0.
168 : else
169 0 : do i = 1,jri(n)
170 : x(i) = ( f(n,i,1,l,1)*f(n,i,1,lp,2)+
171 0 : + f(n,i,2,l,1)*f(n,i,2,lp,2) )*beff(i,lh,n)
172 : enddo
173 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),uju(l,lp,lh))
174 :
175 0 : do i = 1,jri(n)
176 : x(i) = ( g(n,i,1,l,1)*f(n,i,1,lp,2)+
177 0 : + g(n,i,2,l,1)*f(n,i,2,lp,2) )*beff(i,lh,n)
178 : enddo
179 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),dju(l,lp,lh))
180 :
181 0 : do i = 1,jri(n)
182 : x(i) = ( f(n,i,1,l,1)*g(n,i,1,lp,2)+
183 0 : + f(n,i,2,l,1)*g(n,i,2,lp,2) )*beff(i,lh,n)
184 : enddo
185 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujd(l,lp,lh))
186 :
187 0 : do i = 1,jri(n)
188 : x(i) = ( g(n,i,1,l,1)*g(n,i,1,lp,2)+
189 0 : + g(n,i,2,l,1)*g(n,i,2,lp,2) )*beff(i,lh,n)
190 : enddo
191 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),djd(l,lp,lh))
192 : endif
193 : enddo ! l
194 :
195 : c********************************************************************
196 : c...overlaps of the lo's with the apws
197 : c********************************************************************
198 0 : if (nlo(n).GE.1) then
199 0 : do lo = 1,nlo(n)
200 0 : l = llo(lo,n)
201 0 : lmini = abs(lp-l)
202 0 : lmaxi = lp + l
203 : c..gaunt conditions
204 0 : if ((mod(l+lp+lpp,2).eq.1) .OR. (lpp.lt.lmini) .or.
205 0 : + (lpp.gt.lmaxi)) then
206 0 : ujulo(lo,lp,lh) = 0.
207 0 : djulo(lo,lp,lh) = 0.
208 0 : loju(lo,lp,lh) = 0.
209 0 : lojd(lo,lp,lh) = 0.
210 : else
211 0 : do i = 1,jri(n)
212 : x(i) = ( flo(n,i,1,lo,1)*f(n,i,1,lp,2)+
213 0 : + flo(n,i,2,lo,1)*f(n,i,2,lp,2) )*beff(i,lh,n)
214 : enddo
215 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),loju(lo,lp,lh))
216 :
217 0 : do i = 1,jri(n)
218 : x(i) = ( flo(n,i,1,lo,1)*g(n,i,1,lp,2)+
219 0 : + flo(n,i,2,lo,1)*g(n,i,2,lp,2) )*beff(i,lh,n)
220 : enddo
221 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),lojd(lo,lp,lh))
222 :
223 0 : do i = 1,jri(n)
224 : x(i) = ( flo(n,i,1,lo,2)*f(n,i,1,lp,1)+
225 0 : + flo(n,i,2,lo,2)*f(n,i,2,lp,1) )*beff(i,lh,n)
226 : enddo
227 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujulo(lo,lp,lh))
228 :
229 0 : do i = 1,jri(n)
230 : x(i) = ( flo(n,i,1,lo,2)*g(n,i,1,lp,1)+
231 0 : + flo(n,i,2,lo,2)*g(n,i,2,lp,1) )*beff(i,lh,n)
232 : enddo
233 0 : call intgr3(x,rmsh(1:,n),dx(n),jri(n),djulo(lo,lp,lh))
234 :
235 : endif
236 : enddo !lo
237 : endif ! local orbitals
238 : enddo !lp
239 : enddo !lh
240 : c********************************************************************
241 : c multiply with gaunt-coefficient (apw-apw)
242 : c********************************************************************
243 0 : ujug(:,:,n)=cmplx(0.0,0.0)
244 0 : ujdg(:,:,n)=cmplx(0.0,0.0)
245 0 : djug(:,:,n)=cmplx(0.0,0.0)
246 0 : djdg(:,:,n)=cmplx(0.0,0.0)
247 0 : do l = 0,lwn
248 0 : do m = -l,l
249 0 : lm=l*(l+1)+m
250 0 : do lp = 0,lwn
251 0 : do mp = -lp,lp
252 0 : lpmp=lp*(lp+1)+mp
253 0 : do lh = 0,nlh(ns)
254 0 : lpp=llh(lh,ns)
255 0 : llpp = lpp*(lpp+1)
256 0 : mpp = mp - m
257 0 : lmpp = llpp + mpp
258 0 : lmini = abs(l-lpp)
259 0 : lmaxi = l+lpp
260 0 : jmem=0
261 0 : do j=1,nmem(lh,ns)
262 0 : if(mlh(j,lh,ns)==mpp)then
263 : jmem=j
264 : exit
265 : endif
266 : enddo
267 0 : l_doit= (lmini.le.lp)
268 0 : l_doit=l_doit.and.(lp.le.lmaxi)
269 0 : l_doit=l_doit.and.(mod(l+lp+lpp,2).eq.0)
270 0 : l_doit=l_doit.and.(abs(mpp).LE.lpp)
271 0 : l_doit=l_doit.and.(jmem.ne.0)
272 0 : if ( l_doit )then
273 : factor=ic**(l-lp)*
274 : * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)*
275 0 : * clnu(jmem,lh,ns)
276 : ujug(lpmp,lm,n)=ujug(lpmp,lm,n)+
277 0 : + factor*uju(lp,l,lh)
278 : ujdg(lpmp,lm,n)=ujdg(lpmp,lm,n)+
279 0 : + factor*ujd(lp,l,lh)
280 : djug(lpmp,lm,n)=djug(lpmp,lm,n)+
281 0 : + factor*dju(lp,l,lh)
282 : djdg(lpmp,lm,n)=djdg(lpmp,lm,n)+
283 0 : + factor*djd(lp,l,lh)
284 :
285 : endif
286 : enddo ! lh
287 : enddo ! mp
288 : enddo ! lp
289 : enddo ! m
290 : enddo ! l
291 : c******************************************************************
292 : c multiply with the gaunt-coefficient (apw-lo)
293 : c******************************************************************
294 0 : if(nlo(n).ge.1) then
295 0 : ulojug(:,:,:,n)=cmplx(0.0,0.0)
296 0 : ulojdg(:,:,:,n)=cmplx(0.0,0.0)
297 0 : ujulog(:,:,:,n)=cmplx(0.0,0.0)
298 0 : djulog(:,:,:,n)=cmplx(0.0,0.0)
299 :
300 0 : do lo = 1,nlo(n)
301 0 : l = llo(lo,n)
302 0 : do m = -l,l
303 0 : lm=l*(l+1)+m
304 :
305 0 : do lp=0,lwn
306 0 : do mp = -lp,lp
307 0 : lpmp=lp*(lp+1)+mp
308 :
309 0 : do lh = 0,nlh(ns)
310 0 : lpp=llh(lh,ns)
311 0 : llpp = lpp*(lpp+1)
312 0 : mpp = mp - m ! <mp|m+mpp>
313 0 : lmpp = llpp + mpp
314 0 : lmini = abs(l-lpp)
315 0 : lmaxi = l+lpp
316 0 : jmem=0
317 0 : do j=1,nmem(lh,ns)
318 0 : if(mlh(j,lh,ns)==mpp)then
319 : jmem=j
320 : exit
321 : endif
322 : enddo
323 0 : l_doit= (lmini.le.lp)
324 0 : l_doit=l_doit.and.(lp.le.lmaxi)
325 0 : l_doit=l_doit.and.(mod(l+lp+lpp,2).eq.0)
326 0 : l_doit=l_doit.and.(abs(mpp).LE.lpp)
327 0 : l_doit=l_doit.and.(jmem.ne.0)
328 0 : if ( l_doit )then
329 : factor=ic**(l-lp)*
330 : * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)*
331 0 : * clnu(jmem,lh,ns)
332 :
333 : ujulog(lpmp,lo,m,n)=
334 : = ujulog(lpmp,lo,m,n)+
335 0 : * ujulo(lo,lp,lpp)*factor
336 :
337 : djulog(lpmp,lo,m,n)=
338 : = djulog(lpmp,lo,m,n)+
339 0 : * djulo(lo,lp,lpp)*factor
340 : endif
341 :
342 0 : mpp = m - mp ! <m|mp+mpp>
343 0 : lmpp = llpp + mpp
344 0 : lmini = abs(lp-lpp)
345 0 : lmaxi = lp+lpp
346 0 : jmem=0
347 0 : do j=1,nmem(lh,ns)
348 0 : if(mlh(j,lh,ns)==mpp)then
349 : jmem=j
350 : exit
351 : endif
352 : enddo
353 0 : l_doit= (lmini.le.l)
354 0 : l_doit=l_doit.and.(l.le.lmaxi)
355 0 : l_doit=l_doit.and.(mod(l+lp+lpp,2).eq.0)
356 0 : l_doit=l_doit.and.(abs(mpp).LE.lpp)
357 0 : l_doit=l_doit.and.(jmem.ne.0)
358 0 : if ( l_doit )then
359 : factor=ic**(lp-l)*
360 : * gaunt1(l,lpp,lp,m,mpp,mp,lmaxd)*
361 0 : * clnu(jmem,lh,ns)
362 :
363 : ulojug(lpmp,lo,m,n)=
364 : = ulojug(lpmp,lo,m,n)+
365 0 : * loju(lo,lp,lpp)*factor
366 :
367 : ulojdg(lpmp,lo,m,n)=
368 : = ulojdg(lpmp,lo,m,n)+
369 0 : * lojd(lo,lp,lpp)*factor
370 : endif
371 :
372 : enddo ! lh
373 : enddo ! mp
374 : enddo ! lop
375 : enddo ! m
376 : enddo ! lo
377 : endif ! local orbitals on this atom
378 :
379 :
380 : c*************************************************************
381 : c multiply with the gaunt-coefficient (lo-lo)
382 : c*************************************************************
383 0 : if(nlo(n).ge.1) then
384 :
385 0 : ulojulog(:,:,:,:,n)=cmplx(0.0,0.0)
386 0 : do lo = 1,nlo(n)
387 0 : l = llo(lo,n)
388 0 : do m = -l,l
389 0 : lm=l*(l+1)+m
390 :
391 0 : do lop = 1,nlo(n)
392 0 : lp = llo(lop,n)
393 0 : do mp = -lp,lp
394 0 : lpmp=lp*(lp+1)+mp
395 :
396 0 : do lh = 0,nlh(ns)
397 0 : lpp=llh(lh,ns)
398 0 : llpp = lpp*(lpp+1)
399 0 : mpp = mp - m ! <mp|m+mpp>
400 0 : lmpp = llpp + mpp
401 0 : lmini = abs(l-lpp)
402 0 : lmaxi = l+lpp
403 0 : jmem=0
404 0 : do j=1,nmem(lh,ns)
405 0 : if(mlh(j,lh,ns)==mpp)then
406 : jmem=j
407 : exit
408 : endif
409 : enddo
410 0 : l_doit= (lmini.le.lp)
411 0 : l_doit=l_doit.and.(lp.le.lmaxi)
412 0 : l_doit=l_doit.and.(mod(l+lp+lpp,2).eq.0)
413 0 : l_doit=l_doit.and.(abs(mpp).LE.lpp)
414 0 : l_doit=l_doit.and.(jmem.ne.0)
415 0 : if ( l_doit )then
416 : factor=ic**(l-lp)*
417 : * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)*
418 0 : * clnu(jmem,lh,ns)
419 :
420 : ulojulog(lop,mp,lo,m,n)=
421 : = ulojulog(lop,mp,lo,m,n)+
422 0 : * ulojulo(lop,lo,lpp)*factor
423 :
424 : endif
425 : enddo ! lh
426 : enddo ! mp
427 : enddo ! lop
428 : enddo ! m
429 : enddo ! lo
430 : endif ! local orbitals on this atom
431 :
432 :
433 0 : na=na+neq(n)
434 : enddo !ntype
435 0 : deallocate(djd,ujd,uju,dju)
436 :
437 0 : deallocate(ujulo,djulo,ulojulo,loju,lojd)
438 :
439 0 : call timestop("wann_gabeffgagaunt")
440 0 : end subroutine wann_gabeffgagaunt
441 :
442 : end module m_wann_gabeffgagaunt
|
