wann_dipole2.f 4.09 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
c********************************************************
c     Calculate the electric polarization.
c     Frank Freimuth, October 2006
c********************************************************
      module m_wann_dipole2
      use m_juDFT
      contains
      subroutine wann_dipole2(
     >               jspins,pos,omtil,natd,
     >               l_nocosoc)
      use m_wann_readcenters
      implicit none
      integer,intent(in) :: jspins
      real,intent(in)    :: omtil
      integer,intent(in) :: natd
      real,intent(in)    :: pos(3,natd)
      logical,intent(in) :: l_nocosoc

      integer            :: jspin,jspins2
      logical            :: l_file
      integer,allocatable:: ind(:)
      integer            :: nwf,nwfs
      real,allocatable   :: distance(:,:)
      real               :: moment(3,jspins)
      real               :: moment2(3)
      real               :: elemchargmu,bohrtocm
      integer            :: shifted,ishift,ind1,ind2
      real,allocatable   :: centers(:,:)
      real,allocatable   :: spreads(:)

      open(204,file='dipole_out')
      write(204,*)"*****************************************"
      write(204,*)"Calculation of the electric polarization."
      write(204,*)"*****************************************"
      elemchargmu=1.60217646e-13
      bohrtocm=0.529177e-8

      jspins2=jspins
      if(l_nocosoc)jspins2=1
      do jspin=1,jspins2
         inquire(file='proj',exist=l_file)
         IF(.NOT.l_file) CALL juDFT_error("proj not found",calledby
     +        ="wann_dipole2")
         open(203,file='proj',form='formatted')
         read(203,*)nwfs
         print*,"nwfs=",nwfs
         allocate(ind(nwfs))
         allocate(distance(3,nwfs))
         allocate(centers(3,nwfs))
         allocate(spreads(nwfs))
         do nwf=1,nwfs
          read(203,*)ind(nwf)
          read(203,*)
         enddo
         close(203)
         write(204,*)"spin=",jspin
         call wann_readcenters(nwfs,jspin,centers,spreads)
         write(204,*)"*****centers:******"
         do nwf=1,nwfs
           write(204,fmt=1000)nwf,centers(:,nwf),spreads(nwf)
         enddo
         moment(:,jspin)=0.0
         do nwf=1,nwfs
            distance(:,nwf)=pos(:,ind(nwf))-centers(:,nwf)
            moment(:,jspin)=moment(:,jspin)+distance(:,nwf)
         enddo
         if((jspins.eq.1).and..not.l_nocosoc) 
     &                    moment(:,jspin)=moment(:,jspin)*2
         moment(:,jspin)=moment(:,jspin)/omtil
         moment(:,jspin)=moment(:,jspin)*elemchargmu/((bohrtocm)**2)
         write(204,*)"polarization due to shift of wannierorbitals"
         write(204,*)moment(:,jspin),"uC/cm2"
         deallocate(ind,distance,centers,spreads)
      enddo!jspin
1000  format(2x,'WF centre and spread', 
     &       i5,2x,'(',f10.6,',',f10.6,',',f10.6,' )',f15.8)
      inquire(file='dipole',exist=l_file)
      if(.not.l_file)then
            write(204,*)"no file dipole found"
      else
         moment2(:)=0.0
         open(224,file='dipole',form='formatted')
         read(224,*)shifted
         do ishift=1,shifted
               read(224,*)ind1,ind2 !electron taken from 1 and moved to 2
               moment2(:)=moment2(:)+pos(:,ind1)-pos(:,ind2)
         enddo
         close(224)
         moment2(:)=moment2(:)/omtil
         moment2(:)=moment2(:)*elemchargmu/(bohrtocm)**2
         write(204,*)"polarization due to electrons 
     &              being catched by other atoms"
         write(204,*)moment2(:),"uC/cm2"

      endif
      write(204,*)"total moment:"
      if(jspins2.eq.2)moment(:,1)=moment(:,1)+moment(:,2)
      if(l_file)then
         moment(:,1)=moment(:,1)+moment2(:)
      endif
      write(204,*)moment(:,1),"uC/cm2"
      close(204)
      end subroutine wann_dipole2
      end module m_wann_dipole2