Line data Source code
1 : c*************************************************
2 : c Calculate the matrix elements of the
3 : c Pauli matrix in real space from the
4 : c files WF1.chk (and WF1_um.dat) produced
5 : c by wannier90.
6 : c FF, January 2009
7 : c*************************************************
8 : module m_wann_pauliat_rs
9 : USE m_juDFT
10 : contains
11 0 : subroutine wann_pauliat_rs(
12 0 : > rvecnum,rvec,kpoints,
13 : > jspins_in,nkpts,l_bzsym,film,
14 : > l_soc,band_min,band_max,neigd,
15 0 : > l_socmmn0,ntype,neq,l_ndegen,ndegen,
16 : > wan90version,l_unformatted)
17 :
18 : use m_constants, only:pimach
19 : use m_wann_read_umatrix
20 :
21 : implicit none
22 : integer, intent(in) :: rvecnum
23 : integer, intent(in) :: rvec(:,:)
24 : real, intent(in) :: kpoints(:,:)
25 : integer, intent(in) :: jspins_in
26 : integer, intent(in) :: nkpts
27 : logical,intent (in) :: l_bzsym,l_soc
28 : logical,intent(in) :: film
29 : integer,intent(in) :: band_min(2),band_max(2),neigd
30 : logical, intent(in) :: l_socmmn0
31 : integer, intent(in) :: ntype
32 : integer, intent(in) :: neq(:)
33 : logical, intent(in) :: l_ndegen
34 : integer, intent(in) :: ndegen(:)
35 : integer, intent(in) :: wan90version
36 : logical, intent(in) :: l_unformatted
37 :
38 : integer :: ikpt,jspins
39 : integer :: kpts
40 : logical :: l_file
41 : integer :: num_wann,num_kpts,num_nnmax,jspin
42 : integer :: kspin,kkspin
43 : integer :: wann_shift,num_wann2
44 : integer :: i,j,k,m,info,r1,r2,r3,dummy1
45 : integer :: dummy2,dummy3
46 : integer :: hopmin,hopmax,counter,m1,m2
47 : integer :: num_bands2
48 : integer,allocatable :: iwork(:)
49 : real,allocatable :: energy(:,:),ei(:)
50 : real,allocatable :: eigw(:,:),rwork(:)
51 : complex,allocatable :: work(:),vec(:,:)
52 0 : complex,allocatable :: u_matrix(:,:,:,:),hwann(:,:,:,:)
53 0 : complex,allocatable :: hreal(:,:,:,:)
54 : complex,allocatable :: hrealsoc(:,:,:,:,:,:,:)
55 : complex,allocatable :: hwannsoc(:,:,:,:,:)
56 0 : complex,allocatable :: paulimat(:,:,:,:)
57 0 : complex,allocatable :: paulimat_opt(:,:,:,:)
58 0 : complex,allocatable :: paulimat2(:,:,:,:)
59 : complex :: fac,eulav,eulav1
60 : real :: tmp_omi,rdotk,tpi,minenerg,maxenerg
61 : real, allocatable :: minieni(:),maxieni(:)
62 : character :: jobz,uplo
63 : integer :: kpt,band,lee,lwork,lrwork,liwork,n,lda
64 : complex :: value(4)
65 : logical :: um_format
66 : logical :: repro_eig
67 : logical :: l_chk,l_proj
68 : logical :: have_disentangled
69 0 : integer,allocatable :: ndimwin(:,:)
70 0 : logical,allocatable :: lwindow(:,:,:)
71 : integer :: chk_unit,nkp,ntmp,ierr
72 : character(len=33) :: header
73 : character(len=20) :: checkpoint
74 : real :: tmp_latt(3,3), tmp_kpt_latt(3,nkpts)
75 : real :: omega_invariant
76 0 : complex,allocatable :: u_matrix_opt(:,:,:,:)
77 : integer :: num_bands,counter1,counter2
78 : logical :: l_umdat
79 : real,allocatable :: eigval2(:,:)
80 : real,allocatable :: eigval_opt(:,:)
81 : real :: scale,a,b
82 : character(len=2) :: spinspin12(0:2)
83 : character(len=3) :: spin12(2)
84 : character(len=6) :: filename
85 : integer :: jp,mp,kk
86 : complex,parameter :: ci=(0.0,1.0)
87 : integer :: dir,rvecind
88 : integer :: spin1,spin2
89 : character(len=10) :: torquename
90 : character(len=9) :: torquename2
91 : character(len=15) :: torquename2ndegen
92 : integer :: at,itype,nn
93 : integer :: nbnd,fullnkpts,nwfs
94 0 : complex,allocatable :: amn(:,:,:)
95 0 : complex,allocatable :: perpmag(:,:,:,:)
96 0 : complex,allocatable :: hwannmix(:,:,:,:)
97 0 : complex,allocatable :: paulimatmix(:,:,:,:)
98 :
99 : data spinspin12/' ','.1' , '.2'/
100 : data spin12/'WF1','WF2'/
101 :
102 0 : call timestart("wann_pauliat_rs")
103 0 : tpi=2*pimach()
104 :
105 0 : jspins=jspins_in
106 0 : if(l_soc)jspins=1
107 :
108 0 : write(6,*)"nkpts=",nkpts
109 :
110 : c*****************************************************
111 : c get num_bands and num_wann from the proj file
112 : c*****************************************************
113 0 : do j=jspins,0,-1
114 0 : inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
115 0 : if(l_file)then
116 0 : filename='proj'//spinspin12(j)
117 0 : exit
118 : endif
119 : enddo
120 0 : if(l_file)then
121 0 : open (203,file=trim(filename),status='old')
122 0 : rewind (203)
123 : else
124 : CALL judft_error("no proj/proj.1/proj.2",calledby
125 0 : + ="wann_pauli_rs")
126 : endif
127 0 : read (203,*) num_wann,num_bands
128 0 : close (203)
129 0 : write(6,*)'According to proj there are ',num_bands,' bands'
130 0 : write(6,*)"and ",num_wann," wannier functions."
131 :
132 0 : allocate( u_matrix_opt(num_bands,num_wann,nkpts,jspins) )
133 0 : allocate( u_matrix(num_wann,num_wann,nkpts,jspins) )
134 0 : allocate( lwindow(num_bands,nkpts,jspins) )
135 0 : allocate( ndimwin(nkpts,jspins) )
136 :
137 0 : do jspin=1,jspins !spin loop
138 :
139 :
140 : c****************************************************************
141 : c read in chk
142 : c****************************************************************
143 0 : num_kpts=nkpts
144 :
145 : call wann_read_umatrix2(
146 : > nkpts,num_wann,num_bands,
147 : > um_format,jspin,wan90version,
148 : < have_disentangled,
149 : < lwindow(:,:,jspin),
150 : < ndimwin(:,jspin),
151 : < u_matrix_opt(:,:,:,jspin),
152 0 : < u_matrix(:,:,:,jspin))
153 :
154 :
155 : enddo
156 :
157 0 : jspin=1
158 :
159 0 : num_bands2=jspins*num_bands
160 0 : num_wann2=jspins*num_wann
161 0 : write(*,*)"l_unformatted=",l_unformatted
162 : c****************************************************
163 : c Read the file "WF1.socmmn0".
164 : c****************************************************
165 0 : allocate( paulimat(3,num_bands2,num_bands2,nkpts) )
166 0 : at=0
167 0 : do itype=1,ntype
168 0 : do nn=1,neq(itype)
169 0 : at=at+1
170 0 : write(*,*)"atom=",at
171 : write(torquename,
172 0 : & fmt='("mmn0up_",i3.3)')at
173 0 : if(l_unformatted)then
174 :
175 0 : paulimat=0.0
176 : c$$$ if(l_soc)then
177 : c$$$ open(304,file='WF1.socmmn0',form='formatted')
178 : c$$$ else
179 : c$$$ open(304,file='WF1.mmn0',form='formatted')
180 : c$$$ endif
181 : c$$$ read(304,*)
182 : c$$$ read(304,*)
183 : c$$$ do nkp=1,num_kpts
184 : c$$$ do i=1,num_bands
185 : c$$$ do j=1,num_bands
186 : c$$$ read(304,*)dummy1,dummy2,dummy3,a,b
187 : c$$$ paulimat(3,i,j,nkp)=cmplx(a,b)
188 : c$$$ enddo !j
189 : c$$$ enddo !i
190 : c$$$ enddo !nkp
191 : c$$$ close(304)
192 :
193 : c****************************************************
194 : c Read the file "WF2.socmmn0".
195 : c****************************************************
196 : c$$$ if(l_soc)then
197 : c$$$ open(304,file='WF2.socmmn0',form='formatted')
198 : c$$$ else
199 : c$$$ open(304,file='WF2.mmn0',form='formatted')
200 : c$$$ endif
201 : c$$$ read(304,*)
202 : c$$$ read(304,*)
203 : c$$$ do nkp=1,num_kpts
204 : c$$$ do i=1+(jspins-1)*num_bands,num_bands2
205 : c$$$ do j=1+(jspins-1)*num_bands,num_bands2
206 : c$$$ read(304,*)dummy1,dummy2,dummy3,a,b
207 : c$$$ paulimat(3,i,j,nkp)=
208 : c$$$ & paulimat(3,i,j,nkp)-cmplx(a,b)
209 : c$$$ enddo !j
210 : c$$$ enddo !i
211 : c$$$ enddo !nkp
212 : c$$$ close(304)
213 :
214 : c****************************************************
215 : c Read the file "updown.mmn0".
216 : c****************************************************
217 0 : inquire(file=trim(torquename)//'_unf',exist=l_file)
218 0 : if(.not.l_file)goto 12345
219 : open(304,file=trim(torquename)//'_unf',
220 0 : & form='unformatted')
221 0 : read(304)nbnd,fullnkpts,nwfs
222 0 : if(fullnkpts.ne.num_kpts)stop'fullnkpts.ne.num_kpts'
223 0 : if(nbnd.ne.num_bands)stop'nbnd.ne.num_bands'
224 0 : if(nwfs.ne.num_bands)stop'nwfs.ne.num_bands'
225 0 : allocate(amn(nbnd,nwfs,fullnkpts))
226 0 : read(304)amn
227 0 : close(304)
228 0 : do nkp=1,num_kpts
229 0 : do i=1,num_bands
230 0 : do j=1,num_bands
231 0 : paulimat(1,j,i+num_bands*(jspins-1),nkp)=amn(j,i,nkp)
232 : enddo !j
233 : enddo !i
234 0 : paulimat(2,:,:,nkp)=ci*paulimat(1,:,:,nkp)
235 : paulimat(1,:,:,nkp)=paulimat(1,:,:,nkp)+
236 0 : & transpose(conjg( paulimat(1,:,:,nkp) ))
237 : paulimat(2,:,:,nkp)=paulimat(2,:,:,nkp)+
238 0 : & transpose(conjg( paulimat(2,:,:,nkp) ))
239 : enddo !nkp
240 0 : deallocate(amn)
241 : else!unformatted
242 0 : stop 'not yet finished'
243 : endif
244 :
245 : c****************************************************************
246 : c Calculate matrix elements of Pauli in the basis of
247 : c rotated Bloch functions.
248 : c****************************************************************
249 0 : allocate( paulimat2(3,num_wann2,num_wann2,nkpts) )
250 : write(6,*)"calculate matrix elements of spin operator
251 0 : & between wannier orbitals"
252 :
253 0 : if(have_disentangled) then
254 :
255 0 : allocate( paulimat_opt(3,num_bands2,num_bands2,nkpts) )
256 0 : allocate( paulimatmix(3,num_wann2,num_bands2,nkpts))
257 :
258 0 : do nkp=1,num_kpts
259 : counter1=0
260 0 : do m=1,num_bands2
261 0 : spin1=(m-1)/num_bands
262 0 : if(lwindow(m-spin1*num_bands,nkp,spin1+1))then
263 0 : counter1=counter1+1
264 0 : counter2=0
265 0 : do mp=1,num_bands2
266 0 : spin2=(mp-1)/num_bands
267 0 : if(lwindow(mp-spin2*num_bands,nkp,spin2+1))then
268 0 : counter2=counter2+1
269 0 : do k=1,3
270 : paulimat_opt(k,counter2,counter1,nkp)=
271 0 : & paulimat(k,mp,m,nkp)
272 : enddo
273 : endif
274 : enddo !mp
275 : endif
276 : enddo !m
277 : enddo
278 :
279 0 : paulimatmix=0.0
280 0 : do spin1=0,jspins-1
281 0 : do spin2=0,jspins-1
282 0 : do nkp=1,num_kpts
283 0 : do jp=1,num_wann
284 0 : do m=1,ndimwin(nkp,spin2+1)
285 0 : do mp=1,ndimwin(nkp,spin1+1)
286 0 : do k=1,3
287 :
288 : paulimatmix(k,jp+spin1*num_wann,m+spin2*ndimwin(nkp,1),nkp)=
289 : = paulimatmix(k,jp+spin1*num_wann,m+spin2*ndimwin(nkp,1),nkp)+
290 : & conjg(u_matrix_opt(mp,jp,nkp,spin1+1))*
291 : & paulimat_opt(k,mp+spin1*ndimwin(nkp,1),
292 0 : & m+spin2*ndimwin(nkp,1),nkp)
293 : enddo !k
294 : enddo !mp
295 : enddo !m
296 : enddo !jp
297 : enddo !nkp
298 : enddo !spin2
299 : enddo !spin1
300 :
301 0 : paulimat2=0.0
302 0 : do spin1=0,jspins-1
303 0 : do spin2=0,jspins-1
304 0 : do nkp=1,num_kpts
305 0 : do j=1,num_wann
306 0 : do jp=1,num_wann
307 0 : do m=1,ndimwin(nkp,spin2+1)
308 0 : do k=1,3
309 :
310 :
311 : paulimat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)=
312 : = paulimat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)+
313 : & paulimatmix(k,jp+spin1*num_wann,
314 : & m+spin2*ndimwin(nkp,1),nkp)*
315 0 : & u_matrix_opt(m,j,nkp,spin2+1)
316 : enddo !k
317 : enddo !m
318 : enddo !jp
319 : enddo !j
320 : enddo !nkp
321 : enddo !spin2
322 : enddo !spin1
323 :
324 :
325 0 : deallocate(paulimat_opt)
326 0 : deallocate(paulimatmix)
327 : else
328 0 : paulimat2=paulimat
329 : end if !have_disentangled
330 :
331 0 : allocate(hwann(3,num_wann2,num_wann2,num_kpts))
332 0 : hwann=cmplx(0.0,0.0)
333 0 : wann_shift=0
334 :
335 0 : allocate(hwannmix(3,num_wann2,num_wann2,num_kpts))
336 0 : hwannmix=cmplx(0.0,0.0)
337 :
338 0 : do spin1=0,jspins-1
339 0 : do spin2=0,jspins-1
340 0 : do k=1,num_kpts
341 0 : do i=1,num_wann
342 0 : do mp=1,num_wann
343 0 : do j=1,num_wann
344 0 : do kk=1,3
345 : hwannmix(kk,mp+spin1*num_wann,i+spin2*num_wann,k)=
346 : = hwannmix(kk,mp+spin1*num_wann,i+spin2*num_wann,k)+
347 : * conjg(u_matrix(j,mp,k,spin1+1))*
348 : * paulimat2(kk,j+spin1*num_wann,
349 0 : * i+spin2*num_wann,k)
350 : enddo !kk
351 : enddo !j
352 : enddo !mp
353 : enddo !i
354 : enddo !k
355 : enddo !spin2
356 : enddo !spin1
357 :
358 0 : do spin1=0,jspins-1
359 0 : do spin2=0,jspins-1
360 0 : do k=1,num_kpts
361 0 : do m=1,num_wann
362 0 : do i=1,num_wann
363 0 : do mp=1,num_wann
364 0 : do kk=1,3
365 : hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)=
366 : = hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)+
367 : * hwannmix(kk,mp+spin1*num_wann,
368 : * i+spin2*num_wann,k)*
369 0 : * u_matrix(i,m,k,spin2+1)
370 : enddo !kk
371 : enddo !mp
372 : enddo !i
373 : enddo !m
374 : enddo !k
375 : enddo !spin2
376 : enddo !spin1
377 :
378 0 : deallocate(hwannmix)
379 :
380 : c************************************************************
381 : c Calculate matrix elements in real space.
382 : c***********************************************************
383 0 : write(6,*)"calculate pauli-mat in rs"
384 0 : allocate(hreal(3,num_wann2,num_wann2,rvecnum))
385 0 : hreal = cmplx(0.0,0.0)
386 :
387 0 : if(l_ndegen)then
388 0 : do rvecind=1,rvecnum
389 0 : do k=1,nkpts
390 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
391 : & kpoints(2,k)*rvec(2,rvecind)+
392 0 : & kpoints(3,k)*rvec(3,rvecind) )
393 0 : fac=cmplx(cos(rdotk),-sin(rdotk))/real(ndegen(rvecind))
394 0 : do m2=1,num_wann2
395 0 : do m1=1,num_wann2
396 0 : do dir=1,3
397 : hreal(dir,m1,m2,rvecind)=hreal(dir,m1,m2,rvecind)+
398 0 : & fac*hwann(dir,m1,m2,k)
399 : enddo !dir
400 : enddo !m1
401 : enddo !m2
402 : enddo !k
403 : enddo !rvecind
404 : else !l_ndegen
405 0 : do rvecind=1,rvecnum
406 0 : do k=1,nkpts
407 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
408 : & kpoints(2,k)*rvec(2,rvecind)+
409 0 : & kpoints(3,k)*rvec(3,rvecind) )
410 0 : fac=cmplx(cos(rdotk),-sin(rdotk))
411 0 : do m2=1,num_wann2
412 0 : do m1=1,num_wann2
413 0 : do dir=1,3
414 : hreal(dir,m1,m2,rvecind)=
415 : & hreal(dir,m1,m2,rvecind)+
416 0 : & fac*hwann(dir,m1,m2,k)
417 : enddo !dir
418 : enddo !m1
419 : enddo !m2
420 : enddo !k
421 : enddo !rvecind
422 : endif
423 0 : hreal=hreal/cmplx(real(nkpts),0.0)
424 :
425 : if(l_unformatted)then
426 0 : allocate(perpmag(num_wann2,num_wann2,3,rvecnum))
427 0 : if(l_ndegen)then
428 : write(torquename2ndegen,
429 0 : & fmt='("rspaundegen_",i3.3)')at
430 : open(321,file=trim(torquename2ndegen)//'_unf',
431 : & form='unformatted',
432 0 : & convert='BIG_ENDIAN')
433 : else
434 : write(torquename2,
435 0 : & fmt='("rspau_",i3.3)')at
436 : open(321,file=trim(torquename2)//'_unf',
437 : & form='unformatted',
438 0 : & convert='BIG_ENDIAN')
439 : endif
440 0 : do rvecind=1,rvecnum
441 0 : do j=1,num_wann2
442 0 : do i=1,num_wann2
443 0 : do kk=1,3
444 0 : perpmag(i,j,kk,rvecind)=hreal(kk,i,j,rvecind)
445 : enddo !kk
446 : enddo !i
447 : enddo !j
448 : enddo !rvecind
449 0 : write(321)perpmag
450 0 : deallocate(perpmag)
451 : else !l_unformatted
452 : stop 'still need new file names'
453 : open(321,file='rspauli'//spinspin12(1),form='formatted')
454 : do rvecind=1,rvecnum
455 : r3=rvec(3,rvecind)
456 : r2=rvec(2,rvecind)
457 : r1=rvec(1,rvecind)
458 : do j=1,num_wann2
459 : do i=1,num_wann2
460 : do kk=1,3
461 : write(321,'(i3,1x,i3,1x,i3,1x,i3,
462 : & 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
463 : & r1,r2,r3,i,j,kk,hreal(kk,i,j,rvecind)
464 : enddo !kk
465 : enddo !i
466 : enddo !j
467 : enddo !rvecnum
468 : close(321)
469 : endif !l_unformatted
470 :
471 0 : deallocate(paulimat2)
472 0 : deallocate(hwann,hreal)
473 :
474 0 : 12345 continue
475 :
476 : enddo !nn
477 : enddo !itype
478 :
479 0 : deallocate(lwindow,u_matrix_opt,ndimwin)
480 0 : deallocate(u_matrix)
481 :
482 0 : call timestop("wann_pauliat_rs")
483 0 : end subroutine wann_pauliat_rs
484 : end module m_wann_pauliat_rs
|
