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