TABLE OF CONTENTS
ABINIT/m_stress [ Modules ]
NAME
m_stress
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 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 ABI_FREE(vxctotg) 302 end if 303 call timab(551,2,tsec) 304 end if 305 306 !Local ionic potential by method 2 307 if (vloc_method==2) then 308 option=3 309 ABI_MALLOC(dyfr_dum,(3,3,natom)) 310 ABI_MALLOC(gr_dum,(3,natom)) 311 ABI_MALLOC(v_dum,(nfft)) 312 call mklocl_recipspace(dyfr_dum,eei,gmet,gprimd,gr_dum,gsqcut,icutcoul,lpsstr,mgfft,& 313 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,nkpt,ntypat,option,ph1d,qgrid,& 314 & qprtrb_dum,dum_rcut,rhog,rprimd,ucvol,vcutgeo,vlspl,vprtrb_dum,v_dum) 315 ABI_FREE(dyfr_dum) 316 ABI_FREE(gr_dum) 317 ABI_FREE(v_dum) 318 end if 319 320 !Pseudo core electron density by method 2 321 if (coredens_method==2.or.coretau_method==2) then 322 if (n1xccc/=0) then 323 call timab(55,1,tsec) 324 option=3 325 ABI_MALLOC(dyfr_dum,(3,3,natom)) 326 ABI_MALLOC(gr_dum,(3,natom)) 327 ABI_MALLOC(v_dum,(nfft)) 328 if (coredens_method==2) then 329 if (psps%usewvl==0.and.usepaw==0.and.icoulomb==0) then 330 if(opt_hybr==0) then 331 call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),& 332 & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc,& 333 & xcccrc,xccc1d,xccc3d,xred) 334 else 335 call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),& 336 & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc_hf,& 337 & xcccrc,xccc1d,xccc3d,xred) 338 end if 339 else if (psps%usewvl==0.and.(usepaw==1.or.icoulomb==1)) then 340 call mkcore_alt(atindx1,corstr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,& 341 & nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 342 & ucvol,vxc,xcccrc,xccc1d,xccc3d,xred,pawrad,pawtab,usepaw) 343 end if 344 end if 345 if (usekden==1.and.coretau_method==2) then 346 call mkcore_alt(atindx1,taustr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,& 347 & nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 348 & ucvol,vxctau(:,:,1),xcccrc,xccc1d,xcctau3d,xred,pawrad,pawtab,usepaw,& 349 & usekden=.true.) 350 351 end if 352 ABI_FREE(dyfr_dum) 353 ABI_FREE(gr_dum) 354 ABI_FREE(v_dum) 355 call timab(55,2,tsec) 356 else 357 corstr(:)=zero 358 end if 359 end if 360 361 !======================================================================= 362 !======================= Hartree energy contribution =================== 363 !======================================================================= 364 365 call strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd) 366 367 !======================================================================= 368 !======================= Ewald contribution ============================ 369 !======================================================================= 370 371 call timab(38,1,tsec) 372 call ewald2(gmet,natom,ntypat,rmet,rprimd,ewestr,typat,ucvol,xred,zion) 373 374 !======================================================================= 375 !================== VdW DFT-D contribution ============================ 376 !======================================================================= 377 378 if (vdw_xc==5) then 379 call vdw_dftd2(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_tol,& 380 & xred,znucl,str_vdw_dftd2=vdwstr) 381 elseif (vdw_xc==6.or.vdw_xc==7) then 382 call vdw_dftd3(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_xc,& 383 & vdw_tol,vdw_tol_3bt,xred,znucl,str_vdw_dftd3=vdwstr) 384 end if 385 386 call timab(38,2,tsec) 387 388 !HONG no Berry phase contribution if using reduced ebar or d according to 389 !HONG PRL 89, 117602 (2002) [[cite:Souza2002]] 390 !HONG Nature Physics: M. Stengel et.al. (2009)) [[cite:Stengel1999]] 391 !======================================================================= 392 !=================== Berry phase contribution ========================== 393 !======================================================================= 394 395 !if (berryopt==4) then 396 !berrystr_tmp(:,:) = zero 397 !Diagonal: 398 !do mu = 1, 3 399 !do ii = 1, 3 400 !berrystr_tmp(mu,mu) = berrystr_tmp(mu,mu) - & 401 !& efield(mu)*rprimd(mu,ii)*(pel(ii) + pion(ii))/ucvol 402 !end do 403 !end do 404 !Off-diagonal (symmetrized before adding it to strten): 405 !do ii = 1, 3 406 !berrystr_tmp(3,2) = berrystr_tmp(3,2) & 407 !& - efield(3)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol 408 !berrystr_tmp(2,3) = berrystr_tmp(2,3) & 409 !& - efield(2)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol 410 !berrystr_tmp(3,1) = berrystr_tmp(3,1) & 411 !& - efield(3)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol 412 !berrystr_tmp(1,3) = berrystr_tmp(1,3) & 413 !& - efield(1)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol 414 !berrystr_tmp(2,1) = berrystr_tmp(2,1) & 415 !& - efield(2)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol 416 !berrystr_tmp(1,2) = berrystr_tmp(1,2) & 417 !& - efield(1)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol 418 !end do 419 !berrystr(1) = berrystr_tmp(1,1) 420 !berrystr(2) = berrystr_tmp(2,2) 421 !berrystr(3) = berrystr_tmp(3,3) 422 !berrystr(4) = (berrystr_tmp(3,2) + berrystr_tmp(2,3))/two 423 !berrystr(5) = (berrystr_tmp(3,1) + berrystr_tmp(1,3))/two 424 !berrystr(6) = (berrystr_tmp(2,1) + berrystr_tmp(1,2))/two 425 !end if 426 427 !======================================================================= 428 !================= Other (trivial) contributions ======================= 429 !======================================================================= 430 431 !Nonlocal part of stress has already been computed 432 !(in forstrnps(norm-conserving) or pawgrnl(PAW)) 433 434 !Kinetic part of stress has already been computed 435 !(in forstrnps) 436 437 !cDFT part of stress tensor has already been computed in "constrained_residual" 438 439 !XC part of stress tensor has already been computed in "strsxc" 440 441 !ii part of stress (diagonal) is trivial! 442 strsii=-eii/ucvol 443 !qvpotzero is non zero, only when usepotzero=1 444 strsii=strsii+qvpotzero/ucvol 445 446 !====================================================================== 447 !HONG Maxwell stress when electric/displacement field is non-zero===== 448 !====================================================================== 449 efield_flag = (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 450 & berryopt==14 .or. berryopt==16 .or. berryopt==17) 451 calc_epaw3_stress = (efield_flag .and. usepaw == 1) 452 if ( efield_flag ) then 453 ModE=dot_product(efield,efield) 454 do ii=1,3 455 Maxstr(ii)=two*efield(ii)*efield(ii)-ModE 456 end do 457 Maxstr(4)=two*efield(3)*efield(2) 458 Maxstr(5)=two*efield(3)*efield(1) 459 Maxstr(6)=two*efield(2)*efield(1) 460 ! Converting to units of Ha/Bohr^3 461 ! Maxstr(:)=Maxstr(:)*e_Cb*Bohr_Ang*1.0d-10/(Ha_J*8.0d0*pi) 462 463 Maxstr(:)=Maxstr(:)*eps0*Ha_J*Bohr_Ang*1.0d-10/(8.0d0*pi*e_Cb**2) 464 465 write(message, '(a,a)' )ch10,& 466 & ' Cartesian components of Maxwell stress tensor (hartree/bohr^3)' 467 call wrtout(ab_out,message,'COLL') 468 call wrtout(std_out, message,'COLL') 469 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 470 & ' Maxstr(1 1)=',Maxstr(1),' Maxstr(3 2)=',Maxstr(4) 471 call wrtout(ab_out,message,'COLL') 472 call wrtout(std_out, message,'COLL') 473 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 474 & ' Maxstr(2 2)=',Maxstr(2),' Maxstr(3 1)=',Maxstr(5) 475 call wrtout(ab_out,message,'COLL') 476 call wrtout(std_out, message,'COLL') 477 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 478 & ' Maxstr(3 3)=',Maxstr(3),' Maxstr(2 1)=',Maxstr(6) 479 call wrtout(ab_out,message,'COLL') 480 call wrtout(std_out, message,'COLL') 481 write(message, '(a)' ) ' ' 482 call wrtout(ab_out,message,'COLL') 483 call wrtout(std_out, message,'COLL') 484 485 end if 486 487 ! compute additional F3-type stress due to projectors for electric field with PAW 488 if ( efield_flag .and. calc_epaw3_stress ) then 489 do sdir = 1, 6 490 ep3(:) = zero 491 do idir = 1, 3 492 vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir)) 493 do iatom = 1, natom 494 ep3(idir) = ep3(idir) + vol_element*dtefield%epaws3(iatom,idir,sdir) 495 end do ! end loop over atoms 496 end do ! end loop over idir (components of P) 497 ! note no appearance of ucvol here unlike in forces, stress definition includes 498 ! division by ucvol, which cancels the factor in -ucvol e . p 499 epaws3red(sdir) = -dot_product(red_efieldbar(1:3),ep3(1:3)) 500 end do 501 502 ! write(message, '(a,a)' )ch10,& 503 !& ' Cartesian components of PAW sigma_3 stress tensor (hartree/bohr^3)' 504 ! call wrtout(ab_out,message,'COLL') 505 ! call wrtout(std_out, message,'COLL') 506 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 507 !& ' epaws3red(1 1)=',epaws3red(1),' epaws3red(3 2)=',epaws3red(4) 508 ! call wrtout(ab_out,message,'COLL') 509 ! call wrtout(std_out, message,'COLL') 510 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 511 !& ' epaws3red(2 2)=',epaws3red(2),' epaws3red(3 1)=',epaws3red(5) 512 ! call wrtout(ab_out,message,'COLL') 513 ! call wrtout(std_out, message,'COLL') 514 ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 515 !& ' epaws3red(3 3)=',epaws3red(3),' epaws3red(2 1)=',epaws3red(6) 516 ! call wrtout(ab_out,message,'COLL') 517 ! call wrtout(std_out, message,'COLL') 518 ! write(message, '(a)' ) ' ' 519 ! call wrtout(ab_out,message,'COLL') 520 ! call wrtout(std_out, message,'COLL') 521 522 end if 523 524 !======================================================================= 525 !===== Assemble the various contributions to the stress tensor ========= 526 !======================================================================= 527 !In cartesian coordinates (symmetric storage) 528 529 strten(:)=kinstr(:)+ewestr(:)+corstr(:)+strscondft(:)+strsxc(:)+harstr(:)+lpsstr(:)+nlstr(:) 530 531 if (usefock==1 .and. associated(fock)) then 532 if (fock%fock_common%optstr) then 533 strten(:)=strten(:)+fock%fock_common%stress(:) 534 end if 535 end if 536 537 !Add contributions for constant E or D calculation. 538 if ( efield_flag ) then 539 strten(:)=strten(:)+Maxstr(:) 540 if ( calc_epaw3_stress ) strten(:) = strten(:) + epaws3red(:) 541 end if 542 if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)+vdwstr(:) 543 544 !Additional stuff for electron-positron 545 ipositron=0 546 if (present(electronpositron)) then 547 if (associated(electronpositron)) then 548 if (allocated(electronpositron%stress_ep)) ipositron=electronpositron_calctype(electronpositron) 549 end if 550 end if 551 if (abs(ipositron)==1) then 552 strten(:)=strten(:)-harstr(:)-ewestr(:)-corstr(:)-lpsstr(:) 553 harstr(:)=zero;ewestr(:)=zero;corstr(:)=zero;strsii=zero 554 lpsstr(:)=-lpsstr(:);lpsstr(1:3)=lpsstr(1:3)-two*eei/ucvol 555 strten(:)=strten(:)+lpsstr(:) 556 if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)-vdwstr(:) 557 if (vdw_xc>=5.and.vdw_xc<=7) vdwstr(:)=zero 558 end if 559 if (abs(ipositron)==2) then 560 ABI_MALLOC(rhog_ep,(2,nfft)) 561 ABI_MALLOC(dummy,(6)) 562 call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,1,ngfft,0) 563 rhog_ep=-rhog_ep 564 call strhar(electronpositron%e_hartree,gsqcut,dummy,mpi_enreg,nfft,ngfft,rhog_ep,rprimd) 565 strten(:)=strten(:)+dummy(:);harstr(:)=harstr(:)+dummy(:) 566 ABI_FREE(rhog_ep) 567 ABI_FREE(dummy) 568 end if 569 if (ipositron>0) strten(:)=strten(:)+electronpositron%stress_ep(:) 570 571 !Symmetrize resulting tensor if nsym>1 572 if (nsym>1) then 573 call stresssym(gprimd,nsym,strten,symrec) 574 end if 575 576 !Set to zero very small values of stress 577 do mu=1,6 578 if (abs(strten(mu))<tol) strten(mu)=zero 579 end do 580 581 !Include diagonal terms, save uncorrected stress for output 582 do mu=1,3 583 uncorr(mu)=strten(mu)+strsii 584 strten(mu)=uncorr(mu) 585 end do 586 587 !Adding the extfpmd continous contribution to stress tensor 588 if(associated(extfpmd)) then 589 strten(1:3)=strten(1:3)-(2./3.)*extfpmd%e_kinetic/extfpmd%ucvol 590 end if 591 592 !======================================================================= 593 !================ Print out info about stress tensor =================== 594 !======================================================================= 595 if (prtvol>=10.and.ipositron>=0) then 596 write(message, '(a)' ) ' ' 597 call wrtout(std_out,message,'COLL') 598 do mu=1,6 599 write(message, '(a,i5,a,1p,e22.12)' )& 600 & ' stress: component',mu,' of hartree stress is',harstr(mu) 601 call wrtout(std_out,message,'COLL') 602 end do 603 write(message, '(a)' ) ' ' 604 call wrtout(std_out,message,'COLL') 605 do mu=1,6 606 write(message, '(a,i5,a,1p,e22.12)' )& 607 & ' stress: component',mu,' of loc psp stress is',lpsstr(mu) 608 call wrtout(std_out,message,'COLL') 609 end do 610 write(message, '(a)' ) ' ' 611 call wrtout(std_out,message,'COLL') 612 do mu=1,6 613 write(message, '(a,i5,a,1p,e22.12)' )& 614 & ' stress: component',mu,& 615 & ' of kinetic stress is',kinstr(mu) 616 call wrtout(std_out,message,'COLL') 617 end do 618 write(message, '(a)' ) ' ' 619 call wrtout(std_out,message,'COLL') 620 do mu=1,6 621 write(message, '(a,i5,a,1p,e22.12)' )& 622 & ' stress: component',mu,' of nonlocal ps stress is',nlstr(mu) 623 call wrtout(std_out,message,'COLL') 624 end do 625 write(message, '(a)' ) ' ' 626 call wrtout(std_out,message,'COLL') 627 do mu=1,6 628 write(message, '(a,i5,a,1p,e22.12)' )& 629 & ' stress: component',mu,' of core xc stress is',corstr(mu) 630 call wrtout(std_out,message,'COLL') 631 end do 632 write(message, '(a)' ) ' ' 633 call wrtout(std_out,message,'COLL') 634 do mu=1,6 635 write(message, '(a,i5,a,1p,e22.12)' )& 636 & ' stress: component',mu,& 637 & ' of Ewald energ stress is',ewestr(mu) 638 call wrtout(std_out,message,'COLL') 639 end do 640 write(message, '(a)' ) ' ' 641 call wrtout(std_out,message,'COLL') 642 do mu=1,6 643 write(message, '(a,i5,a,1p,e22.12)' ) & 644 & ' stress: component',mu,' of xc stress is',strsxc(mu) 645 call wrtout(std_out,message,'COLL') 646 end do 647 648 if( any( abs(strscondft(:))>tol8 ) )then 649 write(message, '(a)' ) ' ' 650 call wrtout(std_out,message,'COLL') 651 do mu=1,6 652 write(message, '(a,i5,a,1p,e22.12)' ) & 653 & ' stress: component',mu,' of cDFT stress is',strscondft(mu) 654 call wrtout(std_out,message,'COLL') 655 end do 656 endif 657 658 if (vdw_xc>=5.and.vdw_xc<=7) then 659 write(message, '(a)' ) ' ' 660 call wrtout(std_out,message,'COLL') 661 do mu=1,6 662 write(message, '(a,i5,a,1p,e22.12)' )& 663 & ' stress: component',mu,& 664 & ' of VdW DFT-D stress is',vdwstr(mu) 665 call wrtout(std_out,message,'COLL') 666 end do 667 end if 668 write(message, '(a)' ) ' ' 669 call wrtout(std_out,message,'COLL') 670 write(message, '(a,1p,e22.12)' ) & 671 & ' stress: ii (diagonal) part is',strsii 672 call wrtout(std_out,message,'COLL') 673 if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 674 & berryopt==14 .or. berryopt==16 .or. berryopt==17) then !!HONG 675 write(message, '(a)' ) ' ' 676 call wrtout(std_out,message,'COLL') 677 do mu = 1, 6 678 write(message, '(a,i2,a,1p,e22.12)' )& 679 & ' stress: component',mu,' of Maxwell stress is',& 680 & Maxstr(mu) 681 call wrtout(std_out,message,'COLL') 682 end do 683 end if 684 if (ipositron/=0) then 685 write(message, '(a)' ) ' ' 686 call wrtout(std_out,message,'COLL') 687 do mu=1,6 688 write(message, '(a,i5,3a,1p,e22.12)' ) & 689 & ' stress: component',mu,' of ',EPName(abs(ipositron)), & 690 & ' stress is',electronpositron%stress_ep(mu) 691 call wrtout(std_out,message,'COLL') 692 end do 693 end if 694 695 end if ! prtvol 696 if (ipositron>=0) then 697 write(message, '(a,a)' )ch10,& 698 & ' Cartesian components of stress tensor (hartree/bohr^3)' 699 call wrtout(ab_out,message,'COLL') 700 call wrtout(std_out, message,'COLL') 701 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 702 & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4) 703 call wrtout(ab_out,message,'COLL') 704 call wrtout(std_out, message,'COLL') 705 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 706 & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5) 707 call wrtout(ab_out,message,'COLL') 708 call wrtout(std_out, message,'COLL') 709 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 710 & ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6) 711 call wrtout(ab_out,message,'COLL') 712 call wrtout(std_out, message,'COLL') 713 write(message, '(a)' ) ' ' 714 call wrtout(ab_out,message,'COLL') 715 call wrtout(std_out, message,'COLL') 716 end if 717 call timab(37,2,tsec) 718 719 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
749 subroutine strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd,& 750 & rhog2) ! optional argument 751 752 !Arguments ------------------------------------ 753 !scalars 754 integer,intent(in) :: nfft 755 real(dp),intent(in) :: ehart,gsqcut 756 type(MPI_type),intent(in) :: mpi_enreg 757 !arrays 758 integer,intent(in) :: ngfft(18) 759 real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft) 760 real(dp),intent(in),optional :: rhog2(2,nfft) 761 real(dp),intent(out) :: harstr(6) 762 763 !Local variables------------------------------- 764 !scalars 765 integer,parameter :: im=2,re=1 766 integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,me_fft,n1,n2,n3,nproc_fft 767 real(dp) :: cutoff,gsquar,rhogsq,tolfix=1.000000001_dp,ucvol 768 !arrays 769 real(dp) :: gcart(3),gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2) 770 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 771 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 772 773 ! ************************************************************************* 774 775 call timab(568,1,tsec) 776 777 harstr(:)=zero 778 !ehtest=0.0_dp (used for testing) 779 780 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 781 782 irho2=0;if (present(rhog2)) irho2=1 783 784 !Conduct looping over all fft grid points to find G vecs inside gsqcut 785 !Include G**2 on surface of cutoff sphere as well as inside: 786 cutoff=gsqcut*tolfix 787 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 788 me_fft=ngfft(11) 789 nproc_fft=ngfft(10) 790 id1=n1/2+2 791 id2=n2/2+2 792 id3=n3/2+2 793 ii=0 794 795 ! Get the distrib associated with this fft_grid 796 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 797 798 do i3=1,n3 799 ig3=i3-(i3/id3)*n3-1 800 do i2=1,n2 801 ig2=i2-(i2/id2)*n2-1 802 if (fftn2_distrib(i2)==me_fft) then 803 do i1=1,n1 804 ig1=i1-(i1/id1)*n1-1 805 ! ii=ii+1 806 ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1)) 807 ! ** GET RID OF THIS IF STATEMENT LATER for speed if needed 808 ! Avoid G=0: 809 ! if (ii>1) then 810 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 811 ! Compute cartesian components of G 812 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+gprimd(1,3)*dble(ig3) 813 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+gprimd(2,3)*dble(ig3) 814 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+gprimd(3,3)*dble(ig3) 815 ! Compute |G|^2 816 gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2 817 818 ! Keep only G**2 inside larger cutoff (not sure this is needed): 819 if (gsquar<=cutoff) then 820 ! take |rho(G)|^2 for complex rhog 821 if (irho2==0) then 822 rhogsq=rhog(re,ii)**2+rhog(im,ii)**2 823 else 824 rhogsq=rhog(re,ii)*rhog2(re,ii)+rhog(im,ii)*rhog2(im,ii) 825 end if 826 harstr(1)=harstr(1)+(rhogsq/gsquar**2)*gcart(1)*gcart(1) 827 harstr(2)=harstr(2)+(rhogsq/gsquar**2)*gcart(2)*gcart(2) 828 harstr(3)=harstr(3)+(rhogsq/gsquar**2)*gcart(3)*gcart(3) 829 harstr(4)=harstr(4)+(rhogsq/gsquar**2)*gcart(3)*gcart(2) 830 harstr(5)=harstr(5)+(rhogsq/gsquar**2)*gcart(3)*gcart(1) 831 harstr(6)=harstr(6)+(rhogsq/gsquar**2)*gcart(2)*gcart(1) 832 end if 833 ! end if 834 end do 835 end if 836 end do 837 end do 838 839 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel 840 #ifdef FC_IBM 841 write(std_out,*)' strhar : before mpi_comm, harstr=',harstr 842 #endif 843 844 !Init mpi_comm 845 if(mpi_enreg%nproc_fft>1)then 846 call timab(48,1,tsec) 847 call xmpi_sum(harstr,mpi_enreg%comm_fft ,ierr) 848 call timab(48,2,tsec) 849 end if 850 851 #ifdef FC_IBM 852 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel 853 write(std_out,*)' strhar : after mpi_comm, harstr=',harstr 854 write(std_out,*)' strhar : ehart,ucvol=',ehart,ucvol 855 #endif 856 857 !Normalize and add term -ehart/ucvol on diagonal 858 harstr(1)=harstr(1)/pi-ehart/ucvol 859 harstr(2)=harstr(2)/pi-ehart/ucvol 860 harstr(3)=harstr(3)/pi-ehart/ucvol 861 harstr(4)=harstr(4)/pi 862 harstr(5)=harstr(5)/pi 863 harstr(6)=harstr(6)/pi 864 865 call timab(568,2,tsec) 866 867 end subroutine strhar