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_nabla_rs
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_nabla_rs(
11 0 : > rvecnum,rvec,kpoints,
12 : > jspins_in,nkpts,l_bzsym,film,
13 : > l_soc,band_min,band_max,neigd,
14 : > l_socmmn0,wan90version)
15 : c*************************************************
16 : c Calculate the matrix elements of nabla
17 : c in real space from the
18 : c files WF1.chk (and WF1_um.dat) produced
19 : c by wannier90.
20 : c FF, December 2009
21 : c*************************************************
22 : use m_constants
23 : use m_wann_read_umatrix
24 :
25 : implicit none
26 : integer, intent(in) :: rvecnum
27 : integer, intent(in) :: rvec(:,:)
28 : real, intent(in) :: kpoints(:,:)
29 : integer, intent(in) :: jspins_in
30 : integer, intent(in) :: nkpts
31 : logical, intent(in) :: l_bzsym
32 : logical, intent(in) :: film
33 :
34 : logical, intent(in) :: l_soc
35 : integer, intent(in) :: band_min(2),band_max(2),neigd
36 :
37 : logical, intent(in) :: l_socmmn0
38 : integer, intent(in) :: wan90version
39 :
40 : integer :: ikpt,jspins
41 : integer :: kpts
42 : logical :: l_file
43 : c real :: kpoints(3,nkpts)
44 : integer :: num_wann,num_kpts,num_nnmax,jspin
45 : integer :: kspin,kkspin
46 : integer :: wann_shift,num_wann2
47 : integer :: i,j,k,m,info,r1,r2,r3,dummy1
48 : integer :: dummy2,dummy3,dummy4
49 : integer :: hopmin,hopmax,counter,m1,m2
50 : integer :: num_bands2
51 : integer,allocatable :: iwork(:)
52 : real,allocatable :: energy(:,:),ei(:)
53 : real,allocatable :: eigw(:,:),rwork(:)
54 : complex,allocatable :: work(:),vec(:,:)
55 0 : complex,allocatable :: u_matrix(:,:,:),hwann(:,:,:,:)
56 0 : complex,allocatable :: hreal(:,:,:,:)
57 : complex,allocatable :: hrealsoc(:,:,:,:,:,:,:)
58 : complex,allocatable :: hwannsoc(:,:,:,:,:)
59 0 : complex,allocatable :: nablamat(:,:,:,:,:)
60 0 : complex,allocatable :: nablamat2(:,:,:,:)
61 : complex :: fac,eulav,eulav1
62 : real :: tmp_omi,rdotk,tpi,minenerg,maxenerg
63 : real, allocatable :: minieni(:),maxieni(:)
64 : character :: jobz,uplo
65 : integer :: kpt,band,lee,lwork,lrwork,liwork,n,lda
66 : complex :: value(4)
67 : logical :: um_format,l_she
68 : logical :: repro_eig
69 : logical :: l_chk,l_proj
70 : logical :: have_disentangled
71 0 : integer,allocatable :: ndimwin(:)
72 0 : logical,allocatable :: lwindow(:,:)
73 : integer :: chk_unit,nkp,ntmp,ierr
74 : character(len=33) :: header
75 : character(len=20) :: checkpoint
76 : real :: tmp_latt(3,3), tmp_kpt_latt(3,nkpts)
77 : real :: omega_invariant
78 0 : complex,allocatable :: u_matrix_opt(:,:,:)
79 : integer :: num_bands
80 : logical :: l_umdat
81 : real,allocatable :: eigval2(:,:)
82 : real,allocatable :: eigval_opt(:,:)
83 : real :: scale,a,b
84 : character(len=2) :: spinspin12(0:2)
85 : character(len=3) :: spin12(2)
86 : character(len=6) :: filename
87 : integer :: jp,mp,kk
88 : integer :: rvecind
89 :
90 : data spinspin12/' ','.1' , '.2'/
91 : data spin12/'WF1','WF2'/
92 :
93 0 : call timestart("wann_nabla_rs")
94 0 : inquire(file='she',exist=l_she)
95 :
96 0 : tpi=2*pimach()
97 :
98 0 : jspins=jspins_in
99 0 : if(l_soc)jspins=1
100 :
101 0 : write(oUnit,*)"nkpts=",nkpts
102 :
103 : c$$$c***************************************************
104 : c$$$c read in the kpoints from w90kpts or kpts
105 : c$$$c***************************************************
106 : c$$$ if (l_bzsym) then
107 : c$$$ l_file=.false.
108 : c$$$ inquire(file='w90kpts',exist=l_file)
109 : c$$$ if (.not.l_file) stop 'where is w90kpts?'
110 : c$$$ open(177,file='w90kpts',form='formatted')
111 : c$$$ read(177,*)kpts,scale
112 : c$$$ else
113 : c$$$ l_file=.false.
114 : c$$$ inquire(file='kpts',exist=l_file)
115 : c$$$ if(.not.l_file) stop 'where is kpts?'
116 : c$$$ open(177,file='kpts',form='formatted')
117 : c$$$ read(177,*)kpts,scale
118 : c$$$ endif
119 : c$$$ if(kpts.ne.nkpts) stop 'mismatch in number of kpoints'
120 : c$$$
121 : c$$$ do ikpt = 1,nkpts
122 : c$$$ read(177,*)kpoints(:,ikpt)
123 : c$$$ enddo
124 : c$$$
125 : c$$$ if(film)then
126 : c$$$ kpoints(3,:)=0.0
127 : c$$$ endif
128 : c$$$ kpoints=kpoints/scale
129 : c$$$ close(177)
130 :
131 0 : do jspin=1,jspins !spin loop
132 : c*****************************************************
133 : c get num_bands and num_wann from the proj file
134 : c*****************************************************
135 0 : do j=jspin,0,-1
136 0 : inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
137 0 : if(l_file)then
138 0 : filename='proj'//spinspin12(j)
139 0 : exit
140 : endif
141 : enddo
142 0 : if(l_file)then
143 0 : open (203,file=trim(filename),status='old')
144 0 : rewind (203)
145 : else
146 : CALL juDFT_error("no proj/proj.1/proj.2",calledby
147 0 : + ="wann_nabla_rs")
148 : endif
149 0 : read (203,*) num_wann,num_bands
150 0 : close (203)
151 0 : write(oUnit,*)'According to proj there are ',
152 0 : + num_bands,' bands'
153 0 : write(oUnit,*)"and ",num_wann," wannier functions."
154 :
155 : c****************************************************************
156 : c read in chk
157 : c****************************************************************
158 0 : num_kpts=nkpts
159 0 : allocate( u_matrix_opt(num_bands,num_wann,nkpts) )
160 0 : allocate( u_matrix(num_wann,num_wann,nkpts) )
161 0 : allocate( lwindow(num_bands,nkpts) )
162 0 : allocate( ndimwin(nkpts) )
163 : call wann_read_umatrix2(
164 : > nkpts,num_wann,num_bands,
165 : > um_format,jspin,wan90version,
166 : < have_disentangled,
167 : < lwindow,ndimwin,
168 0 : < u_matrix_opt,u_matrix)
169 :
170 0 : num_bands2=num_bands
171 0 : if(l_soc.and.l_socmmn0)then
172 0 : num_bands2=neigd
173 : endif
174 :
175 : c************************************************
176 : c Read the files "WF1.nabl" and "WF2.nabl"
177 : c************************************************
178 0 : allocate( nablamat(3,num_bands2,num_bands2,nkpts,2) )
179 0 : open(304,file=spin12(jspin)//'.nabl',form='formatted')
180 0 : read(304,*)
181 0 : read(304,*)
182 0 : do nkp=1,num_kpts
183 0 : do i=1,num_bands2
184 0 : do j=1,num_bands2
185 0 : do k=1,3
186 0 : read(304,*)dummy1,dummy2,dummy3,dummy4,a,b
187 0 : nablamat(k,j,i,nkp,1)=cmplx(a,-b)
188 : enddo !k
189 : enddo !j
190 : enddo !i
191 : enddo !nkp
192 0 : close(304)
193 :
194 0 : if(l_soc)then
195 0 : open(304,file=spin12(2)//'.nabl',form='formatted')
196 0 : read(304,*)
197 0 : read(304,*)
198 0 : do nkp=1,num_kpts
199 0 : do i=1,num_bands2
200 0 : do j=1,num_bands2
201 0 : do k=1,3
202 0 : read(304,*)dummy1,dummy2,dummy3,dummy4,a,b
203 0 : nablamat(k,j,i,nkp,2)=cmplx(a,-b)
204 : enddo !k
205 : enddo !j
206 : enddo !i
207 : enddo !nkp
208 0 : close(304)
209 : endif
210 :
211 : c**************************************
212 : c Enforce Hermiticity
213 : c**************************************
214 0 : do nkp=1,num_kpts
215 0 : do i=1,num_bands2
216 0 : do j=1,i
217 0 : do k=1,3
218 : nablamat(k,i,j,nkp,1)=( nablamat(k,i,j,nkp,1)+
219 0 : & conjg(nablamat(k,j,i,nkp,1)) )/2.0
220 0 : nablamat(k,j,i,nkp,1)=conjg(nablamat(k,i,j,nkp,1))
221 0 : if(l_soc)then
222 : nablamat(k,i,j,nkp,2)=( nablamat(k,i,j,nkp,2)+
223 0 : & conjg(nablamat(k,j,i,nkp,2)) )/2.0
224 0 : nablamat(k,j,i,nkp,2)=conjg(nablamat(k,i,j,nkp,2))
225 : endif
226 : enddo !k
227 : enddo !j
228 : enddo !i
229 : enddo !nkp
230 :
231 : c**************************************
232 : c Nabla-sum
233 : c**************************************
234 0 : if( l_soc )then
235 0 : if(l_she)then
236 0 : do nkp=1,num_kpts
237 0 : do i=1,num_bands2
238 0 : do j=1,num_bands2
239 0 : do k=1,3
240 : nablamat(k,j,i,nkp,1)=nablamat(k,j,i,nkp,1)
241 0 : & -nablamat(k,j,i,nkp,2)
242 : enddo !k
243 : enddo !j
244 : enddo !i
245 : enddo !nkp
246 : else
247 0 : do nkp=1,num_kpts
248 0 : do i=1,num_bands2
249 0 : do j=1,num_bands2
250 0 : do k=1,3
251 : nablamat(k,j,i,nkp,1)=nablamat(k,j,i,nkp,1)
252 0 : & +nablamat(k,j,i,nkp,2)
253 : enddo !k
254 : enddo !j
255 : enddo !i
256 : enddo !nkp
257 : endif
258 : endif !jspins.eq.2
259 :
260 : c****************************************************************
261 : c Calculate matrix elements of Nabla in the basis of
262 : c rotated Bloch functions.
263 : c****************************************************************
264 0 : allocate( nablamat2(3,num_wann,num_wann,nkpts) )
265 : write(oUnit,*)"calculate matrix elements of momentum operator
266 0 : & between wannier orbitals"
267 :
268 0 : if(have_disentangled) then
269 0 : nablamat2=0.0
270 0 : do nkp=1,num_kpts
271 0 : do j=1,num_wann
272 0 : do jp=1,num_wann
273 0 : do m=1,ndimwin(nkp)
274 0 : do mp=1,ndimwin(nkp)
275 0 : do k=1,3
276 : nablamat2(k,jp,j,nkp)=nablamat2(k,jp,j,nkp)+
277 : & conjg(u_matrix_opt(mp,jp,nkp))*
278 : & nablamat(k,mp,m,nkp,1)*
279 0 : & u_matrix_opt(m,j,nkp)
280 : enddo !k
281 : enddo !mp
282 : enddo !m
283 : enddo !jp
284 : enddo !j
285 : enddo !nkp
286 : else
287 0 : nablamat2(:,:,:,:)=nablamat(:,:,:,:,1)
288 : end if !have_disentangled
289 :
290 0 : allocate(hwann(3,num_wann,num_wann,num_kpts))
291 0 : hwann=cmplx(0.0,0.0)
292 0 : wann_shift=0
293 : if(l_socmmn0)then
294 : wann_shift=band_min(jspin)-1
295 : endif
296 0 : do k=1,num_kpts
297 0 : do m=1,num_wann
298 0 : do i=1,num_wann
299 0 : do j=1,num_wann
300 0 : do mp=1,num_wann
301 0 : do kk=1,3
302 : hwann(kk,mp,m,k)=hwann(kk,mp,m,k)+
303 : * conjg(u_matrix(j,mp,k))*
304 : * nablamat2(kk,j,i,k)*
305 0 : * u_matrix(i,m,k)
306 : enddo !kk
307 : enddo !mp
308 : enddo !j
309 : enddo !i
310 : enddo !m
311 : enddo !k
312 :
313 : c************************************************************
314 : c Calculate matrix elements in real space.
315 : c***********************************************************
316 0 : write(oUnit,*)"calculate nabla-mat in rs"
317 : c$$$ hopmin=-5
318 : c$$$ hopmax=5
319 : c$$$ allocate(hreal(3,num_wann,num_wann,hopmin:hopmax,
320 : c$$$ & hopmin:hopmax,hopmin:hopmax))
321 : c$$$ hreal=cmplx(0.0,0.0)
322 : c$$$ do r3=hopmin,hopmax
323 : c$$$ do r2=hopmin,hopmax
324 : c$$$ do r1=hopmin,hopmax
325 : c$$$ do k=1,nkpts
326 : c$$$ rdotk=tpi*(kpoints(1,k)*r1+kpoints(2,k)*r2+
327 : c$$$ & kpoints(3,k)*r3)
328 : c$$$ fac=cmplx(cos(rdotk),-sin(rdotk))/nkpts
329 : c$$$ do i=1,num_wann
330 : c$$$ do j=1,num_wann
331 : c$$$ do kk=1,3
332 : c$$$ hreal(kk,j,i,r1,r2,r3)=
333 : c$$$ & hreal(kk,j,i,r1,r2,r3)+
334 : c$$$ & fac*hwann(kk,j,i,k)
335 : c$$$ enddo !kk
336 : c$$$ enddo !j
337 : c$$$ enddo !i
338 : c$$$ enddo !k
339 : c$$$ enddo !r1
340 : c$$$ enddo !r2
341 : c$$$ enddo !r3
342 :
343 :
344 0 : allocate(hreal(3,num_wann,num_wann,rvecnum))
345 0 : hreal=cmplx(0.0,0.0)
346 :
347 0 : do rvecind=1,rvecnum
348 0 : do k=1,nkpts
349 : rdotk=tpi*(kpoints(1,k)*rvec(1,rvecind)+
350 : & kpoints(2,k)*rvec(2,rvecind)+
351 0 : & kpoints(3,k)*rvec(3,rvecind) )
352 0 : fac=cmplx(cos(rdotk),-sin(rdotk))
353 0 : do i=1,num_wann
354 0 : do j=1,num_wann
355 0 : do kk=1,3
356 : hreal(kk,j,i,rvecind)=
357 : & hreal(kk,j,i,rvecind)+
358 0 : & fac*hwann(kk,j,i,k)
359 : enddo !kk
360 : enddo !j
361 : enddo !i
362 : enddo !k
363 : enddo !rvecind
364 0 : hreal=hreal/cmplx(real(nkpts),0.0)
365 :
366 0 : if(l_she)then
367 0 : open(321,file='rs_sv'//spinspin12(jspin),form='formatted')
368 : else
369 0 : open(321,file='rsnabla'//spinspin12(jspin),form='formatted')
370 : endif
371 :
372 0 : do rvecind=1,rvecnum
373 0 : r3=rvec(3,rvecind)
374 0 : r2=rvec(2,rvecind)
375 0 : r1=rvec(1,rvecind)
376 0 : do j=1,num_wann
377 0 : do i=1,num_wann
378 0 : do kk=1,3
379 : write(321,'(i3,1x,i3,1x,i3,1x,i3,
380 : & 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
381 0 : & r1,r2,r3,i,j,kk,hreal(kk,i,j,rvecind)
382 : enddo !kk
383 : enddo !j
384 : enddo !i
385 : enddo !rvecind
386 0 : close(321)
387 :
388 0 : deallocate(lwindow,u_matrix_opt,ndimwin)
389 0 : deallocate(u_matrix,hwann,hreal)
390 0 : deallocate(nablamat,nablamat2)
391 : enddo !jspin
392 :
393 0 : call timestop("wann_nabla_rs")
394 0 : end subroutine wann_nabla_rs
395 : end module m_wann_nabla_rs
|
