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
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_dipole(
11 0 : > jspins,omtil,natd,pos,
12 : > amat,
13 0 : > ntype,neq,zatom)
14 : c***************************************
15 : c Calculate electronic polarization.
16 : c Frank Freimuth
17 : c***************************************
18 : use m_wann_readcenters
19 : implicit none
20 : integer,intent(in) :: jspins
21 : real,intent(in) :: omtil
22 : integer,intent(in) :: natd
23 : real,intent(in) :: pos(3,natd)
24 : real,intent(in) :: amat(3,3)
25 : integer,intent(in) :: ntype
26 : integer,intent(in) :: neq(ntype)
27 : real,intent(in) :: zatom(ntype)
28 :
29 : integer :: nwfs,i
30 : integer :: num_atoms,num_wann,num_wann2
31 : integer :: j,k,ind,jj,jspin
32 0 : real,allocatable :: ioncharge(:)
33 0 : character(len=2),allocatable :: namat(:)
34 : character(len=2) :: symbol
35 0 : real,allocatable :: wann_centers1(:,:)
36 0 : real,allocatable :: wann_centers2(:,:)
37 : integer,allocatable :: wann_of_at(:)
38 : real :: polarization,charge
39 : integer :: num_symbols
40 : real :: ioni_polari(3)
41 : real :: shifted_polari(3)
42 : real :: size_polari
43 : real :: smallest_polari
44 : real :: final_polari(3)
45 : real :: electroni_polari(3,2)
46 : real :: elemchargmu,bohrtocm
47 : character*2 :: namat2(0:103)
48 : character(len=2) :: spin12(0:2)
49 : data spin12/' ', '.1', '.2'/
50 : character(len=6) :: filename
51 : logical :: l_file
52 :
53 : DATA namat2/'va',' H','He','Li','Be',
54 : + ' B',' C',' N',' O',' F','Ne',
55 : + 'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca','Sc','Ti',
56 : + ' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se',
57 : + 'Br','Kr','Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd',
58 : + 'Ag','Cd','In','Sn','Sb','Te',' I','Xe','Cs','Ba','La','Ce',
59 : + 'Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb',
60 : + 'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb',
61 : + 'Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa',' U','Np','Pu',
62 : + 'Am','Cm','Bk','Cf','Es','Fm','Md','No','Lw'/
63 :
64 0 : elemchargmu=1.60217646e-13
65 0 : bohrtocm=0.529177e-8
66 :
67 0 : open(666,file='polarization_out')
68 :
69 0 : num_atoms=0
70 0 : do j=1,ntype
71 0 : do k=1,neq(j)
72 0 : num_atoms=num_atoms+1
73 : enddo
74 : enddo
75 :
76 0 : allocate( namat(num_atoms) )
77 0 : ind=0
78 0 : do j=1,ntype
79 0 : do k=1,neq(j)
80 0 : ind=ind+1
81 0 : namat(ind)=namat2(nint(zatom(j)))
82 : enddo
83 : enddo
84 :
85 0 : do j=1,num_atoms
86 0 : print*,namat(j)," pos3=",pos(3,j)
87 : enddo
88 :
89 0 : allocate( ioncharge(num_atoms) )
90 0 : ioncharge=0.0
91 0 : open(400,file='IONS')
92 0 : read(400,*)num_symbols
93 0 : do j=1,num_symbols
94 0 : read(400,fmt=333)symbol,charge
95 0 : write(*,fmt=333)symbol,charge
96 0 : do k=1,num_atoms
97 0 : if(namat(k)==symbol)then
98 0 : ioncharge(k)=charge
99 : endif
100 : enddo
101 : enddo
102 0 : close(400)
103 : 333 format(a2,1x,f10.6)
104 :
105 0 : open(300,file='ioncharge')
106 0 : do j=1,num_atoms
107 0 : write(300,*)ioncharge(j)
108 : enddo
109 0 : close(300)
110 :
111 0 : open(300,file='ioncharge')
112 0 : do j=1,num_atoms
113 0 : read(300,*)ioncharge(j)
114 : enddo
115 0 : close(300)
116 :
117 0 : ioni_polari=0.0
118 0 : do j=1,num_atoms
119 : ioni_polari(:)=
120 0 : & ioni_polari(:)+ioncharge(j)*pos(:,j)
121 : enddo
122 :
123 0 : ioni_polari=ioni_polari/omtil*elemchargmu/((bohrtocm)**2)
124 :
125 0 : print*,"ioni_polari=",ioni_polari,"uC/cm2"
126 0 : write(666,*) "ioni_polari=",ioni_polari,"uC/cm2"
127 :
128 : c*************************************************
129 : c..reading the proj.1 / proj.2 / proj file
130 : c*************************************************
131 0 : do j=1,0,-1
132 0 : inquire(file=trim('proj'//spin12(j)),exist=l_file)
133 0 : if(l_file)then
134 0 : filename='proj'//spin12(j)
135 0 : exit
136 : endif
137 : enddo
138 0 : if(l_file)then
139 0 : open (203,file=trim(filename),status='old')
140 0 : rewind (203)
141 : else
142 : CALL juDFT_error("no proj/proj.1/proj.2",
143 0 : + calledby ="wann_dipole")
144 : endif
145 0 : read(203,*)num_wann
146 0 : close(203)
147 0 : print*,"number of wannier functions= ",num_wann
148 0 : allocate(wann_centers1(3,num_wann))
149 :
150 0 : if(jspins.eq.2)then
151 0 : do j=2,0,-1
152 0 : inquire(file=trim('proj'//spin12(j)),exist=l_file)
153 0 : if(l_file)then
154 0 : filename='proj'//spin12(j)
155 0 : exit
156 : endif
157 : enddo
158 0 : if(l_file)then
159 0 : open (203,file=trim(filename),status='old')
160 0 : rewind (203)
161 : else
162 : CALL juDFT_error("no proj/proj.1/proj.2",calledby
163 0 : + ="wann_dipole")
164 : endif
165 0 : read(203,*)num_wann2
166 0 : close(203)
167 0 : print*,"number of wannier functions spin 2= ",num_wann2
168 0 : allocate(wann_centers2(3,num_wann2))
169 : endif
170 :
171 0 : electroni_polari=0.0
172 : call wann_readcenters(
173 : > num_wann,1,
174 0 : < wann_centers1(:,:))
175 :
176 0 : do k=1,num_wann
177 : electroni_polari(:,1)=
178 0 : & electroni_polari(:,1)-wann_centers1(:,k)
179 : enddo
180 :
181 0 : if(jspins.eq.2)then
182 : call wann_readcenters(
183 : > num_wann2,2,
184 0 : < wann_centers2(:,:))
185 0 : do k=1,num_wann2
186 : electroni_polari(:,2)=
187 0 : & electroni_polari(:,2)-wann_centers2(:,k)
188 : enddo
189 : endif
190 :
191 :
192 0 : if(jspins.eq.1) electroni_polari = electroni_polari*2.0
193 : electroni_polari=
194 0 : & electroni_polari/omtil*elemchargmu/((bohrtocm)**2)
195 :
196 0 : do j=1,jspins
197 0 : print*,"spin ",j," electronic_polarization=",
198 0 : & electroni_polari(:,j)
199 0 : write(666,*)"spin ",j," electronic_polarization=",
200 0 : & electroni_polari(:,j)
201 : enddo
202 :
203 0 : ioni_polari(:)=ioni_polari(:)+electroni_polari(:,1)
204 0 : if(jspins.eq.2)ioni_polari(:)=ioni_polari(:)
205 0 : & +electroni_polari(:,2)
206 :
207 0 : print*,"total polarization=",ioni_polari(:)
208 0 : write(666,*)"total polarization=",ioni_polari(:)
209 : ! Check if the polarization may be reduced by adding
210 : ! contributions due to electrons shifted by primitive
211 : ! lattice translations.
212 0 : final_polari(:)=ioni_polari(:)
213 : size_polari=sqrt( (final_polari(1))**2 +
214 : + (final_polari(2))**2 +
215 0 : + (final_polari(3))**2 )
216 0 : smallest_polari=size_polari
217 0 : do i=-5,5
218 0 : do j=-5,5
219 0 : do k=-5,5
220 : shifted_polari(:)=ioni_polari(:)+
221 : + elemchargmu/((bohrtocm)**2)/omtil *
222 0 : * (amat(:,1)*k+amat(:,2)*j+amat(:,3)*i)
223 : size_polari=sqrt( (shifted_polari(1))**2 +
224 : + (shifted_polari(2))**2 +
225 0 : + (shifted_polari(3))**2 )
226 0 : if(size_polari.lt.smallest_polari)then
227 0 : final_polari=shifted_polari
228 0 : smallest_polari=size_polari
229 : endif
230 : enddo
231 : enddo
232 : enddo
233 0 : print*,"final polarization after shifting=",final_polari(:)
234 0 : write(666,*)"final polarization after shifting=",final_polari(:)
235 0 : close(666)
236 :
237 :
238 :
239 0 : end subroutine wann_dipole
240 : end module m_wann_dipole
|
