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

Generated by: LCOV version 1.13