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