Line data Source code
1 : c****************************************c
2 : c Muffin tin contribution to uHu c
3 : c < u_{k+b1} | H_{k}^{mt} | u_{k+b2} > c
4 : c****************************************c
5 : c Routine to set up T(lm,lmp) for all c
6 : c pairs (b1,b2) and every atom n c
7 : c c
8 : c Includes spherical and non-sph. c
9 : c contributions to Hamiltonian: c
10 : c H^{mt} ~ H^{sph} + H^{non-sph} c
11 : c c
12 : c possible SOC contribution is dealt c
13 : c with in a separate routine c
14 : c****************************************c
15 : c J.-P. Hanke, Dec. 2015 c
16 : c****************************************c
17 : MODULE m_wann_uHu_tlmplm
18 : CONTAINS
19 0 : SUBROUTINE wann_uHu_tlmplm(
20 : > memd,nlhd,ntypsd,ntypd,jmtd,lmaxd,jspd,
21 0 : > ntype,dx,rmsh,jri,lmax,ntypsy,natd,
22 0 : > lnonsph,lmd,lmplmd,clnu,mlh,nmem,llh,nlh,neq,
23 : > irank,mlotot,mlolotot,
24 0 : > vr,nlod,llod,loplod,ello,llo,nlo,lo1l,l_dulo,
25 0 : > ulo_der,f,g,flo,f_b,g_b,flo_b,
26 0 : > kdiff,kdiff2,nntot,nntot2,bmat,bbmat,vr0,epar,
27 0 : > invsat,
28 : > l_skip_sph,l_skip_non,l_skip_loc,
29 0 : < tuu,tud,tdu,tdd,tuulo,tulou,tdulo,tulod,tuloulo)
30 :
31 : USE m_intgr, ONLY : intgr3
32 : USE m_radflo
33 : USE m_radfun
34 : USE m_tlo
35 : USE m_sphbes
36 : USE m_ylm
37 : USE m_gaunt, ONLY: gaunt1
38 : USE m_dujdr
39 : USE m_wann_uHu_radintsra5
40 : USE m_wann_uHu_radintsra3
41 : USE m_wann_uHu_tlo
42 : USE m_constants
43 : #ifdef CPP_MPI
44 : USE mpi
45 : #endif
46 :
47 : IMPLICIT NONE
48 : C ..
49 : C .. Scalar Arguments ..
50 : INTEGER, INTENT (IN) :: memd,nlhd,ntypsd,ntypd,jmtd,lmaxd,jspd
51 : INTEGER, INTENT (IN) :: lmd,lmplmd,ntype,irank
52 : INTEGER, INTENT (IN) :: nlod,llod,loplod,natd,mlotot,mlolotot
53 : INTEGER, INTENT (IN) :: nntot,nntot2
54 : LOGICAL, INTENT (IN) :: l_skip_sph,l_skip_non,l_skip_loc
55 : C ..
56 : C .. Array Arguments ..
57 : COMPLEX, INTENT (OUT)::
58 : > tdd(0:lmd,0:lmd,ntypd,nntot*nntot2),
59 : > tdu(0:lmd,0:lmd,ntypd,nntot*nntot2),
60 : > tud(0:lmd,0:lmd,ntypd,nntot*nntot2),
61 : > tuu(0:lmd,0:lmd,ntypd,nntot*nntot2)
62 : COMPLEX, INTENT (OUT)::
63 : > tdulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
64 : > tuulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
65 : > tulou(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
66 : > tulod(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
67 : > tuloulo(nlod,-llod:llod,nlod,-llod:llod,ntypd,nntot*nntot2)
68 : COMPLEX, INTENT (IN) :: clnu(memd,0:nlhd,ntypsd)
69 : INTEGER, INTENT (IN) :: llo(nlod,ntypd),nlo(ntypd)
70 : INTEGER, INTENT (IN) :: lo1l(0:llod,ntypd),neq(ntypd)
71 : INTEGER, INTENT (IN) :: mlh(memd,0:nlhd,ntypsd),nlh(ntypsd)
72 : INTEGER, INTENT (IN) :: nmem(0:nlhd,ntypsd),llh(0:nlhd,ntypsd)
73 : INTEGER, INTENT (IN) :: jri(ntypd),lmax(ntypd),ntypsy(natd)
74 : INTEGER, INTENT (IN) :: lnonsph(ntypd),ulo_der(nlod,ntypd)
75 : INTEGER, INTENT (IN) :: invsat(natd)
76 : REAL, INTENT (IN) :: dx(ntypd),rmsh(jmtd,ntypd)
77 : REAL, INTENT (IN) :: vr(jmtd,0:nlhd,ntypd)
78 : REAL, INTENT (IN) :: ello(nlod,ntypd)
79 : REAL, INTENT (IN) :: epar(0:lmaxd,ntypd)
80 : REAL, INTENT (IN) :: vr0(jmtd,ntypd)
81 : REAL, INTENT (IN) :: f (ntypd,jmtd,2,0:lmaxd),
82 : > g (ntypd,jmtd,2,0:lmaxd),
83 : > f_b(ntypd,jmtd,2,0:lmaxd),
84 : > g_b(ntypd,jmtd,2,0:lmaxd),
85 : > flo (ntypd,jmtd,2,nlod),
86 : > flo_b(ntypd,jmtd,2,nlod)
87 : REAL, INTENT (IN) :: bmat(3,3),bbmat(3,3)
88 : REAL, INTENT (IN) :: kdiff (3,nntot)
89 : REAL, INTENT (IN) :: kdiff2(3,nntot2)
90 : LOGICAL, INTENT (IN) :: l_dulo(nlod,ntypd)
91 : C ..
92 : C .. Local Scalars ..
93 : COMPLEX cil,cci,cil1,fac1,fac2,fac3,cilp,cciljj1,ciljj2
94 : INTEGER i,l,l2,lamda,lh,lm,lmin,lmp,lp
95 : INTEGER lp1,lpl,m,mem,mems,mp,mu,n,nh,noded,nodeu,nrec,nsym,na,ne
96 : INTEGER ljj1,ljj2,mjj1,mjj2,lo,lop
97 : INTEGER ltil1,ltil2,mtil1,mtil2
98 : INTEGER mlo, mlolo, iu
99 : INTEGER :: lmini,lmini2,lmini3
100 : INTEGER :: lmaxi,lmaxi2,lmaxi3
101 : INTEGER :: lmx,lmx2,lmx3
102 : INTEGER :: ll,llp,lljj1,lljj2,lmjj1,lmjj2
103 : INTEGER :: ikpt_b,ikpt_b2
104 : INTEGER :: lwn,indexx,indexx2
105 : REAL gs,rk1,rk2,sign
106 : REAL gc1,gc2,gc3,invsfct,e,temp_val1,temp_val2
107 : REAL t1,t0,t_sphbes,t_ylm,t_radint,t_gggint,t_sphint,t_lo
108 : REAL tt1,tt2,tt3,sy1,sy2
109 : LOGICAL :: l_nns,l_oldsym
110 : C ..
111 : C .. Local Arrays ..
112 0 : COMPLEX, ALLOCATABLE :: yl1(:),yl2(:)
113 0 : REAL, ALLOCATABLE :: dvd_non(:,:,:,:,:),dvu_non(:,:,:,:,:)
114 0 : REAL, ALLOCATABLE :: uvd_non(:,:,:,:,:),uvu_non(:,:,:,:,:)
115 0 : REAL, ALLOCATABLE :: dvd_sph(:,:,:,:,:),dvu_sph(:,:,:,:,:)
116 0 : REAL, ALLOCATABLE :: uvd_sph(:,:,:,:,:),uvu_sph(:,:,:,:,:)
117 0 : REAL, ALLOCATABLE :: p(:,:,:,:),p_b(:,:,:,:),dp_b(:,:,:,:)
118 0 : REAL, ALLOCATABLE :: q(:,:,:,:),q_b(:,:,:,:),dq_b(:,:,:,:)
119 0 : REAL, ALLOCATABLE :: dp(:,:,:,:),dq(:,:,:,:)
120 0 : REAL, ALLOCATABLE :: x(:)
121 0 : REAL, ALLOCATABLE :: jlpp(:),jj1(:,:),jj2(:,:)
122 : REAL :: bpt(3),bpt2(3),bkrot(3)
123 :
124 : #if (defined(CPP_MPI) && !defined(CPP_T90))
125 : INTEGER ierr(3)
126 : #endif
127 : C .. Intrinsic Functions ..
128 : INTRINSIC abs,cmplx,max,mod
129 : C ..
130 :
131 0 : l_oldsym=.false.
132 :
133 : sy1=1.0
134 : sy2=0.0
135 : if(l_oldsym) then
136 : sy1=0.5
137 : sy2=0.5
138 : endif
139 :
140 0 : t_sphbes = 0.
141 0 : t_ylm = 0.
142 0 : t_radint = 0.
143 0 : t_sphint = 0.
144 0 : t_gggint = 0.
145 0 : t_lo = 0.
146 0 : tt1 = 0.
147 0 : tt2 = 0.
148 0 : tt3 = 0.
149 :
150 0 : cci = cmplx(0.,-1.)
151 0 : tdulo = cmplx(0.,0.); tuulo = cmplx(0.,0.); tuloulo = cmplx(0.,0.)
152 0 : tuu = cmplx(0.,0.); tdu = cmplx(0.,0.)
153 0 : tud = cmplx(0.,0.); tdd = cmplx(0.,0.)
154 :
155 0 : allocate( yl1((lmaxd+1)**2), yl2((lmaxd+1)**2) )
156 0 : allocate( dvd_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
157 0 : allocate( dvu_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
158 0 : allocate( uvd_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
159 0 : allocate( uvu_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
160 0 : allocate( dvd_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
161 0 : allocate( dvu_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
162 0 : allocate( uvd_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
163 0 : allocate( uvu_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
164 0 : allocate( p(jmtd,2,0:lmaxd,0:lmaxd), p_b(jmtd,2,0:lmaxd,0:lmaxd) )
165 0 : allocate( q(jmtd,2,0:lmaxd,0:lmaxd), q_b(jmtd,2,0:lmaxd,0:lmaxd) )
166 0 : allocate( dp(jmtd,2,0:lmaxd,0:lmaxd) )
167 0 : allocate( dp_b(jmtd,2,0:lmaxd,0:lmaxd) )
168 0 : allocate( dq(jmtd,2,0:lmaxd,0:lmaxd) )
169 0 : allocate( dq_b(jmtd,2,0:lmaxd,0:lmaxd) )
170 0 : allocate( x(jmtd) )
171 0 : allocate( jlpp(0:lmaxd), jj1(0:lmaxd,jmtd), jj2(0:lmaxd,jmtd) )
172 :
173 : c***************************************c
174 : c begin 1st neighbor loop: <k+b1| c
175 : c***************************************c
176 0 : do ikpt_b=1,nntot
177 0 : bpt = kdiff(:,ikpt_b)
178 0 : rk1 = sqrt(dot_product(bpt,matmul(bbmat,bpt)))
179 : ! rk1= sqrt(dotirp(bpt,bpt,bbmat))
180 :
181 : c***************************************c
182 : c begin 2nd neighbor loop: |k+b2> c
183 : c***************************************c
184 0 : do ikpt_b2=1,nntot2
185 0 : indexx = ikpt_b2 + (ikpt_b-1)*nntot2
186 0 : bpt2 = kdiff2(:,ikpt_b2)
187 0 : rk2 = sqrt(dot_product(bpt2,matmul(bbmat,bpt2)))
188 : ! rk2= sqrt(dotirp(bpt2,bpt2,bbmat))
189 :
190 0 : na = 1
191 0 : mlo = 1 ; mlolo = 1
192 :
193 : ! loop over atoms
194 0 : DO 210 n = 1,ntype
195 0 : lwn = lmax(n)
196 0 : nsym = ntypsy(na)
197 0 : nh = nlh(nsym)
198 :
199 : l_nns = l_skip_non .OR. (lnonsph(n).LT.0) .OR.
200 0 : > ( (invsat(na).ne.0).and.(invsat(na).ne.1) )
201 0 : if(invsat(na).eq.0) invsfct = 1.0
202 0 : if(invsat(na).eq.1) invsfct = 2.0
203 0 : if(indexx.eq.1 .and. irank.eq.0) write(*,*)'invsfct=',invsfct
204 :
205 : ! set up spherical Bessel functions
206 : ! jj1(l,b1*r) and jj2(l,b2*r)
207 0 : call cpu_time(t0)
208 0 : do i=1,jri(n)
209 0 : gs = rk1*rmsh(i,n)
210 0 : call sphbes(lwn,gs,jlpp)
211 0 : jj1(:,i) = jlpp(:)
212 :
213 0 : gs = rk2*rmsh(i,n)
214 0 : call sphbes(lwn,gs,jlpp)
215 0 : jj2(:,i) = jlpp(:)
216 : enddo
217 0 : call cpu_time(t1)
218 0 : t_sphbes=t_sphbes + t1-t0
219 :
220 : ! set up spherical harmonics
221 : ! yl1(lm) for b1 and yl2(lm) for b2
222 0 : yl1 = cmplx(0.,0.)
223 0 : bkrot=MATMUL(bpt,bmat)
224 0 : call ylm4(lwn,bkrot,yl1)
225 :
226 0 : yl2 = cmplx(0.,0.)
227 0 : bkrot=MATMUL(bpt2,bmat)
228 0 : call ylm4(lwn,bkrot,yl2)
229 :
230 0 : call cpu_time(t0)
231 0 : t_ylm=t_ylm + t0-t1
232 :
233 : ! set up products of radial functions and sph. Bessel
234 0 : DO ljj1 = 0, lwn
235 0 : DO l = 0, lwn
236 0 : DO i=1,jri(n)
237 0 : p(i,:,l,ljj1) = f(n,i,:,l)*jj1(ljj1,i)
238 0 : q(i,:,l,ljj1) = g(n,i,:,l)*jj1(ljj1,i)
239 0 : p_b(i,:,l,ljj1) = f_b(n,i,:,l)*jj2(ljj1,i)
240 0 : q_b(i,:,l,ljj1) = g_b(n,i,:,l)*jj2(ljj1,i)
241 : ENDDO
242 : CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),f(n,:,:,l),
243 0 : > jj1,rk1,ljj1,lmaxd,dp(:,:,l,ljj1))
244 : CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),g(n,:,:,l),
245 0 : > jj1,rk1,ljj1,lmaxd,dq(:,:,l,ljj1))
246 : CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),f_b(n,:,:,l),
247 0 : > jj2,rk2,ljj1,lmaxd,dp_b(:,:,l,ljj1))
248 : CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),g_b(n,:,:,l),
249 0 : > jj2,rk2,ljj1,lmaxd,dq_b(:,:,l,ljj1))
250 : ENDDO
251 : ENDDO
252 :
253 : c****************************************************c
254 : c compute radial integrals c
255 : c <u(l')jj1(ljj1) | v(lamda,nu) | u(l)jj2(ljj2)> c
256 : c where v(lamda,nu) ~ V^{sph} + V^{non-sph} c
257 : c****************************************************c
258 : ! spherical part
259 0 : IF(.not.l_skip_sph) THEN
260 0 : uvu_sph = 0.
261 0 : uvd_sph = 0.
262 0 : dvu_sph = 0.
263 0 : dvd_sph = 0.
264 0 : DO ljj1 = 0, lwn
265 0 : DO lp = 0, lwn
266 0 : lmini = abs(lp-ljj1)
267 0 : lmaxi = lp+ljj1
268 0 : lmx = min(lmaxi,lwn)
269 0 : DO lh=lmini,lmx
270 0 : if(mod(lmaxi+lh,2).eq.1) cycle
271 0 : DO ljj2 = 0, lwn
272 0 : lmini2= abs(lh-ljj2)
273 0 : lmaxi2= lh+ljj2
274 0 : lmx2 = min(lmaxi2,lwn)
275 0 : DO l = lmini2, lmx2
276 0 : if(mod(lmaxi2+l,2).eq.1) cycle
277 0 : e = (epar(l,n)+epar(lp,n))/2.0!epar(lh,n)
278 0 : if(.not.l_oldsym) then
279 : call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
280 : > e,vr0(1,n),p(:,:,lp,ljj1),p_b(:,:,l,ljj2),
281 : > dp(:,:,lp,ljj1),dp_b(:,:,l,ljj2),lmaxd,lh,
282 0 : > uvu_sph(lh,l,lp,ljj2,ljj1),irank)
283 :
284 : call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
285 : > e,vr0(1,n),p(:,:,lp,ljj1),q_b(:,:,l,ljj2),
286 : > dp(:,:,lp,ljj1),dq_b(:,:,l,ljj2),lmaxd,lh,
287 0 : > uvd_sph(lh,l,lp,ljj2,ljj1),irank)
288 :
289 : call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
290 : > e,vr0(1,n),q(:,:,lp,ljj1),p_b(:,:,l,ljj2),
291 : > dq(:,:,lp,ljj1),dp_b(:,:,l,ljj2),lmaxd,lh,
292 0 : > dvu_sph(lh,l,lp,ljj2,ljj1),irank)
293 :
294 : call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
295 : > e,vr0(1,n),q(:,:,lp,ljj1),q_b(:,:,l,ljj2),
296 : > dq(:,:,lp,ljj1),dq_b(:,:,l,ljj2),lmaxd,lh,
297 0 : > dvd_sph(lh,l,lp,ljj2,ljj1),irank)
298 :
299 : else
300 :
301 : e = epar(l,n)!( epar(l,n) + epar(lp,n) )/2.0
302 : call wann_uHu_radintsra3(jmtd,jri(n),rmsh(1,n),dx(n),
303 : > e,vr0(1,n),p(:,:,lp,ljj1),p_b(:,:,l,ljj2),
304 : > dp_b(:,:,l,ljj2),lmaxd,lh,uvu_sph(lh,l,lp,ljj2,ljj1))
305 :
306 : call wann_uHu_radintsra3(jmtd,jri(n),rmsh(1,n),dx(n),
307 : > e,vr0(1,n),p(:,:,lp,ljj1),q_b(:,:,l,ljj2),
308 : > dq_b(:,:,l,ljj2),lmaxd,lh,uvd_sph(lh,l,lp,ljj2,ljj1))
309 :
310 : call wann_uHu_radintsra3(jmtd,jri(n),rmsh(1,n),dx(n),
311 : > e,vr0(1,n),q(:,:,lp,ljj1),p_b(:,:,l,ljj2),
312 : > dp_b(:,:,l,ljj2),lmaxd,lh,dvu_sph(lh,l,lp,ljj2,ljj1))
313 :
314 : call wann_uHu_radintsra3(jmtd,jri(n),rmsh(1,n),dx(n),
315 : > e,vr0(1,n),q(:,:,lp,ljj1),q_b(:,:,l,ljj2),
316 : > dq_b(:,:,l,ljj2),lmaxd,lh,dvd_sph(lh,l,lp,ljj2,ljj1))
317 : endif
318 : ENDDO
319 : ENDDO
320 : ENDDO
321 : ENDDO
322 : ENDDO
323 :
324 : if(.false. .and. (irank.eq.0)) then
325 : DO ljj1=0,lwn
326 : DO ljj2=0,lwn
327 : DO lp=0,lwn
328 : DO l=0,lwn
329 : DO lh=0,lwn
330 : if(abs(uvu_sph(lh,l,lp,ljj2,ljj1)-
331 : > uvu_sph(lh,lp,l,ljj1,ljj2)).gt.1e-10) then
332 : write(*,*)'uvu',ikpt_b,ikpt_b2,n
333 : write(*,*)'bpt',bpt
334 : write(*,*)'bpt2',bpt2
335 : write(*,*)ljj1,ljj2,lp,l,lh
336 : write(*,*)uvu_sph(lh,l,lp,ljj2,ljj1)
337 : write(*,*)uvu_sph(lh,lp,l,ljj1,ljj2)
338 : endif
339 : if(abs(dvd_sph(lh,l,lp,ljj2,ljj1)-
340 : > dvd_sph(lh,lp,l,ljj1,ljj2)).gt.1e-10) then
341 : write(*,*)'dvd',ikpt_b,ikpt_b2,n
342 : write(*,*)'bpt',bpt
343 : write(*,*)'bpt2',bpt2
344 : write(*,*)ljj1,ljj2,lp,l,lh
345 : write(*,*)dvd_sph(lh,l,lp,ljj2,ljj1)
346 : write(*,*)dvd_sph(lh,lp,l,ljj1,ljj2)
347 : endif
348 : if(abs(dvu_sph(lh,l,lp,ljj2,ljj1)-
349 : > uvd_sph(lh,lp,l,ljj1,ljj2)).gt.1e-10) then
350 : write(*,*)'dvu',ikpt_b,ikpt_b2,n
351 : write(*,*)'bpt',bpt
352 : write(*,*)'bpt2',bpt2
353 : write(*,*)ljj1,ljj2,lp,l,lh
354 : write(*,*)dvu_sph(lh,l,lp,ljj2,ljj1)
355 : write(*,*)uvd_sph(lh,lp,l,ljj1,ljj2)
356 : endif
357 : if(abs(uvd_sph(lh,l,lp,ljj2,ljj1)-
358 : > dvu_sph(lh,lp,l,ljj1,ljj2)).gt.1e-10) then
359 : write(*,*)'uvd',ikpt_b,ikpt_b2,n
360 : write(*,*)'bpt',bpt
361 : write(*,*)'bpt2',bpt2
362 : write(*,*)ljj1,ljj2,lp,l,lh
363 : write(*,*)uvd_sph(lh,l,lp,ljj2,ljj1)
364 : write(*,*)dvu_sph(lh,lp,l,ljj1,ljj2)
365 : goto 777
366 : endif
367 : ENDDO
368 : ENDDO
369 : ENDDO
370 : ENDDO
371 : ENDDO
372 :
373 : 777 continue
374 : endif
375 :
376 : ENDIF
377 :
378 :
379 : ! nonspherical part
380 0 : IF(.not.l_nns) then
381 0 : DO ljj1 = 0, lwn
382 0 : DO ljj2 = 0, lwn
383 0 : DO lp = 0, lwn
384 0 : DO l = 0, lp
385 : c DO l = 0, lwn
386 0 : DO lh = 1, nh
387 0 : DO i = 1,jri(n)
388 : x(i) = ( p(i,1,lp,ljj1)*p_b(i,1,l,ljj2)
389 : > + p(i,2,lp,ljj1)*p_b(i,2,l,ljj2) )
390 0 : > * vr(i,lh,n)
391 : ENDDO
392 : CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
393 0 : > uvu_non(lh,l,lp,ljj2,ljj1))
394 0 : DO i = 1,jri(n)
395 : x(i) = ( q(i,1,lp,ljj1)*p_b(i,1,l,ljj2)
396 : > + q(i,2,lp,ljj1)*p_b(i,2,l,ljj2) )
397 0 : > * vr(i,lh,n)
398 : ENDDO
399 : CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
400 0 : > dvu_non(lh,l,lp,ljj2,ljj1))
401 0 : DO i = 1,jri(n)
402 : x(i) = ( p(i,1,lp,ljj1)*q_b(i,1,l,ljj2)
403 : > + p(i,2,lp,ljj1)*q_b(i,2,l,ljj2) )
404 0 : > * vr(i,lh,n)
405 : ENDDO
406 : CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
407 0 : > uvd_non(lh,l,lp,ljj2,ljj1))
408 0 : DO i = 1,jri(n)
409 : x(i) = ( q(i,1,lp,ljj1)*q_b(i,1,l,ljj2)
410 : > + q(i,2,lp,ljj1)*q_b(i,2,l,ljj2) )
411 0 : > * vr(i,lh,n)
412 : ENDDO
413 : CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
414 0 : > dvd_non(lh,l,lp,ljj2,ljj1))
415 : ENDDO
416 : ENDDO
417 : ENDDO
418 : ENDDO
419 : ENDDO
420 :
421 : ! symmetry can be used if f=f_b and g=g_b
422 : ! i.e. if jspin = jspin_b
423 0 : DO lp=0,lwn
424 0 : DO l=lp+1,lwn
425 0 : uvu_non(:,l,lp,:,:) = uvu_non(:,lp,l,:,:)
426 0 : dvd_non(:,l,lp,:,:) = dvd_non(:,lp,l,:,:)
427 0 : uvd_non(:,l,lp,:,:) = dvu_non(:,lp,l,:,:)
428 0 : dvu_non(:,l,lp,:,:) = uvd_non(:,lp,l,:,:)
429 : ENDDO
430 : ENDDO
431 : ENDIF
432 :
433 0 : call cpu_time(t0)
434 0 : t_radint = t_radint + t0-t1
435 :
436 0 : tuu(:,:,n,indexx) = cmplx(0.,0.)
437 0 : tdu(:,:,n,indexx) = cmplx(0.,0.)
438 0 : tud(:,:,n,indexx) = cmplx(0.,0.)
439 0 : tdd(:,:,n,indexx) = cmplx(0.,0.)
440 :
441 0 : if(l_skip_sph) GOTO 444
442 : c************** SPHERICAL CONTRIBUTION *******************c
443 : c compute product of the two Gaunt coefficients c
444 : c with the radial integrals (+prefactors) c
445 : c*********************************************************c
446 : c We deal with two Gaunt coefficients: c
447 : c G1 = G( (ltil1,mtil1), (ljj1,mjj1) , ( lp, mp) ) c
448 : c G2 = G( (l , m) , (ljj2, mjj2), (ltil1, mtil1) ) c
449 : c*********************************************************c
450 : c use Gaunt conditions to reduce number of operations. c
451 : c coefficient G(L1,L2,L3) only nonzero if c
452 : c a) l1 + l2 + l3 = even c
453 : c b) |l2-l3| <= l1 <= l2+l3 c
454 : c c) m1 = m2 + m3 c
455 : c*********************************************************c
456 0 : DO ljj1=0,lwn
457 0 : lljj1 = ljj1 * (ljj1 + 1)
458 0 : cciljj1 = cci**ljj1
459 0 : DO mjj1=-ljj1,ljj1
460 0 : lmjj1 = lljj1 + mjj1
461 0 : fac1 = cciljj1*conjg(yl1(lmjj1 + 1))
462 0 : if(fac1.eq.0.0) cycle
463 :
464 0 : DO lp=0,lwn
465 0 : llp = lp*(lp+1)
466 0 : lmini = abs(ljj1-lp)
467 0 : lmaxi = ljj1+lp
468 0 : lmx = min(lmaxi,lwn)
469 0 : DO mp=-lp,lp
470 0 : lmp = llp+mp
471 0 : mtil1=mjj1+mp
472 0 : IF(abs(mtil1).gt.lmx) cycle
473 :
474 0 : DO ltil1=lmini,lmx
475 : ! Gaunt conditions (G1):
476 0 : if((mod(lmaxi+ltil1,2).eq.1).or.
477 : > (abs(mtil1).gt.ltil1)) cycle
478 :
479 : gc1 = gaunt1(ltil1, ljj1, lp,
480 0 : > mtil1, mjj1, mp, lmaxd)
481 :
482 0 : cil1 = (ImagUnit**lp) * fac1 * gc1
483 :
484 0 : DO ljj2=0,lwn
485 0 : lljj2 = ljj2 * (ljj2 + 1)
486 0 : lmini2 = abs(ltil1-ljj2)
487 0 : lmaxi2 = ltil1+ljj2
488 0 : lmx2 = min(lmaxi2,lwn)
489 0 : DO mjj2=-ljj2,ljj2
490 0 : lmjj2 = lljj2 + mjj2
491 0 : m=mjj2+mtil1
492 0 : IF(abs(m).gt.lmx2) cycle
493 :
494 0 : DO l=lmini2,lmx2
495 : ! Gaunt conditions (G2):
496 0 : if((mod(lmaxi2+l,2).eq.1).or.
497 : > (abs(m).gt.l)) cycle
498 0 : ll = l*(l+1)
499 0 : lm = ll+m
500 :
501 : gc2 = gaunt1( l, ljj2, ltil1,
502 0 : > m, mjj2, mtil1, lmaxd)
503 :
504 : ! set up prefactor
505 : cil = cil1* ( ImagUnit ** (ljj2 - l) )
506 0 : + * gc2 * conjg( yl2(lmjj2 + 1) )
507 :
508 0 : if(cil.eq.0.0) cycle
509 :
510 : ! additional factor from symmetrization (below)
511 : c cil = 0.5 * cil
512 :
513 : ! compute T-coefficients
514 : ! using symmetrized uvu,uvd,dvu,dvd integrals
515 : tuu(lm, lmp, n, indexx)
516 : > = tuu(lm, lmp, n, indexx)
517 : > + cil * ( sy1*uvu_sph(ltil1, l, lp, ljj2, ljj1)
518 0 : > + sy2*uvu_sph(ltil1,lp, l, ljj1, ljj2) )
519 :
520 : tdd(lm, lmp, n, indexx)
521 : > = tdd(lm, lmp, n, indexx)
522 : > + cil * ( sy1*dvd_sph(ltil1, l, lp, ljj2, ljj1)
523 0 : > + sy2*dvd_sph(ltil1,lp, l, ljj1, ljj2) )
524 :
525 : tud(lm, lmp, n, indexx)
526 : > = tud(lm, lmp, n, indexx)
527 : > + cil * ( sy1*uvd_sph(ltil1, l, lp, ljj2, ljj1)
528 0 : > + sy2*dvu_sph(ltil1,lp, l, ljj1, ljj2) )
529 :
530 : tdu(lm, lmp, n, indexx)
531 : > = tdu(lm, lmp, n, indexx)
532 : > + cil * ( sy1*dvu_sph(ltil1, l, lp, ljj2, ljj1)
533 0 : > + sy2*uvd_sph(ltil1,lp, l, ljj1, ljj2) )
534 : ENDDO!ljj2
535 : ENDDO!m
536 : ENDDO!l
537 : ENDDO!ljj1
538 : ENDDO!mp
539 : ENDDO!lp
540 : ENDDO!mtil1
541 : ENDDO!ltil1
542 :
543 : 444 CONTINUE
544 0 : call cpu_time(t1)
545 0 : t_sphint = t_sphint + t1-t0
546 :
547 0 : IF( l_nns ) GOTO 555
548 : c************** NON-SPHERICAL CONTRIBUTION ***************c
549 : c compute product of the three Gaunt coefficients c
550 : c with the radial integrals (+prefactors) c
551 : c*********************************************************c
552 : c We deal with three Gaunt coefficients: c
553 : c G1 = G( (ltil1,mtil1) , (lamda,mu) , (ltil2,mtil2) ) c
554 : c G2 = G( (ltil1,mtil1) , (ljj1,mjj1) , (lp,mp) ) c
555 : c G3 = G( (l,m) , (ljj2,mjj2) , (ltil2,mtil2) ) c
556 : c*********************************************************c
557 : c use Gaunt conditions to reduce number of operations. c
558 : c coefficient G(L1,L2,L3) only nonzero if c
559 : c a) l1 + l2 + l3 = even c
560 : c b) |l2-l3| <= l1 <= l2+l3 c
561 : c c) m1 = m2 + m3 c
562 : c*********************************************************c
563 0 : DO ltil2 = 0, lwn
564 0 : DO mtil2 = -ltil2, ltil2
565 0 : DO lh = 1, nh
566 0 : lamda = llh(lh,nsym)
567 0 : lmini = abs(ltil2 - lamda)
568 0 : lmaxi = ltil2 + lamda
569 0 : lmx = min(lmaxi,lwn)
570 0 : mems = nmem(lh,nsym)
571 0 : DO mem = 1, mems
572 0 : mu = mlh(mem,lh,nsym)
573 0 : mtil1 = mtil2+mu
574 0 : IF(abs(mtil1).GT.lmx) CYCLE
575 :
576 0 : DO ltil1 = lmini, lmx
577 :
578 : ! Gaunt conditions (G1):
579 0 : if((mod(lmaxi+ltil1,2).eq.1).or.
580 : > (abs(mtil1).gt.ltil1)) cycle
581 :
582 : gc1 = gaunt1(ltil1, lamda, ltil2,
583 0 : > mtil1, mu, mtil2, lmaxd)
584 0 : fac1 = conjg(clnu(mem,lh,nsym)) * gc1 * invsfct
585 0 : if(fac1.eq.0.0) CYCLE
586 :
587 0 : DO lp = 0, lnonsph(n)
588 0 : llp = lp * (lp + 1)
589 0 : cilp = ImagUnit**lp
590 0 : DO mp = -lp, lp
591 0 : lmp = llp + mp
592 0 : mjj1 = mtil1 - mp
593 0 : IF(abs(mjj1).GT.lwn) CYCLE
594 :
595 0 : DO ljj1 = 0, lwn
596 0 : lmini2 = abs(lp-ljj1)
597 0 : lmaxi2 = lp+ljj1
598 :
599 : ! Gaunt conditions (G2):
600 : if((lmini2.gt.ltil1).or.(lmaxi2.lt.ltil1).or.
601 0 : > (mod(lmaxi2+ltil1,2).eq.1).or.
602 : > (abs(mjj1).gt.ljj1)) cycle
603 :
604 0 : lljj1 = ljj1 * (ljj1 + 1)
605 0 : lmjj1 = lljj1 + mjj1
606 :
607 : gc2 = gaunt1(ltil1, ljj1, lp,
608 0 : > mtil1, mjj1, mp, lmaxd)
609 : fac2 = fac1 * cilp * (cci**ljj1) * gc2
610 0 : > * conjg( yl1(lmjj1 + 1) )
611 0 : IF(fac2.eq.0.0) CYCLE
612 :
613 0 : DO ljj2 = 0, lwn
614 0 : lljj2 = ljj2 * (ljj2 + 1)
615 0 : ciljj2 = ImagUnit**ljj2
616 0 : lmini3 = abs(ltil2 - ljj2)
617 0 : lmaxi3 = ltil2 + ljj2
618 0 : lmx3 = min(lmaxi3,lnonsph(n))
619 0 : DO mjj2 = -ljj2,ljj2
620 0 : m = mjj2+mtil2
621 0 : IF(abs(m).GT.lmx3) CYCLE
622 0 : lmjj2 = lljj2 + mjj2
623 0 : fac3 = fac2 * ciljj2 * conjg(yl2(lmjj2 + 1))
624 0 : if(fac3.eq.0.0) CYCLE
625 :
626 0 : DO l = lmini3,lmx3
627 :
628 : ! Gaunt conditions (G3):
629 0 : if((mod(lmaxi3+l,2).eq.1).or.
630 : > (abs(m).gt.l)) cycle
631 :
632 0 : ll = l * (l + 1)
633 0 : lm = ll + m
634 :
635 : gc3 = gaunt1(l, ljj2, ltil2,
636 0 : > m, mjj2, mtil2, lmaxd)
637 :
638 : ! set up prefactor
639 0 : cil = fac3 * (cci**l) * gc3
640 :
641 : ! compute T-coefficients
642 : tuu(lm, lmp, n, indexx)
643 : > = tuu(lm, lmp, n, indexx)
644 0 : > + cil * uvu_non(lh, l, lp, ljj2, ljj1)
645 :
646 : tdd(lm, lmp, n, indexx)
647 : > = tdd(lm, lmp, n, indexx)
648 0 : > + cil * dvd_non(lh, l, lp, ljj2, ljj1)
649 :
650 : tud(lm, lmp, n, indexx)
651 : > = tud(lm, lmp, n, indexx)
652 0 : > + cil * uvd_non(lh, l, lp, ljj2, ljj1)
653 :
654 : tdu(lm, lmp, n, indexx)
655 : > = tdu(lm, lmp, n, indexx)
656 0 : > + cil * dvu_non(lh, l, lp, ljj2, ljj1)
657 : ENDDO !l
658 : ENDDO !mjj2
659 : ENDDO !ljj2
660 : ENDDO !lp
661 : ENDDO !mjj1
662 : ENDDO !mjj1
663 : ENDDO !ltil1
664 : ENDDO !mem
665 : ENDDO !lh
666 : ENDDO !mtil2
667 : ENDDO !ltil2
668 :
669 : 555 CONTINUE
670 0 : call cpu_time(t0)
671 0 : t_gggint = t_gggint + t0-t1
672 :
673 : c************** LOCAL ORBITAL CONTRIBUTION ***************c
674 : c determine contributions to T-matrix due to loc. orb. c
675 : c*********************************************************c
676 0 : IF((nlo(n).GE.1).AND.(.NOT.l_skip_loc)) THEN
677 : call wann_uHu_tlo(
678 : > memd,nlhd,ntypsd,jmtd,lmaxd,dx(n),rmsh(1,n),jri(n),
679 : > lmax(n),ntypsy,natd,lnonsph(n),lmd,lmplmd,clnu,
680 : > mlh,nmem,llh,nlh,irank,vr(1,0,n),nlod,llod,
681 : > loplod,ello(1,n),
682 : > llo(1,n),nlo(n),lo1l(0,n),l_dulo(1,n),ulo_der(1,n),
683 : > f(n,1:,1:,0:),g(n,1:,1:,0:),flo(n,1:,1:,1:),
684 : > f_b(n,1:,1:,0:),g_b(n,1:,1:,0:),flo_b(n,1:,1:,1:),
685 : > vr0(1,n),epar(0,n),l_skip_sph,
686 : > l_nns,rk1,rk2,invsfct,yl1,yl2,jj1,jj2,
687 : > p,dp,p_b,dp_b,q,dq,q_b,dq_b,
688 : > tuulo(0,1,-llod,n,indexx),
689 : > tdulo(0,1,-llod,n,indexx),
690 : > tulou(0,1,-llod,n,indexx),
691 : > tulod(0,1,-llod,n,indexx),
692 0 : > tuloulo(1,-llod,1,-llod,n,indexx))
693 0 : call cpu_time(t1)
694 0 : t_lo = t_lo + t1-t0
695 : ENDIF
696 :
697 0 : na = na + neq(n)
698 0 : 210 ENDDO !loop over atoms
699 :
700 : c write(*,*)ikpt_b,ikpt_b2
701 : c write(*,*)'t_radint=',t_radint-tt1
702 : c write(*,*)'t_sphint=',t_sphint-tt2
703 : c write(*,*)'t_gggint=',t_gggint-tt3
704 : c tt1 = t_radint
705 : c tt2 = t_sphint
706 : c tt3 = t_gggint
707 :
708 : enddo !ikpt_b2
709 : enddo !ikpt_b
710 :
711 : if(.false. .and. (irank.eq.0)) then
712 : do ikpt_b=1,nntot
713 : do ikpt_b2=1,nntot2
714 : indexx = ikpt_b2 + (ikpt_b-1)*nntot2
715 : indexx2= ikpt_b + (ikpt_b2-1)*nntot
716 : do n = 1,ntype
717 : do l=0,lwn
718 : do m=-l,l
719 : lm=l*(l+1)+m
720 : do lp=0,lwn
721 : do mp=-lp,lp
722 : lmp=lp*(lp+1)+mp
723 :
724 : if(abs(tuu(lm,lmp,n,indexx)-conjg(tuu(lmp,lm,n,indexx2)))
725 : > .gt.1.e-10) then
726 : write(*,*)'tuu',lm,lmp,n,ikpt_b,ikpt_b2
727 : write(*,*)tuu(lm,lmp,n,indexx)
728 : write(*,*)tuu(lmp,lm,n,indexx2)
729 : endif
730 : if(abs(tdd(lm,lmp,n,indexx)-conjg(tdd(lmp,lm,n,indexx2)))
731 : > .gt.1.e-10) then
732 : write(*,*)'tdd',lm,lmp,n,ikpt_b,ikpt_b2
733 : write(*,*)tdd(lm,lmp,n,indexx)
734 : write(*,*)tdd(lmp,lm,n,indexx2)
735 : endif
736 : if(abs(tud(lm,lmp,n,indexx)-conjg(tdu(lmp,lm,n,indexx2)))
737 : > .gt.1.e-10) then
738 : write(*,*)'tud',lm,lmp,n,ikpt_b,ikpt_b2
739 : write(*,*)tud(lm,lmp,n,indexx)
740 : write(*,*)tdu(lmp,lm,n,indexx2)
741 : endif
742 : if(abs(tdu(lm,lmp,n,indexx)-conjg(tud(lmp,lm,n,indexx2)))
743 : > .gt.1.e-10) then
744 : write(*,*)'tdu',lm,lmp,n,ikpt_b,ikpt_b2
745 : write(*,*)tdu(lm,lmp,n,indexx)
746 : write(*,*)tud(lmp,lm,n,indexx2)
747 : endif
748 :
749 : enddo
750 : enddo
751 : enddo
752 : enddo
753 : enddo
754 : enddo
755 : enddo
756 : endif
757 :
758 : #if (defined(CPP_MPI) && !defined(CPP_T90))
759 0 : CALL MPI_BARRIER(MPI_COMM_WORLD,ierr(1))
760 : #endif
761 :
762 0 : deallocate( yl1, yl2 )
763 0 : deallocate( dvd_sph, dvu_sph, uvd_sph, uvu_sph )
764 0 : deallocate( dvd_non, dvu_non, uvd_non, uvu_non )
765 0 : deallocate( p, p_b, dp_b, q, q_b, dq_b )
766 0 : deallocate( dp, dq )
767 0 : deallocate( x )
768 0 : deallocate( jlpp, jj1, jj2 )
769 :
770 0 : if(irank.eq.0) then
771 0 : write(*,*)'t_sphbes=',t_sphbes
772 0 : write(*,*)'t_ylm =',t_ylm
773 0 : write(*,*)'t_radint=',t_radint
774 0 : write(*,*)'t_sphint=',t_sphint
775 0 : write(*,*)'t_gggint=',t_gggint
776 0 : write(*,*)'t_locorb=',t_lo
777 : endif
778 :
779 0 : END SUBROUTINE wann_uHu_tlmplm
780 : END MODULE m_wann_uHu_tlmplm
|