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_dipole_electronic
       8             :       use m_juDFT
       9             :       contains
      10           0 :       subroutine wann_dipole_electronic(
      11           0 :      >               natd,pos,omtil,
      12             :      >               jspins,l_absolute,num_wann,
      13             :      <               electronic_moment)
      14             : c*********************************************
      15             : c     Calculate the dipole moment due to 
      16             : c     Wannier-electrons.
      17             : c
      18             : c     This subroutine can operate in two modes:
      19             : c
      20             : c     i)  l_absolute=.true.  => Simply evaluate
      21             : c         the center of mass moment of all 
      22             : c         Wannier functions.
      23             : c 
      24             : c     ii) l_absolute=.false. => Sum up the dipoles
      25             : c         that arise from the displacements of the
      26             : c         Wannier functions' centers from their 
      27             : c         host atoms. 
      28             : c
      29             : c     Frank Freimuth
      30             : c*********************************************
      31             : 
      32             :       USE m_constants
      33             :       use m_wann_readcenters
      34             : 
      35             :       implicit none
      36             : 
      37             :       integer, intent(in)  :: natd
      38             :       real,    intent(in)  :: pos(3,natd)
      39             :       real,    intent(in)  :: omtil
      40             :       integer, intent(in)  :: jspins
      41             :       logical, intent(in)  :: l_absolute
      42             :       integer, intent(in)  :: num_wann(2)
      43             :       real,    intent(out) :: electronic_moment(3,2)
      44             : 
      45             :       real                 :: electronic_polarization(3,2)
      46             :       real,parameter       :: elemchargmu=1.60217646e-13
      47             :       real,parameter       :: bohrtocm=0.529177e-8
      48             :       integer              :: jspin,j
      49             :       logical              :: l_file
      50             :       character(len=2)     :: spin12(0:2)
      51             :       character(len=6)     :: filename
      52             :       integer              :: num_wann_dum
      53           0 :       integer, allocatable :: waind(:)
      54             :       integer              :: nwf,k
      55           0 :       real, allocatable    :: wann_centers(:,:)
      56             : 
      57             :       data spin12/'  ', '.1', '.2'/
      58             : 
      59           0 :       call timestart("wann_dipole_electronic")
      60           0 :       write(666,*)"*****************************"
      61           0 :       write(666,*)" Electronic terms            "
      62           0 :       write(666,*)"*****************************"
      63           0 :       write(*,*)  "*****************************"
      64           0 :       write(*,*)  " Electronic terms            "
      65           0 :       write(*,*)  "*****************************"
      66             : 
      67           0 :       electronic_moment=0.0
      68             :       
      69           0 :       do jspin=1,jspins
      70           0 :         if(.not.l_absolute)then 
      71             : c--------reading the proj.1 / proj.2 / proj file
      72           0 :          do j=jspin,0,-1
      73           0 :           inquire(file=trim('proj'//spin12(j)),exist=l_file)
      74           0 :           if(l_file)then
      75           0 :              filename='proj'//spin12(j)
      76           0 :              exit
      77             :           endif
      78             :          enddo
      79             :          
      80           0 :          if(l_file)then
      81           0 :           open (203,file=trim(filename),status='old')
      82           0 :           rewind (203)
      83             :          else
      84             :             CALL juDFT_error("no proj/proj.1/proj.2",calledby
      85           0 :      +           ="wann_dipole_electronic")
      86             :          endif  
      87           0 :          read(203,*)num_wann_dum
      88           0 :          if(num_wann_dum/=num_wann(jspin))  CALL juDFT_error
      89           0 :      +        ("num_wann_dum",calledby ="wann_dipole_electronic")
      90           0 :          allocate( waind(num_wann(jspin)) )
      91           0 :          do nwf=1,num_wann(jspin)
      92           0 :           read(203,*)waind(nwf)
      93           0 :           read(203,*)
      94             :          enddo
      95           0 :          close(203)
      96             :         endif 
      97           0 :         print*,"number of wannier functions= ",num_wann(jspin)
      98           0 :         write(oUnit,*)"number of wannier functions",num_wann(jspin)
      99           0 :         allocate( wann_centers(3,num_wann(jspin)) )
     100             :         call wann_readcenters(
     101             :      >         num_wann(jspin),jspin,
     102           0 :      <         wann_centers(:,:))
     103           0 :         if(l_absolute)then
     104           0 :          do k=1,num_wann(jspin)  
     105             :            electronic_moment(:,jspin)=
     106           0 :      &     electronic_moment(:,jspin)-wann_centers(:,k)
     107             :          enddo  
     108             :         else
     109           0 :          do k=1,num_wann(jspin)
     110             :            electronic_moment(:,jspin)=
     111             :      &     electronic_moment(:,jspin)+
     112             :      +        pos(:,waind(k))-
     113           0 :      -        wann_centers(:,k)
     114             :          enddo
     115           0 :          deallocate(waind)
     116             :         endif 
     117           0 :         deallocate(wann_centers)
     118           0 :         write(*,  fmt=555)jspin,electronic_moment(:,jspin)
     119           0 :         write(666,fmt=555)jspin,electronic_moment(:,jspin)
     120           0 :         if(.not.l_absolute)then
     121             :            electronic_polarization=
     122           0 :      &            electronic_moment/omtil*elemchargmu/((bohrtocm)**2)
     123           0 :            write(*,  fmt=777)jspin,electronic_polarization(:,jspin)
     124           0 :            write(666,fmt=777)jspin,electronic_polarization(:,jspin)
     125             :         endif
     126             :       enddo  
     127             :  555  format("spin ",i1.1," electronic moment = (",3f12.6,") a.u.")
     128             :  777  format("spin ",i1.1,
     129             :      &          " electronic polarization = (",3f12.6,") uC/cm**2")
     130             : 
     131           0 :       call timestop("wann_dipole_electronic")
     132           0 :       end subroutine wann_dipole_electronic
     133             :       end module m_wann_dipole_electronic