Line data Source code
1 : c*****************************************c
2 : c Spin-orbit coupling cont. to uHu c
3 : c < u_{k+b1} | V_{k}^{soc} | 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 J.-P. Hanke, Dec. 2015 c
9 : c*****************************************c
10 : MODULE m_wann_uHu_soc
11 : #if (defined(CPP_MPI) && !defined(CPP_T90))
12 : use mpi
13 : #endif
14 : CONTAINS
15 0 : SUBROUTINE wann_uHu_soc(
16 : > input,atoms,
17 0 : > ntypd,jmtd,lmaxd,jspd,ntype,dx,rmsh,jri,lmax,
18 : > natd,lmd,lmplmd,neq,irank,nlod,llod,loplod,
19 0 : > ello,llo,nlo,lo1l,l_dulo,ulo_der,f,g,flo,f_b,
20 0 : > g_b,flo_b,kdiff,kdiff2,nntot,nntot2,
21 0 : > bmat,bbmat,vr,epar,jspin,
22 0 : > jspin_b,jspins,l_spav,theta,phi,alph,beta,
23 : > l_noco,l_skip_loc,
24 0 : < tuu,tud,tdu,tdd,tuulo,tulou,tdulo,tulod,tuloulo)
25 :
26 : USE m_intgr, ONLY : intgr0
27 : USE m_sphbes
28 : USE m_ylm
29 : USE m_gaunt, ONLY: gaunt1
30 : USE m_sointg
31 : USE m_constants
32 : USE m_wann_uHu_soc_tlo
33 : USE m_types
34 : USE m_anglso
35 :
36 : IMPLICIT NONE
37 :
38 : REAL :: sgml
39 : EXTERNAL sgml
40 : C .. Intrinsic Functions ..
41 : INTRINSIC abs,cmplx,max,mod
42 :
43 : TYPE(t_input),INTENT(IN) :: input
44 : TYPE(t_atoms),INTENT(IN) :: atoms
45 :
46 : C .. Scalar Arguments ..
47 : LOGICAL, INTENT (IN) :: l_noco,l_skip_loc
48 : INTEGER, INTENT (IN) :: ntypd,jmtd,lmaxd,jspd,jspins
49 : INTEGER, INTENT (IN) :: lmd,lmplmd,ntype,irank,jspin,jspin_b
50 : INTEGER, INTENT (IN) :: nlod,llod,loplod,natd
51 : INTEGER, INTENT (IN) :: nntot,nntot2
52 : REAL, INTENT (IN) :: theta,phi
53 : C ..
54 : C .. Array Arguments ..
55 : COMPLEX, INTENT (INOUT)::
56 : > tdd(0:lmd,0:lmd,ntypd,nntot*nntot2),
57 : > tdu(0:lmd,0:lmd,ntypd,nntot*nntot2),
58 : > tud(0:lmd,0:lmd,ntypd,nntot*nntot2),
59 : > tuu(0:lmd,0:lmd,ntypd,nntot*nntot2)
60 : COMPLEX, INTENT (INOUT)::
61 : > tdulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
62 : > tuulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
63 : > tulou(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
64 : > tulod(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
65 : > tuloulo(nlod,-llod:llod,nlod,-llod:llod,ntypd,nntot*nntot2)
66 : INTEGER, INTENT (IN) :: llo(nlod,ntypd),nlo(ntypd)
67 : INTEGER, INTENT (IN) :: lo1l(0:llod,ntypd),neq(ntypd)
68 : INTEGER, INTENT (IN) :: jri(ntypd),lmax(ntypd)
69 : INTEGER, INTENT (IN) :: ulo_der(nlod,ntypd)
70 : REAL, INTENT (IN) :: alph(ntypd),beta(ntypd)
71 : REAL, INTENT (IN) :: dx(ntypd),rmsh(jmtd,ntypd)
72 : REAL, INTENT (IN) :: ello(nlod,ntypd,max(jspd,2))
73 : REAL, INTENT (IN) :: epar(0:lmaxd,ntypd,max(jspd,2))
74 : REAL, INTENT (IN) :: vr(jmtd,ntypd,jspd)
75 : REAL, INTENT (IN) :: f (ntypd,jmtd,2,0:lmaxd),
76 : > g (ntypd,jmtd,2,0:lmaxd),
77 : > f_b(ntypd,jmtd,2,0:lmaxd),
78 : > g_b(ntypd,jmtd,2,0:lmaxd),
79 : > flo (ntypd,jmtd,2,nlod),
80 : > flo_b(ntypd,jmtd,2,nlod)
81 : REAL, INTENT (IN) :: bmat(3,3),bbmat(3,3)
82 : REAL, INTENT (IN) :: kdiff (3,nntot)
83 : REAL, INTENT (IN) :: kdiff2(3,nntot2)
84 : LOGICAL, INTENT (IN) :: l_dulo(nlod,ntypd),l_spav
85 : C ..
86 : C .. Local Scalars ..
87 : REAL gs,rk1,rk2,socscreen_fac
88 : REAL gc1,gc2,c,e,r0
89 : REAL t0,t1,t11,t22,t33,t44,t55,t66
90 : COMPLEX cil,cci,cciljj1,ciljj2,fac1,fac2,fac3,dump1,dump2
91 : INTEGER i,l,l2,lh,lm,lmp,lp
92 : INTEGER lp1,lpl,m,mp,mu,n,nh,noded,nodeu,nrec,na,ne
93 : INTEGER ljj1,ljj2,mjj1,mjj2
94 : INTEGER ltil1,mtil1,mtil2
95 : INTEGER mlo, mlolo, iu, indexx
96 : INTEGER lmini,lmini2, lmaxi,lmaxi2,lmx,lmx2
97 : INTEGER ll,llp,lljj1,lljj2,lmjj1,lmjj2
98 : INTEGER lltil1,lmtil1,lmtil2,sp1,sp2
99 : INTEGER mmin,mmax,lwn, ikpt_b,ikpt_b2
100 : LOGICAL l_socscreen
101 : C ..
102 : C .. Local Arrays ..
103 0 : REAL, ALLOCATABLE :: dvd(:,:,:,:,:),dvu(:,:,:,:,:)
104 0 : REAL, ALLOCATABLE :: uvd(:,:,:,:,:),uvu(:,:,:,:,:)
105 0 : REAL, ALLOCATABLE :: p(:,:,:),p_b(:,:,:),q(:,:,:),q_b(:,:,:)
106 0 : REAL, ALLOCATABLE :: x(:)
107 0 : REAL, ALLOCATABLE :: jlpp(:),jj1(:,:),jj2(:,:)
108 0 : REAL, ALLOCATABLE :: v0(:),vso(:,:)
109 : REAL :: bpt(3),bpt2(3),bkrot(3)
110 0 : COMPLEX, ALLOCATABLE :: yl1(:),yl2(:)
111 0 : COMPLEX, ALLOCATABLE :: angso(:,:,:)
112 : INTEGER :: ispjsp(2)
113 : DATA ispjsp/1,-1/
114 :
115 : #if (defined(CPP_MPI) && !defined(CPP_T90))
116 : INTEGER ierr(3)
117 : #endif
118 :
119 0 : t11 = 0.; t22 = 0.; t33 = 0.; t44 = 0.; t55 = 0.; t66 = 0.
120 0 : cci = cmplx(0.0,-1.0)
121 : c = c_light(1.0)
122 0 : sp1 = ispjsp(jspin)
123 0 : sp2 = ispjsp(jspin_b)
124 :
125 0 : l_socscreen = .false. ; socscreen_fac = 1.0
126 0 : INQUIRE(file='socscreen',exist=l_socscreen)
127 0 : IF(l_socscreen) THEN
128 0 : OPEN(757,file='socscreen')
129 0 : READ(757,*)socscreen_fac
130 0 : CLOSE(757)
131 0 : IF(irank.eq.0) write(*,*)'socscreen with factor',socscreen_fac
132 : ENDIF
133 :
134 0 : allocate( dvd(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
135 0 : allocate( dvu(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
136 0 : allocate( uvd(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
137 0 : allocate( uvu(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
138 0 : allocate( p(jmtd,0:lmaxd,0:lmaxd), p_b(jmtd,0:lmaxd,0:lmaxd) )
139 0 : allocate( q(jmtd,0:lmaxd,0:lmaxd), q_b(jmtd,0:lmaxd,0:lmaxd) )
140 0 : allocate( x(jmtd) )
141 0 : allocate( jlpp(0:lmaxd), jj1(0:lmaxd,jmtd), jj2(0:lmaxd,jmtd) )
142 0 : allocate( v0(jmtd), vso(jmtd,2) )
143 0 : allocate( angso(0:lmd,0:lmd,ntypd) )
144 0 : allocate( yl1((lmaxd+1)**2), yl2((lmaxd+1)**2) )
145 :
146 0 : call cpu_time(t0)
147 0 : angso = cmplx(0.,0.)
148 : c*******************************************c
149 : c set up < l m sp1 | sigma L | l mp sp2 > c
150 : c*******************************************c
151 0 : DO n=1,ntype
152 0 : DO l=1,lmaxd ! TODO: exclude s-states?
153 0 : ll=l*(l+1)
154 0 : DO m=-l,l
155 0 : lm = ll+m
156 0 : DO mp=-l,l
157 0 : lmp = ll+mp
158 0 : if(l_noco) then
159 0 : angso(lm,lmp,n) = anglso(beta(n),alph(n),l,m,sp1,l,mp,sp2)
160 : else
161 0 : angso(lm,lmp,n) = anglso(theta,phi,l,m,sp1,l,mp,sp2)
162 : endif
163 : ENDDO
164 : ENDDO
165 : ENDDO
166 : ENDDO
167 0 : call cpu_time(t1)
168 0 : t11 = t11 + t1-t0
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 :
178 : c***************************************c
179 : c begin 2nd neighbor loop: |k+b2> c
180 : c***************************************c
181 0 : do ikpt_b2=1,nntot2
182 0 : indexx = ikpt_b2 + (ikpt_b-1)*nntot2
183 :
184 0 : bpt2 = kdiff2(:,ikpt_b2)
185 0 : rk2 = sqrt(dot_product(bpt2,matmul(bbmat,bpt2)))
186 : ! rk2= sqrt(dotirp(bpt2,bpt2,bbmat))
187 :
188 : ! loop over atoms
189 0 : DO 210 n = 1,ntype
190 0 : lwn = lmax(n)
191 0 : r0 = rmsh(1,n)
192 :
193 : ! set up spherical Bessel functions
194 : ! jj1(l,b1*r) and jj2(l,b2*r)
195 0 : call cpu_time(t0)
196 0 : do i=1,jri(n)
197 0 : gs = rk1*rmsh(i,n)
198 0 : call sphbes(lwn,gs,jlpp)
199 0 : jj1(:,i) = jlpp(:)
200 :
201 0 : gs = rk2*rmsh(i,n)
202 0 : call sphbes(lwn,gs,jlpp)
203 0 : jj2(:,i) = jlpp(:)
204 : enddo
205 0 : call cpu_time(t1)
206 0 : t22 = t22 + t1-t0
207 :
208 : ! set up spherical harmonics
209 : ! yl1(lm) for b1 and yl2(lm) for b2
210 0 : yl1 = cmplx(0.,0.)
211 0 : bkrot=MATMUL(bpt,bmat)
212 0 : call ylm4(lwn,bkrot,yl1)
213 :
214 0 : yl2 = cmplx(0.,0.)
215 0 : bkrot=MATMUL(bpt2,bmat)
216 0 : call ylm4(lwn,bkrot,yl2)
217 0 : call cpu_time(t0)
218 0 : t33 = t33 + t0-t1
219 :
220 : ! set up products of radial functions and sph. Bessel
221 0 : DO ljj1 = 0, lwn
222 0 : DO l = 0, lwn
223 0 : DO i=1,jri(n)
224 0 : p(i,l,ljj1) = f(n,i,1,l)*jj1(ljj1,i)
225 0 : q(i,l,ljj1) = g(n,i,1,l)*jj1(ljj1,i)
226 0 : p_b(i,l,ljj1) = f_b(n,i,1,l)*jj2(ljj1,i)
227 0 : q_b(i,l,ljj1) = g_b(n,i,1,l)*jj2(ljj1,i)
228 : ENDDO
229 : ENDDO
230 : ENDDO
231 :
232 : c****************************************************c
233 : c compute radial integrals c
234 : c <u(l')jj1(ljj1) | vso(lh) | u(l)jj2(ljj2)> c
235 : c****************************************************c
236 0 : uvu = 0.0; uvd = 0.0; dvu = 0.0; dvd = 0.0
237 0 : DO lp = 0, lwn
238 0 : DO l = 0, lwn
239 :
240 : ! construct vso
241 0 : v0 = 0.
242 0 : IF (jspins.EQ.1) THEN
243 0 : v0(1:jri(n)) = vr(1:jri(n),n,1)
244 0 : e = ( epar(l,n,1) + epar(lp,n,1) )/2.
245 : ELSE
246 0 : DO i = 1,jri(n)
247 0 : v0(i) = (vr(i,n,1)+vr(i,n,jspins))/2.
248 : END DO
249 : e = ( epar(l,n,1)+epar(l,n,jspins)
250 0 : > +epar(lp,n,1)+epar(lp,n,jspins) )/4.
251 : END IF
252 :
253 0 : CALL sointg(n,e,vr(:,n,:),v0,atoms,input,vso)
254 :
255 : ! spin-averaging
256 0 : if(l_spav)then
257 0 : DO i= 1,jmtd
258 0 : vso(i,1)= (vso(i,1)+vso(i,2))/2.
259 0 : vso(i,2)= vso(i,1)
260 : ENDDO
261 : endif
262 :
263 0 : if(l_socscreen) vso(:,:) = vso(:,:)*socscreen_fac
264 :
265 : ! calculate radial integrals
266 0 : DO ljj1 = 0, lwn
267 0 : lmini = abs(lp-ljj1)
268 0 : lmaxi = lp+ljj1
269 0 : lmx = min(lmaxi,lwn)
270 0 : DO ljj2 = 0, lwn
271 0 : DO lh = lmini, lmx
272 0 : lmini2= abs(lh-ljj2)
273 0 : lmaxi2= lh+ljj2
274 : if((lmini2.gt.l) .or. (lmaxi2.lt.l).or.
275 0 : > (mod(lmaxi+lh,2).eq.1).or.
276 : > (mod(lmaxi2+l,2).eq.1) ) CYCLE
277 :
278 0 : DO i=1,jri(n)
279 0 : x(i) = p(i,lp,ljj1)*p_b(i,l,ljj2)*vso(i,jspin)
280 : ENDDO
281 : call intgr0(x(1:jri(n)),r0,dx(n),jri(n),
282 0 : > uvu(lh,l,lp,ljj2,ljj1))
283 :
284 0 : DO i=1,jri(n)
285 0 : x(i) = q(i,lp,ljj1)*p_b(i,l,ljj2)*vso(i,jspin)
286 : ENDDO
287 : call intgr0(x(1:jri(n)),r0,dx(n),jri(n),
288 0 : > dvu(lh,l,lp,ljj2,ljj1))
289 :
290 0 : DO i=1,jri(n)
291 0 : x(i) = p(i,lp,ljj1)*q_b(i,l,ljj2)*vso(i,jspin)
292 : ENDDO
293 : call intgr0(x(1:jri(n)),r0,dx(n),jri(n),
294 0 : > uvd(lh,l,lp,ljj2,ljj1))
295 :
296 0 : DO i=1,jri(n)
297 0 : x(i) = q(i,lp,ljj1)*q_b(i,l,ljj2)*vso(i,jspin)
298 : ENDDO
299 : call intgr0(x(1:jri(n)),r0,dx(n),jri(n),
300 0 : > dvd(lh,l,lp,ljj2,ljj1))
301 : ENDDO
302 : ENDDO
303 : ENDDO
304 : ENDDO!l
305 : ENDDO!lp
306 0 : call cpu_time(t1)
307 0 : t44 = t44 +t1-t0
308 :
309 : c***************** SOC CONTRIBUTION **********************c
310 : c compute product of the two Gaunt coefficients c
311 : c with the radial integrals (+prefactors) c
312 : c*********************************************************c
313 : c We deal with two Gaunt coefficients: c
314 : c G1 = G( (ltil1, mtil1), (ljj1, mjj1), (lp, mp) ) c
315 : c G2 = G( (l , m) , (ljj2, mjj2), (ltil1, mtil2) ) c
316 : c*********************************************************c
317 : c use Gaunt conditions to reduce number of operations. c
318 : c coefficient G(L1,L2,L3) only nonzero if c
319 : c a) l1 + l2 + l3 = even c
320 : c b) |l2-l3| <= l1 <= l2+l3 c
321 : c c) m1 = m2 + m3 c
322 : c*********************************************************c
323 0 : DO ljj1=0,lwn
324 0 : lljj1 = ljj1 * (ljj1 + 1)
325 0 : cciljj1 = cci**ljj1
326 0 : DO mjj1 = -ljj1,ljj1
327 0 : lmjj1 = lljj1 + mjj1
328 0 : fac1 = cciljj1 * conjg(yl1(lmjj1+1))
329 0 : IF(fac1.eq.0.0) CYCLE
330 :
331 0 : DO lp=0,lwn
332 0 : llp = lp*(lp+1)
333 0 : lmini = abs(lp-ljj1)
334 0 : lmaxi = lp+ljj1
335 0 : lmx = min(lmaxi,lwn)
336 0 : DO mp=-lp,lp
337 0 : lmp = llp+mp
338 0 : mtil1=mjj1+mp
339 0 : IF(abs(mtil1).GT.lmx) CYCLE
340 :
341 0 : DO ltil1=lmini,lmx
342 : ! Gaunt conditions (G1):
343 0 : if((mod(lmaxi+ltil1,2).eq.1).or.
344 : > (abs(mtil1).gt.ltil1)) cycle
345 :
346 0 : lltil1 = ltil1*(ltil1+1)
347 0 : lmtil1 = lltil1+mtil1
348 0 : mmin = max(-ltil1,mtil1-1)
349 0 : mmax = min( ltil1,mtil1+1)
350 :
351 : gc1 = gaunt1(ltil1, ljj1, lp,
352 0 : > mtil1, mjj1, mp, lmaxd)
353 0 : fac2 = fac1 * gc1 * (ImagUnit**lp)
354 :
355 0 : DO ljj2=0,lwn
356 0 : lljj2 = ljj2 * (ljj2 + 1)
357 0 : ciljj2 = ImagUnit**ljj2
358 0 : lmini2 = abs(ltil1-ljj2)
359 0 : lmaxi2 = ltil1+ljj2
360 0 : lmx2 = min(lmaxi2,lwn)
361 0 : DO mjj2 = -ljj2,ljj2
362 0 : lmjj2 = lljj2 + mjj2
363 0 : fac3 = fac2 * ciljj2 * conjg(yl2(lmjj2+1))
364 0 : IF(fac3.eq.0.0) CYCLE
365 :
366 0 : DO mtil2=mmin,mmax ! mtil1-1 <= mtil2 <= mtil1+1
367 0 : lmtil2 = lltil1+mtil2
368 0 : m=mjj2+mtil2
369 0 : IF(abs(m).GT.lmx2) cycle
370 :
371 0 : DO l=lmini2,lmx2
372 : ! Gaunt conditions (G2):
373 0 : if((mod(lmaxi2+l,2).eq.1).or.
374 : > (abs(m).gt.l)) cycle
375 0 : ll = l*(l+1)
376 0 : lm = ll+m
377 :
378 : gc2 = gaunt1( l, ljj2, ltil1,
379 0 : > m, mjj2, mtil2, lmaxd)
380 :
381 : ! set up prefactor
382 0 : cil = fac3 * (cci**l) * gc2 * conjg(angso(lmtil1,lmtil2,n))
383 :
384 : ! compute T-coefficients
385 : tuu(lm, lmp, n, indexx)
386 : > = tuu(lm, lmp, n, indexx)
387 0 : > + cil * uvu(ltil1, l, lp, ljj2, ljj1)
388 :
389 : tdd(lm, lmp, n, indexx)
390 : > = tdd(lm, lmp, n, indexx)
391 0 : > + cil * dvd(ltil1, l, lp, ljj2, ljj1)
392 :
393 : tud(lm, lmp, n, indexx)
394 : > = tud(lm, lmp, n, indexx)
395 0 : > + cil * uvd(ltil1, l, lp, ljj2, ljj1)
396 :
397 : tdu(lm, lmp, n, indexx)
398 : > = tdu(lm, lmp, n, indexx)
399 0 : > + cil * dvu(ltil1, l, lp, ljj2, ljj1)
400 : ENDDO!ljj2
401 : ENDDO!m
402 : ENDDO!l
403 : ENDDO!mtil2
404 : ENDDO!ljj1
405 : ENDDO!mp
406 : ENDDO!lp
407 : ENDDO!mtil1
408 : ENDDO!ltil1
409 0 : call cpu_time(t0)
410 0 : t55 = t55 + t0-t1
411 :
412 : c************** LOCAL ORBITAL CONTRIBUTION ***************c
413 : c determine contributions to T-matrix due to loc. orb. c
414 : c*********************************************************c
415 0 : IF((nlo(n).GE.1).AND.(.NOT.l_skip_loc)) THEN
416 :
417 : call wann_uHu_soc_tlo(
418 : > input,atoms,
419 : > ntypd,jmtd,lmaxd,jspd,ntype,dx,rmsh,jri,lmax,
420 : > natd,lmd,irank,nlod,llod,loplod,ello,llo,nlo,
421 : > lo1l,l_dulo,ulo_der,flo,flo_b,kdiff,kdiff2,
422 : > nntot,nntot2,
423 : > vr,epar,jspin,jspin_b,jspins,l_spav,theta,phi,
424 : > yl1,yl2,jj1,jj2,p,p_b,q,q_b,n,angso(0,0,n),
425 : > l_socscreen,socscreen_fac,
426 : > tuulo(0,1,-llod,n,indexx),
427 : > tdulo(0,1,-llod,n,indexx),
428 : > tulou(0,1,-llod,n,indexx),
429 : > tulod(0,1,-llod,n,indexx),
430 0 : > tuloulo(1,-llod,1,-llod,n,indexx))
431 :
432 0 : call cpu_time(t1)
433 0 : t66 = t66 + t1-t0
434 : ENDIF
435 :
436 0 : 210 ENDDO !loop over atoms
437 :
438 : enddo !ikpt_b2
439 : enddo !ikpt_b
440 :
441 : #if (defined(CPP_MPI) && !defined(CPP_T90))
442 0 : CALL MPI_BARRIER(MPI_COMM_WORLD,ierr(1))
443 : #endif
444 :
445 0 : deallocate( dvd, dvu, uvd, uvu )
446 0 : deallocate( p, p_b, q, q_b )
447 0 : deallocate( x )
448 0 : deallocate( jlpp, jj1, jj2 )
449 0 : deallocate( yl1, yl2 )
450 0 : deallocate( v0, vso, angso )
451 :
452 0 : if(irank.eq.0) then
453 0 : write(*,*)'t_angso =',t11
454 0 : write(*,*)'t_sphbes=',t22
455 0 : write(*,*)'t_ylm =',t33
456 0 : write(*,*)'t_radint=',t44
457 0 : write(*,*)'t_gaunt =',t55
458 0 : write(*,*)'t_locorb=',t66
459 : endif
460 :
461 0 : END SUBROUTINE wann_uHu_soc
462 : END MODULE m_wann_uHu_soc
|
