wann_dipole_ionic.f 4.9 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 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
      module m_wann_dipole_ionic
      contains
      subroutine wann_dipole_ionic(
     >               natd,pos,omtil,
     >               amat,taual,ntype,
     >               neq,zatom,l_absolute,
     >               l_invsubtract,
     <               ionic_moment)
      use m_wann_ioncharge_gen
      implicit none
      real,intent(in)              :: omtil
      integer,intent(in)           :: natd
      real,intent(in)              :: pos(3,natd)
      real,intent(in)              :: amat(3,3)
      real,intent(in)              :: taual(3,natd)
      integer,intent(in)           :: ntype
      integer,intent(in)           :: neq(ntype)
      real,intent(in)              :: zatom(ntype)
      logical,intent(in)           :: l_absolute
      logical,intent(in)           :: l_invsubtract
      real,intent(out)             :: ionic_moment(3)

      integer                      :: i
      integer                      :: num_atoms
      integer                      :: j,k,ind,jj,jspin
      real,allocatable             :: ioncharge(:)
      real                         :: polarization
      real                         :: ionic_polarization(3)
      real,parameter               :: elemchargmu=1.60217646e-13
      real,parameter               :: bohrtocm=0.529177e-8
      character(len=6)             :: filename
      logical                      :: l_file
      real                         :: pos_inv(3,natd)
      real                         :: taual_inv(3,natd)
      real                         :: coordinate
      real                         :: ionchargesum

      write(666,*)"*****************************"
      write(666,*)" Ionic terms                 "
      write(666,*)"*****************************"
      write(*,*)  "*****************************"
      write(*,*)  " Ionic terms                 "
      write(*,*)  "*****************************"

c-----count atoms
      num_atoms=0
      do j=1,ntype
         do k=1,neq(j)
            num_atoms=num_atoms+1
         enddo
      enddo

c-----get charges of ions
      allocate( ioncharge(num_atoms) )
      call wann_ioncharge_gen(
     >         num_atoms,ntype,natd,
     >         neq,zatom,pos,
     <         ioncharge)
      ionchargesum=sum(ioncharge(1:num_atoms))
      write(666,*)"ionchargesum=",ionchargesum

c-----Try to find the atomic positions of the inversion symmetric counterpart.
c-----The ferroelectric system is assumed to be connected to an inversion 
c-----symmetric system by an adiabatic path. The polarization of the inversion
c-----symmetric system is zero and consequently ambiguities in the evaluation of
c-----the polarization can be removed by subtracting the calculated polarizations
c-----of the two systems. We do not care here, whether the calculated inversion
c-----symmetric positions are physically reasonable: For example two atoms might
c-----end up at the same location. But this does not matter.
      pos_inv=0.0
      if(l_invsubtract)then
        do j=1,num_atoms
          do k=1,3
               coordinate=taual(k,j)
               if(abs(coordinate).lt.0.125)then
                  taual_inv(k,j)=0.0
               elseif(abs(coordinate).lt.0.375)then
                  taual_inv(k,j)=0.25*coordinate/abs(coordinate)
               elseif(abs(coordinate).lt.0.625)then
                  taual_inv(k,j)=0.5*coordinate/abs(coordinate)
               elseif(abs(coordinate).lt.0.875)then
                  taual_inv(k,j)=0.75*coordinate/abs(coordinate)
               else
                  taual_inv(k,j)=coordinate/abs(coordinate)
               endif   
          enddo
        enddo
        do j=1,num_atoms
           pos_inv(:,j)=matmul(amat,taual_inv(:,j))
        enddo
        do j=1,num_atoms
          print*,"old position:"
          print*,pos(:,j)
          print*,"new position:"
          print*,pos_inv(:,j)
          print*,"************************"
        enddo
      endif

c-----compute the ionic moment
      ionic_moment=0.0
      do j=1,num_atoms
         ionic_moment(:)=
     &        ionic_moment(:)+ioncharge(j)*(pos(:,j)-pos_inv(:,j))
      enddo

      ionic_polarization =
     =   ionic_moment/omtil*elemchargmu/((bohrtocm)**2)

      write(*,  fmt=555)ionic_moment(:)
      write(666,fmt=555)ionic_moment(:)
      if(.not.l_absolute)then
        write(*,  fmt=777)ionic_polarization(:)
        write(666,fmt=777)ionic_polarization(:)
      endif

 555  format("ionic moment = (",3f12.6,")a.u.")
 777  format("ionic polarization = (",3f12.6,")uC/cm**2")

      end subroutine wann_dipole_ionic
      end module m_wann_dipole_ionic