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