LCOV - code coverage report
Current view: top level - wannier - wann_rmat.f (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 0 174 0.0 %
Date: 2024-03-28 04:22:06 Functions: 0 1 0.0 %

          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_rmat
       8             :       USE m_juDFT
       9             :       contains
      10           0 :       subroutine wann_rmat(
      11           0 :      >          bmat,amat, 
      12           0 :      >          rvecnum,rvec,kpoints,
      13             :      >          jspins_in,nkpts,l_bzsym,film,
      14             :      >          l_nocosoc,band_min,band_max,neigd,
      15             :      >          l_socmmn0,wan90version)
      16             : c***********************************************************
      17             : c     Calculate the matrix elements of the position operator
      18             : c     between wannier functions.
      19             : c
      20             : c     FF, 2009
      21             : c***********************************************************
      22             :       use m_constants
      23             :       use m_wann_read_umatrix
      24             : 
      25             :       implicit none
      26             : 
      27             :       real, intent(in) :: bmat(:,:)
      28             :       real, intent(in) :: amat(:,:)
      29             : 
      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
      36             :       logical, intent(in) :: film
      37             :    
      38             :       logical, intent(in) :: l_nocosoc
      39             :       integer, intent(in) :: band_min(2),band_max(2),neigd
      40             : 
      41             :       logical, intent(in) :: l_socmmn0
      42             :       integer, intent(in) :: wan90version
      43             : 
      44           0 :       real,allocatable    :: bpunkt(:,:)
      45             :       real                :: kdiff(3)
      46           0 :       real,allocatable    :: wb(:)
      47           0 :       integer,allocatable :: gb(:,:,:),bpt(:,:)
      48             :       integer             :: ikpt,jspins,ikpt_help
      49             :       integer             :: kpts,nntot,nn
      50             :       logical             :: l_file
      51             :       integer             :: num_wann,num_kpts,num_nnmax,jspin
      52             :       integer             :: kspin,kkspin
      53             :       integer             :: num_wann2
      54             :       integer             :: i,j,k,m,info,r1,r2,r3,dummy1
      55             :       integer             :: dummy2,dummy3,dummy4
      56             :       integer             :: hopmin,hopmax,counter,m1,m2
      57             :       integer,allocatable :: iwork(:)
      58           0 :       real,allocatable    :: energy(:,:),ei(:)
      59             :       real,allocatable    :: eigw(:,:),rwork(:)
      60             :       complex,allocatable :: work(:),vec(:,:)
      61           0 :       complex,allocatable :: u_matrix(:,:,:),m_matrix(:,:,:,:)
      62           0 :       complex,allocatable :: posop(:,:,:,:),posop2(:,:,:,:)
      63             :       complex             :: fac,eulav,eulav1
      64             :       real                :: tmp_omi,rdotk,tpi,minenerg,maxenerg
      65             :       real, allocatable   :: minieni(:),maxieni(:)
      66             :       character           :: jobz,uplo
      67             :       integer             :: kpt,band,lee,lwork,lrwork,liwork,n,lda
      68             :       real                :: rvalue
      69             :       logical             :: um_format
      70             :       logical             :: repro_eig
      71             :       logical             :: l_chk,l_proj
      72             :       logical             :: have_disentangled
      73           0 :       integer,allocatable :: ndimwin(:)
      74           0 :       logical,allocatable :: lwindow(:,:)
      75             :       integer             :: chk_unit,nkp,ntmp,ierr
      76             :       character(len=33)   :: header
      77             :       character(len=20)   :: checkpoint
      78             :       real                :: tmp_latt(3,3), tmp_kpt_latt(3,nkpts)
      79             :       real                :: omega_invariant
      80           0 :       complex,allocatable :: u_matrix_opt(:,:,:)
      81             :       integer             :: num_bands,nn2
      82             :       logical             :: l_umdat
      83           0 :       real,allocatable    :: eigval2(:,:)
      84           0 :       real,allocatable    :: eigval_opt(:,:)
      85             :       real                :: scale,a,b
      86             :       character(len=2)    :: spinspin12(0:2)
      87             :       character(len=3)    :: spin12(2)
      88             :       character(len=6)    :: filename
      89             :       integer             :: jp,mp,kk
      90             :       integer             :: rvecind,rvecind_0
      91             :       real                :: realvec(3)
      92             :       complex,allocatable :: mmnk(:,:,:,:)
      93             :       logical             :: l_worksout,l_she
      94             : 
      95             :       data spinspin12/'  ','.1' , '.2'/
      96             :       data spin12/'WF1','WF2'/
      97             : 
      98           0 :       call timestart("wann_rmat")
      99             : 
     100           0 :       tpi=2*pimach()
     101             : 
     102           0 :       jspins=jspins_in
     103           0 :       if(l_nocosoc)jspins=1
     104             : 
     105           0 :       write(oUnit,*)"nkpts=",nkpts
     106             : 
     107           0 :       do jspin=1,jspins  !spin loop
     108             : c*****************************************************
     109             : c        get num_bands and num_wann from the proj file
     110             : c*****************************************************
     111           0 :          do j=jspin,0,-1
     112           0 :           inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
     113           0 :           if(l_file)then
     114           0 :             filename='proj'//spinspin12(j)
     115           0 :             exit
     116             :           endif
     117             :          enddo
     118           0 :          if(l_file)then
     119           0 :           open (203,file=trim(filename),status='old')
     120           0 :           rewind (203)
     121             :          else
     122             :             CALL judft_error("no proj/proj.1/proj.2",calledby
     123           0 :      +           ="wann_rmat")
     124             :          endif
     125           0 :          read (203,*) num_wann,num_bands
     126           0 :          close (203)
     127           0 :          write(oUnit,*)'According to proj there are ',
     128           0 :      +                 num_bands,' bands'
     129           0 :          write(oUnit,*)"and ",num_wann," wannier functions."
     130             : 
     131             : c****************************************************************
     132             : c        get nntot and bk and wb from bkpts file
     133             : c****************************************************************
     134           0 :          inquire (file='bkpts',exist=l_file)
     135           0 :          if (.not.l_file)  CALL judft_error("need bkpts"
     136           0 :      +      ,calledby ="wann_rmat")
     137           0 :          open (202,file='bkpts',form='formatted',status='old')
     138           0 :          read (202,'(i4)') nntot
     139           0 :          allocate ( gb(3,nntot,nkpts) )
     140           0 :          allocate ( bpt(nntot,nkpts) )
     141           0 :          do ikpt=1,nkpts
     142           0 :           do nn=1,nntot
     143             :            read (202,'(2i6,3x,3i4)')
     144           0 :      &      ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
     145             :           enddo !nn
     146             :          enddo !ikpt
     147           0 :          allocate( wb(nntot) )
     148           0 :          allocate( bpunkt(3,nntot) )
     149           0 :          do nn=1,nntot
     150           0 :             read(202,*)bpunkt(:,nn),wb(nn)
     151             :          enddo
     152           0 :          close (202)
     153             : 
     154           0 :          do ikpt=1,nkpts
     155           0 :           do nn=1,nntot
     156             :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     157           0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     158             :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     159           0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     160             :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     161           0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     162           0 :            kdiff=matmul(transpose(bmat),kdiff)
     163           0 :            l_worksout=.false.
     164           0 :            do nn2=1,nntot
     165           0 :              if (abs(kdiff(1)-bpunkt(1,nn2)).lt.1.e-5)then
     166           0 :                if (abs(kdiff(2)-bpunkt(2,nn2)).lt.1.e-5)then
     167           0 :                  if (abs(kdiff(3)-bpunkt(3,nn2)).lt.1.e-5)then
     168             :                      l_worksout=.true. 
     169             :                      exit
     170             :                  endif          
     171             :                endif
     172             :              endif 
     173             :            enddo !nn2
     174           0 :            if(.not.l_worksout)then
     175           0 :               write(*,*)"ikpt,nn=",ikpt,nn
     176           0 :               write(*,*)"kdiff(1)=",kdiff(1)
     177           0 :               write(*,*)"kdiff(2)=",kdiff(2)
     178           0 :               write(*,*)"kdiff(3)=",kdiff(3)
     179           0 :               stop 'worksout'
     180             :            endif   
     181             :           enddo !nn
     182             :          enddo !ikpt 
     183             : 
     184             : c****************************************************************
     185             : c        read in chk
     186             : c****************************************************************
     187           0 :          num_kpts=nkpts
     188           0 :          allocate( u_matrix_opt(num_bands,num_wann,nkpts) )
     189           0 :          allocate( u_matrix(num_wann,num_wann,nkpts) )
     190           0 :          allocate( lwindow(num_bands,nkpts) )
     191           0 :          allocate( ndimwin(nkpts) )
     192           0 :          allocate( m_matrix(num_wann,num_wann,nntot,num_kpts) )
     193             :          call wann_read_umatrix2(
     194             :      >            nkpts,num_wann,num_bands,
     195             :      >            um_format,jspin,wan90version,
     196             :      <            have_disentangled,
     197             :      <            lwindow,ndimwin,
     198           0 :      <            u_matrix_opt,u_matrix,m_matrix)
     199             : c****************************************************************
     200             : c        read in eig-file
     201             : c****************************************************************
     202           0 :          write(oUnit,*)"read in eig-file"
     203           0 :          allocate(energy(num_bands,num_kpts))
     204           0 :          inquire(file=spin12(jspin)//'.eig',exist=l_umdat)
     205           0 :          IF(.NOT.l_umdat)  CALL judft_error
     206             :      +        ("Thou shall not hide your eig file",calledby
     207           0 :      +        ="wann_hopping")
     208           0 :          open(300,file=spin12(jspin)//'.eig',form='formatted')
     209           0 :          do i=1,num_kpts
     210           0 :            do j=1,num_bands
     211           0 :               read(300,*)band,kpt,energy(j,i)
     212             :            enddo
     213             :          enddo
     214           0 :          close(300)
     215             : 
     216           0 :          minenerg=minval(energy(:,:))
     217           0 :          maxenerg=maxval(energy(:,:))
     218           0 :          write(oUnit,*)"minenerg=",minenerg
     219           0 :          write(oUnit,*)"maxenerg=",maxenerg
     220             : 
     221             : 
     222           0 :          allocate(eigval_opt(num_bands,nkpts))
     223           0 :          allocate(eigval2(num_wann,nkpts))
     224           0 :          eigval_opt=0.0
     225           0 :          eigval2=0.0
     226             : 
     227           0 :          if(have_disentangled) then
     228             : 
     229           0 :            do nkp=1,num_kpts
     230           0 :             counter=0
     231           0 :             do j=1,num_bands
     232           0 :               if(lwindow(j,nkp)) then
     233           0 :                 counter=counter+1
     234           0 :                 eigval_opt(counter,nkp)=energy(j,nkp)
     235             :               end if
     236             :             end do
     237             :            end do
     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 = energy
     250             :          end if                    !have_disentangled
     251             : 
     252           0 :          deallocate(eigval_opt)
     253           0 :          deallocate(energy)
     254             : 
     255             : c****************************************************************
     256             : c        Set up posop.
     257             : c****************************************************************
     258           0 :          write(oUnit,*)"Set up posop."
     259           0 :          allocate( posop2(3,num_wann,num_wann,num_kpts) )
     260           0 :          posop2=cmplx(0.0,0.0)
     261             :          
     262           0 :          do ikpt=1,nkpts
     263           0 :           do nn=1,nntot
     264             :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     265           0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     266             :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     267           0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     268             :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     269           0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     270           0 :            kdiff=matmul(transpose(bmat),kdiff)          
     271           0 :            do i=1,num_wann
     272           0 :             do j=1,num_wann
     273           0 :              if(j.eq.i)then
     274           0 :               do kk=1,3
     275             :                posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
     276           0 :      &             wb(nn)*kdiff(kk)*(m_matrix(j,i,nn,ikpt)-1.0)  
     277             :               enddo !kk
     278             :              else
     279           0 :               do kk=1,3
     280             :                posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
     281           0 :      &             wb(nn)*kdiff(kk)*m_matrix(j,i,nn,ikpt)  
     282             :               enddo !kk
     283             :              endif
     284             :             enddo !j
     285             :            enddo !i
     286             :           enddo !nn   
     287             :          enddo !ikpt
     288             : 
     289           0 :          do ikpt=1,nkpts
     290           0 :           do i=1,num_wann
     291           0 :            do j=1,i
     292           0 :             do kk=1,3
     293             :               posop2(kk,j,i,ikpt)=(posop2(kk,j,i,ikpt)+
     294           0 :      &                        conjg(posop2(kk,i,j,ikpt)))/2.0
     295           0 :               posop2(kk,i,j,ikpt)=posop2(kk,j,i,ikpt)
     296             :             enddo !kk
     297             :            enddo !j
     298             :           enddo !i
     299             :          enddo !ikpt
     300             : 
     301           0 :          allocate( posop(3,num_wann,num_wann,rvecnum) )
     302           0 :          posop=cmplx(0.0,0.0)
     303           0 :          do rvecind=1,rvecnum
     304           0 :             do k=1,nkpts  
     305             :               rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
     306             :      &                    kpoints(2,k)*rvec(2,rvecind)+
     307           0 :      &                    kpoints(3,k)*rvec(3,rvecind) )
     308           0 :               fac=cmplx(cos(rdotk),-sin(rdotk))
     309           0 :               do i=1,num_wann
     310           0 :                do j=1,num_wann
     311           0 :                 do kk=1,3
     312             :                  posop(kk,j,i,rvecind)=
     313             :      &                   posop(kk,j,i,rvecind)+
     314           0 :      &                   fac*posop2(kk,j,i,k)
     315             :                 enddo !kk
     316             :                enddo !j
     317             :               enddo !i
     318             :             enddo !k
     319             :          enddo !rvecind
     320           0 :          posop=posop/cmplx(real(nkpts),0.0)
     321           0 :          deallocate( posop2 )
     322             : 
     323           0 :          do rvecind=1,rvecnum
     324           0 :           if(rvec(1,rvecind).eq.0)then
     325           0 :            if(rvec(2,rvecind).eq.0)then
     326           0 :             if(rvec(3,rvecind).eq.0)then
     327           0 :                rvecind_0=rvecind
     328             :                goto 123
     329             :             endif
     330             :            endif
     331             :           endif
     332             :          enddo !rvecind
     333           0 :          stop 'Ou est ce point-la ?'
     334             :  123     continue
     335             : 
     336           0 :          do i=1,num_wann
     337           0 :           do kk=1,3 
     338           0 :              posop(kk,i,i,rvecind_0)=0.0
     339             :           enddo !kk
     340             :          enddo !i
     341             : 
     342           0 :          do ikpt=1,nkpts
     343           0 :           do nn=1,nntot  
     344             :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     345           0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     346             :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     347           0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     348             :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     349           0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     350           0 :            kdiff=matmul(transpose(bmat),kdiff)/real(nkpts)             
     351           0 :            do i=1,num_wann
     352           0 :             do kk=1,3 
     353             :              posop(kk,i,i,rvecind_0)=posop(kk,i,i,rvecind_0)-
     354           0 :      &        wb(nn)*kdiff(kk)*aimag(log(m_matrix(i,i,nn,ikpt)))
     355             :             enddo !kk
     356             :            enddo !i
     357             :           enddo !nn
     358             :          enddo !ikpt
     359             : 
     360             : c********************************************
     361             : c        Print posop.
     362             : c********************************************
     363           0 :          open(321,file='posop'//spinspin12(jspin),form='formatted')
     364           0 :          do rvecind=1,rvecnum
     365           0 :             r3=rvec(3,rvecind)
     366           0 :             r2=rvec(2,rvecind)
     367           0 :             r1=rvec(1,rvecind)             
     368           0 :             do j=1,num_wann
     369           0 :              do i=1,num_wann
     370           0 :               do kk=1,3  
     371             :                write(321,'(i3,1x,i3,1x,i3,1x,i3,
     372             :      &                 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
     373           0 :      &                 r1,r2,r3,i,j,kk,posop(kk,i,j,rvecind) 
     374             :               enddo !kk  
     375             :              enddo !j
     376             :             enddo !i
     377             :          enddo !rvecind  
     378           0 :          close(321)
     379             : 
     380           0 :          deallocate( lwindow,u_matrix_opt,ndimwin )
     381           0 :          deallocate( u_matrix,m_matrix )
     382           0 :          deallocate( posop,eigval2 )
     383           0 :          deallocate( gb,bpt,wb,bpunkt )
     384             :       enddo !jspin
     385             : 
     386           0 :       call timestop("wann_rmat")
     387           0 :       end subroutine wann_rmat
     388           0 :       end module m_wann_rmat

Generated by: LCOV version 1.14