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