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