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_hopping
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_hopping(
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 hoppings using the information
18 : c in file WF1.chk produced by wannier90.
19 : c
20 : c Frank Freimuth
21 : c****************************************************
22 : use m_constants
23 : use m_wann_read_umatrix
24 : use m_juDFT
25 : c use m_wann_get_mp
26 : c use m_wann_get_kpts
27 : c use m_wann_wigner_seitz
28 :
29 : implicit none
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,l_soc
36 : logical, intent(in) :: film
37 : integer, intent(in) :: band_min(2),band_max(2),neigd
38 : logical, intent(in) :: l_socmmn0
39 : logical, intent(in) :: l_ndegen
40 : integer, intent(in) :: ndegen(:)
41 : integer, intent(in) :: wan90version
42 : logical, intent(in) :: l_unformatted
43 :
44 : integer :: ikpt,jspins
45 : integer :: kpts
46 : logical :: l_file
47 : c real :: kpoints(3,nkpts)
48 : integer :: num_wann,num_kpts,num_nnmax,jspin
49 : integer :: kspin,kkspin
50 : integer :: wann_shift,num_wann2
51 : integer :: i,j,k,m,info,r1,r2,r3,dummy1
52 : integer :: dummy2,dummy3
53 : integer :: counter,m1,m2
54 : integer :: num_bands2
55 0 : integer,allocatable :: iwork(:)
56 0 : real,allocatable :: energy(:,:),ei(:)
57 0 : real,allocatable :: eigw(:,:),rwork(:)
58 0 : complex,allocatable :: work(:),vec(:,:)
59 0 : complex,allocatable :: u_matrix(:,:,:),hwann(:,:,:)
60 0 : complex,allocatable :: hreal(:,:,:)
61 0 : complex,allocatable :: hrealsoc(:,:,:,:,:)
62 0 : complex,allocatable :: hwannsoc(:,:,:,:,:)
63 0 : complex,allocatable :: mmn0(:,:,:,:)
64 : complex :: fac,eulav,eulav1
65 : real :: tmp_omi,rdotk,tpi,minenerg,maxenerg
66 : real, allocatable :: minieni(:),maxieni(:)
67 : character :: jobz,uplo
68 : integer :: kpt,band,lee,lwork,lrwork,liwork,n,lda
69 : complex :: value(4)
70 : logical :: um_format
71 : logical :: repro_eig
72 : logical :: l_chk,l_proj
73 : logical :: have_disentangled
74 0 : integer,allocatable :: ndimwin(:)
75 0 : logical,allocatable :: lwindow(:,:)
76 : integer :: chk_unit,nkp,ntmp,ierr
77 : character(len=33) :: header
78 : character(len=20) :: checkpoint
79 : real :: tmp_latt(3,3), tmp_kpt_latt(3,nkpts)
80 : real :: omega_invariant
81 0 : complex,allocatable :: u_matrix_opt(:,:,:)
82 : integer :: num_bands,num(3)
83 : logical :: l_umdat
84 0 : real,allocatable :: eigval2(:,:)
85 0 : real,allocatable :: eigval_opt(:,:)
86 : real :: scale,a,b
87 : character(len=2) :: spinspin12(0:2)
88 : character(len=3) :: spin12(2)
89 : character(len=6) :: filename
90 : logical :: l_socham
91 : ! integer :: hopmin_x,hopmin_y,hopmin_z
92 : ! integer :: hopmax_x,hopmax_y,hopmax_z
93 : integer :: rvecind!,rvecnum
94 : ! integer,allocatable :: rvec(:,:)
95 : integer :: ii,int_dummy
96 :
97 : data spinspin12/' ','.1' , '.2'/
98 : data spin12/'WF1','WF2'/
99 :
100 0 : call timestart("wann_hopping")
101 0 : um_format=.false. !if you would like to get a formatted um_dat
102 0 : repro_eig=.false. !if you would like to check unitarity of u_matrix
103 :
104 0 : tpi=2*pimach()
105 :
106 0 : jspins=jspins_in
107 0 : if(l_soc)jspins=1
108 0 : l_socham=.false.
109 :
110 0 : write(oUnit,*)"nkpts=",nkpts
111 :
112 0 : do jspin=1,jspins !spin loop
113 : c*****************************************************
114 : c get num_bands and num_wann from the proj file
115 : c*****************************************************
116 0 : do j=jspin,0,-1
117 0 : inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
118 0 : if(l_file)then
119 0 : filename='proj'//spinspin12(j)
120 0 : exit
121 : endif
122 : enddo
123 0 : if(l_file)then
124 0 : open (203,file=trim(filename),status='old')
125 0 : rewind (203)
126 : else
127 : CALL juDFT_error("no proj/proj.1/proj.2",calledby
128 0 : + ="wann_hopping")
129 : endif
130 0 : read (203,*) num_wann,num_bands
131 0 : close (203)
132 0 : write(oUnit,*)'According to proj there are ',
133 0 : + num_bands,' bands'
134 0 : write(oUnit,*)"and ",num_wann," wannier functions."
135 : c****************************************************************
136 : c read in chk
137 : c****************************************************************
138 0 : num_kpts=nkpts
139 0 : allocate( u_matrix_opt(num_bands,num_wann,nkpts) )
140 0 : allocate( u_matrix(num_wann,num_wann,nkpts) )
141 0 : allocate( lwindow(num_bands,nkpts) )
142 0 : allocate( ndimwin(nkpts) )
143 : call wann_read_umatrix2(
144 : > nkpts,num_wann,num_bands,
145 : > um_format,jspin,wan90version,
146 : < have_disentangled,
147 : < lwindow,ndimwin,
148 0 : < u_matrix_opt,u_matrix)
149 : c****************************************************************
150 : c read in eig-file
151 : c****************************************************************
152 0 : num_bands2=num_bands
153 0 : if(l_soc.and.l_socmmn0)then
154 0 : num_bands2=neigd
155 : endif
156 0 : write(oUnit,*)"read in eig-file"
157 0 : allocate(energy(num_bands2,num_kpts))
158 0 : inquire(file=spin12(jspin)//'.eig',exist=l_umdat)
159 0 : IF(.NOT.l_umdat) CALL juDFT_error
160 : + ("Thou shall not hide your eig file",calledby
161 0 : + ="wann_hopping")
162 0 : open(300,file=spin12(jspin)//'.eig',form='formatted')
163 0 : do i=1,num_kpts
164 0 : do j=1,num_bands2
165 0 : read(300,*)band,kpt,energy(j,i)
166 : enddo
167 : enddo
168 0 : close(300)
169 :
170 0 : minenerg=minval(energy(:,:))
171 0 : maxenerg=maxval(energy(:,:))
172 0 : write(oUnit,*)"minenerg=",minenerg
173 0 : write(oUnit,*)"maxenerg=",maxenerg
174 : c*********************************************
175 : c Preparations for spin-orbit coupling
176 : c*********************************************
177 : if(l_socham)then
178 : allocate(mmn0(nkpts,num_bands2,num_bands2,2))
179 : open(304,file='spinup.mmn0',form='formatted')
180 : read(304,*)
181 : read(304,*)
182 : do nkp=1,num_kpts
183 : do i=1,num_bands2
184 : do j=1,num_bands2
185 : read(304,*)dummy1,dummy2,dummy3,a,b
186 : c mmn0(nkp,j,i,1)=cmplx(a,b)
187 : mmn0(nkp,i,j,1)=cmplx(a,-b)
188 : enddo !j
189 : enddo !i
190 : enddo !nkp
191 : close(304)
192 : open(304,file='spindown.mmn0',form='formatted')
193 : read(304,*)
194 : read(304,*)
195 : do nkp=1,num_kpts
196 : do i=1,num_bands2
197 : do j=1,num_bands2
198 : read(304,*)dummy1,dummy2,dummy3,a,b
199 : c mmn0(nkp,j,i,2)=cmplx(a,b)
200 : mmn0(nkp,i,j,2)=cmplx(a,-b)
201 : enddo !j
202 : enddo !i
203 : enddo !nkp
204 : close(304)
205 : endif!l_soc
206 : c****************************************************************
207 : c calculate matrix elements of hamiltonian
208 : c****************************************************************
209 :
210 : write(oUnit,*)"calculate matrix elements of hamiltonian
211 0 : & between wannier orbitals"
212 :
213 :
214 0 : allocate(eigval_opt(num_bands,nkpts))
215 0 : allocate(eigval2(num_bands2,nkpts))
216 0 : eigval_opt=0.0
217 0 : eigval2=0.0
218 :
219 :
220 0 : if(have_disentangled) then
221 :
222 : ! slim down eigval to contain states within the outer window
223 :
224 0 : do nkp=1,num_kpts
225 0 : counter=0
226 0 : do j=1,num_bands
227 0 : if(lwindow(j,nkp)) then
228 0 : counter=counter+1
229 0 : eigval_opt(counter,nkp)=energy(j,nkp)
230 : end if
231 : end do
232 : end do
233 :
234 : ! rotate eigval into the optimal subspace
235 : ! in general eigval would be a matrix at each kpoints
236 : ! but we choose u_matrix_opt such that the Hamiltonian is
237 : ! diagonal at each kpoint. (I guess we should check it here)
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(1:num_bands2,:)=energy(1:num_bands2,:)
250 : end if !have_disentangled
251 :
252 0 : deallocate(eigval_opt)
253 0 : deallocate(energy)
254 :
255 0 : allocate(hwann(num_wann,num_wann,num_kpts))
256 0 : hwann=cmplx(0.0,0.0)
257 0 : wann_shift=0
258 0 : if(l_socmmn0)then
259 0 : wann_shift=band_min(jspin)-1
260 : endif
261 :
262 0 : do k=1,num_kpts
263 0 : do i=1,num_wann
264 0 : do j=1,num_wann
265 0 : do m=1,num_wann
266 : hwann(i,j,k)=hwann(i,j,k)
267 : + +eigval2(m+wann_shift,k)*u_matrix(m,j,k)*
268 0 : * conjg(u_matrix(m,i,k))
269 : enddo
270 : enddo
271 : enddo
272 : enddo
273 :
274 : if(l_socham)then
275 :
276 : num_wann2=num_wann
277 : wann_shift=0
278 : if(l_socmmn0)then
279 : num_wann2=neigd
280 : wann_shift=band_min(jspin)-1
281 : endif
282 : write(oUnit,*)"num_wann2=",num_wann2
283 : write(oUnit,*)"num_wann=",num_wann
284 : write(oUnit,*)"wann_shift=",wann_shift
285 : allocate(hwannsoc(num_wann,num_wann,num_kpts,2,2))
286 : hwannsoc=cmplx(0.0,0.0)
287 : do kspin=1,2
288 : do kkspin=1,2
289 : do i=1,num_wann
290 : do j=1,num_wann
291 : do k=1,num_kpts
292 : do m=1,num_wann2
293 : do m1=1,num_wann
294 : do m2=1,num_wann
295 : hwannsoc(i,j,k,kkspin,kspin)=
296 : = hwannsoc(i,j,k,kkspin,kspin)
297 : + +eigval2(m,k)*u_matrix(m2,j,k)*conjg(u_matrix(m1,i,k))*
298 : * mmn0(k,m1+wann_shift,m,kkspin)*
299 : * mmn0(k,m,m2+wann_shift,kspin)
300 : enddo
301 : enddo
302 : enddo !m
303 : enddo !k
304 : enddo !j
305 : enddo !i
306 : enddo !kkspin
307 : enddo !kspin
308 : deallocate(mmn0)
309 :
310 :
311 : do i=1,num_wann
312 : do j=1,num_wann
313 : do k=1,num_kpts
314 : eulav=hwannsoc(i,j,k,1,1)+hwannsoc(i,j,k,1,2)+
315 : & hwannsoc(i,j,k,2,1)+hwannsoc(i,j,k,2,2)
316 : eulav=eulav-hwann(i,j,k)
317 : if(abs(eulav).gt.1.e-6)then
318 : write(oUnit,*)"soc-hop: something wrong:",eulav
319 : endif
320 : enddo
321 : enddo
322 : enddo
323 :
324 :
325 :
326 : endif !l_soc
327 :
328 :
329 :
330 :
331 : c***************************************************************
332 : c repro_eig
333 : c***************************************************************
334 :
335 : if(repro_eig)then
336 :
337 : write(oUnit,*)
338 : + "As a check, try to reproduce the old eigenvalues"
339 :
340 : c Note, that this check will give positive result,
341 : c if u_matrix is unitary
342 :
343 : allocate (vec(num_wann,num_wann),ei(num_wann))
344 : allocate(eigw(num_kpts,num_wann))
345 :
346 : lee = log( dble(num_wann) )/log(2.d0) + 1
347 : lwork = 1 + 5*num_wann + 2*num_wann*lee + 3*(num_wann**2)
348 : lrwork = 1 + 4*num_wann + 2*num_wann*lee + 3*(num_wann**2)
349 : liwork = 2 + 5*num_wann +1
350 :
351 : allocate (work(lwork),rwork(lrwork),iwork(liwork))
352 :
353 :
354 :
355 : do i=1,num_kpts
356 : do j = 1,num_wann
357 : do k = 1,num_wann
358 : vec(j,k) = hwann(j,k,i)
359 : enddo
360 : enddo
361 :
362 :
363 : jobz = 'V' ; uplo = 'L' ; n = num_wann ; lda = num_wann
364 : call zheevd(jobz,uplo,n,vec,lda,ei,work,lwork,
365 : & rwork,lrwork,iwork,liwork,info)
366 :
367 : if (info.lt.0) write (*,*)
368 : & 'ith argument had an illegal value ',info
369 : IF (info>0) CALL juDFT_error("not converged diagonalization"
370 : + ,calledby ="wann_hopping")
371 :
372 : do j = 1,num_wann
373 : eigw(i,j) = ei(j)
374 : if(abs(eigw(i,j)-eigval2(j,i)).gt.0.0001)then
375 : write(oUnit,*)"found different eigenvalues:"
376 : write(oUnit,*)"kpt=",i
377 : write(oUnit,*)"band=",j
378 : write(oUnit,*)"eig=",eigval2(j,i)
379 : write(oUnit,*)"neweig=",eigw(i,j)
380 : endif
381 : enddo
382 :
383 : enddo
384 : deallocate(work,rwork,iwork)
385 : deallocate(vec,ei)
386 : open(500,file='reeig'//spin12(jspin),form='formatted')
387 : do i=1,num_kpts
388 : do j=1,num_wann
389 : write(500,'(i5,3x,i5,3x,f20.16)')i,j,eigw(i,j)
390 : enddo
391 : enddo
392 : close(500)
393 : deallocate(eigw)
394 : endif !repro_eig
395 : c************************************************************
396 : c Calculate hoppings.
397 : c***********************************************************
398 0 : write(oUnit,*)"calculate hoppings"
399 0 : allocate( hreal(num_wann,num_wann,rvecnum) )
400 0 : hreal=cmplx(0.0,0.0)
401 0 : do rvecind=1,rvecnum
402 0 : do k=1,nkpts
403 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
404 : + kpoints(2,k)*rvec(2,rvecind)+
405 0 : + kpoints(3,k)*rvec(3,rvecind) )
406 0 : fac=cmplx(cos(rdotk),-sin(rdotk))
407 0 : do m2=1,num_wann
408 0 : do m1=1,num_wann
409 : hreal(m1,m2,rvecind)=hreal(m1,m2,rvecind)+
410 0 : & fac*hwann(m1,m2,k)
411 : enddo !m1
412 : enddo !m2
413 : enddo !k
414 : enddo !rvecind
415 0 : hreal=hreal/cmplx(real(nkpts),0.0)
416 :
417 0 : if(l_ndegen)then
418 0 : do rvecind=1,rvecnum
419 0 : do m2=1,num_wann
420 0 : do m1=1,num_wann
421 : hreal(m1,m2,rvecind)=
422 0 : & hreal(m1,m2,rvecind)/cmplx(real(ndegen(rvecind)),0.0)
423 : enddo !m1
424 : enddo !m2
425 : enddo !rvecind
426 : endif !l_ndegen
427 :
428 : c************************************************************
429 : c calculate hoppings for soc
430 : c***********************************************************
431 : if(l_socham)then
432 : allocate(hrealsoc(num_wann,num_wann,rvecnum,
433 : & 2,2))
434 : hrealsoc=cmplx(0.0,0.0)
435 : do rvecind=1,rvecnum
436 : do k=1,nkpts
437 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
438 : + kpoints(2,k)*rvec(2,rvecind)+
439 : + kpoints(3,k)*rvec(3,rvecind) )
440 : fac=cmplx(cos(rdotk),-sin(rdotk))/nkpts
441 : do kspin=1,2
442 : do kkspin=1,2
443 : hrealsoc(:,:,rvecind,kspin,kkspin)=
444 : = hrealsoc(:,:,rvecind,kspin,kkspin)+
445 : + fac*hwannsoc(:,:,k,kspin,kkspin)
446 : enddo
447 : enddo
448 : enddo !k
449 : enddo !rvecnum
450 : endif !l_soc
451 :
452 : c****************************************************************
453 : c make the hoppings real
454 : c****************************************************************
455 : c$$$ if(.false.)then
456 : c$$$ do i=1,num_wann
457 : c$$$ do j=i+1,num_wann
458 : c$$$ eulav=cmplx(0.0,0.0)
459 : c$$$ do r1=hopmin,hopmax
460 : c$$$ do r2=hopmin,hopmax
461 : c$$$ do r3=hopmin,hopmax
462 : c$$$ eulav1=hreal(i,j,r1,r2,r3)
463 : c$$$ if(abs(eulav1).gt.abs(eulav))then
464 : c$$$ eulav=eulav1
465 : c$$$ endif
466 : c$$$ enddo !r1
467 : c$$$ enddo !r2
468 : c$$$ enddo !r3
469 : c$$$ if(abs(eulav).gt.1e-6)then
470 : c$$$ eulav=eulav/abs(eulav)
471 : c$$$ do r1=hopmin,hopmax
472 : c$$$ do r2=hopmin,hopmax
473 : c$$$ do r3=hopmin,hopmax
474 : c$$$ hreal(i,j,r1,r2,r3)=hreal(i,j,r1,r2,r3)/eulav
475 : c$$$ hreal(j,i,-r1,-r2,-r3)=conjg(hreal(i,j,r1,r2,r3))
476 : c$$$! hreal(j,i,r1,r2,r3)=hreal(j,i,r1,r2,r3)/conjg(eulav)
477 : c$$$ if(l_soc)then
478 : c$$$ hrealsoc(i,j,r1,r2,r3,:,:)=
479 : c$$$ & hrealsoc(i,j,r1,r2,r3,:,:)/eulav
480 : c$$$ hrealsoc(j,i,-r1,-r2,-r3,:,:)=
481 : c$$$ & conjg( hrealsoc(i,j,r1,r2,r3,:,:) )
482 : c$$$! hrealsoc(j,i,r1,r2,r3,:,:)=
483 : c$$$! & hrealsoc(j,i,r1,r2,r3,:,:)/conjg(eulav)
484 : c$$$ endif
485 : c$$$ enddo !r3
486 : c$$$ enddo !r2
487 : c$$$ enddo !r1
488 : c$$$ endif !eulav
489 : c$$$ enddo !j
490 : c$$$ enddo !i
491 : c$$$ endif
492 :
493 : if(l_socham)then
494 : do kspin=1,2
495 : do kkspin=1,2
496 : open(321,file='hopping'//spinspin12(kspin)//spinspin12(kkspin),
497 : & form='formatted')
498 : do rvecind=1,rvecnum
499 : r1=rvec(1,rvecind)
500 : r2=rvec(2,rvecind)
501 : r3=rvec(3,rvecind)
502 : do i=1,num_wann
503 : do j=1,num_wann
504 : write(321,'(i3,i3,i3,i3,i3,f20.8,f20.8)')
505 : & r1,r2,r3,i,j,hrealsoc(i,j,rvecind,kspin,kkspin)
506 : enddo !j
507 : enddo !i
508 : enddo !rvecind
509 : close(321)
510 : enddo !kkspin
511 : enddo !kspin
512 : deallocate(hrealsoc,hwannsoc)
513 : endif !l_soc
514 :
515 0 : if(l_unformatted)then
516 0 : if(l_soc)then
517 0 : if(l_ndegen)then
518 : open(321,file='hunfndegen',form='unformatted'
519 : #ifdef CPP_INTEL
520 : & ,convert='BIG_ENDIAN'
521 : #endif
522 0 : & )
523 : else
524 : open(321,file='hunf',form='unformatted'
525 : #ifdef CPP_INTEL
526 : & ,convert='BIG_ENDIAN'
527 : #endif
528 0 : & )
529 : endif
530 : else
531 0 : if(l_ndegen)then
532 : open(321,file='hunfndegen'//spinspin12(jspin),
533 : & form='unformatted'
534 : #ifdef CPP_INTEL
535 : & ,convert='BIG_ENDIAN'
536 : #endif
537 0 : & )
538 : else
539 : open(321,file='hunf'//spinspin12(jspin),
540 : & form='unformatted'
541 : #ifdef CPP_INTEL
542 : & ,convert='BIG_ENDIAN'
543 : #endif
544 0 : & )
545 : endif
546 : endif
547 0 : write(321)num_wann
548 0 : write(321)rvecnum
549 0 : write(321)rvec
550 0 : write(321)hreal
551 : else
552 0 : if(l_ndegen)then
553 : open(321,file='hopping_ndegen'//spinspin12(jspin),
554 0 : & form='formatted')
555 : else
556 : open(321,file='hopping'//spinspin12(jspin),
557 0 : & form='formatted')
558 : endif
559 0 : do rvecind=1,rvecnum
560 0 : r3=rvec(3,rvecind)
561 0 : r2=rvec(2,rvecind)
562 0 : r1=rvec(1,rvecind)
563 0 : do j=1,num_wann
564 0 : do i=1,num_wann
565 : write(321,'(i3,1x,i3,1x,i3,1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
566 0 : & r1,r2,r3,i,j,hreal(i,j,rvecind)
567 : enddo
568 : enddo
569 : enddo !rvecnum
570 : endif !l_unformatted
571 0 : close(321)
572 :
573 0 : deallocate(lwindow,u_matrix_opt,ndimwin)
574 0 : deallocate(eigval2,u_matrix,hwann,hreal)
575 : ! deallocate(rvec)
576 : enddo !jspin
577 :
578 :
579 0 : call timestop("wann_hopping")
580 :
581 0 : end subroutine wann_hopping
582 : end module m_wann_hopping
|
