Line data Source code
1 : c***************************************c
2 : c Vacuum contribution to uHu matrix c
3 : c in FLAPW film mode c
4 : c***************************************c
5 : c < u_{k+b1} | H_{k} | u_{k+b2} > c
6 : c***************************************c
7 : c J.-P. Hanke, Dec. 2015 c
8 : c***************************************c
9 : MODULE m_wann_uHu_vac
10 : USE m_juDFT
11 :
12 : CONTAINS
13 0 : SUBROUTINE wann_uHu_vac(
14 : > chi,l_noco,l_soc,jspins,nlotot,qss,nbnd,z1,nmzxyd,
15 0 : > nmzd,n2d,nv2d,k1d,k2d,k3d,n3d,nvac,ig,rgphs,nmzxy,
16 0 : > nmz,delz,ig2,nq2,kv2,area,bmat,bbmat,evac,evac_b,
17 0 : > bkpt,bkpt_b,vzxy,vz,nslibd,nslibd_b,jspin,jspin_b,
18 0 : > ico,k1,k2,k3,k1_b,k2_b,k3_b,jspd,nvd,nbasfcn,neigd,
19 0 : > zMat,zMat_b,nv,nv_b,omtil,gb,gb2,sign2,uHu)
20 :
21 : USE m_types
22 : use m_constants, only : pimach
23 : use m_intgr, only : intgz0
24 : use m_d2fdz2cmplx
25 : USE m_vacuz
26 : USE m_vacudz
27 :
28 : implicit none
29 :
30 : TYPE(t_mat), INTENT(IN) :: zMat, zMat_b
31 :
32 : c .. scalar Arguments..
33 : logical, intent (in) :: l_noco,l_soc
34 : integer, intent (in) :: nlotot,jspin_b,n2d,jspins,ico
35 : integer, intent (in) :: nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n3d,nbnd
36 : integer, intent (in) :: nmzxy,nmz,nslibd,nslibd_b,nvac
37 : integer, intent (in) :: jspin,jspd,nvd,nq2
38 : integer, intent (in) :: nbasfcn,neigd
39 : real, intent (in) :: delz,z1,omtil,area
40 : complex, intent (in) :: chi
41 :
42 : c ..array arguments..
43 : real, intent (in) :: bkpt(3),bkpt_b(3),evac(2),evac_b(2)
44 : real, intent (in) :: qss(3)
45 : real, intent (in) :: vz(nmzd,2,4)
46 : real, intent (in) :: bbmat(3,3),bmat(3,3)
47 : ! complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
48 : complex, intent (in) :: rgphs(-k1d:,-k2d:,-k3d:)
49 : ! integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
50 : integer, intent (in) :: ig(-k1d:,-k2d:,-k3d:)
51 : integer, intent (in) :: ig2(n3d),nv(jspd),nv_b(jspd)
52 : integer, intent (in) :: kv2(2,n2d)
53 : integer, intent (in) :: gb(3),gb2(3)
54 : ! integer, intent (in) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
55 : integer, intent (in) :: k1(:,:),k2(:,:),k3(:,:)
56 : ! integer, intent (in) :: k1_b(nvd,jspd),k2_b(nvd,jspd)
57 : ! integer, intent (in) :: k3_b(nvd,jspd)
58 : integer, intent (in) :: k1_b(:,:),k2_b(:,:)
59 : integer, intent (in) :: k3_b(:,:)
60 : complex, intent (in) :: vzxy(nmzxyd,n2d-1,2)
61 : complex, intent (inout) :: uHu(nbnd,nbnd)
62 :
63 : c ..local arrays..
64 0 : complex, allocatable :: ac(:,:),bc(:,:),ac_b(:,:),bc_b(:,:)
65 0 : complex, allocatable :: vv(:),fac(:),fac_b(:),fac1(:)
66 0 : complex, allocatable :: du(:),due(:),du_b(:),due_b(:)
67 0 : complex, allocatable :: vxy(:,:),vxy_help(:)
68 0 : real, allocatable :: dt(:),dte(:)
69 0 : real, allocatable :: t(:),te(:),tei(:)
70 0 : real, allocatable :: u(:,:),ue(:,:),v(:)
71 0 : real, allocatable :: dt_b(:),dte_b(:)
72 0 : real, allocatable :: t_b(:),te_b(:),tei_b(:)
73 0 : real, allocatable :: u_b(:,:),ue_b(:,:)
74 0 : real, allocatable :: zmsh(:)
75 :
76 : c ..local scalars..
77 : logical tail
78 : real wronk,arg,zks,zks_b,tpi,vz0(2),vz0_b(2),scale,evacp,ev,const
79 0 : real xx(nmz),xximag(nmz),phase,phase2,rk,rk_b
80 : real :: qss1,qss2,xv,yv,s(3)
81 : integer i,m,l,j,k,n,nv2,nv2_b,ivac,n2,n2_b,sign,ik,sign2
82 : integer :: lp,np0,np1,addnoco,addnoco2,jspin2,jspin2_b,jspin2H
83 : complex :: av,bv,ic,c_1
84 : c complex :: tuu,tud,tdu,tdd
85 0 : complex, allocatable :: tuu(:,:),tud(:,:),tdu(:,:),tdd(:,:)
86 0 : complex, allocatable :: mat(:,:)
87 0 : integer, allocatable :: kvac1(:),kvac2(:),map2(:)
88 0 : integer, allocatable :: kvac1_b(:),kvac2_b(:),map2_b(:)
89 : integer symvac,symvacvac,igvm2,igvm2i
90 : integer :: j1,j2,i1,i2,ind3,ind2
91 : character(len=20) :: fname
92 :
93 : allocate ( ac(nv2d,nslibd),bc(nv2d,nslibd),
94 : + ac_b(nv2d,nslibd_b),bc_b(nv2d,nslibd_b),
95 : + dt(nv2d),dte(nv2d),t(nv2d),te(nv2d),
96 : + tei(nv2d),u(nmzd,nv2d),ue(nmzd,nv2d),
97 : + dt_b(nv2d),dte_b(nv2d),t_b(nv2d),te_b(nv2d),
98 : + tei_b(nv2d),u_b(nmzd,nv2d),ue_b(nmzd,nv2d),
99 : + v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd),
100 : + kvac1_b(nv2d),kvac2_b(nv2d),map2_b(nvd),
101 : + vv(nmzd),zmsh(nmzd),vxy(nmzxyd,n2d-1),vxy_help(n2d-1),
102 : + fac1(nmzd), fac(nmzd), fac_b(nmzd),
103 : + du(nmzd), due(nmzd), du_b(nmzd), due_b(nmzd),
104 : + tuu(nv2d,nv2d),tud(nv2d,nv2d),tdu(nv2d,nv2d),
105 0 : + tdd(nv2d,nv2d),mat(nslibd_b,nv2d) )
106 :
107 0 : tpi = 2 * pimach() ; ic = cmplx(0.,1.)
108 0 : tail = .true.
109 0 : np0 = nmzxy + 1
110 0 : np1 = nmz + 1
111 :
112 0 : jspin2 = jspin
113 0 : jspin2_b = jspin_b
114 0 : jspin2H = ico
115 0 : if(l_soc.and.jspins.eq.1) then
116 0 : jspin2 = 1
117 0 : jspin2_b = 1
118 0 : jspin2H = 1
119 : endif
120 :
121 : c.. determining the indexing array (in-plane stars)
122 : c.. for the k-point
123 :
124 0 : wronk = 2.0
125 0 : const = 1.0 / ( sqrt(omtil)*wronk )
126 :
127 0 : do ivac = 1,2
128 0 : vz0(ivac) = vz(nmz,ivac,jspin2)
129 0 : vz0_b(ivac) = vz(nmz,ivac,jspin2_b)
130 : c write(*,*)'vz0',ivac,vz0(ivac),vz0_b(ivac)
131 : enddo
132 :
133 0 : n2 = 0 ; n2_b = 0
134 :
135 0 : addnoco=0
136 0 : addnoco2=0
137 0 : if(l_noco.and.jspin.eq.2)then
138 0 : addnoco=nv(1)+nlotot
139 : endif
140 0 : if(l_noco.and.jspin_b.eq.2)then
141 0 : addnoco2=nv_b(1)+nlotot
142 : endif
143 :
144 0 : do 40 k = 1,nv(jspin)
145 0 : do 30 j = 1,n2
146 0 : if ( k1(k,jspin).eq.kvac1(j) .and.
147 : + k2(k,jspin).eq.kvac2(j) ) then
148 0 : map2(k) = j
149 0 : goto 40
150 : endif
151 0 : 30 continue
152 0 : n2 = n2 + 1
153 :
154 0 : IF(n2>nv2d) then
155 0 : write(*,*)n2,nv2d,'jspin',jspin
156 : endif
157 :
158 0 : IF (n2>nv2d) CALL juDFT_error("wannier uHu vac",calledby
159 0 : + ="wann_uHu_vac")
160 :
161 0 : kvac1(n2) = k1(k,jspin)
162 0 : kvac2(n2) = k2(k,jspin)
163 0 : map2(k) = n2
164 0 : 40 continue
165 :
166 :
167 : c.. and for the b-point
168 :
169 0 : do 41 k = 1,nv_b(jspin_b)
170 0 : do 31 j = 1,n2_b
171 0 : if ( k1_b(k,jspin_b).eq.kvac1_b(j) .and.
172 : + k2_b(k,jspin_b).eq.kvac2_b(j) ) then
173 0 : map2_b(k) = j
174 0 : goto 41
175 : endif
176 0 : 31 continue
177 0 : n2_b = n2_b + 1
178 :
179 0 : IF(n2_b>nv2d) then
180 0 : write(*,*)n2_b,nv2d,'jspin_b',jspin_b
181 : endif
182 :
183 0 : IF (n2_b>nv2d) CALL juDFT_error("wannier uHu vac",calledby
184 0 : + ="wann_uHu_vac")
185 :
186 0 : kvac1_b(n2_b) = k1_b(k,jspin_b)
187 0 : kvac2_b(n2_b) = k2_b(k,jspin_b)
188 0 : map2_b(k) = n2_b
189 0 : 41 continue
190 :
191 :
192 :
193 :
194 : c...cycle by the vacua
195 0 : do 140 ivac = 1,2!nvac
196 :
197 0 : sign = 3. - 2.*ivac
198 0 : evacp = evac(ivac)
199 :
200 0 : vxy(:,:) = vzxy(:,:,ivac)
201 : IF(l_noco) THEN ! symmetrization
202 : !IF (nvac.EQ.1 .and. ivac.EQ.2 .AND.(.NOT.zrfs) ) THEN
203 : ! DO i=1,nmzxy
204 : ! DO igvm2 = 2,nq2
205 : ! igvm2i = ig2(ig(-kv2(1,igvm2),-kv2(2,igvm2),0))
206 : ! vxy_help(igvm2-1) = vxy(i,igvm2i-1)
207 : ! ENDDO
208 : ! DO igvm2 = 2,nq2
209 : ! vxy(i,igvm2-1) = vxy_help(igvm2-1)
210 : ! ENDDO
211 : ! ENDDO
212 : !ENDIF
213 : ENDIF
214 :
215 0 : tuu = cmplx(0.,0.)
216 0 : tud = cmplx(0.,0.)
217 0 : tdu = cmplx(0.,0.)
218 0 : tdd = cmplx(0.,0.)
219 :
220 0 : nv2 = n2 ; nv2_b = n2_b
221 :
222 : c.. the body of the routine
223 :
224 0 : qss1=0.0
225 0 : qss2=0.0
226 0 : if(l_noco.and.jspin.eq.1)then
227 0 : qss1=-qss(1)/2.0
228 0 : qss2=-qss(2)/2.0
229 0 : elseif(l_noco.and.jspin.eq.2)then
230 0 : qss1=qss(1)/2.0
231 0 : qss2=qss(2)/2.0
232 : endif
233 :
234 0 : do ik = 1,nv2
235 0 : v(1) = bkpt(1) + kvac1(ik) + qss1
236 0 : v(2) = bkpt(2) + kvac2(ik) + qss2
237 0 : v(3) = 0.
238 0 : ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
239 : ! ev = evacp - 0.5*dotirp(v,v,bbmat)
240 : call vacuz(ev,vz(1,ivac,jspin2),vz0(ivac),nmz,delz,t(ik),
241 0 : + dt(ik),u(1,ik))
242 : call vacudz(ev,vz(1,ivac,jspin2),vz0(ivac),nmz,delz,te(ik),
243 : + dte(ik),tei(ik),ue(1,ik),dt(ik),
244 0 : + u(1,ik))
245 0 : scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
246 0 : te(ik) = scale*te(ik)
247 0 : dte(ik) = scale*dte(ik)
248 0 : tei(ik) = scale*tei(ik)
249 0 : do j = 1,nmz
250 0 : ue(j,ik) = scale*ue(j,ik)
251 : enddo
252 : enddo
253 : c-----> construct a and b coefficients for the k-point
254 : symvacvac=1
255 : c if (nvac==1) symvacvac=2
256 0 : do symvac=1,symvacvac
257 0 : do i = 1,nv2d
258 0 : do n = 1,nslibd
259 0 : ac(i,n) = cmplx(0.0,0.0)
260 0 : bc(i,n) = cmplx(0.0,0.0)
261 : enddo
262 0 : do n = 1,nslibd_b
263 0 : ac_b(i,n) = cmplx(0.0,0.0)
264 0 : bc_b(i,n) = cmplx(0.0,0.0)
265 : enddo
266 : enddo
267 :
268 : if (symvac==2) sign=-1.0
269 :
270 0 : do k = 1,nv(jspin)
271 0 : l = map2(k)
272 0 : zks = k3(k,jspin)*bmat(3,3)*sign
273 0 : arg = zks*z1
274 0 : c_1 = cmplx(cos(arg),sin(arg)) * const
275 0 : av = -c_1 * cmplx( dte(l),zks*te(l) )
276 0 : bv = c_1 * cmplx( dt(l),zks* t(l) )
277 : c-----> loop over basis functions
278 :
279 0 : IF(zMat%l_real) THEN
280 0 : do n = 1,nslibd
281 0 : ac(l,n) = ac(l,n) + zMat%data_r(k+addnoco,n)*av
282 0 : bc(l,n) = bc(l,n) + zMat%data_r(k+addnoco,n)*bv
283 : enddo
284 : ELSE
285 0 : do n = 1,nslibd
286 0 : ac(l,n) = ac(l,n) + zMat%data_c(k+addnoco,n)*av
287 0 : bc(l,n) = bc(l,n) + zMat%data_c(k+addnoco,n)*bv
288 : enddo
289 : END IF
290 : enddo
291 :
292 :
293 : c...now for the k+b point
294 :
295 0 : evacp = evac_b(ivac)
296 0 : do ik = 1,nv2_b
297 0 : v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
298 0 : v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
299 0 : v(3) = 0.
300 0 : ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
301 : ! ev = evacp - 0.5*dotirp(v,v,bbmat)
302 : call vacuz(ev,vz(1,ivac,jspin2_b),vz0_b(ivac),nmz,delz,
303 0 : + t_b(ik),dt_b(ik),u_b(1,ik))
304 : call vacudz(ev,vz(1,ivac,jspin2_b),vz0_b(ivac),nmz,delz,
305 : + te_b(ik),dte_b(ik),tei_b(ik),ue_b(1,ik),dt_b(ik),
306 0 : + u_b(1,ik))
307 0 : scale = wronk/ (te_b(ik)*dt_b(ik)-dte_b(ik)*t_b(ik))
308 0 : te_b(ik) = scale*te_b(ik)
309 0 : dte_b(ik) = scale*dte_b(ik)
310 0 : tei_b(ik) = scale*tei_b(ik)
311 0 : do j = 1,nmz
312 0 : ue_b(j,ik) = scale*ue_b(j,ik)
313 : enddo
314 : enddo
315 : c-----> construct a and b coefficients for the k+b point
316 :
317 0 : do k = 1,nv_b(jspin_b)
318 0 : l = map2_b(k)
319 0 : zks = k3_b(k,jspin_b)*bmat(3,3)*sign
320 0 : arg = zks*z1
321 0 : c_1 = cmplx(cos(arg),sin(arg)) * const
322 0 : av = -c_1 * cmplx( dte_b(l),zks*te_b(l) )
323 0 : bv = c_1 * cmplx( dt_b(l),zks*t_b(l) )
324 : c-----> loop over basis functions
325 :
326 0 : IF(zMat_b%l_real) THEN
327 0 : do n = 1,nslibd_b
328 0 : ac_b(l,n) = ac_b(l,n) + zMat_b%data_r(k+addnoco2,n)*av
329 0 : bc_b(l,n) = bc_b(l,n) + zMat_b%data_r(k+addnoco2,n)*bv
330 : enddo
331 : ELSE
332 0 : do n = 1,nslibd_b
333 0 : ac_b(l,n) = ac_b(l,n) + zMat_b%data_c(k+addnoco2,n)*av
334 0 : bc_b(l,n) = bc_b(l,n) + zMat_b%data_c(k+addnoco2,n)*bv
335 : enddo
336 : END IF
337 : enddo
338 :
339 :
340 : ! set up z-mesh and plane-wave factors
341 0 : zks = gb(3) *bmat(3,3) ! TODO: different sign before
342 0 : zks_b= -gb2(3)*bmat(3,3) ! TODO: different sign before
343 0 : do i=1,nmz
344 0 : zmsh(i)= ( z1 + (i-1)*delz )*sign
345 0 : fac(i) = cmplx( cos(zmsh(i)*zks), sin(zmsh(i)*zks) )
346 0 : fac_b(i)= cmplx( cos(zmsh(i)*zks_b),sin(zmsh(i)*zks_b))
347 0 : fac1(i) = fac(i) * fac_b(i)
348 : enddo
349 :
350 :
351 : ! calculate uHu matrix elements
352 0 : do l = 1,nv2
353 :
354 0 : j1 = kvac1(l) - gb(1)
355 0 : j2 = kvac2(l) - gb(2)
356 :
357 0 : do lp = 1,nv2_b
358 0 : i1 = j1 - kvac1_b(lp) + gb2(1)
359 0 : i2 = j2 - kvac2_b(lp) + gb2(2)
360 :
361 : c i1 = -i1
362 : c i2 = -i2
363 :
364 0 : ind3 = ig(sign2*i1,sign2*i2,0)
365 0 : IF (ind3.EQ.0) CYCLE
366 :
367 0 : phase = rgphs(i1,i2,0) ! TODO: sign also here?
368 0 : phase2= rgphs(sign2*i1,sign2*i2,0)
369 0 : if(phase.ne.phase2) then
370 0 : stop 'rgphs in wann_uHu_vac'
371 : endif
372 :
373 0 : ind2 = ig2(ind3)
374 0 : IF (ind2.EQ.0) CALL juDFT_error("error in map2 for 2-d stars",
375 0 : > calledby="wann_uHu_vac")
376 0 : ind2 = ind2 - 1
377 0 : IF (ind2.NE.0) THEN ! warping components G=/=0
378 : ! need integral with vxy and plane-waves
379 : ! tuu
380 0 : DO i=1,nmzxy
381 : xx(np0-i) = u(i,l)*u_b(i,lp)
382 0 : > * real( fac1(i)*vxy(i,ind2) )
383 : xximag(np0-i) = u(i,l)*u_b(i,lp)
384 0 : > * aimag( fac1(i)*vxy(i,ind2) )
385 : ENDDO
386 0 : CALL intgz0(xx,delz,nmzxy,xv,tail)
387 0 : CALL intgz0(xximag,delz,nmzxy,yv,tail)
388 0 : tuu(lp,l) = phase*cmplx(xv,yv)
389 :
390 : ! tud
391 0 : DO i=1,nmzxy
392 : xx(np0-i) = u(i,l)*ue_b(i,lp)
393 0 : > * real( fac1(i)*vxy(i,ind2) )
394 : xximag(np0-i) = u(i,l)*ue_b(i,lp)
395 0 : > * aimag( fac1(i)*vxy(i,ind2) )
396 : ENDDO
397 0 : CALL intgz0(xx,delz,nmzxy,xv,tail)
398 0 : CALL intgz0(xximag,delz,nmzxy,yv,tail)
399 0 : tud(lp,l) = phase*cmplx(xv,yv)
400 :
401 : ! tdu
402 0 : DO i=1,nmzxy
403 : xx(np0-i) = ue(i,l)*u_b(i,lp)
404 0 : > * real( fac1(i)*vxy(i,ind2) )
405 : xximag(np0-i) = ue(i,l)*u_b(i,lp)
406 0 : > * aimag( fac1(i)*vxy(i,ind2) )
407 : ENDDO
408 0 : CALL intgz0(xx,delz,nmzxy,xv,tail)
409 0 : CALL intgz0(xximag,delz,nmzxy,yv,tail)
410 0 : tdu(lp,l) = phase*cmplx(xv,yv)
411 :
412 : ! tdd
413 0 : DO i=1,nmzxy
414 : xx(np0-i) = ue(i,l)*ue_b(i,lp)
415 0 : > * real( fac1(i)*vxy(i,ind2) )
416 : xximag(np0-i) = ue(i,l)*ue_b(i,lp)
417 0 : > * aimag( fac1(i)*vxy(i,ind2) )
418 : ENDDO
419 0 : CALL intgz0(xx,delz,nmzxy,xv,tail)
420 0 : CALL intgz0(xximag,delz,nmzxy,yv,tail)
421 0 : tdd(lp,l) = phase*cmplx(xv,yv)
422 :
423 : ELSE ! non-warping components G==0
424 : ! need integral with H(z) = -1/2 d2/dz2 + vz + (G+k)^2
425 :
426 0 : IF ( (ico.EQ.1) .OR. (ico.EQ.2) ) THEN ! spin-diagonal
427 :
428 : ! determine second derivative of (u*fac) etc.
429 : CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
430 0 : > u(:,l),fac,du)
431 : CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
432 0 : > ue(:,l),fac,due)
433 : CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
434 0 : > u_b(:,lp),fac_b,du_b)
435 : CALL d2fdz2cmplx(nmzd,nmz,zmsh,delz,
436 0 : > ue_b(:,lp),fac_b,due_b)
437 :
438 : ! determine |G+k+b1|^2 and |G'+k+b2|^2
439 0 : s(1) = bkpt(1) + kvac1(l)
440 0 : s(2) = bkpt(2) + kvac2(l)
441 0 : s(3) = 0.0
442 0 : rk = dot_product(s,matmul(bbmat,s))
443 : ! rk = dotirp(s,s,bbmat)
444 :
445 0 : s(1) = bkpt_b(1) + kvac1_b(lp)
446 0 : s(2) = bkpt_b(2) + kvac2_b(lp)
447 : s(3) = 0.0
448 0 : rk_b = dot_product(s,matmul(bbmat,s))
449 : ! rk_b = dotirp(s,s,bbmat)
450 :
451 : c ! no kinetic energy if jspin =/= jspin_b
452 : c if(jspin.ne.jspin_b) then
453 : c rk = 0.0
454 : c rk_b = 0.0
455 : c du=cmplx(0.0,0.0)
456 : c due=cmplx(0.0,0.0)
457 : c du_b=cmplx(0.0,0.0)
458 : c due_b=cmplx(0.0,0.0)
459 : c endif
460 :
461 : c rk_b=0.0
462 : c du_b=cmplx(0.,0.)
463 : c due_b=cmplx(0.,0.)
464 :
465 : ! tuu
466 0 : DO i=1,nmz
467 : xx(np1-i) = u(i,l)*u_b(i,lp)
468 : > * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
469 : > - 0.25*real( fac(i)*u(i,l)*du_b(i)
470 0 : > + du(i)*u_b(i,lp)*fac_b(i) )
471 :
472 : xximag(np1-i) = u(i,l)*u_b(i,lp)
473 : > * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
474 : > - 0.25*aimag( fac(i)*u(i,l)*du_b(i)
475 0 : > + du(i)*u_b(i,lp)*fac_b(i) )
476 : ENDDO
477 0 : CALL intgz0(xx,delz,nmz,xv,tail)
478 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
479 0 : tuu(lp,l) = cmplx(xv,yv)
480 :
481 : ! tud
482 0 : DO i=1,nmz
483 : xx(np1-i) = u(i,l)*ue_b(i,lp)
484 : > * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
485 : > - 0.25*real( fac(i)*u(i,l)*due_b(i)
486 0 : > + du(i)*ue_b(i,lp)*fac_b(i) )
487 :
488 : xximag(np1-i) = u(i,l)*ue_b(i,lp)
489 : > * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
490 : > - 0.25*aimag( fac(i)*u(i,l)*due_b(i)
491 0 : > + du(i)*ue_b(i,lp)*fac_b(i) )
492 : ENDDO
493 0 : CALL intgz0(xx,delz,nmz,xv,tail)
494 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
495 0 : tud(lp,l) = cmplx(xv,yv)
496 :
497 : ! tdu
498 0 : DO i=1,nmz
499 : xx(np1-i) = ue(i,l)*u_b(i,lp)
500 : > * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
501 : > - 0.25*real( fac(i)*ue(i,l)*du_b(i)
502 0 : > + due(i)*u_b(i,lp)*fac_b(i) )
503 :
504 : xximag(np1-i) = ue(i,l)*u_b(i,lp)
505 : > * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
506 : > - 0.25*aimag( fac(i)*ue(i,l)*du_b(i)
507 0 : > + due(i)*u_b(i,lp)*fac_b(i) )
508 : ENDDO
509 0 : CALL intgz0(xx,delz,nmz,xv,tail)
510 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
511 0 : tdu(lp,l) = cmplx(xv,yv)
512 :
513 : ! tdd
514 0 : DO i=1,nmz
515 : xx(np1-i) = ue(i,l)*ue_b(i,lp)
516 : > * real(fac1(i)) * ( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
517 : > - 0.25*real( fac(i)*ue(i,l)*due_b(i)
518 0 : > + due(i)*ue_b(i,lp)*fac_b(i) )
519 :
520 : xximag(np1-i) = ue(i,l)*ue_b(i,lp)
521 : > * aimag(fac1(i)) *( vz(i,ivac,jspin2H)+0.25*(rk+rk_b) )
522 : > - 0.25*aimag( fac(i)*ue(i,l)*due_b(i)
523 0 : > + due(i)*ue_b(i,lp)*fac_b(i) )
524 : ENDDO
525 0 : CALL intgz0(xx,delz,nmz,xv,tail)
526 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
527 0 : tdd(lp,l) = cmplx(xv,yv)
528 :
529 : ELSE ! spin-off-diagonal
530 :
531 0 : DO i=1,nmz
532 0 : vv(i) = cmplx( vz(i,ivac,3), vz(i,ivac,4) )
533 : ENDDO
534 0 : IF ( ico.EQ.4 ) vv = conjg(vv)
535 :
536 : ! tuu
537 0 : DO i=1,nmz
538 : xx(np1-i) = u(i,l)*u_b(i,lp)
539 0 : > * real( fac1(i)*vv(i) )
540 : xximag(np1-i) = u(i,l)*u_b(i,lp)
541 0 : > * aimag( fac1(i)*vv(i) )
542 : ENDDO
543 0 : CALL intgz0(xx,delz,nmz,xv,tail)
544 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
545 0 : tuu(lp,l) = cmplx(xv,yv)
546 :
547 : ! tud
548 0 : DO i=1,nmz
549 : xx(np1-i) = u(i,l)*ue_b(i,lp)
550 0 : > * real( fac1(i)*vv(i) )
551 : xximag(np1-i) = u(i,l)*ue_b(i,lp)
552 0 : > * aimag( fac1(i)*vv(i) )
553 : ENDDO
554 0 : CALL intgz0(xx,delz,nmz,xv,tail)
555 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
556 0 : tud(lp,l) = cmplx(xv,yv)
557 :
558 : ! tdu
559 0 : DO i=1,nmz
560 : xx(np1-i) = ue(i,l)*u_b(i,lp)
561 0 : > * real( fac1(i)*vv(i) )
562 : xximag(np1-i) = ue(i,l)*u_b(i,lp)
563 0 : > * aimag( fac1(i)*vv(i) )
564 : ENDDO
565 0 : CALL intgz0(xx,delz,nmz,xv,tail)
566 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
567 0 : tdu(lp,l) = cmplx(xv,yv)
568 :
569 : ! tdd
570 0 : DO i=1,nmz
571 : xx(np1-i) = ue(i,l)*ue_b(i,lp)
572 0 : > * real( fac1(i)*vv(i) )
573 : xximag(np1-i) = ue(i,l)*ue_b(i,lp)
574 0 : > * aimag( fac1(i)*vv(i) )
575 : ENDDO
576 0 : CALL intgz0(xx,delz,nmz,xv,tail)
577 0 : CALL intgz0(xximag,delz,nmz,yv,tail)
578 0 : tdd(lp,l) = cmplx(xv,yv)
579 :
580 : ENDIF !(ico.eq.1).or.(ico.eq.2)
581 :
582 : ENDIF !(ind2.ne.0)
583 :
584 :
585 : ! contraction of integrals with a,b coefficients
586 : ! yields contribution to uHu matrix
587 : c do i = 1,nslibd
588 : c do j = 1,nslibd_b
589 : c uHu(i,j) = uHu(i,j) + area*
590 : c * ( ac(l,i)*conjg(ac_b(lp,j))*tuu +
591 : c + ac(l,i)*conjg(bc_b(lp,j))*tud +
592 : c * bc(l,i)*conjg(ac_b(lp,j))*tdu +
593 : c + bc(l,i)*conjg(bc_b(lp,j))*tdd )
594 : c enddo
595 : c enddo
596 :
597 :
598 : enddo ! lp
599 : enddo ! l
600 :
601 : if(.false.) then
602 : call zgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
603 : > ac_b(1,1),nv2d,tuu(1,1),nv2d,cmplx(0.0),
604 : > mat(1,1),nslibd_b)
605 : call zgemm('T','T',nslibd,nslibd_b,nv2d,chi,
606 : > ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
607 : > uHu,nbnd)
608 :
609 : call zgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
610 : > bc_b(1,1),nv2d,tud(1,1),nv2d,cmplx(0.0),
611 : > mat(1,1),nslibd_b)
612 : call zgemm('T','T',nslibd,nslibd_b,nv2d,chi,
613 : > ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
614 : > uHu,nbnd)
615 :
616 : call zgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
617 : > ac_b(1,1),nv2d,tdu(1,1),nv2d,cmplx(0.0),
618 : > mat(1,1),nslibd_b)
619 : call zgemm('T','T',nslibd,nslibd_b,nv2d,chi,
620 : > bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
621 : > uHu,nbnd)
622 :
623 : call zgemm('C','N',nslibd_b,nv2d,nv2d,cmplx(area),
624 : > bc_b(1,1),nv2d,tdd(1,1),nv2d,cmplx(0.0),
625 : > mat(1,1),nslibd_b)
626 : call zgemm('T','T',nslibd,nslibd_b,nv2d,chi,
627 : > bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
628 : > uHu,nbnd)
629 : endif
630 :
631 0 : if(.true.)then ! standard so far
632 : call zgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
633 : > ac_b(1,1),nv2d,tuu(1,1),nv2d,cmplx(0.0),
634 0 : > mat(1,1),nslibd_b)
635 : call zgemm('T','C',nslibd,nslibd_b,nv2d,chi,
636 : > ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
637 0 : > uHu,nbnd)
638 :
639 : call zgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
640 : > bc_b(1,1),nv2d,tud(1,1),nv2d,cmplx(0.0),
641 0 : > mat(1,1),nslibd_b)
642 : call zgemm('T','C',nslibd,nslibd_b,nv2d,chi,
643 : > ac(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
644 0 : > uHu,nbnd)
645 :
646 : call zgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
647 : > ac_b(1,1),nv2d,tdu(1,1),nv2d,cmplx(0.0),
648 0 : > mat(1,1),nslibd_b)
649 : call zgemm('T','C',nslibd,nslibd_b,nv2d,chi,
650 : > bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
651 0 : > uHu,nbnd)
652 :
653 : call zgemm('T','N',nslibd_b,nv2d,nv2d,cmplx(area),
654 : > bc_b(1,1),nv2d,tdd(1,1),nv2d,cmplx(0.0),
655 0 : > mat(1,1),nslibd_b)
656 : call zgemm('T','C',nslibd,nslibd_b,nv2d,chi,
657 : > bc(1,1),nv2d,mat(1,1),nslibd_b,cmplx(1.0),
658 0 : > uHu,nbnd)
659 : endif
660 :
661 : enddo !symvac
662 : c... cycle by the vacua finishes
663 0 : 140 enddo
664 :
665 0 : deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
666 0 : + v,kvac1,kvac2,map2 )
667 0 : deallocate ( ac_b,bc_b,dt_b,dte_b,t_b,te_b,tei_b,u_b,ue_b,
668 0 : + kvac1_b,kvac2_b,map2_b )
669 0 : deallocate ( vxy, vxy_help, vv, zmsh, fac1, fac, fac_b )
670 0 : deallocate ( du,due,du_b,due_b )
671 0 : deallocate ( tuu,tud,tdu,tdd,mat )
672 :
673 : c call fleur_end("stop")
674 :
675 0 : END SUBROUTINE wann_uHu_vac
676 : END MODULE m_wann_uHu_vac
|
