LCOV - code coverage report
Current view: top level - wannier - wann_nabla_rs.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 130 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_nabla_rs
       8             :       use m_juDFT
       9             :       contains 
      10           0 :       subroutine wann_nabla_rs(
      11           0 :      >          rvecnum,rvec,kpoints,
      12             :      >          jspins_in,nkpts,l_bzsym,film,l_onedimens,
      13             :      >          l_soc,band_min,band_max,neigd,
      14             :      >          l_socmmn0,wan90version)
      15             : c*************************************************
      16             : c     Calculate the matrix elements of nabla
      17             : c     in real space from the
      18             : c     files WF1.chk (and WF1_um.dat) produced
      19             : c     by wannier90.
      20             : c     FF, December 2009
      21             : c*************************************************
      22             :       use m_constants, only:pimach
      23             :       use m_wann_read_umatrix
      24             : 
      25             :       implicit none
      26             :       integer, intent(in) :: rvecnum
      27             :       integer, intent(in) :: rvec(:,:)
      28             :       real,    intent(in) :: kpoints(:,:)
      29             :       integer, intent(in) :: jspins_in
      30             :       integer, intent(in) :: nkpts
      31             :       logical, intent(in) :: l_bzsym
      32             :       logical, intent(in) :: film
      33             :       logical, intent(in) :: l_onedimens
      34             : 
      35             :       logical, intent(in) :: l_soc
      36             :       integer, intent(in) :: band_min(2),band_max(2),neigd
      37             : 
      38             :       logical, intent(in) :: l_socmmn0
      39             :       integer, intent(in) :: wan90version
      40             : 
      41             :       integer             :: ikpt,jspins
      42             :       integer             :: kpts
      43             :       logical             :: l_file
      44             : c      real                :: kpoints(3,nkpts)
      45             :       integer             :: num_wann,num_kpts,num_nnmax,jspin
      46             :       integer             :: kspin,kkspin
      47             :       integer             :: wann_shift,num_wann2
      48             :       integer             :: i,j,k,m,info,r1,r2,r3,dummy1
      49             :       integer             :: dummy2,dummy3,dummy4
      50             :       integer             :: hopmin,hopmax,counter,m1,m2
      51             :       integer             :: num_bands2
      52             :       integer,allocatable :: iwork(:)
      53             :       real,allocatable    :: energy(:,:),ei(:)
      54             :       real,allocatable    :: eigw(:,:),rwork(:)
      55             :       complex,allocatable :: work(:),vec(:,:)
      56           0 :       complex,allocatable :: u_matrix(:,:,:),hwann(:,:,:,:)
      57           0 :       complex,allocatable :: hreal(:,:,:,:)
      58             :       complex,allocatable :: hrealsoc(:,:,:,:,:,:,:)
      59             :       complex,allocatable :: hwannsoc(:,:,:,:,:)
      60           0 :       complex,allocatable :: nablamat(:,:,:,:,:)
      61           0 :       complex,allocatable :: nablamat2(:,:,:,:)
      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,l_she
      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
      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             :: rvecind
      90             : 
      91             :       data spinspin12/'  ','.1' , '.2'/
      92             :       data spin12/'WF1','WF2'/
      93             : 
      94           0 :       inquire(file='she',exist=l_she)
      95             : 
      96           0 :       tpi=2*pimach()
      97             : 
      98           0 :       jspins=jspins_in
      99           0 :       if(l_soc)jspins=1
     100             : 
     101           0 :       write(6,*)"nkpts=",nkpts
     102             : 
     103             : c$$$c***************************************************
     104             : c$$$c     read in the kpoints from w90kpts or kpts
     105             : c$$$c***************************************************      
     106             : c$$$      if (l_bzsym) then
     107             : c$$$         l_file=.false.
     108             : c$$$         inquire(file='w90kpts',exist=l_file)
     109             : c$$$         if (.not.l_file) stop 'where is w90kpts?'
     110             : c$$$         open(177,file='w90kpts',form='formatted')
     111             : c$$$         read(177,*)kpts,scale
     112             : c$$$      else
     113             : c$$$         l_file=.false.
     114             : c$$$         inquire(file='kpts',exist=l_file)
     115             : c$$$         if(.not.l_file) stop 'where is kpts?'
     116             : c$$$         open(177,file='kpts',form='formatted')
     117             : c$$$         read(177,*)kpts,scale
     118             : c$$$      endif   
     119             : c$$$      if(kpts.ne.nkpts) stop 'mismatch in number of kpoints'
     120             : c$$$
     121             : c$$$      do ikpt = 1,nkpts 
     122             : c$$$         read(177,*)kpoints(:,ikpt)
     123             : c$$$      enddo 
     124             : c$$$
     125             : c$$$      if(film.and..not.l_onedimens)then
     126             : c$$$         kpoints(3,:)=0.0
     127             : c$$$      endif   
     128             : c$$$      kpoints=kpoints/scale
     129             : c$$$      close(177)
     130             : 
     131           0 :       do jspin=1,jspins  !spin loop
     132             : c*****************************************************
     133             : c     get num_bands and num_wann from the proj file
     134             : c*****************************************************
     135           0 :          do j=jspin,0,-1
     136           0 :           inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
     137           0 :           if(l_file)then
     138           0 :             filename='proj'//spinspin12(j)
     139           0 :             exit
     140             :           endif
     141             :          enddo
     142           0 :          if(l_file)then
     143           0 :           open (203,file=trim(filename),status='old')
     144           0 :           rewind (203)
     145             :          else
     146             :             CALL juDFT_error("no proj/proj.1/proj.2",calledby
     147           0 :      +           ="wann_nabla_rs")
     148             :          endif
     149           0 :          read (203,*) num_wann,num_bands
     150           0 :          close (203)
     151           0 :          write(6,*)'According to proj there are ',num_bands,' bands'
     152           0 :          write(6,*)"and ",num_wann," wannier functions."
     153             : 
     154             : c****************************************************************
     155             : c        read in chk
     156             : c****************************************************************
     157           0 :          num_kpts=nkpts
     158           0 :          allocate( u_matrix_opt(num_bands,num_wann,nkpts) )
     159           0 :          allocate( u_matrix(num_wann,num_wann,nkpts) )
     160           0 :          allocate( lwindow(num_bands,nkpts) )
     161           0 :          allocate( ndimwin(nkpts) )
     162             :          call wann_read_umatrix2(
     163             :      >            nkpts,num_wann,num_bands,
     164             :      >            um_format,jspin,wan90version,
     165             :      <            have_disentangled,
     166             :      <            lwindow,ndimwin,
     167           0 :      <            u_matrix_opt,u_matrix)
     168             : 
     169           0 :          num_bands2=num_bands
     170           0 :          if(l_soc.and.l_socmmn0)then
     171           0 :           num_bands2=neigd
     172             :          endif
     173             :  
     174             : c************************************************
     175             : c        Read the files "WF1.nabl" and "WF2.nabl"
     176             : c************************************************
     177           0 :          allocate( nablamat(3,num_bands2,num_bands2,nkpts,2) )
     178           0 :          open(304,file=spin12(jspin)//'.nabl',form='formatted')
     179           0 :          read(304,*)
     180           0 :          read(304,*)
     181           0 :          do nkp=1,num_kpts
     182           0 :           do i=1,num_bands2
     183           0 :            do j=1,num_bands2
     184           0 :             do k=1,3  
     185           0 :              read(304,*)dummy1,dummy2,dummy3,dummy4,a,b
     186           0 :              nablamat(k,j,i,nkp,1)=cmplx(a,-b)
     187             :             enddo !k 
     188             :            enddo !j
     189             :           enddo !i
     190             :          enddo !nkp
     191           0 :          close(304)
     192             : 
     193           0 :          if(l_soc)then
     194           0 :           open(304,file=spin12(2)//'.nabl',form='formatted')
     195           0 :           read(304,*)
     196           0 :           read(304,*)
     197           0 :           do nkp=1,num_kpts
     198           0 :            do i=1,num_bands2
     199           0 :             do j=1,num_bands2
     200           0 :              do k=1,3  
     201           0 :               read(304,*)dummy1,dummy2,dummy3,dummy4,a,b
     202           0 :               nablamat(k,j,i,nkp,2)=cmplx(a,-b)
     203             :              enddo !k 
     204             :             enddo !j
     205             :            enddo !i
     206             :           enddo !nkp
     207           0 :           close(304)
     208             :          endif
     209             : 
     210             : c**************************************
     211             : c        Enforce Hermiticity
     212             : c**************************************
     213           0 :          do nkp=1,num_kpts
     214           0 :           do i=1,num_bands2
     215           0 :            do j=1,i
     216           0 :             do k=1,3
     217             :              nablamat(k,i,j,nkp,1)=( nablamat(k,i,j,nkp,1)+
     218           0 :      &                         conjg(nablamat(k,j,i,nkp,1)) )/2.0
     219           0 :              nablamat(k,j,i,nkp,1)=conjg(nablamat(k,i,j,nkp,1))
     220           0 :              if(l_soc)then
     221             :                nablamat(k,i,j,nkp,2)=( nablamat(k,i,j,nkp,2)+
     222           0 :      &                         conjg(nablamat(k,j,i,nkp,2)) )/2.0
     223           0 :                nablamat(k,j,i,nkp,2)=conjg(nablamat(k,i,j,nkp,2))
     224             :              endif
     225             :             enddo !k
     226             :            enddo !j
     227             :           enddo !i
     228             :          enddo !nkp
     229             : 
     230             : c**************************************
     231             : c        Nabla-sum
     232             : c**************************************
     233           0 :          if( l_soc )then
     234           0 :           if(l_she)then
     235           0 :            do nkp=1,num_kpts
     236           0 :             do i=1,num_bands2
     237           0 :              do j=1,num_bands2
     238           0 :               do k=1,3  
     239             :                nablamat(k,j,i,nkp,1)=nablamat(k,j,i,nkp,1)
     240           0 :      &                            -nablamat(k,j,i,nkp,2)
     241             :               enddo !k 
     242             :              enddo !j
     243             :             enddo !i
     244             :            enddo !nkp
     245             :           else   
     246           0 :            do nkp=1,num_kpts
     247           0 :             do i=1,num_bands2
     248           0 :              do j=1,num_bands2
     249           0 :               do k=1,3  
     250             :                nablamat(k,j,i,nkp,1)=nablamat(k,j,i,nkp,1)
     251           0 :      &                            +nablamat(k,j,i,nkp,2)
     252             :               enddo !k 
     253             :              enddo !j
     254             :             enddo !i
     255             :            enddo !nkp
     256             :           endif
     257             :          endif !jspins.eq.2
     258             : 
     259             : c****************************************************************
     260             : c        Calculate matrix elements of Nabla in the basis of
     261             : c        rotated Bloch functions.
     262             : c****************************************************************
     263           0 :          allocate( nablamat2(3,num_wann,num_wann,nkpts) )
     264             :          write(6,*)"calculate matrix elements of momentum operator
     265           0 :      &   between wannier orbitals"
     266             : 
     267           0 :          if(have_disentangled) then       
     268           0 :           nablamat2=0.0  
     269           0 :           do nkp=1,num_kpts
     270           0 :            do j=1,num_wann
     271           0 :             do jp=1,num_wann  
     272           0 :              do m=1,ndimwin(nkp)
     273           0 :               do mp=1,ndimwin(nkp)  
     274           0 :                do k=1,3  
     275             :                 nablamat2(k,jp,j,nkp)=nablamat2(k,jp,j,nkp)+ 
     276             :      &                 conjg(u_matrix_opt(mp,jp,nkp))*
     277             :      &                        nablamat(k,mp,m,nkp,1)*
     278           0 :      &                       u_matrix_opt(m,j,nkp)
     279             :                enddo !k 
     280             :               enddo !mp  
     281             :              enddo !m
     282             :             enddo !jp 
     283             :            enddo !j
     284             :           enddo !nkp
     285             :          else
     286           0 :           nablamat2(:,:,:,:)=nablamat(:,:,:,:,1)
     287             :          end if                    !have_disentangled
     288             : 
     289           0 :          allocate(hwann(3,num_wann,num_wann,num_kpts))
     290           0 :          hwann=cmplx(0.0,0.0)
     291           0 :          wann_shift=0
     292             :          if(l_socmmn0)then
     293             :             wann_shift=band_min(jspin)-1
     294             :          endif
     295           0 :          do k=1,num_kpts
     296           0 :           do m=1,num_wann
     297           0 :            do i=1,num_wann
     298           0 :             do j=1,num_wann
     299           0 :              do mp=1,num_wann
     300           0 :               do kk=1,3  
     301             :                 hwann(kk,mp,m,k)=hwann(kk,mp,m,k)+
     302             :      *              conjg(u_matrix(j,mp,k))*
     303             :      *                nablamat2(kk,j,i,k)*
     304           0 :      *                    u_matrix(i,m,k)
     305             :               enddo !kk
     306             :              enddo !mp   
     307             :             enddo !j
     308             :            enddo !i
     309             :           enddo !m
     310             :          enddo !k
     311             : 
     312             : c************************************************************
     313             : c        Calculate matrix elements in real space.
     314             : c***********************************************************      
     315           0 :          write(6,*)"calculate nabla-mat in rs"
     316             : c$$$         hopmin=-5
     317             : c$$$         hopmax=5
     318             : c$$$         allocate(hreal(3,num_wann,num_wann,hopmin:hopmax,
     319             : c$$$     &         hopmin:hopmax,hopmin:hopmax))
     320             : c$$$         hreal=cmplx(0.0,0.0)
     321             : c$$$         do r3=hopmin,hopmax
     322             : c$$$          do r2=hopmin,hopmax
     323             : c$$$           do r1=hopmin,hopmax
     324             : c$$$            do k=1,nkpts  
     325             : c$$$              rdotk=tpi*(kpoints(1,k)*r1+kpoints(2,k)*r2+
     326             : c$$$     &                                   kpoints(3,k)*r3)
     327             : c$$$              fac=cmplx(cos(rdotk),-sin(rdotk))/nkpts
     328             : c$$$              do i=1,num_wann
     329             : c$$$               do j=1,num_wann
     330             : c$$$                do kk=1,3
     331             : c$$$                 hreal(kk,j,i,r1,r2,r3)=
     332             : c$$$     &                   hreal(kk,j,i,r1,r2,r3)+
     333             : c$$$     &                   fac*hwann(kk,j,i,k)
     334             : c$$$                enddo !kk
     335             : c$$$               enddo !j
     336             : c$$$              enddo !i
     337             : c$$$            enddo !k
     338             : c$$$           enddo !r1
     339             : c$$$          enddo !r2
     340             : c$$$         enddo !r3
     341             : 
     342             : 
     343           0 :          allocate(hreal(3,num_wann,num_wann,rvecnum))
     344           0 :          hreal=cmplx(0.0,0.0)
     345             : 
     346           0 :          do rvecind=1,rvecnum
     347           0 :             do k=1,nkpts  
     348             :               rdotk=tpi*(kpoints(1,k)*rvec(1,rvecind)+
     349             :      &                   kpoints(2,k)*rvec(2,rvecind)+
     350           0 :      &                   kpoints(3,k)*rvec(3,rvecind) )
     351           0 :               fac=cmplx(cos(rdotk),-sin(rdotk))
     352           0 :               do i=1,num_wann
     353           0 :                do j=1,num_wann
     354           0 :                 do kk=1,3
     355             :                  hreal(kk,j,i,rvecind)=
     356             :      &                   hreal(kk,j,i,rvecind)+
     357           0 :      &                   fac*hwann(kk,j,i,k)
     358             :                 enddo !kk
     359             :                enddo !j
     360             :               enddo !i
     361             :             enddo !k
     362             :          enddo !rvecind
     363           0 :          hreal=hreal/cmplx(real(nkpts),0.0)
     364             : 
     365           0 :          if(l_she)then
     366           0 :           open(321,file='rs_sv'//spinspin12(jspin),form='formatted')
     367             :          else   
     368           0 :           open(321,file='rsnabla'//spinspin12(jspin),form='formatted')
     369             :          endif
     370             : 
     371           0 :          do rvecind=1,rvecnum
     372           0 :             r3=rvec(3,rvecind)
     373           0 :             r2=rvec(2,rvecind)
     374           0 :             r1=rvec(1,rvecind)             
     375           0 :             do j=1,num_wann
     376           0 :              do i=1,num_wann
     377           0 :               do kk=1,3  
     378             :                write(321,'(i3,1x,i3,1x,i3,1x,i3,
     379             :      &                 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
     380           0 :      &                 r1,r2,r3,i,j,kk,hreal(kk,i,j,rvecind) 
     381             :               enddo !kk  
     382             :              enddo !j
     383             :             enddo !i
     384             :          enddo !rvecind  
     385           0 :          close(321)
     386             : 
     387           0 :          deallocate(lwindow,u_matrix_opt,ndimwin)
     388           0 :          deallocate(u_matrix,hwann,hreal)
     389           0 :          deallocate(nablamat,nablamat2)
     390             :       enddo !jspin
     391             : 
     392           0 :       end subroutine wann_nabla_rs
     393             :       end module m_wann_nabla_rs

Generated by: LCOV version 1.13