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_ionic
8 : use m_judft
9 : contains
10 0 : subroutine wann_dipole_ionic(
11 0 : > natd,pos,omtil,
12 0 : > amat,taual,ntype,
13 0 : > neq,zatom,l_absolute,
14 : > l_invsubtract,
15 : < ionic_moment)
16 : use m_wann_ioncharge_gen
17 : implicit none
18 : real,intent(in) :: omtil
19 : integer,intent(in) :: natd
20 : real,intent(in) :: pos(3,natd)
21 : real,intent(in) :: amat(3,3)
22 : real,intent(in) :: taual(3,natd)
23 : integer,intent(in) :: ntype
24 : integer,intent(in) :: neq(ntype)
25 : real,intent(in) :: zatom(ntype)
26 : logical,intent(in) :: l_absolute
27 : logical,intent(in) :: l_invsubtract
28 : real,intent(out) :: ionic_moment(3)
29 :
30 : integer :: i
31 : integer :: num_atoms
32 : integer :: j,k,ind,jj,jspin
33 0 : real,allocatable :: ioncharge(:)
34 : real :: polarization
35 : real :: ionic_polarization(3)
36 : real,parameter :: elemchargmu=1.60217646e-13
37 : real,parameter :: bohrtocm=0.529177e-8
38 : character(len=6) :: filename
39 : logical :: l_file
40 0 : real :: pos_inv(3,natd)
41 0 : real :: taual_inv(3,natd)
42 : real :: coordinate
43 : real :: ionchargesum
44 :
45 0 : call timestart("wann_dipole_ionic")
46 0 : write(666,*)"*****************************"
47 0 : write(666,*)" Ionic terms "
48 0 : write(666,*)"*****************************"
49 0 : write(*,*) "*****************************"
50 0 : write(*,*) " Ionic terms "
51 0 : write(*,*) "*****************************"
52 :
53 : c-----count atoms
54 0 : num_atoms=0
55 0 : do j=1,ntype
56 0 : do k=1,neq(j)
57 0 : num_atoms=num_atoms+1
58 : enddo
59 : enddo
60 :
61 : c-----get charges of ions
62 0 : allocate( ioncharge(num_atoms) )
63 : call wann_ioncharge_gen(
64 : > num_atoms,ntype,natd,
65 : > neq,zatom,pos,
66 0 : < ioncharge)
67 0 : ionchargesum=sum(ioncharge(1:num_atoms))
68 0 : write(666,*)"ionchargesum=",ionchargesum
69 :
70 : c-----Try to find the atomic positions of the inversion symmetric counterpart.
71 : c-----The ferroelectric system is assumed to be connected to an inversion
72 : c-----symmetric system by an adiabatic path. The polarization of the inversion
73 : c-----symmetric system is zero and consequently ambiguities in the evaluation of
74 : c-----the polarization can be removed by subtracting the calculated polarizations
75 : c-----of the two systems. We do not care here, whether the calculated inversion
76 : c-----symmetric positions are physically reasonable: For example two atoms might
77 : c-----end up at the same location. But this does not matter.
78 0 : pos_inv=0.0
79 0 : if(l_invsubtract)then
80 0 : do j=1,num_atoms
81 0 : do k=1,3
82 0 : coordinate=taual(k,j)
83 0 : if(abs(coordinate).lt.0.125)then
84 0 : taual_inv(k,j)=0.0
85 0 : elseif(abs(coordinate).lt.0.375)then
86 0 : taual_inv(k,j)=0.25*coordinate/abs(coordinate)
87 0 : elseif(abs(coordinate).lt.0.625)then
88 0 : taual_inv(k,j)=0.5*coordinate/abs(coordinate)
89 0 : elseif(abs(coordinate).lt.0.875)then
90 0 : taual_inv(k,j)=0.75*coordinate/abs(coordinate)
91 : else
92 0 : taual_inv(k,j)=coordinate/abs(coordinate)
93 : endif
94 : enddo
95 : enddo
96 0 : do j=1,num_atoms
97 0 : pos_inv(:,j)=matmul(amat,taual_inv(:,j))
98 : enddo
99 0 : do j=1,num_atoms
100 0 : print*,"old position:"
101 0 : print*,pos(:,j)
102 0 : print*,"new position:"
103 0 : print*,pos_inv(:,j)
104 0 : print*,"************************"
105 : enddo
106 : endif
107 :
108 : c-----compute the ionic moment
109 0 : ionic_moment=0.0
110 0 : do j=1,num_atoms
111 : ionic_moment(:)=
112 0 : & ionic_moment(:)+ioncharge(j)*(pos(:,j)-pos_inv(:,j))
113 : enddo
114 :
115 : ionic_polarization =
116 0 : = ionic_moment/omtil*elemchargmu/((bohrtocm)**2)
117 :
118 0 : write(*, fmt=555)ionic_moment(:)
119 0 : write(666,fmt=555)ionic_moment(:)
120 0 : if(.not.l_absolute)then
121 0 : write(*, fmt=777)ionic_polarization(:)
122 0 : write(666,fmt=777)ionic_polarization(:)
123 : endif
124 :
125 : 555 format("ionic moment = (",3f12.6,")a.u.")
126 : 777 format("ionic polarization = (",3f12.6,")uC/cm**2")
127 :
128 0 : call timestop("wann_dipole_ionic")
129 0 : end subroutine wann_dipole_ionic
130 : end module m_wann_dipole_ionic
131 :
132 :
133 :
134 :
135 :
136 :
137 :
138 :
139 :
140 :
|
