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,
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
23 : use m_wann_read_umatrix
24 :
25 : implicit none
26 :
27 : real, intent(in) :: bmat(:,:)
28 : real, intent(in) :: amat(:,:)
29 :
30 : integer, intent(in) :: rvecnum
31 : integer, intent(in) :: rvec(:,:)
32 : real, intent(in) :: kpoints(:,:)
33 : integer, intent(in) :: jspins_in
34 : integer, intent(in) :: nkpts
35 : logical, intent(in) :: l_bzsym
36 : logical, intent(in) :: film
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 : call timestart("wann_rmat")
99 :
100 0 : tpi=2*pimach()
101 :
102 0 : jspins=jspins_in
103 0 : if(l_nocosoc)jspins=1
104 :
105 0 : write(oUnit,*)"nkpts=",nkpts
106 :
107 0 : do jspin=1,jspins !spin loop
108 : c*****************************************************
109 : c get num_bands and num_wann from the proj file
110 : c*****************************************************
111 0 : do j=jspin,0,-1
112 0 : inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
113 0 : if(l_file)then
114 0 : filename='proj'//spinspin12(j)
115 0 : exit
116 : endif
117 : enddo
118 0 : if(l_file)then
119 0 : open (203,file=trim(filename),status='old')
120 0 : rewind (203)
121 : else
122 : CALL judft_error("no proj/proj.1/proj.2",calledby
123 0 : + ="wann_rmat")
124 : endif
125 0 : read (203,*) num_wann,num_bands
126 0 : close (203)
127 0 : write(oUnit,*)'According to proj there are ',
128 0 : + num_bands,' bands'
129 0 : write(oUnit,*)"and ",num_wann," wannier functions."
130 :
131 : c****************************************************************
132 : c get nntot and bk and wb from bkpts file
133 : c****************************************************************
134 0 : inquire (file='bkpts',exist=l_file)
135 0 : if (.not.l_file) CALL judft_error("need bkpts"
136 0 : + ,calledby ="wann_rmat")
137 0 : open (202,file='bkpts',form='formatted',status='old')
138 0 : read (202,'(i4)') nntot
139 0 : allocate ( gb(3,nntot,nkpts) )
140 0 : allocate ( bpt(nntot,nkpts) )
141 0 : do ikpt=1,nkpts
142 0 : do nn=1,nntot
143 : read (202,'(2i6,3x,3i4)')
144 0 : & ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
145 : enddo !nn
146 : enddo !ikpt
147 0 : allocate( wb(nntot) )
148 0 : allocate( bpunkt(3,nntot) )
149 0 : do nn=1,nntot
150 0 : read(202,*)bpunkt(:,nn),wb(nn)
151 : enddo
152 0 : close (202)
153 :
154 0 : do ikpt=1,nkpts
155 0 : do nn=1,nntot
156 : kdiff(1)=kpoints(1,bpt(nn,ikpt))+
157 0 : + gb(1,nn,ikpt) - kpoints(1,ikpt)
158 : kdiff(2)=kpoints(2,bpt(nn,ikpt))+
159 0 : + gb(2,nn,ikpt) - kpoints(2,ikpt)
160 : kdiff(3)=kpoints(3,bpt(nn,ikpt))+
161 0 : + gb(3,nn,ikpt) - kpoints(3,ikpt)
162 0 : kdiff=matmul(transpose(bmat),kdiff)
163 0 : l_worksout=.false.
164 0 : do nn2=1,nntot
165 0 : if (abs(kdiff(1)-bpunkt(1,nn2)).lt.1.e-5)then
166 0 : if (abs(kdiff(2)-bpunkt(2,nn2)).lt.1.e-5)then
167 0 : if (abs(kdiff(3)-bpunkt(3,nn2)).lt.1.e-5)then
168 : l_worksout=.true.
169 : exit
170 : endif
171 : endif
172 : endif
173 : enddo !nn2
174 0 : if(.not.l_worksout)then
175 0 : write(*,*)"ikpt,nn=",ikpt,nn
176 0 : write(*,*)"kdiff(1)=",kdiff(1)
177 0 : write(*,*)"kdiff(2)=",kdiff(2)
178 0 : write(*,*)"kdiff(3)=",kdiff(3)
179 0 : stop 'worksout'
180 : endif
181 : enddo !nn
182 : enddo !ikpt
183 :
184 : c****************************************************************
185 : c read in chk
186 : c****************************************************************
187 0 : num_kpts=nkpts
188 0 : allocate( u_matrix_opt(num_bands,num_wann,nkpts) )
189 0 : allocate( u_matrix(num_wann,num_wann,nkpts) )
190 0 : allocate( lwindow(num_bands,nkpts) )
191 0 : allocate( ndimwin(nkpts) )
192 0 : allocate( m_matrix(num_wann,num_wann,nntot,num_kpts) )
193 : call wann_read_umatrix2(
194 : > nkpts,num_wann,num_bands,
195 : > um_format,jspin,wan90version,
196 : < have_disentangled,
197 : < lwindow,ndimwin,
198 0 : < u_matrix_opt,u_matrix,m_matrix)
199 : c****************************************************************
200 : c read in eig-file
201 : c****************************************************************
202 0 : write(oUnit,*)"read in eig-file"
203 0 : allocate(energy(num_bands,num_kpts))
204 0 : inquire(file=spin12(jspin)//'.eig',exist=l_umdat)
205 0 : IF(.NOT.l_umdat) CALL judft_error
206 : + ("Thou shall not hide your eig file",calledby
207 0 : + ="wann_hopping")
208 0 : open(300,file=spin12(jspin)//'.eig',form='formatted')
209 0 : do i=1,num_kpts
210 0 : do j=1,num_bands
211 0 : read(300,*)band,kpt,energy(j,i)
212 : enddo
213 : enddo
214 0 : close(300)
215 :
216 0 : minenerg=minval(energy(:,:))
217 0 : maxenerg=maxval(energy(:,:))
218 0 : write(oUnit,*)"minenerg=",minenerg
219 0 : write(oUnit,*)"maxenerg=",maxenerg
220 :
221 :
222 0 : allocate(eigval_opt(num_bands,nkpts))
223 0 : allocate(eigval2(num_wann,nkpts))
224 0 : eigval_opt=0.0
225 0 : eigval2=0.0
226 :
227 0 : if(have_disentangled) then
228 :
229 0 : do nkp=1,num_kpts
230 0 : counter=0
231 0 : do j=1,num_bands
232 0 : if(lwindow(j,nkp)) then
233 0 : counter=counter+1
234 0 : eigval_opt(counter,nkp)=energy(j,nkp)
235 : end if
236 : end do
237 : end do
238 :
239 0 : do nkp=1,num_kpts
240 0 : do j=1,num_wann
241 0 : do m=1,ndimwin(nkp)
242 : eigval2(j,nkp)=eigval2(j,nkp)+eigval_opt(m,nkp)*
243 0 : & real(conjg(u_matrix_opt(m,j,nkp))*u_matrix_opt(m,j,nkp))
244 : enddo
245 : enddo
246 : enddo
247 :
248 : else
249 0 : eigval2 = energy
250 : end if !have_disentangled
251 :
252 0 : deallocate(eigval_opt)
253 0 : deallocate(energy)
254 :
255 : c****************************************************************
256 : c Set up posop.
257 : c****************************************************************
258 0 : write(oUnit,*)"Set up posop."
259 0 : allocate( posop2(3,num_wann,num_wann,num_kpts) )
260 0 : posop2=cmplx(0.0,0.0)
261 :
262 0 : do ikpt=1,nkpts
263 0 : do nn=1,nntot
264 : kdiff(1)=kpoints(1,bpt(nn,ikpt))+
265 0 : + gb(1,nn,ikpt) - kpoints(1,ikpt)
266 : kdiff(2)=kpoints(2,bpt(nn,ikpt))+
267 0 : + gb(2,nn,ikpt) - kpoints(2,ikpt)
268 : kdiff(3)=kpoints(3,bpt(nn,ikpt))+
269 0 : + gb(3,nn,ikpt) - kpoints(3,ikpt)
270 0 : kdiff=matmul(transpose(bmat),kdiff)
271 0 : do i=1,num_wann
272 0 : do j=1,num_wann
273 0 : if(j.eq.i)then
274 0 : do kk=1,3
275 : posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
276 0 : & wb(nn)*kdiff(kk)*(m_matrix(j,i,nn,ikpt)-1.0)
277 : enddo !kk
278 : else
279 0 : do kk=1,3
280 : posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
281 0 : & wb(nn)*kdiff(kk)*m_matrix(j,i,nn,ikpt)
282 : enddo !kk
283 : endif
284 : enddo !j
285 : enddo !i
286 : enddo !nn
287 : enddo !ikpt
288 :
289 0 : do ikpt=1,nkpts
290 0 : do i=1,num_wann
291 0 : do j=1,i
292 0 : do kk=1,3
293 : posop2(kk,j,i,ikpt)=(posop2(kk,j,i,ikpt)+
294 0 : & conjg(posop2(kk,i,j,ikpt)))/2.0
295 0 : posop2(kk,i,j,ikpt)=posop2(kk,j,i,ikpt)
296 : enddo !kk
297 : enddo !j
298 : enddo !i
299 : enddo !ikpt
300 :
301 0 : allocate( posop(3,num_wann,num_wann,rvecnum) )
302 0 : posop=cmplx(0.0,0.0)
303 0 : do rvecind=1,rvecnum
304 0 : do k=1,nkpts
305 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
306 : & kpoints(2,k)*rvec(2,rvecind)+
307 0 : & kpoints(3,k)*rvec(3,rvecind) )
308 0 : fac=cmplx(cos(rdotk),-sin(rdotk))
309 0 : do i=1,num_wann
310 0 : do j=1,num_wann
311 0 : do kk=1,3
312 : posop(kk,j,i,rvecind)=
313 : & posop(kk,j,i,rvecind)+
314 0 : & fac*posop2(kk,j,i,k)
315 : enddo !kk
316 : enddo !j
317 : enddo !i
318 : enddo !k
319 : enddo !rvecind
320 0 : posop=posop/cmplx(real(nkpts),0.0)
321 0 : deallocate( posop2 )
322 :
323 0 : do rvecind=1,rvecnum
324 0 : if(rvec(1,rvecind).eq.0)then
325 0 : if(rvec(2,rvecind).eq.0)then
326 0 : if(rvec(3,rvecind).eq.0)then
327 0 : rvecind_0=rvecind
328 : goto 123
329 : endif
330 : endif
331 : endif
332 : enddo !rvecind
333 0 : stop 'Ou est ce point-la ?'
334 : 123 continue
335 :
336 0 : do i=1,num_wann
337 0 : do kk=1,3
338 0 : posop(kk,i,i,rvecind_0)=0.0
339 : enddo !kk
340 : enddo !i
341 :
342 0 : do ikpt=1,nkpts
343 0 : do nn=1,nntot
344 : kdiff(1)=kpoints(1,bpt(nn,ikpt))+
345 0 : + gb(1,nn,ikpt) - kpoints(1,ikpt)
346 : kdiff(2)=kpoints(2,bpt(nn,ikpt))+
347 0 : + gb(2,nn,ikpt) - kpoints(2,ikpt)
348 : kdiff(3)=kpoints(3,bpt(nn,ikpt))+
349 0 : + gb(3,nn,ikpt) - kpoints(3,ikpt)
350 0 : kdiff=matmul(transpose(bmat),kdiff)/real(nkpts)
351 0 : do i=1,num_wann
352 0 : do kk=1,3
353 : posop(kk,i,i,rvecind_0)=posop(kk,i,i,rvecind_0)-
354 0 : & wb(nn)*kdiff(kk)*aimag(log(m_matrix(i,i,nn,ikpt)))
355 : enddo !kk
356 : enddo !i
357 : enddo !nn
358 : enddo !ikpt
359 :
360 : c********************************************
361 : c Print posop.
362 : c********************************************
363 0 : open(321,file='posop'//spinspin12(jspin),form='formatted')
364 0 : do rvecind=1,rvecnum
365 0 : r3=rvec(3,rvecind)
366 0 : r2=rvec(2,rvecind)
367 0 : r1=rvec(1,rvecind)
368 0 : do j=1,num_wann
369 0 : do i=1,num_wann
370 0 : do kk=1,3
371 : write(321,'(i3,1x,i3,1x,i3,1x,i3,
372 : & 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
373 0 : & r1,r2,r3,i,j,kk,posop(kk,i,j,rvecind)
374 : enddo !kk
375 : enddo !j
376 : enddo !i
377 : enddo !rvecind
378 0 : close(321)
379 :
380 0 : deallocate( lwindow,u_matrix_opt,ndimwin )
381 0 : deallocate( u_matrix,m_matrix )
382 0 : deallocate( posop,eigval2 )
383 0 : deallocate( gb,bpt,wb,bpunkt )
384 : enddo !jspin
385 :
386 0 : call timestop("wann_rmat")
387 0 : end subroutine wann_rmat
388 0 : end module m_wann_rmat
|
