TABLE OF CONTENTS
ABINIT/elpolariz [ Functions ]
NAME
elpolariz
FUNCTION
Calculate corrections to total energy from polarising electric field with or without Berry phases (berryopt keyword)
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) cg(2,mcg)=planewave coefficients of wavefunctions cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables in this dataset gprimd(3,3)=reciprocal space dimensional primitive translations hdr <type(hdr_type)>=the header of wf, den and pot files kg(3,mpw*mkmem)=reduced planewave coordinates mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mkmem=number of k points treated by this node. mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw my_natom=number of atoms treated by current processor natom=number of atoms in cell nattyp(ntypat)= # atoms of each type. nkpt=number of k points npwarr(nkpt)=number of planewaves in basis at this k point nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms in unit cell nkpt=number of k-points option = 1: compute Berryphase polarization 2: compute finite difference expression of the ddk 3: compute polarization & ddk pawrhoij(my_natom*usepaw) <type(pawrhoij_type)> atomic occupancies pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data pel_cg(3) = reduced coordinates of the electronic polarization (a. u.) computed in the SCF loop pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates psps <type(pseudopotential_type)>=variables related to pseudopotentials pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points (see initberry.f) pwind_alloc = first dimension of pwind pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations rprimd(3,3)=dimensional real space primitive translations (bohr) ucvol=unit cell volume in bohr**3. usecprj=1 if cprj datastructure has been allocated xred(3,natom)=reduced atomic coordinates
OUTPUT
(see side effects)
SIDE EFFECTS
dtefield <type(efield_type)> = variables related to Berry phase and electric field calculations (see initberry.f). In case berryopt = 4/6/7/14/16/17, the overlap matrices computed in this routine are stored in dtefield%smat in order to be used in the electric field calculation. enefield=field energy etotal=total energy, might be correct by improved polarization computation pel(3) = reduced coordinates of the electronic polarization (a. u.) pion(3)= reduced coordinates of the ionic polarization (a. u.)
PARENTS
afterscfloop
CHILDREN
berryphase,berryphase_new,metric,uderiv,wrtout
SOURCE
79 #if defined HAVE_CONFIG_H 80 #include "config.h" 81 #endif 82 83 #include "abi_common.h" 84 85 86 subroutine elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,enefield,gprimd,hdr,& 87 & kg,mband,mcg,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,nkpt,& 88 & npwarr,nsppol,ntypat,pawrhoij,pawtab,& 89 & pel,pel_cg,pelev,pion,psps,pwind,pwind_alloc,& 90 & pwnsfac,rprimd,ucvol,usecprj,xred) 91 92 use defs_basis 93 use defs_datatypes 94 use defs_abitypes 95 use m_profiling_abi 96 use m_errors 97 use m_efield 98 99 use m_geometry, only : metric 100 use m_pawtab, only : pawtab_type 101 use m_pawrhoij, only : pawrhoij_type 102 use m_pawcprj, only : pawcprj_type 103 104 !This section has been created automatically by the script Abilint (TD). 105 !Do not modify the following lines by hand. 106 #undef ABI_FUNC 107 #define ABI_FUNC 'elpolariz' 108 use interfaces_14_hidewrite 109 use interfaces_67_common 110 !End of the abilint section 111 112 implicit none 113 114 !Arguments ------------------------------------ 115 !scalars 116 integer,intent(in) :: mband,mcg,mcprj,mkmem,mpw,my_natom,natom,nkpt,nsppol,ntypat 117 integer,intent(in) :: pwind_alloc,usecprj 118 real(dp),intent(in) :: ucvol 119 real(dp),intent(inout) :: enefield,etotal 120 type(MPI_type),intent(in) :: mpi_enreg 121 type(datafiles_type),intent(in) :: dtfil 122 type(dataset_type),intent(inout) :: dtset 123 type(efield_type),intent(inout) :: dtefield 124 type(hdr_type),intent(inout) :: hdr 125 type(pseudopotential_type),intent(in) :: psps 126 !arrays 127 integer,intent(in) :: atindx1(natom),kg(3,mpw*mkmem),nattyp(ntypat) 128 integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3) 129 real(dp),intent(in) :: cg(2,mcg),gprimd(3,3) 130 real(dp),intent(in) :: pel_cg(3),pwnsfac(2,pwind_alloc),rprimd(3,3) 131 real(dp),intent(inout) :: pel(3),pelev(3),pion(3),xred(3,natom) 132 type(pawcprj_type),intent(in) :: cprj(natom,mcprj*usecprj) 133 type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*psps%usepaw) 134 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 135 136 !Local variables------------------------------- 137 !scalars 138 integer :: my_nspinor,option,unit_out,iir,jjr,kkr 139 real(dp) :: pdif_mod,eenth,ucvol_local 140 character(len=500) :: message 141 !arrays 142 real(dp) :: gmet(3,3),gprimdlc(3,3),pdif(3),ptot(3),red_ptot(3),rmet(3,3) 143 !! ptot(3) = total polarization (not reduced) REC 144 !! red_ptot(3) = internal reduced total polarization REC 145 real(dp) :: A(3,3),A1(3,3),A_new(3,3),efield_new(3) 146 147 ! ************************************************************************* 148 149 DBG_ENTER("COLL") 150 151 if (usecprj==0.and.psps%usepaw==1) then 152 write (message,'(3a)')& 153 & 'cprj datastructure must be allocated !',ch10,& 154 & 'Action: change pawusecp input keyword.' 155 MSG_ERROR(message) 156 end if 157 158 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 159 160 if(dtset%berryopt>0 .and. dtset%berryopt/=4 .and. dtset%berryopt/=6 .and. dtset%berryopt/=7 .and. & 161 & dtset%berryopt/=14 .and. dtset%berryopt/=16 .and. dtset%berryopt/=17)then !!HONG 162 163 if (dtset%berryopt==1 .or. dtset%berryopt==3) then 164 call berryphase(atindx1,dtset%bdberry,cg,gprimd,dtset%istwfk,& 165 & dtset%kberry,kg,dtset%kptns,dtset%kptopt,dtset%kptrlatt,& 166 & mband,mcg,mkmem,mpw,natom,nattyp,dtset%nband,dtset%nberry,npwarr,& 167 & my_nspinor,nsppol,psps%ntypat,nkpt,rprimd,ucvol,& 168 & xred,psps%ziontypat) 169 end if 170 171 if (dtset%berryopt==2 .or. dtset%berryopt==3) then 172 call uderiv(dtset%bdberry,cg,gprimd,hdr,dtset%istwfk,& 173 & dtset%kberry,kg,dtset%kptns,dtset%kptopt,& 174 & dtset%kptrlatt,mband,mcg,mkmem,mpi_enreg,mpw,& 175 & natom,dtset%nband,dtset%nberry,npwarr,my_nspinor,nsppol,& 176 & nkpt,dtfil%unddk,dtfil%fnameabo_1wf) 177 end if 178 179 else if(dtset%berryopt<0 .or. dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. & 180 & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17)then !!HONG 181 182 select case (dtset%berryopt) 183 case (-5) 184 option = 2 185 case (-3) 186 option = 3 187 case (-2) 188 option = 2 189 case (-1) 190 option = 1 191 case (4) 192 option = 1 193 pel(:) = zero 194 pelev(:) = zero 195 case (6) !!HONG 196 option = 1 197 pel(:) = zero 198 pelev(:) = zero 199 case (7) !!HONG 200 option = 1 201 pel(:) = zero 202 pelev(:) = zero 203 case (14) !!HONG 204 option = 1 205 pel(:) = zero 206 pelev(:) = zero 207 case (16) !!HONG 208 option = 1 209 pel(:) = zero 210 pelev(:) = zero 211 case (17) !!HONG 212 option = 1 213 pel(:) = zero 214 pelev(:) = zero 215 end select 216 217 unit_out = ab_out 218 call berryphase_new(atindx1,cg,cprj,dtefield,dtfil,dtset,psps,& 219 & gprimd,hdr,psps%indlmn,kg,& 220 & psps%lmnmax,mband,mcg,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,npwarr,& 221 & nsppol,psps%ntypat,nkpt,option,pawrhoij,& 222 & pawtab,pel,pelev,pion,ptot,red_ptot,pwind,& !!REC 223 & pwind_alloc,pwnsfac,rprimd,dtset%typat,ucvol,& 224 & unit_out,usecprj,psps%usepaw,xred,psps%ziontypat) 225 226 dtefield%red_ptot1(:)=red_ptot(:) 227 228 if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. & 229 & dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt == 17 ) then !!HONG 230 231 ! Check if pel has the same value as pel_cg 232 ! if (psps%usepaw == 1) pel(:) = pel(:) + pelev(:) ! add on-site term for PAW 233 ! if (psps%usepaw == 1) red_ptot(:) = red_ptot(:) + pelev(:) ! add on-site term for PAW !! REC 234 ! above line suppressed because in the PAW case, pel already includes all on-site 235 ! terms and pelev should not be added in additionally. We are computing pelev separately for 236 ! reporting purposes only. 237 ! 13 June 2012 J Zwanziger 238 239 pdif(:) = pel_cg(:) - pel(:) 240 pdif_mod = pdif(1)**2 + pdif(2)**2 + pdif(3)**2 241 242 if (pdif_mod > tol8) then 243 write(message,'(11(a),e16.9)')ch10,& 244 & ' scfcv (electric field calculation) : WARNING -',ch10,& 245 & ' The difference between pel (electronic Berry phase updated ',ch10,& 246 & ' at each SCF cycle)',ch10,& 247 & ' and pel_cg (electronic Berryphase computed using the ',& 248 & 'berryphase routine) is',ch10,& 249 & ' pdif_mod = ',pdif_mod 250 call wrtout(std_out,message,'COLL') 251 write(message,'(a,6(a,e16.9,a))') ch10,& 252 & 'pel_cg(1) = ',pel_cg(1),ch10,& 253 & 'pel_cg(2) = ',pel_cg(2),ch10,& 254 & 'pel_cg(3) = ',pel_cg(3),ch10,& 255 & 'pel(1) = ',pel(1),ch10,& 256 & 'pel(2) = ',pel(2),ch10,& 257 & 'pel(3) = ',pel(3),ch10 258 MSG_ERROR(message) 259 end if 260 261 ! Use this (more accurate) value of P to recompute enefield 262 if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then !!HONG 263 etotal = etotal - enefield 264 265 enefield = -dot_product(dtset%red_efieldbar,red_ptot) 266 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 267 eenth = zero 268 do iir=1,3 269 do jjr=1,3 270 eenth= eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j 271 end do 272 end do 273 eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth 274 enefield=enefield+eenth 275 276 etotal = etotal + enefield 277 278 write(message,'(a,a)')ch10,& 279 & ' Stress tensor under a constant electric field:' 280 call wrtout(std_out,message,'COLL') 281 call wrtout(ab_out,message,'COLL') 282 283 end if 284 285 ! ! In finite D-field case, turn it into internal energy !!HONG 286 if (dtset%berryopt == 6 .or. dtset%berryopt == 16 ) then 287 etotal = etotal - enefield 288 289 enefield=zero 290 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 291 do iir=1,3 292 do jjr=1,3 293 enefield= enefield+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j 294 end do 295 end do 296 enefield= ucvol_local/(8.0d0*pi)*enefield 297 298 etotal = etotal + enefield 299 300 write(message,'(a,a)')ch10,& 301 & ' Stress tensor under a constant electric displacement field:' 302 call wrtout(std_out,message,'COLL') 303 call wrtout(ab_out,message,'COLL') 304 305 end if 306 307 ! HONG calculate internal energy and electric enthalpy for mixed BC case. 308 if ( dtset%berryopt == 17 ) then 309 etotal = etotal - enefield 310 enefield = zero 311 312 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 313 A(:,:)=(4*pi/ucvol_local)*rmet(:,:) 314 A1(:,:)=A(:,:) 315 A_new(:,:)=A(:,:) 316 efield_new(:)=dtset%red_efield(:) 317 eenth = zero 318 319 do kkr=1,3 320 if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction 321 ! step 1 add -ebar*p 322 eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr) 323 324 ! step 2 chang to e_new (change e to ebar) 325 efield_new(kkr)=dtset%red_efieldbar(kkr) 326 327 ! step 3 chang matrix A to A1 328 329 do iir=1,3 330 do jjr=1,3 331 if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr) 332 if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) & 333 & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr) 334 if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr) 335 end do 336 end do 337 338 A(:,:)=A1(:,:) 339 A_new(:,:)=A1(:,:) 340 end if 341 342 end do ! end fo kkr 343 344 345 do iir=1,3 346 do jjr=1,3 347 eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr) 348 end do 349 end do 350 351 enefield=eenth 352 etotal = etotal + enefield 353 354 write(message,'(a,a)')ch10,& 355 & ' Stress tensor under a constant (mixed) electric and electric displacement field:' 356 call wrtout(std_out,message,'COLL') 357 call wrtout(ab_out,message,'COLL') 358 359 end if ! berryopt==17 360 361 362 ! MVeithen: to clarify 363 ! Which stress tensor should be used in structural optimizations? 364 ! The one at constant electric field or at constant potential drop. 365 ! write(message,'(a,a)')ch10,& 366 ! & ' Stress tensor imposing a constant electric field:' 367 ! call wrtout(std_out,message,'COLL') 368 ! call wrtout(ab_out,message,'COLL') 369 370 end if ! dtset%berryopt == 4/6/7/14/16/17 371 372 end if ! dtset%berryopt>0 or dtset%berryopt/=4/6/7/14/16/17 373 374 DBG_EXIT("COLL") 375 376 end subroutine elpolariz