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_mmk0_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_mmk0_vac(
18 0 : > l_noco,nlotot,qss,
19 : > z1,nmzd,nv2d,k1d,k2d,k3d,n3d,nvac,
20 : > ig,nmz,delz,ig2,area,bmat,
21 0 : > bbmat,evac,bkpt,vz,
22 0 : > nslibd,jspin,k1,k2,k3,jspd,nvd,
23 0 : > nbasfcn,neigd,zMat,nv,omtil,
24 0 : < mmn)
25 :
26 : USE m_types
27 : use m_constants
28 : USE m_vacuz
29 : USE m_vacudz
30 :
31 : implicit none
32 :
33 : TYPE(t_mat), INTENT(IN) :: zMat
34 :
35 : c .. scalar Arguments..
36 : logical, intent (in) :: l_noco
37 : integer, intent (in) :: nlotot
38 : real, intent (in) :: qss(:) !qss(3)
39 : integer, intent (in) :: nmzd,nv2d,k1d,k2d,k3d,n3d
40 : integer, intent (in) :: nmz,nslibd,nvac
41 : integer, intent (in) :: jspin,jspd,nvd
42 : integer, intent (in) :: nbasfcn,neigd
43 : real, intent (in) :: delz,z1,omtil,area
44 :
45 : c ..array arguments..
46 : real, intent (in) :: bkpt(:) !bkpt(3)
47 : real, intent (in) :: evac(:) !evac(2)
48 : integer, intent (in) :: ig(-k1d:,-k2d:,-k3d:) !ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
49 : integer, intent (in) :: ig2(:) !ig2(n3d)
50 : integer, intent (in) :: nv(:) !nv(jspd)
51 : real, intent (in) :: vz(:,:) !vz(nmzd,2)
52 : real, intent (in) :: bbmat(3,3),bmat(3,3)
53 : integer, intent (in) :: k1(:,:) !k1(nvd,jspd)
54 : integer, intent (in) :: k2(:,:) !k2(nvd,jspd)
55 : integer, intent (in) :: k3(:,:) !k3(nvd,jspd)
56 : complex, intent (inout) :: mmn(:,:) !mmn(nslibd,nslibd)
57 :
58 : c ..basis wavefunctions in the vacuum
59 0 : complex, allocatable :: ac(:,:),bc(:,:)
60 0 : real, allocatable :: dt(:),dte(:)
61 0 : real, allocatable :: t(:),te(:),tei(:)
62 0 : real, allocatable :: u(:,:),ue(:,:),v(:)
63 :
64 : c ..local scalars..
65 : real wronk,arg,zks,tpi,vz0(2),scale,evacp,ev,const
66 : real :: qss1,qss2
67 : integer i,m,l,j,k,n,nv2,ivac,n2,sign,ik,symvac,addnoco
68 : integer symvacvac
69 : complex av,bv,ic,c_1
70 0 : integer, allocatable :: kvac1(:),kvac2(:),map2(:)
71 :
72 0 : call timestart("wann_mmk0_vac")
73 :
74 : allocate ( ac(nv2d,neigd),bc(nv2d,neigd),dt(nv2d),
75 : + dte(nv2d),t(nv2d),te(nv2d),tei(nv2d),
76 : + u(nmzd,nv2d),ue(nmzd,nv2d),
77 0 : + v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd) )
78 :
79 0 : tpi = 2 * pimach() ; ic = cmplx(0.,1.)
80 :
81 : c.. determining the indexing array (in-plane stars)
82 :
83 0 : wronk = 2.0
84 0 : const = 1.0 / ( sqrt(omtil)*wronk )
85 :
86 0 : do ivac = 1,2
87 0 : vz0(ivac) = vz(nmz,ivac)
88 : enddo
89 :
90 0 : addnoco=0
91 0 : if(l_noco .and. jspin.eq.2)then
92 0 : addnoco=nv(1)+nlotot
93 : endif
94 :
95 0 : n2 = 0
96 0 : do 40 k = 1,nv(jspin)
97 0 : do 30 j = 1,n2
98 0 : if ( k1(k,jspin).eq.kvac1(j) .and.
99 : + k2(k,jspin).eq.kvac2(j) ) then
100 0 : map2(k) = j
101 0 : goto 40
102 : endif
103 0 : 30 continue
104 0 : n2 = n2 + 1
105 :
106 0 : IF (n2>nv2d) CALL juDFT_error("wannier Mmn vac",calledby
107 0 : + ="wann_mmk0_vac")
108 :
109 0 : kvac1(n2) = k1(k,jspin)
110 0 : kvac2(n2) = k2(k,jspin)
111 0 : map2(k) = n2
112 0 : 40 continue
113 :
114 : c...cycle by the vacua
115 0 : do 140 ivac = 1,nvac
116 :
117 :
118 0 : sign = 3. - 2.*ivac
119 0 : evacp = evac(ivac)
120 :
121 0 : nv2 = n2
122 :
123 : c.. the body of the routine
124 :
125 0 : qss1=0.0
126 0 : qss2=0.0
127 0 : if(l_noco.and.jspin.eq.1)then
128 0 : qss1=-qss(1)/2.0
129 0 : qss2=-qss(2)/2.0
130 0 : elseif(l_noco.and.jspin.eq.2)then
131 0 : qss1=qss(1)/2.0
132 0 : qss2=qss(2)/2.0
133 : endif
134 :
135 0 : do ik = 1,nv2
136 0 : v(1) = bkpt(1) + kvac1(ik) + qss1
137 0 : v(2) = bkpt(2) + kvac2(ik) + qss2
138 0 : v(3) = 0.
139 0 : ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
140 : call vacuz(ev,vz(1:,ivac),vz0(ivac),nmz,delz,t(ik),dt(ik),
141 0 : + u(1,ik))
142 : call vacudz(ev,vz(1:,ivac),vz0(ivac),nmz,delz,te(ik),
143 : + dte(ik),tei(ik),ue(1,ik),dt(ik),
144 0 : + u(1,ik))
145 0 : scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
146 0 : te(ik) = scale*te(ik)
147 0 : dte(ik) = scale*dte(ik)
148 0 : tei(ik) = scale*tei(ik)
149 0 : do j = 1,nmz
150 0 : ue(j,ik) = scale*ue(j,ik)
151 : enddo
152 : enddo
153 : c-----> construct a and b coefficients
154 :
155 0 : symvacvac=1
156 0 : if (nvac==1) symvacvac=2
157 0 : do symvac=1,symvacvac
158 0 : do 60 n = 1,nslibd
159 0 : do 50 i = 1,nv2d
160 0 : ac(i,n) = cmplx(0.0,0.0)
161 0 : bc(i,n) = cmplx(0.0,0.0)
162 0 : 50 continue
163 0 : 60 continue
164 :
165 0 : if (symvac==2) sign=-1.0
166 :
167 0 : do k = 1,nv(jspin)
168 0 : l = map2(k)
169 0 : zks = k3(k,jspin)*bmat(3,3)*sign
170 0 : arg = zks*z1
171 0 : c_1 = cmplx(cos(arg),sin(arg)) * const
172 0 : av = -c_1 * cmplx( dte(l),zks*te(l) )
173 0 : bv = c_1 * cmplx( dt(l),zks* t(l) )
174 : c-----> loop over basis functions
175 0 : IF(zMat%l_real) THEN
176 0 : do n = 1,nslibd
177 0 : ac(l,n) = ac(l,n) + zMat%data_r(k+addnoco,n)*av
178 0 : bc(l,n) = bc(l,n) + zMat%data_r(k+addnoco,n)*bv
179 : enddo
180 : ELSE
181 0 : do n = 1,nslibd
182 0 : ac(l,n) = ac(l,n) + zMat%data_c(k+addnoco,n)*av
183 0 : bc(l,n) = bc(l,n) + zMat%data_c(k+addnoco,n)*bv
184 : enddo
185 : END IF
186 : enddo
187 :
188 0 : do l = 1,nv2
189 0 : do i = 1,nslibd
190 0 : do j = 1,nslibd
191 : mmn(i,j) = mmn(i,j) +
192 : + area*(ac(l,i)*conjg( ac(l,j))
193 0 : + + tei(l)*bc(l,i)*conjg( bc(l,j)))
194 : enddo
195 : enddo
196 : enddo
197 : enddo !symvac
198 :
199 : c... cycle by the vacua finishes
200 0 : 140 enddo
201 :
202 0 : deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
203 0 : + v,kvac1,kvac2,map2 )
204 :
205 0 : call timestop("wann_mmk0_vac")
206 0 : end subroutine wann_mmk0_vac
207 : end module m_wann_mmk0_vac
|
