Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_wann_mmkb_vac
8 : use m_juDFT
9 : c**************************************************************
10 : c Determines the overlap matrix Mmn(k) in the vacuum
11 : c in the film case for the wannier functions.
12 : c For more details see routine wannier.F
13 : c
14 : c Y. Mokrousov, F. Freimuth
15 : c***************************************************************
16 : CONTAINS
17 0 : SUBROUTINE wann_mmkb_vac(
18 : > vacchi,l_noco,nlotot,qss,
19 : > nbnd,z1,nmzd,nv2d,k1d,k2d,k3d,n3d,nvac,
20 : > ig,nmz,delz,ig2,area,bmat,
21 0 : > bbmat,evac,evac_b,bkpt,bkpt_b,vz,vz_b,
22 : > nslibd,nslibd_b,jspin,jspin_b,
23 0 : > k1,k2,k3,k1_b,k2_b,k3_b,jspd,nvd,
24 0 : > nbasfcn,neigd,zMat,zMat_b,nv,nv_b,omtil,gb,
25 0 : < mmn)
26 :
27 : USE m_types
28 : use m_constants
29 : use m_intgr, only : intgz0
30 : USE m_vacuz
31 : USE m_vacudz
32 :
33 : implicit none
34 :
35 : TYPE(t_mat), INTENT(IN) :: zMat, zMat_b
36 :
37 : c .. scalar Arguments..
38 : logical, intent (in) :: l_noco
39 : integer, intent (in) :: nlotot,jspin_b
40 : real, intent (in) :: qss(3)
41 : integer, intent (in) :: nmzd,nv2d,k1d,k2d,k3d,n3d,nbnd
42 : integer, intent (in) :: nmz,nslibd,nslibd_b,nvac
43 : integer, intent (in) :: jspin,jspd,nvd
44 : integer, intent (in) :: nbasfcn,neigd
45 : integer, intent (in) :: gb(3)
46 : real, intent (in) :: delz,z1,omtil,area
47 : complex, intent (in) :: vacchi
48 :
49 : c ..array arguments..
50 : real, intent (in) :: bkpt(3),bkpt_b(3),evac(2),evac_b(2)
51 : ! integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
52 : integer, intent (in) :: ig(-k1d:,-k2d:,-k3d:)
53 : ! integer, intent (in) :: ig2(n3d)
54 : integer, intent (in) :: ig2(:)!ig2(n3d)
55 : ! integer, intent (in) :: nv(jspd),nv_b(jspd)
56 : integer, intent (in) :: nv(:),nv_b(:)
57 : real, intent (in) :: vz(nmzd,2),vz_b(nmzd,2)
58 : real, intent (in) :: bbmat(3,3),bmat(3,3)
59 : ! integer, intent (in) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
60 : integer, intent (in) :: k1(:,:),k2(:,:),k3(:,:)
61 : ! integer,intent(in)::k1_b(nvd,jspd),k2_b(nvd,jspd),k3_b(nvd,jspd)
62 : integer,intent(in)::k1_b(:,:),k2_b(:,:),k3_b(:,:)
63 : ! complex, intent (inout) :: mmn(nbnd,nbnd)
64 : complex, intent (inout) :: mmn(:,:)
65 :
66 : c ..basis wavefunctions in the vacuum
67 0 : complex, allocatable :: ac(:,:),bc(:,:),ac_b(:,:),bc_b(:,:)
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 :
75 : c ..local scalars..
76 : logical tail
77 : real wronk,arg,zks,tpi,vz0(2),vz0_b(2),scale,evacp,ev,const
78 0 : real uu,ud,du,dd,xx(nmz),xximag(nmz)
79 : real :: uuimag,udimag,duimag,ddimag,qss1,qss2
80 : integer i,m,l,j,k,n,nv2,nv2_b,ivac,n2,n2_b,sign,ik
81 : integer :: lprime,np1,addnoco,addnoco2
82 : complex :: av,bv,ic,c_1
83 0 : integer, allocatable :: kvac1(:),kvac2(:),map2(:)
84 0 : integer, allocatable :: kvac1_b(:),kvac2_b(:),map2_b(:)
85 : integer symvac,symvacvac
86 0 : call timestart("wann_mmkb_vac")
87 : allocate ( ac(nv2d,nslibd),bc(nv2d,nslibd),
88 : + ac_b(nv2d,nslibd_b),bc_b(nv2d,nslibd_b),
89 : + dt(nv2d),dte(nv2d),t(nv2d),te(nv2d),
90 : + tei(nv2d),u(nmzd,nv2d),ue(nmzd,nv2d),
91 : + dt_b(nv2d),dte_b(nv2d),t_b(nv2d),te_b(nv2d),
92 : + tei_b(nv2d),u_b(nmzd,nv2d),ue_b(nmzd,nv2d),
93 : + v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd),
94 0 : + kvac1_b(nv2d),kvac2_b(nv2d),map2_b(nvd) )
95 :
96 0 : tpi = 2 * pimach() ; ic = cmplx(0.,1.)
97 :
98 0 : tail = .true.
99 0 : np1 = nmz + 1
100 :
101 : c.. determining the indexing array (in-plane stars)
102 : c.. for the k-point
103 :
104 0 : wronk = 2.0
105 0 : const = 1.0 / ( sqrt(omtil)*wronk )
106 :
107 0 : do ivac = 1,2
108 0 : vz0(ivac) = vz(nmz,ivac)
109 0 : vz0_b(ivac) = vz_b(nmz,ivac)
110 : enddo
111 :
112 0 : n2 = 0 ; n2_b = 0
113 :
114 0 : addnoco=0
115 0 : addnoco2=0
116 0 : if(l_noco.and.jspin.eq.2)then
117 0 : addnoco=nv(1)+nlotot
118 : endif
119 : if(l_noco.and.jspin_b.eq.2)then
120 : addnoco2=nv_b(1)+nlotot
121 : endif
122 :
123 0 : do 40 k = 1,nv(jspin)
124 0 : do 30 j = 1,n2
125 0 : if ( k1(k,jspin).eq.kvac1(j) .and.
126 : + k2(k,jspin).eq.kvac2(j) ) then
127 0 : map2(k) = j
128 0 : goto 40
129 : endif
130 0 : 30 continue
131 0 : n2 = n2 + 1
132 :
133 0 : IF(n2>nv2d) then
134 0 : write(*,*)n2,nv2d,'jspin',jspin
135 : endif
136 :
137 0 : IF (n2>nv2d) CALL juDFT_error("wannier Mmn vac",calledby
138 0 : + ="wann_mmkb_vac")
139 :
140 0 : kvac1(n2) = k1(k,jspin)
141 0 : kvac2(n2) = k2(k,jspin)
142 0 : map2(k) = n2
143 0 : 40 continue
144 : !write(*,*)'ok',n2,nv2d,'jspin',jspin
145 :
146 : c.. and for the b-point
147 :
148 0 : do 41 k = 1,nv_b(jspin_b)
149 0 : do 31 j = 1,n2_b
150 0 : if ( k1_b(k,jspin_b).eq.kvac1_b(j) .and.
151 : + k2_b(k,jspin_b).eq.kvac2_b(j) ) then
152 0 : map2_b(k) = j
153 0 : goto 41
154 : endif
155 0 : 31 continue
156 0 : n2_b = n2_b + 1
157 :
158 0 : IF(n2_b>nv2d) then
159 0 : write(*,*)n2_b,nv2d,'jspin_b',jspin_b
160 : endif
161 :
162 0 : IF (n2_b>nv2d) CALL juDFT_error("wannier Mmn vac",calledby
163 0 : + ="wann_mmkb_vac")
164 :
165 0 : kvac1_b(n2_b) = k1_b(k,jspin_b)
166 0 : kvac2_b(n2_b) = k2_b(k,jspin_b)
167 0 : map2_b(k) = n2_b
168 0 : 41 continue
169 : !write(*,*)'ok',n2_b,nv2d,'jspin_b',jspin_b
170 :
171 : c...cycle by the vacua
172 0 : do 140 ivac = 1,nvac
173 :
174 :
175 :
176 0 : sign = 3. - 2.*ivac
177 0 : evacp = evac(ivac)
178 :
179 0 : nv2 = n2 ; nv2_b = n2_b
180 :
181 : c.. the body of the routine
182 :
183 0 : qss1=0.0
184 0 : qss2=0.0
185 0 : if(l_noco.and.jspin.eq.1)then
186 0 : qss1=-qss(1)/2.0
187 0 : qss2=-qss(2)/2.0
188 0 : elseif(l_noco.and.jspin.eq.2)then
189 0 : qss1=qss(1)/2.0
190 0 : qss2=qss(2)/2.0
191 : endif
192 :
193 0 : do ik = 1,nv2
194 0 : v(1) = bkpt(1) + kvac1(ik) + qss1
195 0 : v(2) = bkpt(2) + kvac2(ik) + qss2
196 0 : v(3) = 0.
197 0 : ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
198 : call vacuz(ev,vz(1,ivac),vz0(ivac),nmz,delz,t(ik),dt(ik),
199 0 : + u(1,ik))
200 : call vacudz(ev,vz(1,ivac),vz0(ivac),nmz,delz,te(ik),
201 : + dte(ik),tei(ik),ue(1,ik),dt(ik),
202 0 : + u(1,ik))
203 0 : scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
204 0 : te(ik) = scale*te(ik)
205 0 : dte(ik) = scale*dte(ik)
206 0 : tei(ik) = scale*tei(ik)
207 0 : do j = 1,nmz
208 0 : ue(j,ik) = scale*ue(j,ik)
209 : enddo
210 : enddo
211 : c-----> construct a and b coefficients for the k-point
212 0 : symvacvac=1
213 0 : if (nvac==1) symvacvac=2
214 0 : do symvac=1,symvacvac
215 0 : do i = 1,nv2d
216 0 : do n = 1,nslibd
217 0 : ac(i,n) = cmplx(0.0,0.0)
218 0 : bc(i,n) = cmplx(0.0,0.0)
219 : enddo
220 0 : do n = 1,nslibd_b
221 0 : ac_b(i,n) = cmplx(0.0,0.0)
222 0 : bc_b(i,n) = cmplx(0.0,0.0)
223 : enddo
224 : enddo
225 :
226 0 : if (symvac==2) sign=-1.0
227 :
228 0 : do k = 1,nv(jspin)
229 0 : l = map2(k)
230 0 : zks = k3(k,jspin)*bmat(3,3)*sign
231 0 : arg = zks*z1
232 0 : c_1 = cmplx(cos(arg),sin(arg)) * const
233 0 : av = -c_1 * cmplx( dte(l),zks*te(l) )
234 0 : bv = c_1 * cmplx( dt(l),zks* t(l) )
235 : c-----> loop over basis functions
236 0 : IF(zMat%l_real) THEN
237 0 : do n = 1,nslibd
238 0 : ac(l,n) = ac(l,n) + zMat%data_r(k+addnoco,n)*av
239 0 : bc(l,n) = bc(l,n) + zMat%data_r(k+addnoco,n)*bv
240 : enddo
241 : ELSE
242 0 : do n = 1,nslibd
243 0 : ac(l,n) = ac(l,n) + zMat%data_c(k+addnoco,n)*av
244 0 : bc(l,n) = bc(l,n) + zMat%data_c(k+addnoco,n)*bv
245 : enddo
246 : END IF
247 : enddo
248 :
249 : c...now for the k+b point
250 :
251 0 : evacp = evac_b(ivac)
252 0 : do ik = 1,nv2_b
253 0 : v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
254 0 : v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
255 0 : v(3) = 0.
256 0 : ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
257 : call vacuz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,t_b(ik),
258 0 : + dt_b(ik),u_b(1,ik))
259 : call vacudz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,te_b(ik),
260 : + dte_b(ik),tei_b(ik),ue_b(1,ik),dt_b(ik),
261 0 : + u_b(1,ik))
262 0 : scale = wronk/ (te_b(ik)*dt_b(ik)-dte_b(ik)*t_b(ik))
263 0 : te_b(ik) = scale*te_b(ik)
264 0 : dte_b(ik) = scale*dte_b(ik)
265 0 : tei_b(ik) = scale*tei_b(ik)
266 0 : do j = 1,nmz
267 0 : ue_b(j,ik) = scale*ue_b(j,ik)
268 : enddo
269 : enddo
270 : c-----> construct a and b coefficients for the k+b point
271 :
272 0 : do k = 1,nv_b(jspin_b)
273 0 : l = map2_b(k)
274 0 : zks = k3_b(k,jspin_b)*bmat(3,3)*sign
275 0 : arg = zks*z1
276 0 : c_1 = cmplx(cos(arg),sin(arg)) * const
277 0 : av = -c_1 * cmplx( dte_b(l),zks*te_b(l) )
278 0 : bv = c_1 * cmplx( dt_b(l),zks*t_b(l) )
279 : c-----> loop over basis functions
280 0 : IF(zMat_b%l_real) THEN
281 0 : do n = 1,nslibd_b
282 0 : ac_b(l,n) = ac_b(l,n) + zMat_b%data_r(k+addnoco,n)*av
283 0 : bc_b(l,n) = bc_b(l,n) + zMat_b%data_r(k+addnoco,n)*bv
284 : enddo
285 : ELSE
286 0 : do n = 1,nslibd_b
287 0 : ac_b(l,n) = ac_b(l,n) + zMat_b%data_c(k+addnoco,n)*av
288 0 : bc_b(l,n) = bc_b(l,n) + zMat_b%data_c(k+addnoco,n)*bv
289 : enddo
290 : END IF
291 : enddo
292 :
293 :
294 0 : do l = 1,nv2
295 0 : do lprime = 1,nv2_b
296 : if (kvac1(l).eq.(kvac1_b(lprime)-gb(1))
297 0 : & .and. kvac2(l).eq.(kvac2_b(lprime)-gb(2)))then
298 0 : zks = gb(3)*bmat(3,3)*sign
299 :
300 0 : do i = 1,nmz
301 : xx(np1-i) = u(i,l)*u_b(i,lprime)*
302 0 : * cos((z1+(i-1)*delz)*zks)
303 : xximag(np1-i) = u(i,l)*u_b(i,lprime)*
304 0 : * sin((z1+(i-1)*delz)*zks)
305 : enddo
306 0 : call intgz0(xx,delz,nmz,uu,tail)
307 0 : call intgz0(xximag,delz,nmz,uuimag,tail)
308 :
309 0 : do i = 1,nmz
310 : xx(np1-i) = u(i,l)*ue_b(i,lprime)*
311 0 : * cos((z1+(i-1)*delz)*zks)
312 : xximag(np1-i) = u(i,l)*ue_b(i,lprime)*
313 0 : * sin((z1+(i-1)*delz)*zks)
314 : enddo
315 0 : call intgz0(xx,delz,nmz,ud,tail)
316 0 : call intgz0(xximag,delz,nmz,udimag,tail)
317 :
318 0 : do i = 1,nmz
319 : xx(np1-i) = ue(i,l)*u_b(i,lprime)*
320 0 : * cos((z1+(i-1)*delz)*zks)
321 : xximag(np1-i) = ue(i,l)*u_b(i,lprime)*
322 0 : * sin((z1+(i-1)*delz)*zks)
323 : enddo
324 0 : call intgz0(xx,delz,nmz,du,tail)
325 0 : call intgz0(xximag,delz,nmz,duimag,tail)
326 0 : do i = 1,nmz
327 : xx(np1-i) = ue(i,l)*ue_b(i,lprime)*
328 0 : * cos((z1+(i-1)*delz)*zks)
329 : xximag(np1-i) = ue(i,l)*ue_b(i,lprime)*
330 0 : * sin((z1+(i-1)*delz)*zks)
331 : enddo
332 0 : call intgz0(xx,delz,nmz,dd,tail)
333 0 : call intgz0(xximag,delz,nmz,ddimag,tail)
334 :
335 0 : do i = 1,nslibd
336 0 : do j = 1,nslibd_b
337 : mmn(i,j) = mmn(i,j) + area*(
338 : * ac(l,i)*conjg(ac_b(lprime,j))*cmplx(uu,uuimag) +
339 : + ac(l,i)*conjg(bc_b(lprime,j))*cmplx(ud,udimag) +
340 : * bc(l,i)*conjg(ac_b(lprime,j))*cmplx(du,duimag) +
341 0 : + bc(l,i)*conjg(bc_b(lprime,j))*cmplx(dd,ddimag) )*vacchi
342 : enddo
343 : enddo
344 :
345 :
346 :
347 :
348 :
349 : endif ! kvac1=kvac1_b and kvac2=kvac2_b
350 : enddo ! lprime
351 : enddo ! l
352 :
353 : enddo !symvac
354 : c... cycle by the vacua finishes
355 0 : 140 enddo
356 :
357 0 : deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
358 0 : + v,kvac1,kvac2,map2 )
359 0 : deallocate ( ac_b,bc_b,dt_b,dte_b,t_b,te_b,tei_b,u_b,ue_b,
360 0 : + kvac1_b,kvac2_b,map2_b )
361 :
362 0 : call timestop("wann_mmkb_vac")
363 0 : END SUBROUTINE wann_mmkb_vac
364 : END MODULE m_wann_mmkb_vac
|
