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 due to local orbitals c
8 : c*****************************************c
9 : c J.-P. Hanke, Dec. 2015 c
10 : c*****************************************c
11 : MODULE m_wann_uHu_soc_tlo
12 : CONTAINS
13 0 : SUBROUTINE wann_uHu_soc_tlo(
14 : > input,atoms,
15 : > ntypd,jmtd,lmaxd,jspd,
16 0 : > ntype,dx,rmsh,jri,lmax,natd,
17 : > lmd,irank,nlod,llod,
18 0 : > loplod,ello,llo,nlo,lo1l,l_dulo,ulo_der,
19 0 : > flo,flo_b,kdiff,kdiff2,nntot,nntot2,
20 0 : > vr,epar,jspin,jspin_b,jspins,l_spav,theta,phi,
21 0 : > yl1,yl2,jj1,jj2,p,p_b,q,q_b,ntyp,angso,
22 : > l_socscreen,socscreen_fac,
23 0 : < tuulo,tdulo,tulou,tulod,tuloulo)
24 :
25 : USE m_intgr, ONLY : intgr0
26 : USE m_sphbes
27 : USE m_ylm
28 : USE m_gaunt, ONLY: gaunt1
29 : USE m_sointg
30 : USE m_constants
31 : USE m_types
32 : #if (defined(CPP_MPI) && !defined(CPP_T90))
33 : use mpi
34 : #endif
35 :
36 : IMPLICIT NONE
37 : C .. Intrinsic Functions ..
38 : INTRINSIC abs,cmplx,max,mod
39 :
40 : TYPE(t_input),INTENT(IN) :: input
41 : TYPE(t_atoms),INTENT(IN) :: atoms
42 :
43 : C ..
44 : C .. Scalar Arguments ..
45 : INTEGER, INTENT (IN) :: ntypd,jmtd,lmaxd,jspd,jspins
46 : INTEGER, INTENT (IN) :: lmd,ntype,irank,jspin,jspin_b
47 : INTEGER, INTENT (IN) :: nlod,llod,loplod,natd
48 : INTEGER, INTENT (IN) :: nntot,nntot2,ntyp
49 : REAL, INTENT (IN) :: theta,phi,socscreen_fac
50 : LOGICAL, INTENT (IN) :: l_socscreen
51 : C ..
52 : C .. Array Arguments ..
53 : COMPLEX, INTENT (INOUT)::
54 : > tdulo(0:lmd,nlod,-llod:llod),
55 : > tuulo(0:lmd,nlod,-llod:llod),
56 : > tulou(0:lmd,nlod,-llod:llod),
57 : > tulod(0:lmd,nlod,-llod:llod),
58 : > tuloulo(nlod,-llod:llod,nlod,-llod:llod)
59 : COMPLEX, INTENT (IN) :: angso(0:lmaxd,0:lmaxd)
60 : COMPLEX, INTENT (IN) :: yl1((lmaxd+1)**2),yl2((lmaxd+1)**2)
61 : INTEGER, INTENT (IN) :: llo(nlod,ntypd),nlo(ntypd)
62 : INTEGER, INTENT (IN) :: lo1l(0:llod,ntypd)
63 : INTEGER, INTENT (IN) :: jri(ntypd),lmax(ntypd)
64 : INTEGER, INTENT (IN) :: ulo_der(nlod,ntypd)
65 : REAL, INTENT (IN) :: dx(ntypd),rmsh(jmtd,ntypd)
66 : REAL, INTENT (IN) :: ello(nlod,ntypd,max(jspd,2))
67 : REAL, INTENT (IN) :: epar(0:lmaxd,ntypd,max(jspd,2))
68 : REAL, INTENT (IN) :: vr(jmtd,ntypd,jspd)
69 : REAL, INTENT (IN) :: flo (ntypd,jmtd,2,nlod),
70 : > flo_b(ntypd,jmtd,2,nlod)
71 : REAL, INTENT (IN) :: p(jmtd,0:lmaxd,0:lmaxd),
72 : > q(jmtd,0:lmaxd,0:lmaxd),
73 : > p_b(jmtd,0:lmaxd,0:lmaxd),
74 : > q_b(jmtd,0:lmaxd,0:lmaxd)
75 : REAL, INTENT (IN) :: kdiff (3,nntot)
76 : REAL, INTENT (IN) :: kdiff2(3,nntot2)
77 : REAL, INTENT (IN) :: jj1(0:lmaxd,jmtd),jj2(0:lmaxd,jmtd)
78 : LOGICAL, INTENT (IN) :: l_dulo(nlod,ntypd),l_spav
79 : C ..
80 : C .. Local Scalars ..
81 : REAL gc1,gc2,c,e,r0
82 : COMPLEX cil
83 : INTEGER i,l,l2,lh,lm,lmp,lp,lo,lop
84 : INTEGER lp1,lpl,m,mp,mu
85 : INTEGER ljj1,ljj2,mjj1,mjj2
86 : INTEGER ltil1,mtil1,mtil2
87 : INTEGER lmini,lmini2, lmaxi,lmaxi2
88 : INTEGER ll,llp,lljj1,lljj2,lmjj1,lmjj2
89 : INTEGER lltil1,lmtil1,lmtil2,sp1,sp2
90 : INTEGER mmin,mmax,lwn, ikpt_b,ikpt_b2
91 : C ..
92 : C .. Local Arrays ..
93 0 : REAL, ALLOCATABLE :: dvulo(:,:,:,:,:),ulovd(:,:,:,:,:)
94 0 : REAL, ALLOCATABLE :: uvulo(:,:,:,:,:),ulovu(:,:,:,:,:)
95 0 : REAL, ALLOCATABLE :: ulovulo(:,:,:,:,:)
96 0 : REAL, ALLOCATABLE :: plo(:,:,:),plo_b(:,:,:)
97 0 : REAL, ALLOCATABLE :: x(:)
98 0 : REAL, ALLOCATABLE :: v0(:),vso(:,:)
99 : INTEGER :: ispjsp(2)
100 : DATA ispjsp/1,-1/
101 :
102 : #if (defined(CPP_MPI) && !defined(CPP_T90))
103 : INTEGER ierr(3)
104 : #endif
105 :
106 : c = c_light(1.0)
107 0 : sp1 = ispjsp(jspin)
108 0 : sp2 = ispjsp(jspin_b)
109 :
110 0 : allocate( ulovulo(0:lmaxd,0:lmaxd,0:lmaxd,nlod,nlod) )
111 0 : allocate( ulovu(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
112 0 : allocate( ulovd(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
113 0 : allocate( uvulo(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
114 0 : allocate( dvulo(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
115 0 : allocate( x(jmtd) )
116 0 : allocate( v0(jmtd), vso(jmtd,2) )
117 0 : allocate( plo(jmtd,0:lmaxd,nlod), plo_b(jmtd,0:lmaxd,nlod) )
118 :
119 0 : lwn = lmax(ntyp)
120 0 : r0 = rmsh(1,ntyp)
121 :
122 0 : DO lo = 1, nlo(ntyp)
123 0 : DO ljj1 = 0, lwn
124 0 : DO i=1,jri(ntyp)
125 0 : plo(i,ljj1,lo) = flo(ntyp,i,1,lo)*jj1(ljj1,i)
126 0 : plo_b(i,ljj1,lo) = flo_b(ntyp,i,1,lo)*jj2(ljj1,i)
127 : ENDDO
128 : ENDDO
129 : ENDDO
130 :
131 : c****************************************************c
132 : c compute radial integrals c
133 : c <u(l')jj1(ljj1) | vso(lh) | u(l)jj2(ljj2)> c
134 : c****************************************************c
135 :
136 : ! lo-lo
137 0 : DO lop = 1, nlo(ntyp)
138 0 : lp = llo(lop,ntyp)
139 0 : DO lo = 1, nlo(ntyp)
140 0 : l = llo(lop,ntyp)
141 :
142 : ! construct vso
143 0 : v0 = 0.
144 0 : IF (jspins.EQ.1) THEN
145 0 : v0(1:jri(ntyp)) = vr(1:jri(ntyp),ntyp,1)
146 0 : e = ( ello(lo,ntyp,1) + ello(lop,ntyp,1) )/2.
147 : ELSE
148 0 : DO i = 1,jri(ntyp)
149 0 : v0(i) = (vr(i,ntyp,1)+vr(i,ntyp,jspins))/2.
150 : END DO
151 : e = ( ello(lo,ntyp,1) + ello(lo,ntyp,jspins)
152 0 : > +ello(lop,ntyp,1) + ello(lop,ntyp,jspins) )/4.
153 : END IF
154 :
155 0 : CALL sointg(ntyp,e,vr(:,ntyp,:),v0,atoms,input,vso)
156 :
157 : ! spin-averaging
158 0 : if(l_spav)then
159 0 : DO i= 1,jmtd
160 0 : vso(i,1)= (vso(i,1)+vso(i,2))/2.
161 0 : vso(i,2)= vso(i,1)
162 : ENDDO
163 : endif
164 :
165 0 : if(l_socscreen) vso(:,:) = vso(:,:)*socscreen_fac
166 :
167 : ! calculate radial integrals
168 0 : DO ljj1 = 0, lwn
169 0 : DO ljj2 = 0, lwn
170 0 : DO lh=0,lwn
171 0 : lmini = abs(lp-ljj1)
172 0 : lmini2= abs(lh-ljj2)
173 0 : lmaxi = lp+ljj1
174 0 : lmaxi2= lh+ljj2
175 : if((lmini.gt.lh) .or. (lmaxi.lt.lh).or.
176 : > (lmini2.gt.l) .or. (lmaxi2.lt.l).or.
177 0 : > (mod(lp+lh+ljj1,2).ne.0).or.
178 0 : > (mod( l+lh+ljj2,2).ne.0) ) then
179 :
180 0 : ulovulo(lh,ljj2,ljj1,lo,lop) = 0.
181 :
182 : else
183 :
184 0 : DO i=1,jri(ntyp)
185 0 : x(i) = plo(i,ljj1,lop)*plo_b(i,ljj2,lo)*vso(i,jspin)
186 : ENDDO
187 : call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
188 0 : > ulovulo(lh,ljj2,ljj1,lo,lop))
189 :
190 : endif
191 : ENDDO
192 : ENDDO
193 : ENDDO
194 : ENDDO
195 : ENDDO!lh
196 :
197 : ! apw-lo
198 0 : DO lop = 1, nlo(ntyp)
199 0 : lp = llo(lop,ntyp)
200 0 : DO l = 0, lwn
201 :
202 : ! construct vso
203 0 : v0 = 0.
204 0 : IF (jspins.EQ.1) THEN
205 0 : v0(1:jri(ntyp)) = vr(1:jri(ntyp),ntyp,1)
206 0 : e = ( epar(l,ntyp,1) + ello(lop,ntyp,1) )/2.
207 : ELSE
208 0 : DO i = 1,jri(ntyp)
209 0 : v0(i) = (vr(i,ntyp,1)+vr(i,ntyp,jspins))/2.
210 : END DO
211 : e = ( epar(l,ntyp,1) + epar(l,ntyp,jspins)
212 0 : > +ello(lop,ntyp,1) + ello(lop,ntyp,jspins) )/4.
213 : END IF
214 :
215 0 : CALL sointg(ntyp,e,vr(:,ntyp,:),v0,atoms,input,vso)
216 :
217 : ! spin-averaging
218 0 : if(l_spav)then
219 0 : DO i= 1,jmtd
220 0 : vso(i,1)= (vso(i,1)+vso(i,2))/2.
221 0 : vso(i,2)= vso(i,1)
222 : ENDDO
223 : endif
224 :
225 : ! calculate radial integrals
226 0 : DO ljj1 = 0, lwn
227 0 : DO ljj2 = 0, lwn
228 0 : DO lh=0,lwn
229 0 : lmini = abs(lp-ljj1)
230 0 : lmini2= abs(lh-ljj2)
231 0 : lmaxi = lp+ljj1
232 0 : lmaxi2= lh+ljj2
233 : if((lmini.gt.lh) .or. (lmaxi.lt.lh).or.
234 : > (lmini2.gt.l) .or. (lmaxi2.lt.l).or.
235 0 : > (mod(lp+lh+ljj1,2).ne.0).or.
236 : > (mod( l+lh+ljj2,2).ne.0) ) then
237 :
238 0 : ulovu(lh,l,ljj2,ljj1,lop) = 0.
239 0 : ulovd(lh,l,ljj2,ljj1,lop) = 0.
240 :
241 : else
242 :
243 0 : DO i=1,jri(ntyp)
244 0 : x(i) = plo(i,ljj1,lop)*p_b(i,l,ljj2)*vso(i,jspin)
245 : ENDDO
246 : call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
247 0 : > ulovu(lh,l,ljj2,ljj1,lop))
248 0 : DO i=1,jri(ntyp)
249 0 : x(i) = plo(i,ljj1,lop)*q_b(i,l,ljj2)*vso(i,jspin)
250 : ENDDO
251 : call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
252 0 : > ulovd(lh,l,ljj2,ljj1,lop))
253 :
254 : endif
255 :
256 0 : lmini = abs(l-ljj1)
257 0 : lmini2= abs(lh-ljj2)
258 0 : lmaxi = l+ljj1
259 0 : lmaxi2= lh+ljj2
260 : if((lmini.gt.lh) .or. (lmaxi.lt.lh).or.
261 : > (lmini2.gt.lp) .or. (lmaxi2.lt.lp).or.
262 0 : > (mod( l+lh+ljj1,2).ne.0).or.
263 0 : > (mod(lp+lh+ljj2,2).ne.0) ) then
264 :
265 0 : uvulo(lh,l,ljj2,ljj1,lop) = 0.
266 0 : dvulo(lh,l,ljj2,ljj1,lop) = 0.
267 :
268 : else
269 :
270 0 : DO i=1,jri(ntyp)
271 0 : x(i) = p(i,l,ljj1)*plo_b(i,ljj2,lop)*vso(i,jspin)
272 : ENDDO
273 : call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
274 0 : > uvulo(lh,l,ljj2,ljj1,lop))
275 0 : DO i=1,jri(ntyp)
276 0 : x(i) = q(i,l,ljj1)*plo_b(i,ljj2,lop)*vso(i,jspin)
277 : ENDDO
278 : call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
279 0 : > dvulo(lh,l,ljj2,ljj1,lop))
280 :
281 : endif
282 : ENDDO
283 : ENDDO
284 : ENDDO
285 : ENDDO
286 : ENDDO!lh
287 :
288 : c***************** SOC CONTRIBUTION **********************c
289 : c compute product of the two Gaunt coefficients c
290 : c with the radial integrals (+prefactors) c
291 : c*********************************************************c
292 : c We deal with two Gaunt coefficients: c
293 : c G1 = G( (ltil1, mtil1), (ljj1, mjj1), (lp, mp) ) c
294 : c G2 = G( (l , m) , (ljj2, mjj2), (ltil1, mtil2) ) c
295 : c*********************************************************c
296 : c use Gaunt conditions to reduce number of operations. c
297 : c coefficient G(L1,L2,L3) only nonzero if c
298 : c a) l1 + l2 + l3 = even c
299 : c b) |l2-l3| <= l1 <= l2+l3 c
300 : c c) m1 = m2 + m3 c
301 : c*********************************************************c
302 :
303 : ! lo-lo
304 0 : DO lop=1,nlo(ntyp)
305 0 : lp = llo(lop,ntyp)
306 0 : DO mp=-lp,lp
307 0 : DO ltil1=0,lwn
308 0 : lltil1 = ltil1*(ltil1+1)
309 0 : DO mtil1=-ltil1,ltil1
310 0 : lmtil1 = lltil1+mtil1
311 0 : mmin = max(-ltil1,mtil1-1)
312 0 : mmax = min( ltil1,mtil1+1)
313 0 : mjj1=mtil1-mp
314 0 : DO ljj1=0,lwn
315 0 : lljj1 = ljj1 * (ljj1 + 1)
316 0 : lmjj1 = lljj1 + mjj1
317 0 : lmini = abs(lp-ljj1)
318 0 : lmaxi = lp+ljj1
319 : ! Gaunt conditions (G1):
320 : if((lmini.gt.ltil1).or.(lmaxi.lt.ltil1).or.
321 0 : > (mod(lp+ltil1+ljj1,2).ne.0).or.
322 : > (abs(mjj1).gt.ljj1)) cycle
323 :
324 : gc1 = gaunt1(ltil1, ljj1, lp,
325 0 : > mtil1, mjj1, mp, lmaxd)
326 :
327 0 : DO lo=1,nlo(ntyp)
328 0 : l = llo(lo,ntyp)
329 0 : DO m=-l,l
330 0 : DO mtil2=mmin,mmax ! mtil1-1 <= mtil2 <= mtil1+1
331 0 : lmtil2 = lltil1+mtil2
332 0 : mjj2=m-mtil2
333 0 : DO ljj2=0,lwn
334 0 : lljj2 = ljj2 * (ljj2 + 1)
335 0 : lmjj2 = lljj2 + mjj2
336 0 : lmini2 = abs(ltil1-ljj2)
337 0 : lmaxi2 = ltil1+ljj2
338 : ! Gaunt conditions (G2):
339 : if((lmini2.gt.l).or.(lmaxi2.lt.l).or.
340 0 : > (mod(l+ltil1+ljj2,2).ne.0).or.
341 : > (abs(mjj2).gt.ljj2)) cycle
342 :
343 : gc2 = gaunt1( l, ljj2, ltil1,
344 0 : > m, mjj2, mtil2, lmaxd)
345 :
346 : ! set up prefactor
347 : cil = ( ImagUnit ** (lp - l + ljj1 + ljj2) )
348 : + *( (-1.0) ** ljj1 )
349 : + * gc1 * gc2
350 : + * conjg( yl1(lmjj1 + 1) )
351 : + * conjg( yl2(lmjj2 + 1) )
352 0 : + * conjg( angso(lmtil1,lmtil2) )
353 :
354 : ! compute T-coefficients
355 : tuloulo(lo, m, lop, mp) = tuloulo(lo, m, lop, mp)
356 0 : > + cil * ulovulo(ltil1, ljj2, ljj1, lo, lop)
357 :
358 : ENDDO!ljj2
359 : ENDDO!m
360 : ENDDO!l
361 : ENDDO!mtil2
362 : ENDDO!ljj1
363 : ENDDO!mp
364 : ENDDO!lp
365 : ENDDO!mtil1
366 : ENDDO!ltil1
367 :
368 : ! apw-lo
369 0 : DO lop=1,nlo(ntyp)
370 0 : lp = llo(lop,ntyp)
371 0 : DO mp=-lp,lp
372 0 : DO ltil1=0,lwn
373 0 : lltil1 = ltil1*(ltil1+1)
374 0 : DO mtil1=-ltil1,ltil1
375 0 : lmtil1 = lltil1+mtil1
376 0 : mmin = max(-ltil1,mtil1-1)
377 0 : mmax = min( ltil1,mtil1+1)
378 0 : mjj1=mtil1-mp
379 0 : DO ljj1=0,lwn
380 0 : lljj1 = ljj1 * (ljj1 + 1)
381 0 : lmjj1 = lljj1 + mjj1
382 0 : lmini = abs(lp-ljj1)
383 0 : lmaxi = lp+ljj1
384 : ! Gaunt conditions (G1):
385 : if((lmini.gt.ltil1).or.(lmaxi.lt.ltil1).or.
386 0 : > (mod(lp+ltil1+ljj1,2).ne.0).or.
387 : > (abs(mjj1).gt.ljj1)) cycle
388 :
389 : gc1 = gaunt1(ltil1, ljj1, lp,
390 0 : > mtil1, mjj1, mp, lmaxd)
391 :
392 0 : DO mtil2=mmin,mmax ! mtil1-1 <= mtil2 <= mtil1+1
393 0 : lmtil2 = lltil1+mtil2
394 0 : DO l=0,lwn
395 0 : ll = l*(l+1)
396 0 : DO m=-l,l
397 0 : lm = ll+m
398 0 : mjj2=m-mtil2
399 0 : DO ljj2=0,lwn
400 0 : lljj2 = ljj2 * (ljj2 + 1)
401 0 : lmjj2 = lljj2 + mjj2
402 0 : lmini2 = abs(ltil1-ljj2)
403 0 : lmaxi2 = ltil1+ljj2
404 : ! Gaunt conditions (G2):
405 : if((lmini2.gt.l).or.(lmaxi2.lt.l).or.
406 0 : > (mod(l+ltil1+ljj2,2).ne.0).or.
407 : > (abs(mjj2).gt.ljj2)) cycle
408 :
409 : gc2 = gaunt1( l, ljj2, ltil1,
410 0 : > m, mjj2, mtil2, lmaxd)
411 :
412 : ! set up prefactor
413 : cil = ( ImagUnit ** (lp - l + ljj1 + ljj2) )
414 : + *( (-1.0) ** ljj1 )
415 : + * gc1 * gc2
416 : + * conjg( yl1(lmjj1 + 1) )
417 : + * conjg( yl2(lmjj2 + 1) )
418 0 : + * conjg( angso(lmtil1,lmtil2) )
419 :
420 : ! compute T-coefficients
421 : tulou(lm, lop, mp) = tulou(lm, lop, mp)
422 0 : > + cil * ulovu(ltil1, l, ljj2, ljj1, lop)
423 :
424 : tulod(lm, lop, mp) = tulod(lm, lop, mp)
425 0 : > + cil * ulovd(ltil1, l, ljj2, ljj1, lop)
426 :
427 : ENDDO!ljj2
428 : ENDDO!m
429 : ENDDO!l
430 : ENDDO!mtil2
431 : ENDDO!ljj1
432 : ENDDO!mp
433 : ENDDO!lp
434 : ENDDO!mtil1
435 : ENDDO!ltil1
436 :
437 :
438 0 : DO l=0,lwn
439 0 : ll = l*(l+1)
440 0 : DO m=-l,l
441 0 : lm = ll+m
442 0 : DO ltil1=0,lwn
443 0 : lltil1 = ltil1*(ltil1+1)
444 0 : DO mtil1=-ltil1,ltil1
445 0 : lmtil1 = lltil1+mtil1
446 0 : mmin = max(-ltil1,mtil1-1)
447 0 : mmax = min( ltil1,mtil1+1)
448 0 : mjj1=mtil1-m
449 0 : DO ljj1=0,lwn
450 0 : lljj1 = ljj1 * (ljj1 + 1)
451 0 : lmjj1 = lljj1 + mjj1
452 0 : lmini = abs(l-ljj1)
453 0 : lmaxi = l+ljj1
454 : ! Gaunt conditions (G1):
455 : if((lmini.gt.ltil1).or.(lmaxi.lt.ltil1).or.
456 0 : > (mod(l+ltil1+ljj1,2).ne.0).or.
457 : > (abs(mjj1).gt.ljj1)) cycle
458 :
459 : gc1 = gaunt1(ltil1, ljj1, l,
460 0 : > mtil1, mjj1, m, lmaxd)
461 :
462 0 : DO mtil2=mmin,mmax ! mtil1-1 <= mtil2 <= mtil1+1
463 0 : lmtil2 = lltil1+mtil2
464 0 : DO lop=1,nlo(ntyp)
465 0 : lp = llo(lop,ntyp)
466 0 : DO mp=-lp,lp
467 0 : mjj2=mp-mtil2
468 0 : DO ljj2=0,lwn
469 0 : lljj2 = ljj2 * (ljj2 + 1)
470 0 : lmjj2 = lljj2 + mjj2
471 0 : lmini2 = abs(ltil1-ljj2)
472 0 : lmaxi2 = ltil1+ljj2
473 : ! Gaunt conditions (G2):
474 : if((lmini2.gt.lp).or.(lmaxi2.lt.lp).or.
475 0 : > (mod(lp+ltil1+ljj2,2).ne.0).or.
476 : > (abs(mjj2).gt.ljj2)) cycle
477 :
478 : gc2 = gaunt1( lp, ljj2, ltil1,
479 0 : > mp, mjj2, mtil2, lmaxd)
480 :
481 : ! set up prefactor
482 : cil = ( ImagUnit ** (l - lp + ljj1 + ljj2) )
483 : + *( (-1.0) ** ljj1 )
484 : + * gc1 * gc2
485 : + * conjg( yl1(lmjj1 + 1) )
486 : + * conjg( yl2(lmjj2 + 1) )
487 0 : + * conjg( angso(lmtil1,lmtil2) )
488 :
489 : ! compute T-coefficients
490 : tuulo(lm, lop, mp) = tuulo(lm, lop, mp)
491 0 : > + cil * uvulo(ltil1, l, ljj2, ljj1, lop)
492 :
493 : tdulo(lm, lop, mp) = tdulo(lm, lop, mp)
494 0 : > + cil * dvulo(ltil1, l, ljj2, ljj1, lop)
495 :
496 : ENDDO!ljj2
497 : ENDDO!m
498 : ENDDO!l
499 : ENDDO!mtil2
500 : ENDDO!ljj1
501 : ENDDO!mp
502 : ENDDO!lp
503 : ENDDO!mtil1
504 : ENDDO!ltil1
505 :
506 : #if (defined(CPP_MPI) && !defined(CPP_T90))
507 0 : CALL MPI_BARRIER(MPI_COMM_WORLD,ierr(1))
508 : #endif
509 :
510 0 : deallocate( ulovulo )
511 0 : deallocate( uvulo, dvulo, ulovu, ulovd )
512 0 : deallocate( plo,plo_b )
513 0 : deallocate( x )
514 0 : deallocate( v0, vso )
515 :
516 0 : END SUBROUTINE wann_uHu_soc_tlo
517 : END MODULE m_wann_uHu_soc_tlo
|
