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