LCOV - code coverage report
Current view: top level - wannier - wann_rmat.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 170 0.0 %
Date: 2019-09-08 04:53:50 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,l_onedimens,
      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, only:pimach, ImagUnit
      23             :       use m_wann_read_umatrix
      24             : 
      25             :       implicit none
      26             :       real, intent(in) :: bmat(:,:)
      27             :       real, intent(in) :: amat(:,:)
      28             : 
      29             :       integer, intent(in) :: rvecnum
      30             :       integer, intent(in) :: rvec(:,:)
      31             :       real,    intent(in) :: kpoints(:,:)
      32             :       integer, intent(in) :: jspins_in
      33             :       integer, intent(in) :: nkpts
      34             :       logical, intent(in) :: l_bzsym
      35             :       logical, intent(in) :: film
      36             :       logical, intent(in) :: l_onedimens
      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 :       tpi=2*pimach()
      99             : 
     100           0 :       jspins=jspins_in
     101           0 :       if(l_nocosoc)jspins=1
     102             : 
     103           0 :       write(6,*)"nkpts=",nkpts
     104             : 
     105           0 :       do jspin=1,jspins  !spin loop
     106             : c*****************************************************
     107             : c        get num_bands and num_wann from the proj file
     108             : c*****************************************************
     109           0 :          do j=jspin,0,-1
     110           0 :           inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
     111           0 :           if(l_file)then
     112           0 :             filename='proj'//spinspin12(j)
     113           0 :             exit
     114             :           endif
     115             :          enddo
     116           0 :          if(l_file)then
     117           0 :           open (203,file=trim(filename),status='old')
     118           0 :           rewind (203)
     119             :          else
     120             :             CALL judft_error("no proj/proj.1/proj.2",calledby
     121           0 :      +           ="wann_rmat")
     122             :          endif
     123           0 :          read (203,*) num_wann,num_bands
     124           0 :          close (203)
     125           0 :          write(6,*)'According to proj there are ',num_bands,' bands'
     126           0 :          write(6,*)"and ",num_wann," wannier functions."
     127             : 
     128             : c****************************************************************
     129             : c        get nntot and bk and wb from bkpts file
     130             : c****************************************************************
     131           0 :          inquire (file='bkpts',exist=l_file)
     132           0 :          if (.not.l_file)  CALL judft_error("need bkpts"
     133           0 :      +      ,calledby ="wann_rmat")
     134           0 :          open (202,file='bkpts',form='formatted',status='old')
     135           0 :          read (202,'(i4)') nntot
     136           0 :          allocate ( gb(3,nntot,nkpts) )
     137           0 :          allocate ( bpt(nntot,nkpts) )
     138           0 :          do ikpt=1,nkpts
     139           0 :           do nn=1,nntot
     140             :            read (202,'(2i6,3x,3i4)')
     141           0 :      &      ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
     142             :           enddo !nn
     143             :          enddo !ikpt
     144           0 :          allocate( wb(nntot) )
     145           0 :          allocate( bpunkt(3,nntot) )
     146           0 :          do nn=1,nntot
     147           0 :             read(202,*)bpunkt(:,nn),wb(nn)
     148             :          enddo
     149           0 :          close (202)
     150             : 
     151           0 :          do ikpt=1,nkpts
     152           0 :           do nn=1,nntot
     153             :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     154           0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     155             :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     156           0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     157             :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     158           0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     159           0 :            kdiff=matmul(transpose(bmat),kdiff)
     160           0 :            l_worksout=.false.
     161           0 :            do nn2=1,nntot
     162           0 :              if (abs(kdiff(1)-bpunkt(1,nn2)).lt.1.e-5)then
     163           0 :                if (abs(kdiff(2)-bpunkt(2,nn2)).lt.1.e-5)then
     164           0 :                  if (abs(kdiff(3)-bpunkt(3,nn2)).lt.1.e-5)then
     165             :                      l_worksout=.true. 
     166             :                      exit
     167             :                  endif          
     168             :                endif
     169             :              endif 
     170             :            enddo !nn2
     171           0 :            if(.not.l_worksout)then
     172           0 :               write(*,*)"ikpt,nn=",ikpt,nn
     173           0 :               write(*,*)"kdiff(1)=",kdiff(1)
     174           0 :               write(*,*)"kdiff(2)=",kdiff(2)
     175           0 :               write(*,*)"kdiff(3)=",kdiff(3)
     176           0 :               stop 'worksout'
     177             :            endif   
     178             :           enddo !nn
     179             :          enddo !ikpt 
     180             : 
     181             : c****************************************************************
     182             : c        read in chk
     183             : c****************************************************************
     184           0 :          num_kpts=nkpts
     185           0 :          allocate( u_matrix_opt(num_bands,num_wann,nkpts) )
     186           0 :          allocate( u_matrix(num_wann,num_wann,nkpts) )
     187           0 :          allocate( lwindow(num_bands,nkpts) )
     188           0 :          allocate( ndimwin(nkpts) )
     189           0 :          allocate( m_matrix(num_wann,num_wann,nntot,num_kpts) )
     190             :          call wann_read_umatrix2(
     191             :      >            nkpts,num_wann,num_bands,
     192             :      >            um_format,jspin,wan90version,
     193             :      <            have_disentangled,
     194             :      <            lwindow,ndimwin,
     195           0 :      <            u_matrix_opt,u_matrix,m_matrix)
     196             : c****************************************************************
     197             : c        read in eig-file
     198             : c****************************************************************
     199           0 :          write(6,*)"read in eig-file"
     200           0 :          allocate(energy(num_bands,num_kpts))
     201           0 :          inquire(file=spin12(jspin)//'.eig',exist=l_umdat)
     202           0 :          IF(.NOT.l_umdat)  CALL judft_error
     203             :      +        ("Thou shall not hide your eig file",calledby
     204           0 :      +        ="wann_hopping")
     205           0 :          open(300,file=spin12(jspin)//'.eig',form='formatted')
     206           0 :          do i=1,num_kpts
     207           0 :            do j=1,num_bands
     208           0 :               read(300,*)band,kpt,energy(j,i)
     209             :            enddo
     210             :          enddo
     211           0 :          close(300)
     212             : 
     213           0 :          minenerg=minval(energy(:,:))
     214           0 :          maxenerg=maxval(energy(:,:))
     215           0 :          write(6,*)"minenerg=",minenerg
     216           0 :          write(6,*)"maxenerg=",maxenerg
     217             : 
     218             : 
     219           0 :          allocate(eigval_opt(num_bands,nkpts))
     220           0 :          allocate(eigval2(num_wann,nkpts))
     221           0 :          eigval_opt=0.0
     222           0 :          eigval2=0.0
     223             : 
     224           0 :          if(have_disentangled) then
     225             : 
     226           0 :            do nkp=1,num_kpts
     227           0 :             counter=0
     228           0 :             do j=1,num_bands
     229           0 :               if(lwindow(j,nkp)) then
     230           0 :                 counter=counter+1
     231           0 :                 eigval_opt(counter,nkp)=energy(j,nkp)
     232             :               end if
     233             :             end do
     234             :            end do
     235             :        
     236           0 :            do nkp=1,num_kpts
     237           0 :             do j=1,num_wann
     238           0 :              do m=1,ndimwin(nkp)
     239             :                 eigval2(j,nkp)=eigval2(j,nkp)+eigval_opt(m,nkp)* 
     240           0 :      &    real(conjg(u_matrix_opt(m,j,nkp))*u_matrix_opt(m,j,nkp))
     241             :              enddo
     242             :             enddo
     243             :            enddo
     244             : 
     245             :          else
     246           0 :            eigval2 = energy
     247             :          end if                    !have_disentangled
     248             : 
     249           0 :          deallocate(eigval_opt)
     250           0 :          deallocate(energy)
     251             : 
     252             : c****************************************************************
     253             : c        Set up posop.
     254             : c****************************************************************
     255           0 :          write(6,*)"Set up posop."
     256           0 :          allocate( posop2(3,num_wann,num_wann,num_kpts) )
     257           0 :          posop2=cmplx(0.0,0.0)
     258             :          
     259           0 :          do ikpt=1,nkpts
     260           0 :           do nn=1,nntot
     261             :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     262           0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     263             :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     264           0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     265             :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     266           0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     267           0 :            kdiff=matmul(transpose(bmat),kdiff)          
     268           0 :            do i=1,num_wann
     269           0 :             do j=1,num_wann
     270           0 :              if(j.eq.i)then
     271           0 :               do kk=1,3
     272             :                posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
     273           0 :      &             wb(nn)*kdiff(kk)*(m_matrix(j,i,nn,ikpt)-1.0)  
     274             :               enddo !kk
     275             :              else
     276           0 :               do kk=1,3
     277             :                posop2(kk,j,i,ikpt)=posop2(kk,j,i,ikpt)+ImagUnit*
     278           0 :      &             wb(nn)*kdiff(kk)*m_matrix(j,i,nn,ikpt)  
     279             :               enddo !kk
     280             :              endif
     281             :             enddo !j
     282             :            enddo !i
     283             :           enddo !nn   
     284             :          enddo !ikpt
     285             : 
     286           0 :          do ikpt=1,nkpts
     287           0 :           do i=1,num_wann
     288           0 :            do j=1,i
     289           0 :             do kk=1,3
     290             :               posop2(kk,j,i,ikpt)=(posop2(kk,j,i,ikpt)+
     291           0 :      &                        conjg(posop2(kk,i,j,ikpt)))/2.0
     292           0 :               posop2(kk,i,j,ikpt)=posop2(kk,j,i,ikpt)
     293             :             enddo !kk
     294             :            enddo !j
     295             :           enddo !i
     296             :          enddo !ikpt
     297             : 
     298           0 :          allocate( posop(3,num_wann,num_wann,rvecnum) )
     299           0 :          posop=cmplx(0.0,0.0)
     300           0 :          do rvecind=1,rvecnum
     301           0 :             do k=1,nkpts  
     302             :               rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
     303             :      &                    kpoints(2,k)*rvec(2,rvecind)+
     304           0 :      &                    kpoints(3,k)*rvec(3,rvecind) )
     305           0 :               fac=cmplx(cos(rdotk),-sin(rdotk))
     306           0 :               do i=1,num_wann
     307           0 :                do j=1,num_wann
     308           0 :                 do kk=1,3
     309             :                  posop(kk,j,i,rvecind)=
     310             :      &                   posop(kk,j,i,rvecind)+
     311           0 :      &                   fac*posop2(kk,j,i,k)
     312             :                 enddo !kk
     313             :                enddo !j
     314             :               enddo !i
     315             :             enddo !k
     316             :          enddo !rvecind
     317           0 :          posop=posop/cmplx(real(nkpts),0.0)
     318           0 :          deallocate( posop2 )
     319             : 
     320           0 :          do rvecind=1,rvecnum
     321           0 :           if(rvec(1,rvecind).eq.0)then
     322           0 :            if(rvec(2,rvecind).eq.0)then
     323           0 :             if(rvec(3,rvecind).eq.0)then
     324           0 :                rvecind_0=rvecind
     325             :                goto 123
     326             :             endif
     327             :            endif
     328             :           endif
     329             :          enddo !rvecind
     330           0 :          stop 'Ou est ce point-la ?'
     331             :  123     continue
     332             : 
     333           0 :          do i=1,num_wann
     334           0 :           do kk=1,3 
     335           0 :              posop(kk,i,i,rvecind_0)=0.0
     336             :           enddo !kk
     337             :          enddo !i
     338             : 
     339           0 :          do ikpt=1,nkpts
     340           0 :           do nn=1,nntot  
     341             :            kdiff(1)=kpoints(1,bpt(nn,ikpt))+
     342           0 :      +            gb(1,nn,ikpt) - kpoints(1,ikpt)
     343             :            kdiff(2)=kpoints(2,bpt(nn,ikpt))+
     344           0 :      +            gb(2,nn,ikpt) - kpoints(2,ikpt)
     345             :            kdiff(3)=kpoints(3,bpt(nn,ikpt))+
     346           0 :      +            gb(3,nn,ikpt) - kpoints(3,ikpt)
     347           0 :            kdiff=matmul(transpose(bmat),kdiff)/real(nkpts)             
     348           0 :            do i=1,num_wann
     349           0 :             do kk=1,3 
     350             :              posop(kk,i,i,rvecind_0)=posop(kk,i,i,rvecind_0)-
     351           0 :      &        wb(nn)*kdiff(kk)*aimag(log(m_matrix(i,i,nn,ikpt)))
     352             :             enddo !kk
     353             :            enddo !i
     354             :           enddo !nn
     355             :          enddo !ikpt
     356             : 
     357             : c********************************************
     358             : c        Print posop.
     359             : c********************************************
     360           0 :          open(321,file='posop'//spinspin12(jspin),form='formatted')
     361           0 :          do rvecind=1,rvecnum
     362           0 :             r3=rvec(3,rvecind)
     363           0 :             r2=rvec(2,rvecind)
     364           0 :             r1=rvec(1,rvecind)             
     365           0 :             do j=1,num_wann
     366           0 :              do i=1,num_wann
     367           0 :               do kk=1,3  
     368             :                write(321,'(i3,1x,i3,1x,i3,1x,i3,
     369             :      &                 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
     370           0 :      &                 r1,r2,r3,i,j,kk,posop(kk,i,j,rvecind) 
     371             :               enddo !kk  
     372             :              enddo !j
     373             :             enddo !i
     374             :          enddo !rvecind  
     375           0 :          close(321)
     376             : 
     377           0 :          deallocate( lwindow,u_matrix_opt,ndimwin )
     378           0 :          deallocate( u_matrix,m_matrix )
     379           0 :          deallocate( posop,eigval2 )
     380           0 :          deallocate( gb,bpt,wb,bpunkt )
     381             :       enddo !jspin
     382             : 
     383           0 :       end subroutine wann_rmat
     384             :       end module m_wann_rmat

Generated by: LCOV version 1.13