TABLE OF CONTENTS


ABINIT/m_stress [ Modules ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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