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,l_onedimens,
20 : > l_soc,band_min,band_max,neigd,
21 : > l_socmmn0,wan90version)
22 :
23 : use m_constants, only: pimach, ImagUnit
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,l_onedimens
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 : tpi=2*pimach()
98 :
99 0 : jspins=jspins_in
100 0 : if(l_soc)jspins=1
101 :
102 0 : write(6,*)"nkpts=",nkpts
103 :
104 : c$$$c***************************************************
105 : c$$$c read in the kpoints from w90kpts or kpts
106 : c$$$c***************************************************
107 : c$$$ if (l_bzsym) then
108 : c$$$ l_file=.false.
109 : c$$$ inquire(file='w90kpts',exist=l_file)
110 : c$$$ if (.not.l_file) stop 'where is w90kpts?'
111 : c$$$ open(177,file='w90kpts',form='formatted')
112 : c$$$ read(177,*)kpts,scale
113 : c$$$ else
114 : c$$$ l_file=.false.
115 : c$$$ inquire(file='kpts',exist=l_file)
116 : c$$$ if(.not.l_file) stop 'where is kpts?'
117 : c$$$ open(177,file='kpts',form='formatted')
118 : c$$$ read(177,*)kpts,scale
119 : c$$$ endif
120 : c$$$ if(kpts.ne.nkpts) stop 'mismatch in number of kpoints'
121 : c$$$
122 : c$$$ do ikpt = 1,nkpts
123 : c$$$ read(177,*)kpoints(:,ikpt)
124 : c$$$ enddo
125 : c$$$
126 : c$$$ if(film.and..not.l_onedimens)then
127 : c$$$ kpoints(3,:)=0.0
128 : c$$$ endif
129 : c$$$ kpoints=kpoints/scale
130 : c$$$ close(177)
131 :
132 0 : do jspin=1,jspins !spin loop
133 : c*****************************************************
134 : c get num_bands and num_wann from the proj file
135 : c*****************************************************
136 0 : do j=jspin,0,-1
137 0 : inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
138 0 : if(l_file)then
139 0 : filename='proj'//spinspin12(j)
140 0 : exit
141 : endif
142 : enddo
143 0 : if(l_file)then
144 0 : open (203,file=trim(filename),status='old')
145 0 : rewind (203)
146 : else
147 : CALL juDFT_error("no proj/proj.1/proj.2",calledby
148 0 : + ="wann_nabla_pauli_rs")
149 : endif
150 0 : read (203,*) num_wann,num_bands
151 0 : close (203)
152 0 : write(6,*)'According to proj there are ',num_bands,' bands'
153 0 : write(6,*)"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 file "nabl.updown".
177 : c****************************************************
178 0 : allocate( paulimat(3,2,num_bands2,num_bands2,nkpts) )
179 0 : open(304,file='nabl.updown',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 : paulimat(k,1,i,j,nkp)=cmplx(a,b)
188 : enddo !k
189 : enddo !j
190 : enddo !i
191 0 : do k=1,3
192 0 : paulimat(k,2,:,:,nkp)=ImagUnit*paulimat(k,1,:,:,nkp)
193 : paulimat(k,1,:,:,nkp)=paulimat(k,1,:,:,nkp)+
194 0 : & transpose(conjg( paulimat(k,1,:,:,nkp) ))
195 : paulimat(k,2,:,:,nkp)=paulimat(k,2,:,:,nkp)+
196 0 : & transpose(conjg( paulimat(k,2,:,:,nkp) ))
197 : enddo !k
198 : enddo !nkp
199 0 : close(304)
200 :
201 :
202 : c****************************************************************
203 : c Calculate matrix elements of Pauli in the basis of
204 : c rotated Bloch functions.
205 : c****************************************************************
206 0 : allocate( paulimat2(3,2,num_wann,num_wann,nkpts) )
207 : write(6,*)"calculate matrix elements of momentum operator
208 0 : & between wannier orbitals"
209 :
210 0 : if(have_disentangled) then
211 0 : paulimat2=0.0
212 0 : do nkp=1,num_kpts
213 0 : do j=1,num_wann
214 0 : do jp=1,num_wann
215 0 : do m=1,ndimwin(nkp)
216 0 : do mp=1,ndimwin(nkp)
217 0 : do k=1,2
218 0 : do kk=1,3
219 : paulimat2(kk,k,jp,j,nkp)=paulimat2(kk,k,jp,j,nkp)+
220 : & conjg(u_matrix_opt(mp,jp,nkp))*
221 : & paulimat(kk,k,mp,m,nkp)*
222 0 : & u_matrix_opt(m,j,nkp)
223 : enddo !kk
224 : enddo !k
225 : enddo !mp
226 : enddo !m
227 : enddo !jp
228 : enddo !j
229 : enddo !nkp
230 : else
231 0 : paulimat2(:,:,:,:,:)=paulimat(:,:,:,:,:)
232 : end if !have_disentangled
233 :
234 0 : allocate(hwann(3,2,num_wann,num_wann,num_kpts))
235 0 : hwann=cmplx(0.0,0.0)
236 0 : wann_shift=0
237 : if(l_socmmn0)then
238 : wann_shift=band_min(jspin)-1
239 : endif
240 0 : do k=1,num_kpts
241 0 : do m=1,num_wann
242 0 : do mp=1,num_wann
243 0 : do i=1,num_wann
244 0 : do j=1,num_wann
245 0 : do kk=1,2
246 0 : do kkk=1,3
247 : hwann(kkk,kk,mp,m,k)=hwann(kkk,kk,mp,m,k)+
248 : * conjg(u_matrix(j,mp,k))*
249 : * paulimat2(kkk,kk,j,i,k)*
250 0 : * u_matrix(i,m,k)
251 : enddo !kkk
252 : enddo !kk
253 : enddo !j
254 : enddo !i
255 : enddo !mp
256 : enddo !m
257 : enddo !k
258 :
259 : c************************************************************
260 : c Calculate matrix elements in real space.
261 : c************************************************************
262 0 : write(6,*)"calculate pauli-mat in rs"
263 :
264 : c$$$ if(.false.)then !specify r-mesh by its boundaries
265 : c$$$ hopmin_z=-5;hopmax_z=5
266 : c$$$ hopmin_x=0;hopmax_x=0
267 : c$$$ hopmin_y=0;hopmax_y=0
268 : c$$$ rvecnum=(hopmax_z-hopmin_z+1)
269 : c$$$ if(.not.l_onedimens.and.film)then
270 : c$$$ hopmin_x=-5;hopmax_x=5
271 : c$$$ hopmin_y=-5;hopmax_y=5
272 : c$$$ hopmin_z=0; hopmax_z=0
273 : c$$$ else
274 : c$$$ hopmin_x=-5;hopmax_x=5
275 : c$$$ hopmin_y=-5;hopmax_y=5
276 : c$$$ endif
277 : c$$$ rvecnum= (hopmax_z-hopmin_z+1)
278 : c$$$ rvecnum=rvecnum*(hopmax_y-hopmin_y+1)
279 : c$$$ rvecnum=rvecnum*(hopmax_x-hopmin_x+1)
280 : c$$$
281 : c$$$ allocate(rvec(3,rvecnum))
282 : c$$$ rvecind=0
283 : c$$$ do r3=hopmin_z,hopmax_z
284 : c$$$ do r2=hopmin_y,hopmax_y
285 : c$$$ do r1=hopmin_x,hopmax_x
286 : c$$$ rvecind=rvecind+1
287 : c$$$ if(rvecind.gt.rvecnum)stop 'wann_hopping:1'
288 : c$$$ rvec(1,rvecind)=r1
289 : c$$$ rvec(2,rvecind)=r2
290 : c$$$ rvec(3,rvecind)=r3
291 : c$$$ enddo !r1
292 : c$$$ enddo !r2
293 : c$$$ enddo !r3
294 : c$$$ else !determine optimal r-mesh
295 : c$$$ call wann_wigner_seitz(
296 : c$$$ > .true.,num,amat,0,
297 : c$$$ < rvecnum,rvec)
298 : c$$$ allocate(rvec(3,rvecnum))
299 : c$$$ call wann_wigner_seitz(
300 : c$$$ > .false.,num,amat,rvecnum,
301 : c$$$ < int_dummy,rvec)
302 : c$$$
303 : c$$$ open(333,file='wig_vectors',recl=1000)
304 : c$$$ do ii=1,rvecnum
305 : c$$$ write(333,*)ii,rvec(1,ii),rvec(2,ii),rvec(3,ii)
306 : c$$$ enddo
307 : c$$$ close(333)
308 : c$$$
309 : c$$$ endif
310 :
311 0 : allocate(hreal(3,2,num_wann,num_wann,rvecnum))
312 0 : hreal=cmplx(0.0,0.0)
313 :
314 0 : do rvecind=1,rvecnum
315 0 : do k=1,nkpts
316 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
317 : & kpoints(2,k)*rvec(2,rvecind)+
318 0 : & kpoints(3,k)*rvec(3,rvecind) )
319 0 : fac=cmplx(cos(rdotk),-sin(rdotk))
320 0 : do m2=1,num_wann
321 0 : do m1=1,num_wann
322 0 : do spin=1,2
323 0 : do dir=1,3
324 : hreal(dir,spin,m1,m2,rvecind)=
325 : & hreal(dir,spin,m1,m2,rvecind)+
326 0 : & fac*hwann(dir,spin,m1,m2,k)
327 : enddo !dir
328 : enddo !spin
329 : enddo !m1
330 : enddo !m2
331 : enddo !k
332 : enddo !rvecind
333 0 : hreal=hreal/cmplx(real(nkpts),0.0)
334 :
335 0 : open(321,file='rspaulisvx'//spinspin12(jspin),form='formatted')
336 0 : do rvecind=1,rvecnum
337 0 : r3=rvec(3,rvecind)
338 0 : r2=rvec(2,rvecind)
339 0 : r1=rvec(1,rvecind)
340 0 : do j=1,num_wann
341 0 : do i=1,num_wann
342 0 : do kk=1,3
343 : write(321,'(i3,1x,i3,1x,i3,1x,i3,
344 : & 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
345 0 : & r1,r2,r3,i,j,kk,hreal(kk,1,i,j,rvecind)
346 : enddo !kk
347 : enddo !i
348 : enddo !j
349 : enddo !rvecnum
350 0 : close(321)
351 :
352 0 : open(321,file='rspaulisvy'//spinspin12(jspin),form='formatted')
353 0 : do rvecind=1,rvecnum
354 0 : r3=rvec(3,rvecind)
355 0 : r2=rvec(2,rvecind)
356 0 : r1=rvec(1,rvecind)
357 0 : do j=1,num_wann
358 0 : do i=1,num_wann
359 0 : do kk=1,3
360 : write(321,'(i3,1x,i3,1x,i3,1x,i3,
361 : & 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
362 0 : & r1,r2,r3,i,j,kk,hreal(kk,2,i,j,rvecind)
363 : enddo !kk
364 : enddo !i
365 : enddo !j
366 : enddo !rvecnum
367 0 : close(321)
368 :
369 :
370 0 : deallocate(lwindow,u_matrix_opt,ndimwin)
371 0 : deallocate(u_matrix,hwann,hreal)
372 : enddo !jspin
373 :
374 0 : end subroutine wann_nabla_pauli_rs
375 : end module m_wann_nabla_pauli_rs
|
