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 electric polarization.
       9             : c     Frank Freimuth, October 2006
      10             : c********************************************************
      11             :       module m_wann_dipole2
      12             :       use m_juDFT
      13             :       contains
      14           0 :       subroutine wann_dipole2(
      15           0 :      >               jspins,pos,omtil,natd,
      16             :      >               l_nocosoc)
      17             :       use m_wann_readcenters
      18             :       implicit none
      19             :       integer,intent(in) :: jspins
      20             :       real,intent(in)    :: omtil
      21             :       integer,intent(in) :: natd
      22             :       real,intent(in)    :: pos(3,natd)
      23             :       logical,intent(in) :: l_nocosoc
      24             : 
      25             :       integer            :: jspin,jspins2
      26             :       logical            :: l_file
      27           0 :       integer,allocatable:: ind(:)
      28             :       integer            :: nwf,nwfs
      29           0 :       real,allocatable   :: distance(:,:)
      30           0 :       real               :: moment(3,jspins)
      31             :       real               :: moment2(3)
      32             :       real               :: elemchargmu,bohrtocm
      33             :       integer            :: shifted,ishift,ind1,ind2
      34           0 :       real,allocatable   :: centers(:,:)
      35           0 :       real,allocatable   :: spreads(:)
      36             : 
      37           0 :       call timestart("wann_dipole2")
      38           0 :       open(204,file='dipole_out')
      39           0 :       write(204,*)"*****************************************"
      40           0 :       write(204,*)"Calculation of the electric polarization."
      41           0 :       write(204,*)"*****************************************"
      42           0 :       elemchargmu=1.60217646e-13
      43           0 :       bohrtocm=0.529177e-8
      44             : 
      45           0 :       jspins2=jspins
      46           0 :       if(l_nocosoc)jspins2=1
      47           0 :       do jspin=1,jspins2
      48           0 :          inquire(file='proj',exist=l_file)
      49           0 :          IF(.NOT.l_file) CALL juDFT_error("proj not found",calledby
      50           0 :      +        ="wann_dipole2")
      51           0 :          open(203,file='proj',form='formatted')
      52           0 :          read(203,*)nwfs
      53           0 :          print*,"nwfs=",nwfs
      54           0 :          allocate(ind(nwfs))
      55           0 :          allocate(distance(3,nwfs))
      56           0 :          allocate(centers(3,nwfs))
      57           0 :          allocate(spreads(nwfs))
      58           0 :          do nwf=1,nwfs
      59           0 :           read(203,*)ind(nwf)
      60           0 :           read(203,*)
      61             :          enddo
      62           0 :          close(203)
      63           0 :          write(204,*)"spin=",jspin
      64           0 :          call wann_readcenters(nwfs,jspin,centers,spreads)
      65           0 :          write(204,*)"*****centers:******"
      66           0 :          do nwf=1,nwfs
      67           0 :            write(204,fmt=1000)nwf,centers(:,nwf),spreads(nwf)
      68             :          enddo
      69           0 :          moment(:,jspin)=0.0
      70           0 :          do nwf=1,nwfs
      71           0 :             distance(:,nwf)=pos(:,ind(nwf))-centers(:,nwf)
      72           0 :             moment(:,jspin)=moment(:,jspin)+distance(:,nwf)
      73             :          enddo
      74           0 :          if((jspins.eq.1).and..not.l_nocosoc) 
      75           0 :      &                    moment(:,jspin)=moment(:,jspin)*2
      76           0 :          moment(:,jspin)=moment(:,jspin)/omtil
      77           0 :          moment(:,jspin)=moment(:,jspin)*elemchargmu/((bohrtocm)**2)
      78           0 :          write(204,*)"polarization due to shift of wannierorbitals"
      79           0 :          write(204,*)moment(:,jspin),"uC/cm2"
      80           0 :          deallocate(ind,distance,centers,spreads)
      81             :       enddo!jspin
      82             : 1000  format(2x,'WF centre and spread', 
      83             :      &       i5,2x,'(',f10.6,',',f10.6,',',f10.6,' )',f15.8)
      84           0 :       inquire(file='dipole',exist=l_file)
      85           0 :       if(.not.l_file)then
      86           0 :             write(204,*)"no file dipole found"
      87             :       else
      88           0 :          moment2(:)=0.0
      89           0 :          open(224,file='dipole',form='formatted')
      90           0 :          read(224,*)shifted
      91           0 :          do ishift=1,shifted
      92           0 :                read(224,*)ind1,ind2 !electron taken from 1 and moved to 2
      93           0 :                moment2(:)=moment2(:)+pos(:,ind1)-pos(:,ind2)
      94             :          enddo
      95           0 :          close(224)
      96           0 :          moment2(:)=moment2(:)/omtil
      97           0 :          moment2(:)=moment2(:)*elemchargmu/(bohrtocm)**2
      98             :          write(204,*)"polarization due to electrons 
      99           0 :      &              being catched by other atoms"
     100           0 :          write(204,*)moment2(:),"uC/cm2"
     101             : 
     102             :       endif
     103           0 :       write(204,*)"total moment:"
     104           0 :       if(jspins2.eq.2)moment(:,1)=moment(:,1)+moment(:,2)
     105           0 :       if(l_file)then
     106           0 :          moment(:,1)=moment(:,1)+moment2(:)
     107             :       endif
     108           0 :       write(204,*)moment(:,1),"uC/cm2"
     109           0 :       close(204)
     110             : 
     111           0 :       call timestop("wann_dipole2")
     112           0 :       end subroutine wann_dipole2
     113             :       end module m_wann_dipole2