TABLE OF CONTENTS
ABINIT/m_stress [ Modules ]
NAME
m_stress
FUNCTION
COPYRIGHT
Copyright (C) 1998-2022 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_stress 23 24 use defs_basis 25 use m_efield 26 use m_abicore 27 use m_errors 28 use m_xmpi 29 use m_extfpmd 30 31 use defs_abitypes, only : MPI_type 32 use m_time, only : timab 33 use m_geometry, only : metric, stresssym 34 use m_fock, only : fock_type 35 use m_ewald, only : ewald2 36 use defs_datatypes, only : pseudopotential_type 37 use m_pawrad, only : pawrad_type 38 use m_pawtab, only : pawtab_type 39 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 40 use m_fft, only : zerosym, fourdp 41 use m_mpinfo, only : ptabs_fourdp 42 use m_vdw_dftd2, only : vdw_dftd2 43 use m_vdw_dftd3, only : vdw_dftd3 44 use m_atm2fft, only : atm2fft 45 use m_mklocl, only : mklocl_recipspace 46 use m_mkcore, only : mkcore, mkcore_alt 47 48 implicit none 49 50 private
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)].
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=information 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) strscondft(6)=cDFT correction to stress 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. usekden=1 is kinetic energy density has to be taken into account, 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 vxctau(nfft,nspden,4*usekden)=(only for meta-GGA) derivative of XC energy density wrt kinetic energy density (depsxcdtau) 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 xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction 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) [[cite:Nielsen1985a]]. 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) [[cite:Bylander1990]]. 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.
SOURCE
168 subroutine stress(atindx1,berryopt,dtefield,eei,efield,ehart,eii,fock,gsqcut,extfpmd,& 169 & ixc,kinstr,mgfft,mpi_enreg,mqgrid,n1xccc,n3xccc,natom,nattyp,& 170 & nfft,ngfft,nlstr,nspden,nsym,ntypat,psps,pawrad,pawtab,ph1d,& 171 & prtvol,qgrid,red_efieldbar,rhog,rprimd,strten,strscondft,strsxc,symrec,& 172 & typat,usefock,usekden,usepaw,vdw_tol,vdw_tol_3bt,vdw_xc,& 173 & vlspl,vxc,vxctau,vxc_hf,xccc1d,xccc3d,xcctau3d,xcccrc,xred,zion,znucl,qvpotzero,& 174 & electronpositron) ! optional argument 175 176 !Arguments ------------------------------------ 177 !scalars 178 integer,intent(in) :: berryopt,ixc,mgfft,mqgrid,n1xccc,n3xccc,natom,nfft,nspden 179 integer,intent(in) :: nsym,ntypat,prtvol,usefock,usekden,usepaw,vdw_xc 180 real(dp),intent(in) :: eei,ehart,eii,gsqcut,vdw_tol,vdw_tol_3bt,qvpotzero 181 type(efield_type),intent(in) :: dtefield 182 type(extfpmd_type),pointer,intent(inout) :: extfpmd 183 type(pseudopotential_type),intent(in) :: psps 184 type(electronpositron_type),pointer,optional :: electronpositron 185 type(MPI_type),intent(in) :: mpi_enreg 186 type(fock_type),pointer, intent(inout) :: fock 187 !arrays 188 integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),symrec(3,3,nsym) 189 integer,intent(in) :: typat(natom) 190 real(dp),intent(in) :: efield(3),kinstr(6),nlstr(6) 191 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid) 192 real(dp),intent(in) :: red_efieldbar(3),rhog(2,nfft),strscondft(6),strsxc(6) 193 real(dp),intent(in) :: vlspl(mqgrid,2,ntypat),vxc(nfft,nspden),vxctau(nfft,nspden,4*usekden) 194 real(dp),allocatable,intent(in) :: vxc_hf(:,:) 195 real(dp),intent(in) :: xccc1d(n1xccc*(1-usepaw),6,ntypat),xcccrc(ntypat) 196 real(dp),intent(in) :: xred(3,natom),zion(ntypat),znucl(ntypat) 197 real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*usekden),rprimd(3,3) 198 real(dp),intent(out) :: strten(6) 199 type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw) 200 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 201 202 !Local variables------------------------------- 203 !scalars 204 integer :: coredens_method,coretau_method,iatom,icoulomb,idir,ii,ipositron,mu,nkpt=1 205 integer :: optatm,optdyfr,opteltfr,opt_hybr,optgr,option,optn,optn2,optstr,optv,sdir,vloc_method 206 real(dp),parameter :: tol=1.0d-15 207 real(dp) :: e_dum,dum_rcut=zero,strsii,ucvol,vol_element 208 character(len=500) :: message 209 logical :: calc_epaw3_stress, efield_flag 210 !arrays 211 integer :: qprtrb_dum(3),icutcoul=3 212 real(dp) :: corstr(6),dumstr(6),ep3(3),epaws3red(6),ewestr(6),gmet(3,3),vcutgeo(3) 213 real(dp) :: gprimd(3,3),harstr(6),lpsstr(6),rmet(3,3),taustr(6),tsec(2),uncorr(3) 214 real(dp) :: vdwstr(6),vprtrb_dum(2) 215 real(dp) :: Maxstr(6),ModE !Maxwell-stress constribution, and magnitude of efield 216 real(dp) :: dummy_in(0) 217 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0) 218 real(dp),allocatable :: dummy(:),dyfr_dum(:,:,:),gr_dum(:,:),rhog_ep(:,:),v_dum(:) 219 real(dp),allocatable :: vxctotg(:,:) 220 character(len=10) :: EPName(1:2)=(/"Electronic","Positronic"/) 221 222 ! ************************************************************************* 223 224 call timab(37,1,tsec) 225 226 !Compute different geometric tensor, as well as ucvol, from rprimd 227 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 228 229 opt_hybr=0;if (allocated(vxc_hf)) opt_hybr=1 230 icoulomb=0 ! not yet compatible with icoulomb 231 232 !======================================================================= 233 !========= Local pseudopotential and core charge contributions ========= 234 !======================================================================= 235 236 !Determine by which method the local ionic potential and/or the pseudo core charge density 237 ! contributions have to be computed 238 !Local ionic potential: 239 ! Method 1: PAW 240 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 241 vloc_method=1;if (usepaw==0) vloc_method=2 242 if (psps%usewvl==1) vloc_method=2 243 !Pseudo core charge density: 244 ! Method 1: PAW, nc_xccc_gspace 245 ! Method 2: Norm-conserving PP, wavelets 246 coredens_method=1;if (usepaw==0) coredens_method=2 247 if (psps%nc_xccc_gspace==1) coredens_method=1 248 if (psps%nc_xccc_gspace==0) coredens_method=2 249 if (psps%usewvl==1) coredens_method=2 250 coretau_method=0 251 if (usekden==1.and.psps%usepaw==1) then 252 coretau_method=1;if (psps%nc_xccc_gspace==0) coretau_method=2 253 end if 254 255 !Local ionic potential and/or pseudo core charge by method 1 256 if (vloc_method==1.or.coredens_method==1.or.coretau_method==1) then 257 call timab(551,1,tsec) 258 ! Compute Vxc in reciprocal space 259 if (coredens_method==1.and.n3xccc>0) then 260 ABI_MALLOC(v_dum,(nfft)) 261 ABI_MALLOC(vxctotg,(2,nfft)) 262 v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2)) 263 call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0) 264 call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),& 265 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 266 ABI_FREE(v_dum) 267 else 268 ABI_MALLOC(vxctotg,(0,0)) 269 end if 270 ! Compute contribution to stresses from Vloc and/or pseudo core density 271 optv=0;if (vloc_method==1) optv=1 272 optn=0;if (coredens_method==1) optn=n3xccc/nfft 273 optatm=0;optdyfr=0;opteltfr=0;optgr=0;optstr=1;optn2=1 274 if (vloc_method==1.or.coredens_method==1) then 275 call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,& 276 & dummy_out5,dummy_in,gmet,gprimd,dummy_out6,dummy_out7,gsqcut,& 277 & mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,& 278 & psps,pawtab,ph1d,qgrid,qprtrb_dum,dum_rcut,rhog,rprimd,corstr,lpsstr,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,vlspl,& 279 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 280 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 281 end if 282 if (n3xccc==0.and.coredens_method==1) corstr=zero 283 ABI_FREE(vxctotg) 284 if (usekden==1.and.coretau_method==1..and.n3xccc>0) then 285 ! Compute contribution to stresses from pseudo kinetic energy core density 286 optv=0;optn=1;optn2=4 287 ABI_MALLOC(v_dum,(nfft)) 288 ABI_MALLOC(vxctotg,(2,nfft)) 289 v_dum(:)=vxctau(:,1,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxctau(:,2,1)) 290 call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0) 291 call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),& 292 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 293 ABI_FREE(v_dum) 294 call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,& 295 & dummy_out5,dummy_in,gmet,gprimd,dummy_out6,dummy_out7,gsqcut,& 296 & mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,& 297 & psps,pawtab,ph1d,qgrid,qprtrb_dum,dum_rcut,rhog,rprimd,taustr,dumstr,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,vlspl,& 298 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 299 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 300 corstr(1:6)=corstr(1:6)+taustr(1:6) 301 end if 302 call timab(551,2,tsec) 303 end if 304 305 !Local ionic potential by method 2 306 if (vloc_method==2) then 307 option=3 308 ABI_MALLOC(dyfr_dum,(3,3,natom)) 309 ABI_MALLOC(gr_dum,(3,natom)) 310 ABI_MALLOC(v_dum,(nfft)) 311 call mklocl_recipspace(dyfr_dum,eei,gmet,gprimd,gr_dum,gsqcut,icutcoul,lpsstr,mgfft,& 312 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,nkpt,ntypat,option,ph1d,qgrid,& 313 & qprtrb_dum,dum_rcut,rhog,rprimd,ucvol,vcutgeo,vlspl,vprtrb_dum,v_dum) 314 ABI_FREE(dyfr_dum) 315 ABI_FREE(gr_dum) 316 ABI_FREE(v_dum) 317 end if 318 319 !Pseudo core electron density by method 2 320 if (coredens_method==2.or.coretau_method==2) then 321 if (n1xccc/=0) then 322 call timab(55,1,tsec) 323 option=3 324 ABI_MALLOC(dyfr_dum,(3,3,natom)) 325 ABI_MALLOC(gr_dum,(3,natom)) 326 ABI_MALLOC(v_dum,(nfft)) 327 if (coredens_method==2) then 328 if (psps%usewvl==0.and.usepaw==0.and.icoulomb==0) then 329 if(opt_hybr==0) then 330 call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),& 331 & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc,& 332 & xcccrc,xccc1d,xccc3d,xred) 333 else 334 call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),& 335 & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc_hf,& 336 & xcccrc,xccc1d,xccc3d,xred) 337 end if 338 else if (psps%usewvl==0.and.(usepaw==1.or.icoulomb==1)) then 339 call mkcore_alt(atindx1,corstr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,& 340 & nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 341 & ucvol,vxc,xcccrc,xccc1d,xccc3d,xred,pawrad,pawtab,usepaw) 342 end if 343 end if 344 if (usekden==1.and.coretau_method==2) then 345 call mkcore_alt(atindx1,taustr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,& 346 & nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 347 & ucvol,vxctau(:,:,1),xcccrc,xccc1d,xcctau3d,xred,pawrad,pawtab,usepaw,& 348 & usekden=.true.) 349 350 end if 351 ABI_FREE(dyfr_dum) 352 ABI_FREE(gr_dum) 353 ABI_FREE(v_dum) 354 call timab(55,2,tsec) 355 else 356 corstr(:)=zero 357 end if 358 end if 359 360 !======================================================================= 361 !======================= Hartree energy contribution =================== 362 !======================================================================= 363 364 call strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd) 365 366 !======================================================================= 367 !======================= Ewald contribution ============================ 368 !======================================================================= 369 370 call timab(38,1,tsec) 371 call ewald2(gmet,natom,ntypat,rmet,rprimd,ewestr,typat,ucvol,xred,zion) 372 373 !======================================================================= 374 !================== VdW DFT-D contribution ============================ 375 !======================================================================= 376 377 if (vdw_xc==5) then 378 call vdw_dftd2(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_tol,& 379 & xred,znucl,str_vdw_dftd2=vdwstr) 380 elseif (vdw_xc==6.or.vdw_xc==7) then 381 call vdw_dftd3(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_xc,& 382 & vdw_tol,vdw_tol_3bt,xred,znucl,str_vdw_dftd3=vdwstr) 383 end if 384 385 call timab(38,2,tsec) 386 387 !HONG no Berry phase contribution if using reduced ebar or d according to 388 !HONG PRL 89, 117602 (2002) [[cite:Souza2002]] 389 !HONG Nature Physics: M. Stengel et.al. (2009)) [[cite:Stengel1999]] 390 !======================================================================= 391 !=================== Berry phase contribution ========================== 392 !======================================================================= 393 394 !if (berryopt==4) then 395 !berrystr_tmp(:,:) = zero 396 !Diagonal: 397 !do mu = 1, 3 398 !do ii = 1, 3 399 !berrystr_tmp(mu,mu) = berrystr_tmp(mu,mu) - & 400 !& efield(mu)*rprimd(mu,ii)*(pel(ii) + pion(ii))/ucvol 401 !end do 402 !end do 403 !Off-diagonal (symmetrized before adding it to strten): 404 !do ii = 1, 3 405 !berrystr_tmp(3,2) = berrystr_tmp(3,2) & 406 !& - efield(3)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol 407 !berrystr_tmp(2,3) = berrystr_tmp(2,3) & 408 !& - efield(2)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol 409 !berrystr_tmp(3,1) = berrystr_tmp(3,1) & 410 !& - efield(3)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol 411 !berrystr_tmp(1,3) = berrystr_tmp(1,3) & 412 !& - efield(1)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol 413 !berrystr_tmp(2,1) = berrystr_tmp(2,1) & 414 !& - efield(2)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol 415 !berrystr_tmp(1,2) = berrystr_tmp(1,2) & 416 !& - efield(1)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol 417 !end do 418 !berrystr(1) = berrystr_tmp(1,1) 419 !berrystr(2) = berrystr_tmp(2,2) 420 !berrystr(3) = berrystr_tmp(3,3) 421 !berrystr(4) = (berrystr_tmp(3,2) + berrystr_tmp(2,3))/two 422 !berrystr(5) = (berrystr_tmp(3,1) + berrystr_tmp(1,3))/two 423 !berrystr(6) = (berrystr_tmp(2,1) + berrystr_tmp(1,2))/two 424 !end if 425 426 !======================================================================= 427 !================= Other (trivial) contributions ======================= 428 !======================================================================= 429 430 !Nonlocal part of stress has already been computed 431 !(in forstrnps(norm-conserving) or pawgrnl(PAW)) 432 433 !Kinetic part of stress has already been computed 434 !(in forstrnps) 435 436 !cDFT part of stress tensor has already been computed in "constrained_residual" 437 438 !XC part of stress tensor has already been computed in "strsxc" 439 440 !ii part of stress (diagonal) is trivial! 441 strsii=-eii/ucvol 442 !qvpotzero is non zero, only when usepotzero=1 443 strsii=strsii+qvpotzero/ucvol 444 445 !====================================================================== 446 !HONG Maxwell stress when electric/displacement field is non-zero===== 447 !====================================================================== 448 efield_flag = (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 449 & berryopt==14 .or. berryopt==16 .or. berryopt==17) 450 calc_epaw3_stress = (efield_flag .and. usepaw == 1) 451 if ( efield_flag ) then 452 ModE=dot_product(efield,efield) 453 do ii=1,3 454 Maxstr(ii)=two*efield(ii)*efield(ii)-ModE 455 end do 456 Maxstr(4)=two*efield(3)*efield(2) 457 Maxstr(5)=two*efield(3)*efield(1) 458 Maxstr(6)=two*efield(2)*efield(1) 459 ! Converting to units of Ha/Bohr^3 460 ! Maxstr(:)=Maxstr(:)*e_Cb*Bohr_Ang*1.0d-10/(Ha_J*8.0d0*pi) 461 462 Maxstr(:)=Maxstr(:)*eps0*Ha_J*Bohr_Ang*1.0d-10/(8.0d0*pi*e_Cb**2) 463 464 write(message, '(a,a)' )ch10,& 465 & ' Cartesian components of Maxwell stress tensor (hartree/bohr^3)' 466 call wrtout(ab_out,message,'COLL') 467 call wrtout(std_out, message,'COLL') 468 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 469 & ' Maxstr(1 1)=',Maxstr(1),' Maxstr(3 2)=',Maxstr(4) 470 call wrtout(ab_out,message,'COLL') 471 call wrtout(std_out, message,'COLL') 472 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 473 & ' Maxstr(2 2)=',Maxstr(2),' Maxstr(3 1)=',Maxstr(5) 474 call wrtout(ab_out,message,'COLL') 475 call wrtout(std_out, message,'COLL') 476 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 477 & ' Maxstr(3 3)=',Maxstr(3),' Maxstr(2 1)=',Maxstr(6) 478 call wrtout(ab_out,message,'COLL') 479 call wrtout(std_out, message,'COLL') 480 write(message, '(a)' ) ' ' 481 call wrtout(ab_out,message,'COLL') 482 call wrtout(std_out, message,'COLL') 483 484 end if 485 486 ! compute additional F3-type stress due to projectors for electric field with PAW 487 if ( efield_flag .and. calc_epaw3_stress ) then 488 do sdir = 1, 6 489 ep3(:) = zero 490 do idir = 1, 3 491 vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir)) 492 do iatom = 1, natom 493 ep3(idir) = ep3(idir) + vol_element*dtefield%epaws3(iatom,idir,sdir) 494 end do ! end loop over atoms 495 end do ! end loop over idir (components of P) 496 ! note no appearance of ucvol here unlike in forces, stress definition includes 497 ! division by ucvol, which cancels the factor in -ucvol e . p 498 epaws3red(sdir) = -dot_product(red_efieldbar(1:3),ep3(1:3)) 499 end do 500 501 ! write(message, '(a,a)' )ch10,& 502 !& ' Cartesian components of PAW sigma_3 stress tensor (hartree/bohr^3)' 503 ! call wrtout(ab_out,message,'COLL') 504 ! call wrtout(std_out, message,'COLL') 505 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 506 !& ' epaws3red(1 1)=',epaws3red(1),' epaws3red(3 2)=',epaws3red(4) 507 ! call wrtout(ab_out,message,'COLL') 508 ! call wrtout(std_out, message,'COLL') 509 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 510 !& ' epaws3red(2 2)=',epaws3red(2),' epaws3red(3 1)=',epaws3red(5) 511 ! call wrtout(ab_out,message,'COLL') 512 ! call wrtout(std_out, message,'COLL') 513 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 514 !& ' epaws3red(3 3)=',epaws3red(3),' epaws3red(2 1)=',epaws3red(6) 515 ! call wrtout(ab_out,message,'COLL') 516 ! call wrtout(std_out, message,'COLL') 517 ! write(message, '(a)' ) ' ' 518 ! call wrtout(ab_out,message,'COLL') 519 ! call wrtout(std_out, message,'COLL') 520 521 end if 522 523 !======================================================================= 524 !===== Assemble the various contributions to the stress tensor ========= 525 !======================================================================= 526 !In cartesian coordinates (symmetric storage) 527 528 strten(:)=kinstr(:)+ewestr(:)+corstr(:)+strscondft(:)+strsxc(:)+harstr(:)+lpsstr(:)+nlstr(:) 529 530 if (usefock==1 .and. associated(fock)) then 531 if (fock%fock_common%optstr) then 532 strten(:)=strten(:)+fock%fock_common%stress(:) 533 end if 534 end if 535 536 !Add contributions for constant E or D calculation. 537 if ( efield_flag ) then 538 strten(:)=strten(:)+Maxstr(:) 539 if ( calc_epaw3_stress ) strten(:) = strten(:) + epaws3red(:) 540 end if 541 if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)+vdwstr(:) 542 543 !Additional stuff for electron-positron 544 ipositron=0 545 if (present(electronpositron)) then 546 if (associated(electronpositron)) then 547 if (allocated(electronpositron%stress_ep)) ipositron=electronpositron_calctype(electronpositron) 548 end if 549 end if 550 if (abs(ipositron)==1) then 551 strten(:)=strten(:)-harstr(:)-ewestr(:)-corstr(:)-lpsstr(:) 552 harstr(:)=zero;ewestr(:)=zero;corstr(:)=zero;strsii=zero 553 lpsstr(:)=-lpsstr(:);lpsstr(1:3)=lpsstr(1:3)-two*eei/ucvol 554 strten(:)=strten(:)+lpsstr(:) 555 if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)-vdwstr(:) 556 if (vdw_xc>=5.and.vdw_xc<=7) vdwstr(:)=zero 557 end if 558 if (abs(ipositron)==2) then 559 ABI_MALLOC(rhog_ep,(2,nfft)) 560 ABI_MALLOC(dummy,(6)) 561 call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,1,ngfft,0) 562 rhog_ep=-rhog_ep 563 call strhar(electronpositron%e_hartree,gsqcut,dummy,mpi_enreg,nfft,ngfft,rhog_ep,rprimd) 564 strten(:)=strten(:)+dummy(:);harstr(:)=harstr(:)+dummy(:) 565 ABI_FREE(rhog_ep) 566 ABI_FREE(dummy) 567 end if 568 if (ipositron>0) strten(:)=strten(:)+electronpositron%stress_ep(:) 569 570 !Symmetrize resulting tensor if nsym>1 571 if (nsym>1) then 572 call stresssym(gprimd,nsym,strten,symrec) 573 end if 574 575 !Set to zero very small values of stress 576 do mu=1,6 577 if (abs(strten(mu))<tol) strten(mu)=zero 578 end do 579 580 !Include diagonal terms, save uncorrected stress for output 581 do mu=1,3 582 uncorr(mu)=strten(mu)+strsii 583 strten(mu)=uncorr(mu) 584 end do 585 586 !Adding the extfpmd continous contribution to stress tensor 587 if(associated(extfpmd)) then 588 strten(1:3)=strten(1:3)-(2./3.)*extfpmd%e_kinetic/extfpmd%ucvol 589 end if 590 591 !======================================================================= 592 !================ Print out info about stress tensor =================== 593 !======================================================================= 594 if (prtvol>=10.and.ipositron>=0) then 595 write(message, '(a)' ) ' ' 596 call wrtout(std_out,message,'COLL') 597 do mu=1,6 598 write(message, '(a,i5,a,1p,e22.12)' )& 599 & ' stress: component',mu,' of hartree stress is',harstr(mu) 600 call wrtout(std_out,message,'COLL') 601 end do 602 write(message, '(a)' ) ' ' 603 call wrtout(std_out,message,'COLL') 604 do mu=1,6 605 write(message, '(a,i5,a,1p,e22.12)' )& 606 & ' stress: component',mu,' of loc psp stress is',lpsstr(mu) 607 call wrtout(std_out,message,'COLL') 608 end do 609 write(message, '(a)' ) ' ' 610 call wrtout(std_out,message,'COLL') 611 do mu=1,6 612 write(message, '(a,i5,a,1p,e22.12)' )& 613 & ' stress: component',mu,& 614 & ' of kinetic stress is',kinstr(mu) 615 call wrtout(std_out,message,'COLL') 616 end do 617 write(message, '(a)' ) ' ' 618 call wrtout(std_out,message,'COLL') 619 do mu=1,6 620 write(message, '(a,i5,a,1p,e22.12)' )& 621 & ' stress: component',mu,' of nonlocal ps stress is',nlstr(mu) 622 call wrtout(std_out,message,'COLL') 623 end do 624 write(message, '(a)' ) ' ' 625 call wrtout(std_out,message,'COLL') 626 do mu=1,6 627 write(message, '(a,i5,a,1p,e22.12)' )& 628 & ' stress: component',mu,' of core xc stress is',corstr(mu) 629 call wrtout(std_out,message,'COLL') 630 end do 631 write(message, '(a)' ) ' ' 632 call wrtout(std_out,message,'COLL') 633 do mu=1,6 634 write(message, '(a,i5,a,1p,e22.12)' )& 635 & ' stress: component',mu,& 636 & ' of Ewald energ stress is',ewestr(mu) 637 call wrtout(std_out,message,'COLL') 638 end do 639 write(message, '(a)' ) ' ' 640 call wrtout(std_out,message,'COLL') 641 do mu=1,6 642 write(message, '(a,i5,a,1p,e22.12)' ) & 643 & ' stress: component',mu,' of xc stress is',strsxc(mu) 644 call wrtout(std_out,message,'COLL') 645 end do 646 647 if( any( abs(strscondft(:))>tol8 ) )then 648 write(message, '(a)' ) ' ' 649 call wrtout(std_out,message,'COLL') 650 do mu=1,6 651 write(message, '(a,i5,a,1p,e22.12)' ) & 652 & ' stress: component',mu,' of cDFT stress is',strscondft(mu) 653 call wrtout(std_out,message,'COLL') 654 end do 655 endif 656 657 if (vdw_xc>=5.and.vdw_xc<=7) then 658 write(message, '(a)' ) ' ' 659 call wrtout(std_out,message,'COLL') 660 do mu=1,6 661 write(message, '(a,i5,a,1p,e22.12)' )& 662 & ' stress: component',mu,& 663 & ' of VdW DFT-D stress is',vdwstr(mu) 664 call wrtout(std_out,message,'COLL') 665 end do 666 end if 667 write(message, '(a)' ) ' ' 668 call wrtout(std_out,message,'COLL') 669 write(message, '(a,1p,e22.12)' ) & 670 & ' stress: ii (diagonal) part is',strsii 671 call wrtout(std_out,message,'COLL') 672 if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 673 & berryopt==14 .or. berryopt==16 .or. berryopt==17) then !!HONG 674 write(message, '(a)' ) ' ' 675 call wrtout(std_out,message,'COLL') 676 do mu = 1, 6 677 write(message, '(a,i2,a,1p,e22.12)' )& 678 & ' stress: component',mu,' of Maxwell stress is',& 679 & Maxstr(mu) 680 call wrtout(std_out,message,'COLL') 681 end do 682 end if 683 if (ipositron/=0) then 684 write(message, '(a)' ) ' ' 685 call wrtout(std_out,message,'COLL') 686 do mu=1,6 687 write(message, '(a,i5,3a,1p,e22.12)' ) & 688 & ' stress: component',mu,' of ',EPName(abs(ipositron)), & 689 & ' stress is',electronpositron%stress_ep(mu) 690 call wrtout(std_out,message,'COLL') 691 end do 692 end if 693 694 end if ! prtvol 695 if (ipositron>=0) then 696 write(message, '(a,a)' )ch10,& 697 & ' Cartesian components of stress tensor (hartree/bohr^3)' 698 call wrtout(ab_out,message,'COLL') 699 call wrtout(std_out, message,'COLL') 700 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 701 & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4) 702 call wrtout(ab_out,message,'COLL') 703 call wrtout(std_out, message,'COLL') 704 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 705 & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5) 706 call wrtout(ab_out,message,'COLL') 707 call wrtout(std_out, message,'COLL') 708 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 709 & ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6) 710 call wrtout(ab_out,message,'COLL') 711 call wrtout(std_out, message,'COLL') 712 write(message, '(a)' ) ' ' 713 call wrtout(ab_out,message,'COLL') 714 call wrtout(std_out, message,'COLL') 715 end if 716 call timab(37,2,tsec) 717 718 end subroutine stress
ABINIT/strhar [ Functions ]
NAME
strhar
FUNCTION
Compute Hartree energy contribution to stress tensor (Cartesian coordinates).
INPUTS
ehart=Hartree energy (hartree) gsqcut=cutoff value on $G^2$ for (large) sphere inside fft box. $gsqcut=(boxcut^2)*ecut/(2._dp*(\pi^2))$ mpi_enreg=information about MPI parallelization 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 rhog(2,nfft)=Fourier transform of charge density (bohr^-3) rhog(2,nfft)= optional argument: Fourier transform of a second charge density (bohr^-3) rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
harstr(6)=components of Hartree part of stress tensor (Cartesian coordinates, symmetric tensor) in hartree/bohr^3 Definition of symmetric tensor storage: store 6 unique components in the order 11, 22, 33, 32, 31, 21 (suggested by Xavier Gonze).
SOURCE
748 subroutine strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd,& 749 & rhog2) ! optional argument 750 751 !Arguments ------------------------------------ 752 !scalars 753 integer,intent(in) :: nfft 754 real(dp),intent(in) :: ehart,gsqcut 755 type(MPI_type),intent(in) :: mpi_enreg 756 !arrays 757 integer,intent(in) :: ngfft(18) 758 real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft) 759 real(dp),intent(in),optional :: rhog2(2,nfft) 760 real(dp),intent(out) :: harstr(6) 761 762 !Local variables------------------------------- 763 !scalars 764 integer,parameter :: im=2,re=1 765 integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,me_fft,n1,n2,n3,nproc_fft 766 real(dp) :: cutoff,gsquar,rhogsq,tolfix=1.000000001_dp,ucvol 767 !arrays 768 real(dp) :: gcart(3),gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2) 769 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 770 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 771 772 ! ************************************************************************* 773 774 call timab(568,1,tsec) 775 776 harstr(:)=zero 777 !ehtest=0.0_dp (used for testing) 778 779 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 780 781 irho2=0;if (present(rhog2)) irho2=1 782 783 !Conduct looping over all fft grid points to find G vecs inside gsqcut 784 !Include G**2 on surface of cutoff sphere as well as inside: 785 cutoff=gsqcut*tolfix 786 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 787 me_fft=ngfft(11) 788 nproc_fft=ngfft(10) 789 id1=n1/2+2 790 id2=n2/2+2 791 id3=n3/2+2 792 ii=0 793 794 ! Get the distrib associated with this fft_grid 795 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 796 797 do i3=1,n3 798 ig3=i3-(i3/id3)*n3-1 799 do i2=1,n2 800 ig2=i2-(i2/id2)*n2-1 801 if (fftn2_distrib(i2)==me_fft) then 802 do i1=1,n1 803 ig1=i1-(i1/id1)*n1-1 804 ! ii=ii+1 805 ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1)) 806 ! ** GET RID OF THIS IF STATEMENT LATER for speed if needed 807 ! Avoid G=0: 808 ! if (ii>1) then 809 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 810 ! Compute cartesian components of G 811 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+gprimd(1,3)*dble(ig3) 812 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+gprimd(2,3)*dble(ig3) 813 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+gprimd(3,3)*dble(ig3) 814 ! Compute |G|^2 815 gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2 816 817 ! Keep only G**2 inside larger cutoff (not sure this is needed): 818 if (gsquar<=cutoff) then 819 ! take |rho(G)|^2 for complex rhog 820 if (irho2==0) then 821 rhogsq=rhog(re,ii)**2+rhog(im,ii)**2 822 else 823 rhogsq=rhog(re,ii)*rhog2(re,ii)+rhog(im,ii)*rhog2(im,ii) 824 end if 825 harstr(1)=harstr(1)+(rhogsq/gsquar**2)*gcart(1)*gcart(1) 826 harstr(2)=harstr(2)+(rhogsq/gsquar**2)*gcart(2)*gcart(2) 827 harstr(3)=harstr(3)+(rhogsq/gsquar**2)*gcart(3)*gcart(3) 828 harstr(4)=harstr(4)+(rhogsq/gsquar**2)*gcart(3)*gcart(2) 829 harstr(5)=harstr(5)+(rhogsq/gsquar**2)*gcart(3)*gcart(1) 830 harstr(6)=harstr(6)+(rhogsq/gsquar**2)*gcart(2)*gcart(1) 831 end if 832 ! end if 833 end do 834 end if 835 end do 836 end do 837 838 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel 839 #ifdef FC_IBM 840 write(std_out,*)' strhar : before mpi_comm, harstr=',harstr 841 #endif 842 843 !Init mpi_comm 844 if(mpi_enreg%nproc_fft>1)then 845 call timab(48,1,tsec) 846 call xmpi_sum(harstr,mpi_enreg%comm_fft ,ierr) 847 call timab(48,2,tsec) 848 end if 849 850 #ifdef FC_IBM 851 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel 852 write(std_out,*)' strhar : after mpi_comm, harstr=',harstr 853 write(std_out,*)' strhar : ehart,ucvol=',ehart,ucvol 854 #endif 855 856 !Normalize and add term -ehart/ucvol on diagonal 857 harstr(1)=harstr(1)/pi-ehart/ucvol 858 harstr(2)=harstr(2)/pi-ehart/ucvol 859 harstr(3)=harstr(3)/pi-ehart/ucvol 860 harstr(4)=harstr(4)/pi 861 harstr(5)=harstr(5)/pi 862 harstr(6)=harstr(6)/pi 863 864 call timab(568,2,tsec) 865 866 end subroutine strhar