TABLE OF CONTENTS
ABINIT/stress [ Functions ]
NAME
stress
FUNCTION
Compute the stress tensor strten(i,j) = (1/ucvol)*d(Etot)/(d(eps(i,j))) where Etot is energy per unit cell, ucvol is the unstrained unit cell volume, r(i,iat) is the ith position of atom iat, and eps(i,j) is an infinitesimal strain which maps each point r to r(i) -> r(i) + Sum(j) [eps(i,j)*r(j)].
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, FJ, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx berryopt = 4/14: electric field is on -> add the contribution of the -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy = 6/16, or 7/17: electric displacement field is on -> add the contribution of the Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy from Etot(npw) data (at fixed geometry), used for making Pulay correction to stress tensor (hartree). Should be <=0. dtefield <type(efield_type)> = variables related to Berry phase eei=local pseudopotential part of Etot (hartree) efield = cartesian coordinates of the electric field in atomic units ehart=Hartree energy (hartree) eii=pseudoion core correction energy part of Etot (hartree) fock <type(fock_type)>= quantities to calculate Fock exact exchange gsqcut=cutoff value on G**2 for (large) sphere inside FFT box. gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2) ixc = choice of exchange-correlation functional kinstr(6)=kinetic energy part of stress tensor mgfft=maximum size of 1D FFTs mpi_enreg=informations about MPI parallelization mqgrid=dimensioned number of q grid points for local psp spline n1xccc=dimension of xccc1d ; 0 if no XC core correction is used n3xccc=dimension of the xccc3d array (0 or nfft). natom=number of atoms in cell nattyp(ntypat)=number of atoms of each type nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nlstr(6)=nonlocal part of stress tensor nspden=number of spin-density components nsym=number of symmetries in space group ntypat=number of types of atoms psps <type(pseudopotential_type)>=variables related to pseudopotentials pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) array prtvol=integer controlling volume of printed output qgrid(mqgrid)=q point array for local psp spline fits red_efieldbar(3) = efield in reduced units relative to reciprocal lattice rhog(2,nfft)=Fourier transform of charge density (bohr^-3) rprimd(3,3)=dimensional primitive translations in real space (bohr) strsxc(6)=xc correction to stress symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates typat(natom)=type integer for each atom in cell usefock=1 if fock operator is used; 0 otherwise. usepaw= 0 for non paw calculation; =1 for paw calculation vdw_tol= Van der Waals tolerance vdw_tol_3bt= Van der Waals tolerance on the 3-body term (only effective vdw_xc=6) vdw_xc= Van der Waals correction flag vlspl(mqgrid,2,ntypat)=local psp spline vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space vxc_hf(nfft,nspden)=exchange-correlation potential (hartree) in real space for Hartree-Fock corrections xccc1d(n1xccc*(1-usepaw),6,ntypat)=1D core charge function and five derivatives, for each type of atom, from psp (used in Norm-conserving only) xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type xred(3,natom)=reduced dimensionless atomic coordinates zion(ntypat)=valence charge of each type of atom znucl(ntypat)=atomic number of atom type
OUTPUT
strten(6)=components of the stress tensor (hartree/bohr^3) for the
6 unique components of this symmetric 3x3 tensor:
Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
The diagonal components of the returned stress tensor are
CORRECTED for the Pulay stress.
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
NOTES
* Concerning the stress tensor: See O. H. Nielsen and R. M. Martin, PRB 32, 3792 (1985). Note that first term in equation (2) should have minus sign (for kinetic energy contribution to stress tensor). Normalizations in this code differ somewhat from those employed by Nielsen and Martin. For the stress tensor contribution from the nonlocal Kleinman-Bylander separable pseudopotential, see D. M. Bylander, L. Kleinman, and S. Lee, PRB 42, 1394 (1990). Again normalization conventions differ somewhat. See Doug Allan s notes starting page 795 (13 Jan 1992). * This subroutine calls different subroutines to compute the stress tensor contributions from the following parts of the total energy: (1) kinetic energy, (2) exchange-correlation energy, (3) Hartree energy, (4) local pseudopotential energy, (5) pseudoion core correction energy, (6) nonlocal pseudopotential energy, (7) Ewald energy.
PARENTS
forstr
CHILDREN
atm2fft,ewald2,fourdp,metric,mkcore,mkcore_alt,mklocl_recipspace stresssym,strhar,timab,vdw_dftd2,vdw_dftd3,wrtout,zerosym
SOURCE
120 #if defined HAVE_CONFIG_H 121 #include "config.h" 122 #endif 123 124 #include "abi_common.h" 125 126 subroutine stress(atindx1,berryopt,dtefield,eei,efield,ehart,eii,fock,gsqcut,ixc,kinstr,& 127 & mgfft,mpi_enreg,mqgrid,n1xccc,n3xccc,natom,nattyp,& 128 & nfft,ngfft,nlstr,nspden,nsym,ntypat,paral_kgb,psps,pawrad,pawtab,ph1d,& 129 & prtvol,qgrid,red_efieldbar,rhog,rprimd,strten,strsxc,symrec,& 130 & typat,usefock,usepaw,vdw_tol,vdw_tol_3bt,vdw_xc,& 131 & vlspl,vxc,vxc_hf,xccc1d,xccc3d,xcccrc,xred,zion,znucl,qvpotzero,& 132 & electronpositron) ! optional argument 133 134 use defs_basis 135 use defs_abitypes 136 use m_efield 137 use m_profiling_abi 138 use m_errors 139 140 use m_geometry, only : metric 141 use m_fock, only : fock_type 142 use m_ewald, only : ewald2 143 use defs_datatypes, only : pseudopotential_type 144 use m_pawrad, only : pawrad_type 145 use m_pawtab, only : pawtab_type 146 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 147 use m_fft, only : zerosym 148 149 !This section has been created automatically by the script Abilint (TD). 150 !Do not modify the following lines by hand. 151 #undef ABI_FUNC 152 #define ABI_FUNC 'stress' 153 use interfaces_14_hidewrite 154 use interfaces_18_timing 155 use interfaces_41_geometry 156 use interfaces_53_ffts 157 use interfaces_56_xc 158 use interfaces_64_psp 159 use interfaces_67_common, except_this_one => stress 160 !End of the abilint section 161 162 implicit none 163 164 !Arguments ------------------------------------ 165 !scalars 166 integer,intent(in) :: berryopt,ixc,mgfft,mqgrid,n1xccc,n3xccc,natom,nfft,nspden 167 integer,intent(in) :: nsym,ntypat,paral_kgb,prtvol,usefock,usepaw,vdw_xc 168 real(dp),intent(in) :: eei,ehart,eii,gsqcut,vdw_tol,vdw_tol_3bt,qvpotzero 169 type(efield_type),intent(in) :: dtefield 170 type(pseudopotential_type),intent(in) :: psps 171 type(electronpositron_type),pointer,optional :: electronpositron 172 type(MPI_type),intent(in) :: mpi_enreg 173 type(fock_type),pointer, intent(inout) :: fock 174 !arrays 175 integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),symrec(3,3,nsym) 176 integer,intent(in) :: typat(natom) 177 real(dp),intent(in) :: efield(3),kinstr(6),nlstr(6) 178 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid) 179 real(dp),intent(in) :: red_efieldbar(3),rhog(2,nfft),rprimd(3,3),strsxc(6) 180 real(dp),intent(in) :: vlspl(mqgrid,2,ntypat),vxc(nfft,nspden) 181 real(dp),allocatable,intent(in) :: vxc_hf(:,:) 182 real(dp),intent(in) :: xccc1d(n1xccc*(1-usepaw),6,ntypat),xcccrc(ntypat) 183 real(dp),intent(in) :: xred(3,natom),zion(ntypat),znucl(ntypat) 184 real(dp),intent(inout) :: xccc3d(n3xccc) 185 real(dp),intent(out) :: strten(6) 186 type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw) 187 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 188 189 !Local variables------------------------------- 190 !scalars 191 integer :: coredens_method,iatom,icoulomb,idir,ii,ipositron,mu,optatm,optdyfr,opteltfr,opt_hybr,optgr,option 192 integer :: optn,optn2,optstr,optv,sdir,vloc_method 193 real(dp),parameter :: tol=1.0d-15 194 real(dp) :: e_dum,strsii,ucvol,vol_element 195 character(len=500) :: message 196 logical :: calc_epaw3_stress, efield_flag 197 !arrays 198 integer :: qprtrb_dum(3) 199 real(dp) :: corstr(6),ep3(3),epaws3red(6),ewestr(6),gmet(3,3) 200 !Maxwell-stress constribution, and magnitude of efield 201 real(dp) :: Maxstr(6),ModE 202 real(dp) :: gprimd(3,3),harstr(6),lpsstr(6),rmet(3,3),tsec(2),uncorr(3) 203 real(dp) :: vdwstr(6),vprtrb_dum(2) 204 real(dp) :: dummy_in(0) 205 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0) 206 real(dp),allocatable :: dummy(:),dyfr_dum(:,:,:),gr_dum(:,:),rhog_ep(:,:),v_dum(:) 207 real(dp),allocatable :: vxctotg(:,:) 208 character(len=10) :: EPName(1:2)=(/"Electronic","Positronic"/) 209 210 ! ************************************************************************* 211 212 call timab(37,1,tsec) 213 214 !Compute different geometric tensor, as well as ucvol, from rprimd 215 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 216 opt_hybr=0 217 if (allocated(vxc_hf)) opt_hybr=1 218 !======================================================================= 219 !========= Local pseudopotential and core charge contributions ========= 220 !======================================================================= 221 222 !Determine by which method the local ionic potential and/or the pseudo core charge density 223 ! contributions have to be computed 224 !Local ionic potential: 225 ! Method 1: PAW 226 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 227 vloc_method=1;if (usepaw==0) vloc_method=2 228 if (psps%usewvl==1) vloc_method=2 229 !Pseudo core charge density: 230 ! Method 1: PAW, nc_xccc_gspace 231 ! Method 2: Norm-conserving PP, wavelets 232 coredens_method=1;if (usepaw==0) coredens_method=2 233 if (psps%nc_xccc_gspace==1) coredens_method=1 234 if (psps%nc_xccc_gspace==0) coredens_method=2 235 if (psps%usewvl==1) coredens_method=2 236 237 !Local ionic potential and/or pseudo core charge by method 1 238 if (vloc_method==1.or.coredens_method==1) then 239 call timab(551,1,tsec) 240 optv=0;if (vloc_method==1) optv=1 241 optn=0;if (coredens_method==1) optn=n3xccc/nfft 242 optatm=0;optdyfr=0;opteltfr=0;optgr=0;optstr=1;optn2=1 243 if (coredens_method==1.and.n3xccc>0) then 244 ABI_ALLOCATE(v_dum,(nfft)) 245 ABI_ALLOCATE(vxctotg,(2,nfft)) 246 v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2)) 247 call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 248 call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),& 249 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 250 ABI_DEALLOCATE(v_dum) 251 else 252 ABI_ALLOCATE(vxctotg,(0,0)) 253 end if 254 call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,& 255 & dummy_out5,dummy_in,gmet,gprimd,dummy_out6,dummy_out7,gsqcut,& 256 & mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,& 257 & psps,pawtab,ph1d,qgrid,qprtrb_dum,rhog,corstr,lpsstr,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,vlspl,& 258 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 259 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 260 ABI_DEALLOCATE(vxctotg) 261 if (n3xccc==0.and.coredens_method==1) corstr=zero 262 call timab(551,2,tsec) 263 end if 264 265 !Local ionic potential by method 2 266 if (vloc_method==2) then 267 option=3 268 ABI_ALLOCATE(dyfr_dum,(3,3,natom)) 269 ABI_ALLOCATE(gr_dum,(3,natom)) 270 ABI_ALLOCATE(v_dum,(nfft)) 271 call mklocl_recipspace(dyfr_dum,eei,gmet,gprimd,gr_dum,gsqcut,lpsstr,mgfft,& 272 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,option,paral_kgb,ph1d,qgrid,& 273 & qprtrb_dum,rhog,ucvol,vlspl,vprtrb_dum,v_dum) 274 ABI_DEALLOCATE(dyfr_dum) 275 ABI_DEALLOCATE(gr_dum) 276 ABI_DEALLOCATE(v_dum) 277 end if 278 279 !Pseudo core electron density by method 2 280 if (coredens_method==2) then 281 if (n1xccc/=0) then 282 call timab(55,1,tsec) 283 option=3 284 ABI_ALLOCATE(dyfr_dum,(3,3,natom)) 285 ABI_ALLOCATE(gr_dum,(3,natom)) 286 ABI_ALLOCATE(v_dum,(nfft)) 287 icoulomb=0 ! not yet compatible with icoulomb 288 if (psps%usewvl==0.and.usepaw==0.and.icoulomb==0) then 289 if(opt_hybr==0) then 290 call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),& 291 & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc,& 292 & xcccrc,xccc1d,xccc3d,xred) 293 else 294 call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),& 295 & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc_hf,& 296 & xcccrc,xccc1d,xccc3d,xred) 297 end if 298 else if (psps%usewvl==0.and.(usepaw==1.or.icoulomb==1)) then 299 call mkcore_alt(atindx1,corstr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,& 300 & nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 301 & ucvol,vxc,xcccrc,xccc1d,xccc3d,xred,pawrad,pawtab,usepaw) 302 end if 303 ABI_DEALLOCATE(dyfr_dum) 304 ABI_DEALLOCATE(gr_dum) 305 ABI_DEALLOCATE(v_dum) 306 call timab(55,2,tsec) 307 else 308 corstr(:)=zero 309 end if 310 end if 311 312 !======================================================================= 313 !======================= Hartree energy contribution =================== 314 !======================================================================= 315 316 call strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd) 317 318 !======================================================================= 319 !======================= Ewald contribution ============================ 320 !======================================================================= 321 322 call timab(38,1,tsec) 323 call ewald2(gmet,natom,ntypat,rmet,rprimd,ewestr,typat,ucvol,xred,zion) 324 325 !======================================================================= 326 !================== VdW DFT-D contribution ============================ 327 !======================================================================= 328 329 if (vdw_xc==5) then 330 call vdw_dftd2(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_tol,& 331 & xred,znucl,str_vdw_dftd2=vdwstr) 332 elseif (vdw_xc==6.or.vdw_xc==7) then 333 call vdw_dftd3(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_xc,& 334 & vdw_tol,vdw_tol_3bt,xred,znucl,str_vdw_dftd3=vdwstr) 335 end if 336 337 call timab(38,2,tsec) 338 339 !HONG no Berry phase contribution if using reduced ebar or d according to 340 !HONG (PRL 89, 117602 (2002) Nature Physics: M. Stengel et.al. (2009)) 341 !======================================================================= 342 !=================== Berry phase contribution ========================== 343 !======================================================================= 344 345 !if (berryopt==4) then 346 !berrystr_tmp(:,:) = zero 347 !Diagonal: 348 !do mu = 1, 3 349 !do ii = 1, 3 350 !berrystr_tmp(mu,mu) = berrystr_tmp(mu,mu) - & 351 !& efield(mu)*rprimd(mu,ii)*(pel(ii) + pion(ii))/ucvol 352 !end do 353 !end do 354 !Off-diagonal (symmetrized before adding it to strten): 355 !do ii = 1, 3 356 !berrystr_tmp(3,2) = berrystr_tmp(3,2) & 357 !& - efield(3)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol 358 !berrystr_tmp(2,3) = berrystr_tmp(2,3) & 359 !& - efield(2)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol 360 !berrystr_tmp(3,1) = berrystr_tmp(3,1) & 361 !& - efield(3)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol 362 !berrystr_tmp(1,3) = berrystr_tmp(1,3) & 363 !& - efield(1)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol 364 !berrystr_tmp(2,1) = berrystr_tmp(2,1) & 365 !& - efield(2)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol 366 !berrystr_tmp(1,2) = berrystr_tmp(1,2) & 367 !& - efield(1)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol 368 !end do 369 !berrystr(1) = berrystr_tmp(1,1) 370 !berrystr(2) = berrystr_tmp(2,2) 371 !berrystr(3) = berrystr_tmp(3,3) 372 !berrystr(4) = (berrystr_tmp(3,2) + berrystr_tmp(2,3))/two 373 !berrystr(5) = (berrystr_tmp(3,1) + berrystr_tmp(1,3))/two 374 !berrystr(6) = (berrystr_tmp(2,1) + berrystr_tmp(1,2))/two 375 !end if 376 377 !======================================================================= 378 !================= Other (trivial) contributions ======================= 379 !======================================================================= 380 381 !Nonlocal part of stress has already been computed 382 !(in forstrnps(norm-conserving) or pawgrnl(PAW)) 383 384 !Kinetic part of stress has already been computed 385 !(in forstrnps) 386 387 !XC part of stress tensor has already been computed in "strsxc" 388 389 !ii part of stress (diagonal) is trivial! 390 strsii=-eii/ucvol 391 !qvpotzero is non zero, only when usepotzero=1 392 strsii=strsii+qvpotzero/ucvol 393 394 !====================================================================== 395 !HONG Maxwell stress when electric/displacement field is non-zero===== 396 !====================================================================== 397 efield_flag = (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 398 & berryopt==14 .or. berryopt==16 .or. berryopt==17) 399 calc_epaw3_stress = (efield_flag .and. usepaw == 1) 400 if ( efield_flag ) then 401 ModE=dot_product(efield,efield) 402 do ii=1,3 403 Maxstr(ii)=two*efield(ii)*efield(ii)-ModE 404 end do 405 Maxstr(4)=two*efield(3)*efield(2) 406 Maxstr(5)=two*efield(3)*efield(1) 407 Maxstr(6)=two*efield(2)*efield(1) 408 ! Converting to units of Ha/Bohr^3 409 ! Maxstr(:)=Maxstr(:)*e_Cb*Bohr_Ang*1.0d-10/(Ha_J*8.0d0*pi) 410 411 Maxstr(:)=Maxstr(:)*eps0*Ha_J*Bohr_Ang*1.0d-10/(8.0d0*pi*e_Cb**2) 412 413 write(message, '(a,a)' )ch10,& 414 & ' Cartesian components of Maxwell stress tensor (hartree/bohr^3)' 415 call wrtout(ab_out,message,'COLL') 416 call wrtout(std_out, message,'COLL') 417 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 418 & ' Maxstr(1 1)=',Maxstr(1),' Maxstr(3 2)=',Maxstr(4) 419 call wrtout(ab_out,message,'COLL') 420 call wrtout(std_out, message,'COLL') 421 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 422 & ' Maxstr(2 2)=',Maxstr(2),' Maxstr(3 1)=',Maxstr(5) 423 call wrtout(ab_out,message,'COLL') 424 call wrtout(std_out, message,'COLL') 425 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 426 & ' Maxstr(3 3)=',Maxstr(3),' Maxstr(2 1)=',Maxstr(6) 427 call wrtout(ab_out,message,'COLL') 428 call wrtout(std_out, message,'COLL') 429 write(message, '(a)' ) ' ' 430 call wrtout(ab_out,message,'COLL') 431 call wrtout(std_out, message,'COLL') 432 433 end if 434 435 ! compute additional F3-type stress due to projectors for electric field with PAW 436 if ( efield_flag .and. calc_epaw3_stress ) then 437 do sdir = 1, 6 438 ep3(:) = zero 439 do idir = 1, 3 440 vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir)) 441 do iatom = 1, natom 442 ep3(idir) = ep3(idir) + vol_element*dtefield%epaws3(iatom,idir,sdir) 443 end do ! end loop over atoms 444 end do ! end loop over idir (components of P) 445 ! note no appearance of ucvol here unlike in forces, stress definition includes 446 ! division by ucvol, which cancels the factor in -ucvol e . p 447 epaws3red(sdir) = -dot_product(red_efieldbar(1:3),ep3(1:3)) 448 end do 449 450 ! write(message, '(a,a)' )ch10,& 451 !& ' Cartesian components of PAW sigma_3 stress tensor (hartree/bohr^3)' 452 ! call wrtout(ab_out,message,'COLL') 453 ! call wrtout(std_out, message,'COLL') 454 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 455 !& ' epaws3red(1 1)=',epaws3red(1),' epaws3red(3 2)=',epaws3red(4) 456 ! call wrtout(ab_out,message,'COLL') 457 ! call wrtout(std_out, message,'COLL') 458 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 459 !& ' epaws3red(2 2)=',epaws3red(2),' epaws3red(3 1)=',epaws3red(5) 460 ! call wrtout(ab_out,message,'COLL') 461 ! call wrtout(std_out, message,'COLL') 462 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 463 !& ' epaws3red(3 3)=',epaws3red(3),' epaws3red(2 1)=',epaws3red(6) 464 ! call wrtout(ab_out,message,'COLL') 465 ! call wrtout(std_out, message,'COLL') 466 ! write(message, '(a)' ) ' ' 467 ! call wrtout(ab_out,message,'COLL') 468 ! call wrtout(std_out, message,'COLL') 469 470 end if 471 472 !======================================================================= 473 !===== Assemble the various contributions to the stress tensor ========= 474 !======================================================================= 475 !In cartesian coordinates (symmetric storage) 476 477 strten(:)=kinstr(:)+ewestr(:)+corstr(:)+strsxc(:)+harstr(:)+lpsstr(:)+nlstr(:) 478 479 if (usefock==1 .and. associated(fock).and.fock%fock_common%optstr) then 480 strten(:)=strten(:)+fock%fock_common%stress(:) 481 end if 482 483 !Add contributions for constant E or D calculation. 484 if ( efield_flag ) then 485 strten(:)=strten(:)+Maxstr(:) 486 if ( calc_epaw3_stress ) strten(:) = strten(:) + epaws3red(:) 487 end if 488 if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)+vdwstr(:) 489 490 !Additional stuff for electron-positron 491 ipositron=0 492 if (present(electronpositron)) then 493 if (associated(electronpositron)) then 494 if (allocated(electronpositron%stress_ep)) ipositron=electronpositron_calctype(electronpositron) 495 end if 496 end if 497 if (abs(ipositron)==1) then 498 strten(:)=strten(:)-harstr(:)-ewestr(:)-corstr(:)-lpsstr(:) 499 harstr(:)=zero;ewestr(:)=zero;corstr(:)=zero;strsii=zero 500 lpsstr(:)=-lpsstr(:);lpsstr(1:3)=lpsstr(1:3)-two*eei/ucvol 501 strten(:)=strten(:)+lpsstr(:) 502 if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)-vdwstr(:) 503 if (vdw_xc>=5.and.vdw_xc<=7) vdwstr(:)=zero 504 end if 505 if (abs(ipositron)==2) then 506 ABI_ALLOCATE(rhog_ep,(2,nfft)) 507 ABI_ALLOCATE(dummy,(6)) 508 call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 509 rhog_ep=-rhog_ep 510 call strhar(electronpositron%e_hartree,gsqcut,dummy,mpi_enreg,nfft,ngfft,rhog_ep,rprimd) 511 strten(:)=strten(:)+dummy(:);harstr(:)=harstr(:)+dummy(:) 512 ABI_DEALLOCATE(rhog_ep) 513 ABI_DEALLOCATE(dummy) 514 end if 515 if (ipositron>0) strten(:)=strten(:)+electronpositron%stress_ep(:) 516 517 !Symmetrize resulting tensor if nsym>1 518 if (nsym>1) then 519 call stresssym(gprimd,nsym,strten,symrec) 520 end if 521 522 !Set to zero very small values of stress 523 do mu=1,6 524 if (abs(strten(mu))<tol) strten(mu)=zero 525 end do 526 527 !Include diagonal terms, save uncorrected stress for output 528 do mu=1,3 529 uncorr(mu)=strten(mu)+strsii 530 strten(mu)=uncorr(mu) 531 end do 532 533 !======================================================================= 534 !================ Print out info about stress tensor =================== 535 !======================================================================= 536 if (prtvol>=10.and.ipositron>=0) then 537 write(message, '(a)' ) ' ' 538 call wrtout(std_out,message,'COLL') 539 do mu=1,6 540 write(message, '(a,i5,a,1p,e22.12)' )& 541 & ' stress: component',mu,' of hartree stress is',harstr(mu) 542 call wrtout(std_out,message,'COLL') 543 end do 544 write(message, '(a)' ) ' ' 545 call wrtout(std_out,message,'COLL') 546 do mu=1,6 547 write(message, '(a,i5,a,1p,e22.12)' )& 548 & ' stress: component',mu,' of loc psp stress is',lpsstr(mu) 549 call wrtout(std_out,message,'COLL') 550 end do 551 write(message, '(a)' ) ' ' 552 call wrtout(std_out,message,'COLL') 553 do mu=1,6 554 write(message, '(a,i5,a,1p,e22.12)' )& 555 & ' stress: component',mu,& 556 & ' of kinetic stress is',kinstr(mu) 557 call wrtout(std_out,message,'COLL') 558 end do 559 write(message, '(a)' ) ' ' 560 call wrtout(std_out,message,'COLL') 561 do mu=1,6 562 write(message, '(a,i5,a,1p,e22.12)' )& 563 & ' stress: component',mu,' of nonlocal ps stress is',nlstr(mu) 564 call wrtout(std_out,message,'COLL') 565 end do 566 write(message, '(a)' ) ' ' 567 call wrtout(std_out,message,'COLL') 568 do mu=1,6 569 write(message, '(a,i5,a,1p,e22.12)' )& 570 & ' stress: component',mu,' of core xc stress is',corstr(mu) 571 call wrtout(std_out,message,'COLL') 572 end do 573 write(message, '(a)' ) ' ' 574 call wrtout(std_out,message,'COLL') 575 do mu=1,6 576 write(message, '(a,i5,a,1p,e22.12)' )& 577 & ' stress: component',mu,& 578 & ' of Ewald energ stress is',ewestr(mu) 579 call wrtout(std_out,message,'COLL') 580 end do 581 write(message, '(a)' ) ' ' 582 call wrtout(std_out,message,'COLL') 583 do mu=1,6 584 write(message, '(a,i5,a,1p,e22.12)' ) & 585 & ' stress: component',mu,' of xc stress is',strsxc(mu) 586 call wrtout(std_out,message,'COLL') 587 end do 588 if (vdw_xc>=5.and.vdw_xc<=7) then 589 write(message, '(a)' ) ' ' 590 call wrtout(std_out,message,'COLL') 591 do mu=1,6 592 write(message, '(a,i5,a,1p,e22.12)' )& 593 & ' stress: component',mu,& 594 & ' of VdW DFT-D stress is',vdwstr(mu) 595 call wrtout(std_out,message,'COLL') 596 end do 597 end if 598 write(message, '(a)' ) ' ' 599 call wrtout(std_out,message,'COLL') 600 write(message, '(a,1p,e22.12)' ) & 601 & ' stress: ii (diagonal) part is',strsii 602 call wrtout(std_out,message,'COLL') 603 if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 604 & berryopt==14 .or. berryopt==16 .or. berryopt==17) then !!HONG 605 write(message, '(a)' ) ' ' 606 call wrtout(std_out,message,'COLL') 607 do mu = 1, 6 608 write(message, '(a,i2,a,1p,e22.12)' )& 609 & ' stress: component',mu,' of Maxwell stress is',& 610 & Maxstr(mu) 611 call wrtout(std_out,message,'COLL') 612 end do 613 end if 614 if (ipositron/=0) then 615 write(message, '(a)' ) ' ' 616 call wrtout(std_out,message,'COLL') 617 do mu=1,6 618 write(message, '(a,i5,3a,1p,e22.12)' ) & 619 & ' stress: component',mu,' of ',EPName(abs(ipositron)), & 620 & ' stress is',electronpositron%stress_ep(mu) 621 call wrtout(std_out,message,'COLL') 622 end do 623 end if 624 625 end if ! prtvol 626 if (ipositron>=0) then 627 write(message, '(a,a)' )ch10,& 628 & ' Cartesian components of stress tensor (hartree/bohr^3)' 629 call wrtout(ab_out,message,'COLL') 630 call wrtout(std_out, message,'COLL') 631 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 632 & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4) 633 call wrtout(ab_out,message,'COLL') 634 call wrtout(std_out, message,'COLL') 635 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 636 & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5) 637 call wrtout(ab_out,message,'COLL') 638 call wrtout(std_out, message,'COLL') 639 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 640 & ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6) 641 call wrtout(ab_out,message,'COLL') 642 call wrtout(std_out, message,'COLL') 643 write(message, '(a)' ) ' ' 644 call wrtout(ab_out,message,'COLL') 645 call wrtout(std_out, message,'COLL') 646 end if 647 call timab(37,2,tsec) 648 649 end subroutine stress