TABLE OF CONTENTS


ABINIT/m_stress [ Modules ]

[ Top ] [ Modules ]

NAME

  m_stress

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, FJ, MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_stress
28 
29  use defs_basis
30  use defs_abitypes
31  use m_efield
32  use m_abicore
33  use m_errors
34  use m_xmpi
35 
36  use m_time,             only : timab
37  use m_geometry,         only : metric, stresssym
38  use m_fock,             only : fock_type
39  use m_ewald,            only : ewald2
40  use defs_datatypes,     only : pseudopotential_type
41  use m_pawrad,           only : pawrad_type
42  use m_pawtab,           only : pawtab_type
43  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
44  use m_fft,              only : zerosym, fourdp
45  use m_mpinfo,           only : ptabs_fourdp
46  use m_vdw_dftd2,        only : vdw_dftd2
47  use m_vdw_dftd3,        only : vdw_dftd3
48  use m_atm2fft,          only : atm2fft
49  use m_mklocl,           only : mklocl_recipspace
50  use m_mkcore,           only : mkcore, mkcore_alt
51 
52  implicit none
53 
54  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=informations about MPI parallelization
  mqgrid=dimensioned number of q grid points for local psp spline
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  n3xccc=dimension of the xccc3d array (0 or nfft).
  natom=number of atoms in cell
  nattyp(ntypat)=number of atoms of each type
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nlstr(6)=nonlocal part of stress tensor
  nspden=number of spin-density components
  nsym=number of symmetries in space group
  ntypat=number of types of atoms
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) array
  prtvol=integer controlling volume of printed output
  qgrid(mqgrid)=q point array for local psp spline fits
  red_efieldbar(3) = efield in reduced units relative to reciprocal lattice
  rhog(2,nfft)=Fourier transform of charge density (bohr^-3)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  strsxc(6)=xc correction to stress
  symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
  typat(natom)=type integer for each atom in cell
  usefock=1 if fock operator is used; 0 otherwise.
  usepaw= 0 for non paw calculation; =1 for paw calculation
  vdw_tol= Van der Waals tolerance
  vdw_tol_3bt= Van der Waals tolerance on the 3-body term (only effective
               vdw_xc=6)
  vdw_xc= Van der Waals correction flag
  vlspl(mqgrid,2,ntypat)=local psp spline
  vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space
  vxc_hf(nfft,nspden)=exchange-correlation potential (hartree) in real space for Hartree-Fock corrections
  xccc1d(n1xccc*(1-usepaw),6,ntypat)=1D core charge function and five derivatives,
                          for each type of atom, from psp (used in Norm-conserving only)
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
  xred(3,natom)=reduced dimensionless atomic coordinates
  zion(ntypat)=valence charge of each type of atom
  znucl(ntypat)=atomic number of atom type

OUTPUT

  strten(6)=components of the stress tensor (hartree/bohr^3) for the
    6 unique components of this symmetric 3x3 tensor:
    Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
    The diagonal components of the returned stress tensor are
    CORRECTED for the Pulay stress.

SIDE EFFECTS

  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)

NOTES

 * Concerning the stress tensor:
   See O. H. Nielsen and R. M. Martin, PRB 32, 3792 (1985) [[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.

PARENTS

      forstr

CHILDREN

      atm2fft,ewald2,fourdp,metric,mkcore,mkcore_alt,mklocl_recipspace
      stresssym,strhar,timab,vdw_dftd2,vdw_dftd3,wrtout,zerosym

SOURCE

174  subroutine stress(atindx1,berryopt,dtefield,eei,efield,ehart,eii,fock,gsqcut,ixc,kinstr,&
175 &                  mgfft,mpi_enreg,mqgrid,n1xccc,n3xccc,natom,nattyp,&
176 &                  nfft,ngfft,nlstr,nspden,nsym,ntypat,paral_kgb,psps,pawrad,pawtab,ph1d,&
177 &                  prtvol,qgrid,red_efieldbar,rhog,rprimd,strten,strsxc,symrec,&
178 &                  typat,usefock,usepaw,vdw_tol,vdw_tol_3bt,vdw_xc,&
179 &                  vlspl,vxc,vxc_hf,xccc1d,xccc3d,xcccrc,xred,zion,znucl,qvpotzero,&
180 &                  electronpositron) ! optional argument
181 
182 
183 !This section has been created automatically by the script Abilint (TD).
184 !Do not modify the following lines by hand.
185 #undef ABI_FUNC
186 #define ABI_FUNC 'stress'
187 !End of the abilint section
188 
189  implicit none
190 
191 !Arguments ------------------------------------
192 !scalars
193  integer,intent(in) :: berryopt,ixc,mgfft,mqgrid,n1xccc,n3xccc,natom,nfft,nspden
194  integer,intent(in) :: nsym,ntypat,paral_kgb,prtvol,usefock,usepaw,vdw_xc
195  real(dp),intent(in) :: eei,ehart,eii,gsqcut,vdw_tol,vdw_tol_3bt,qvpotzero
196  type(efield_type),intent(in) :: dtefield
197  type(pseudopotential_type),intent(in) :: psps
198  type(electronpositron_type),pointer,optional :: electronpositron
199  type(MPI_type),intent(in) :: mpi_enreg
200  type(fock_type),pointer, intent(inout) :: fock
201 !arrays
202  integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),symrec(3,3,nsym)
203  integer,intent(in) :: typat(natom)
204  real(dp),intent(in) :: efield(3),kinstr(6),nlstr(6)
205  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid)
206  real(dp),intent(in) :: red_efieldbar(3),rhog(2,nfft),rprimd(3,3),strsxc(6)
207  real(dp),intent(in) :: vlspl(mqgrid,2,ntypat),vxc(nfft,nspden)
208  real(dp),allocatable,intent(in) :: vxc_hf(:,:)
209  real(dp),intent(in) :: xccc1d(n1xccc*(1-usepaw),6,ntypat),xcccrc(ntypat)
210  real(dp),intent(in) :: xred(3,natom),zion(ntypat),znucl(ntypat)
211  real(dp),intent(inout) :: xccc3d(n3xccc)
212  real(dp),intent(out) :: strten(6)
213  type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw)
214  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
215 
216 !Local variables-------------------------------
217 !scalars
218  integer :: coredens_method,iatom,icoulomb,idir,ii,ipositron,mu,optatm,optdyfr,opteltfr,opt_hybr,optgr,option
219  integer :: optn,optn2,optstr,optv,sdir,vloc_method
220  real(dp),parameter :: tol=1.0d-15
221  real(dp) :: e_dum,strsii,ucvol,vol_element
222  character(len=500) :: message
223  logical :: calc_epaw3_stress, efield_flag
224 !arrays
225  integer :: qprtrb_dum(3)
226  real(dp) :: corstr(6),ep3(3),epaws3red(6),ewestr(6),gmet(3,3)
227 !Maxwell-stress constribution, and magnitude of efield
228  real(dp) :: Maxstr(6),ModE
229  real(dp) :: gprimd(3,3),harstr(6),lpsstr(6),rmet(3,3),tsec(2),uncorr(3)
230  real(dp) :: vdwstr(6),vprtrb_dum(2)
231  real(dp) :: dummy_in(0)
232  real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0)
233  real(dp),allocatable :: dummy(:),dyfr_dum(:,:,:),gr_dum(:,:),rhog_ep(:,:),v_dum(:)
234  real(dp),allocatable :: vxctotg(:,:)
235  character(len=10) :: EPName(1:2)=(/"Electronic","Positronic"/)
236 
237 ! *************************************************************************
238 
239  call timab(37,1,tsec)
240 
241 !Compute different geometric tensor, as well as ucvol, from rprimd
242  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
243  opt_hybr=0
244  if (allocated(vxc_hf)) opt_hybr=1
245 !=======================================================================
246 !========= Local pseudopotential and core charge contributions =========
247 !=======================================================================
248 
249 !Determine by which method the local ionic potential and/or the pseudo core charge density
250 ! contributions have to be computed
251 !Local ionic potential:
252 ! Method 1: PAW
253 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
254  vloc_method=1;if (usepaw==0) vloc_method=2
255  if (psps%usewvl==1) vloc_method=2
256 !Pseudo core charge density:
257 ! Method 1: PAW, nc_xccc_gspace
258 ! Method 2: Norm-conserving PP, wavelets
259  coredens_method=1;if (usepaw==0) coredens_method=2
260  if (psps%nc_xccc_gspace==1) coredens_method=1
261  if (psps%nc_xccc_gspace==0) coredens_method=2
262  if (psps%usewvl==1) coredens_method=2
263 
264 !Local ionic potential and/or pseudo core charge by method 1
265  if (vloc_method==1.or.coredens_method==1) then
266    call timab(551,1,tsec)
267    optv=0;if (vloc_method==1) optv=1
268    optn=0;if (coredens_method==1) optn=n3xccc/nfft
269    optatm=0;optdyfr=0;opteltfr=0;optgr=0;optstr=1;optn2=1
270    if (coredens_method==1.and.n3xccc>0) then
271      ABI_ALLOCATE(v_dum,(nfft))
272      ABI_ALLOCATE(vxctotg,(2,nfft))
273      v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
274      call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
275      call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
276 &     comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
277      ABI_DEALLOCATE(v_dum)
278    else
279      ABI_ALLOCATE(vxctotg,(0,0))
280    end if
281    call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,&
282 &   dummy_out5,dummy_in,gmet,gprimd,dummy_out6,dummy_out7,gsqcut,&
283 &   mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
284 &   psps,pawtab,ph1d,qgrid,qprtrb_dum,rhog,corstr,lpsstr,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,vlspl,&
285 &   comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
286 &   paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
287    ABI_DEALLOCATE(vxctotg)
288    if (n3xccc==0.and.coredens_method==1) corstr=zero
289    call timab(551,2,tsec)
290  end if
291 
292 !Local ionic potential by method 2
293  if (vloc_method==2) then
294    option=3
295    ABI_ALLOCATE(dyfr_dum,(3,3,natom))
296    ABI_ALLOCATE(gr_dum,(3,natom))
297    ABI_ALLOCATE(v_dum,(nfft))
298    call mklocl_recipspace(dyfr_dum,eei,gmet,gprimd,gr_dum,gsqcut,lpsstr,mgfft,&
299 &   mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,option,paral_kgb,ph1d,qgrid,&
300 &   qprtrb_dum,rhog,ucvol,vlspl,vprtrb_dum,v_dum)
301    ABI_DEALLOCATE(dyfr_dum)
302    ABI_DEALLOCATE(gr_dum)
303    ABI_DEALLOCATE(v_dum)
304  end if
305 
306 !Pseudo core electron density by method 2
307  if (coredens_method==2) then
308    if (n1xccc/=0) then
309      call timab(55,1,tsec)
310      option=3
311      ABI_ALLOCATE(dyfr_dum,(3,3,natom))
312      ABI_ALLOCATE(gr_dum,(3,natom))
313      ABI_ALLOCATE(v_dum,(nfft))
314      icoulomb=0 ! not yet compatible with icoulomb
315      if (psps%usewvl==0.and.usepaw==0.and.icoulomb==0) then
316        if(opt_hybr==0) then
317          call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),&
318 &         n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc,&
319 &         xcccrc,xccc1d,xccc3d,xred)
320        else
321          call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),&
322 &         n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc_hf,&
323 &         xcccrc,xccc1d,xccc3d,xred)
324        end if
325      else if (psps%usewvl==0.and.(usepaw==1.or.icoulomb==1)) then
326        call mkcore_alt(atindx1,corstr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,&
327 &       nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
328 &       ucvol,vxc,xcccrc,xccc1d,xccc3d,xred,pawrad,pawtab,usepaw)
329      end if
330      ABI_DEALLOCATE(dyfr_dum)
331      ABI_DEALLOCATE(gr_dum)
332      ABI_DEALLOCATE(v_dum)
333      call timab(55,2,tsec)
334    else
335      corstr(:)=zero
336    end if
337  end if
338 
339 !=======================================================================
340 !======================= Hartree energy contribution ===================
341 !=======================================================================
342 
343  call strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd)
344 
345 !=======================================================================
346 !======================= Ewald contribution ============================
347 !=======================================================================
348 
349  call timab(38,1,tsec)
350  call ewald2(gmet,natom,ntypat,rmet,rprimd,ewestr,typat,ucvol,xred,zion)
351 
352 !=======================================================================
353 !================== VdW DFT-D contribution ============================
354 !=======================================================================
355 
356  if (vdw_xc==5) then
357    call vdw_dftd2(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_tol,&
358 &   xred,znucl,str_vdw_dftd2=vdwstr)
359  elseif (vdw_xc==6.or.vdw_xc==7) then
360    call vdw_dftd3(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_xc,&
361 &   vdw_tol,vdw_tol_3bt,xred,znucl,str_vdw_dftd3=vdwstr)
362  end if
363 
364  call timab(38,2,tsec)
365 
366 !HONG  no Berry phase contribution if using reduced ebar or d according to
367 !HONG  PRL 89, 117602 (2002) [[cite:Souza2002]]
368 !HONG  Nature Physics: M. Stengel et.al. (2009)) [[cite:Stengel1999]]
369 !=======================================================================
370 !=================== Berry phase contribution ==========================
371 !=======================================================================
372 
373 !if (berryopt==4) then
374 !berrystr_tmp(:,:) = zero
375 !Diagonal:
376 !do mu = 1, 3
377 !do ii = 1, 3
378 !berrystr_tmp(mu,mu) = berrystr_tmp(mu,mu) - &
379 !&       efield(mu)*rprimd(mu,ii)*(pel(ii) + pion(ii))/ucvol
380 !end do
381 !end do
382 !Off-diagonal (symmetrized before adding it to strten):
383 !do ii = 1, 3
384 !berrystr_tmp(3,2) = berrystr_tmp(3,2) &
385 !&     - efield(3)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol
386 !berrystr_tmp(2,3) = berrystr_tmp(2,3) &
387 !&     - efield(2)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol
388 !berrystr_tmp(3,1) = berrystr_tmp(3,1) &
389 !&     - efield(3)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol
390 !berrystr_tmp(1,3) = berrystr_tmp(1,3) &
391 !&     - efield(1)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol
392 !berrystr_tmp(2,1) = berrystr_tmp(2,1) &
393 !&     - efield(2)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol
394 !berrystr_tmp(1,2) = berrystr_tmp(1,2) &
395 !&     - efield(1)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol
396 !end do
397 !berrystr(1) = berrystr_tmp(1,1)
398 !berrystr(2) = berrystr_tmp(2,2)
399 !berrystr(3) = berrystr_tmp(3,3)
400 !berrystr(4) = (berrystr_tmp(3,2) + berrystr_tmp(2,3))/two
401 !berrystr(5) = (berrystr_tmp(3,1) + berrystr_tmp(1,3))/two
402 !berrystr(6) = (berrystr_tmp(2,1) + berrystr_tmp(1,2))/two
403 !end if
404 
405 !=======================================================================
406 !================= Other (trivial) contributions =======================
407 !=======================================================================
408 
409 !Nonlocal part of stress has already been computed
410 !(in forstrnps(norm-conserving) or pawgrnl(PAW))
411 
412 !Kinetic part of stress has already been computed
413 !(in forstrnps)
414 
415 !XC part of stress tensor has already been computed in "strsxc"
416 
417 !ii part of stress (diagonal) is trivial!
418  strsii=-eii/ucvol
419 !qvpotzero is non zero, only when usepotzero=1
420  strsii=strsii+qvpotzero/ucvol
421 
422 !======================================================================
423 !HONG  Maxwell stress when electric/displacement field is non-zero=====
424 !======================================================================
425  efield_flag = (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. &
426 & berryopt==14 .or. berryopt==16 .or. berryopt==17)
427  calc_epaw3_stress = (efield_flag .and. usepaw == 1)
428  if ( efield_flag ) then
429    ModE=dot_product(efield,efield)
430    do ii=1,3
431      Maxstr(ii)=two*efield(ii)*efield(ii)-ModE
432    end do
433    Maxstr(4)=two*efield(3)*efield(2)
434    Maxstr(5)=two*efield(3)*efield(1)
435    Maxstr(6)=two*efield(2)*efield(1)
436 !  Converting to units of Ha/Bohr^3
437 !  Maxstr(:)=Maxstr(:)*e_Cb*Bohr_Ang*1.0d-10/(Ha_J*8.0d0*pi)
438 
439    Maxstr(:)=Maxstr(:)*eps0*Ha_J*Bohr_Ang*1.0d-10/(8.0d0*pi*e_Cb**2)
440 
441    write(message, '(a,a)' )ch10,&
442 &   ' Cartesian components of Maxwell stress tensor (hartree/bohr^3)'
443    call wrtout(ab_out,message,'COLL')
444    call wrtout(std_out,  message,'COLL')
445    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
446 &   ' Maxstr(1 1)=',Maxstr(1),' Maxstr(3 2)=',Maxstr(4)
447    call wrtout(ab_out,message,'COLL')
448    call wrtout(std_out,  message,'COLL')
449    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
450 &   ' Maxstr(2 2)=',Maxstr(2),' Maxstr(3 1)=',Maxstr(5)
451    call wrtout(ab_out,message,'COLL')
452    call wrtout(std_out,  message,'COLL')
453    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
454 &   ' Maxstr(3 3)=',Maxstr(3),' Maxstr(2 1)=',Maxstr(6)
455    call wrtout(ab_out,message,'COLL')
456    call wrtout(std_out,  message,'COLL')
457    write(message, '(a)' ) ' '
458    call wrtout(ab_out,message,'COLL')
459    call wrtout(std_out,  message,'COLL')
460 
461  end if
462 
463 ! compute additional F3-type stress due to projectors for electric field with PAW
464  if ( efield_flag .and. calc_epaw3_stress ) then
465    do sdir = 1, 6
466      ep3(:) = zero
467      do idir = 1, 3
468        vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir))
469        do iatom = 1, natom
470          ep3(idir) = ep3(idir) + vol_element*dtefield%epaws3(iatom,idir,sdir)
471        end do ! end loop over atoms
472      end do ! end loop over idir (components of P)
473 ! note no appearance of ucvol here unlike in forces, stress definition includes
474 ! division by ucvol, which cancels the factor in -ucvol e . p
475      epaws3red(sdir) = -dot_product(red_efieldbar(1:3),ep3(1:3))
476    end do
477 
478 !   write(message, '(a,a)' )ch10,&
479 !&   ' Cartesian components of PAW sigma_3 stress tensor (hartree/bohr^3)'
480 !   call wrtout(ab_out,message,'COLL')
481 !   call wrtout(std_out,  message,'COLL')
482 !   write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
483 !&   ' epaws3red(1 1)=',epaws3red(1),' epaws3red(3 2)=',epaws3red(4)
484 !   call wrtout(ab_out,message,'COLL')
485 !   call wrtout(std_out,  message,'COLL')
486 !   write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
487 !&   ' epaws3red(2 2)=',epaws3red(2),' epaws3red(3 1)=',epaws3red(5)
488 !   call wrtout(ab_out,message,'COLL')
489 !   call wrtout(std_out,  message,'COLL')
490 !   write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
491 !&   ' epaws3red(3 3)=',epaws3red(3),' epaws3red(2 1)=',epaws3red(6)
492 !   call wrtout(ab_out,message,'COLL')
493 !   call wrtout(std_out,  message,'COLL')
494 !   write(message, '(a)' ) ' '
495 !   call wrtout(ab_out,message,'COLL')
496 !   call wrtout(std_out,  message,'COLL')
497 
498  end if
499 
500 !=======================================================================
501 !===== Assemble the various contributions to the stress tensor =========
502 !=======================================================================
503 !In cartesian coordinates (symmetric storage)
504 
505  strten(:)=kinstr(:)+ewestr(:)+corstr(:)+strsxc(:)+harstr(:)+lpsstr(:)+nlstr(:)
506 
507  if (usefock==1 .and. associated(fock).and.fock%fock_common%optstr) then
508    strten(:)=strten(:)+fock%fock_common%stress(:)
509  end if
510 
511 !Add contributions for constant E or D calculation.
512  if ( efield_flag ) then
513    strten(:)=strten(:)+Maxstr(:)
514    if ( calc_epaw3_stress ) strten(:) = strten(:) + epaws3red(:)
515  end if
516  if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)+vdwstr(:)
517 
518 !Additional stuff for electron-positron
519  ipositron=0
520  if (present(electronpositron)) then
521    if (associated(electronpositron)) then
522      if (allocated(electronpositron%stress_ep)) ipositron=electronpositron_calctype(electronpositron)
523    end if
524  end if
525  if (abs(ipositron)==1) then
526    strten(:)=strten(:)-harstr(:)-ewestr(:)-corstr(:)-lpsstr(:)
527    harstr(:)=zero;ewestr(:)=zero;corstr(:)=zero;strsii=zero
528    lpsstr(:)=-lpsstr(:);lpsstr(1:3)=lpsstr(1:3)-two*eei/ucvol
529    strten(:)=strten(:)+lpsstr(:)
530    if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)-vdwstr(:)
531    if (vdw_xc>=5.and.vdw_xc<=7) vdwstr(:)=zero
532  end if
533  if (abs(ipositron)==2) then
534    ABI_ALLOCATE(rhog_ep,(2,nfft))
535    ABI_ALLOCATE(dummy,(6))
536    call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
537    rhog_ep=-rhog_ep
538    call strhar(electronpositron%e_hartree,gsqcut,dummy,mpi_enreg,nfft,ngfft,rhog_ep,rprimd)
539    strten(:)=strten(:)+dummy(:);harstr(:)=harstr(:)+dummy(:)
540    ABI_DEALLOCATE(rhog_ep)
541    ABI_DEALLOCATE(dummy)
542  end if
543  if (ipositron>0) strten(:)=strten(:)+electronpositron%stress_ep(:)
544 
545 !Symmetrize resulting tensor if nsym>1
546  if (nsym>1) then
547    call stresssym(gprimd,nsym,strten,symrec)
548  end if
549 
550 !Set to zero very small values of stress
551  do mu=1,6
552    if (abs(strten(mu))<tol) strten(mu)=zero
553  end do
554 
555 !Include diagonal terms, save uncorrected stress for output
556  do mu=1,3
557    uncorr(mu)=strten(mu)+strsii
558    strten(mu)=uncorr(mu)
559  end do
560 
561 !=======================================================================
562 !================ Print out info about stress tensor ===================
563 !=======================================================================
564  if (prtvol>=10.and.ipositron>=0) then
565    write(message, '(a)' ) ' '
566    call wrtout(std_out,message,'COLL')
567    do mu=1,6
568      write(message, '(a,i5,a,1p,e22.12)' )&
569 &     ' stress: component',mu,' of hartree stress is',harstr(mu)
570      call wrtout(std_out,message,'COLL')
571    end do
572    write(message, '(a)' ) ' '
573    call wrtout(std_out,message,'COLL')
574    do mu=1,6
575      write(message, '(a,i5,a,1p,e22.12)' )&
576 &     ' stress: component',mu,' of loc psp stress is',lpsstr(mu)
577      call wrtout(std_out,message,'COLL')
578    end do
579    write(message, '(a)' ) ' '
580    call wrtout(std_out,message,'COLL')
581    do mu=1,6
582      write(message, '(a,i5,a,1p,e22.12)' )&
583 &     ' stress: component',mu,&
584 &     ' of kinetic stress is',kinstr(mu)
585      call wrtout(std_out,message,'COLL')
586    end do
587    write(message, '(a)' ) ' '
588    call wrtout(std_out,message,'COLL')
589    do mu=1,6
590      write(message, '(a,i5,a,1p,e22.12)' )&
591 &     ' stress: component',mu,' of nonlocal ps stress is',nlstr(mu)
592      call wrtout(std_out,message,'COLL')
593    end do
594    write(message, '(a)' ) ' '
595    call wrtout(std_out,message,'COLL')
596    do mu=1,6
597      write(message, '(a,i5,a,1p,e22.12)' )&
598 &     ' stress: component',mu,' of     core xc stress is',corstr(mu)
599      call wrtout(std_out,message,'COLL')
600    end do
601    write(message, '(a)' ) ' '
602    call wrtout(std_out,message,'COLL')
603    do mu=1,6
604      write(message, '(a,i5,a,1p,e22.12)' )&
605 &     ' stress: component',mu,&
606 &     ' of Ewald energ stress is',ewestr(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,' of xc stress is',strsxc(mu)
614      call wrtout(std_out,message,'COLL')
615    end do
616    if (vdw_xc>=5.and.vdw_xc<=7) then
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,&
622 &       ' of VdW DFT-D stress is',vdwstr(mu)
623        call wrtout(std_out,message,'COLL')
624      end do
625    end if
626    write(message, '(a)' ) ' '
627    call wrtout(std_out,message,'COLL')
628    write(message, '(a,1p,e22.12)' ) &
629 &   ' stress: ii (diagonal) part is',strsii
630    call wrtout(std_out,message,'COLL')
631    if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or.  &
632 &   berryopt==14 .or. berryopt==16 .or. berryopt==17) then  !!HONG
633      write(message, '(a)' ) ' '
634      call wrtout(std_out,message,'COLL')
635      do mu = 1, 6
636        write(message, '(a,i2,a,1p,e22.12)' )&
637 &       ' stress: component',mu,' of Maxwell stress is',&
638 &       Maxstr(mu)
639        call wrtout(std_out,message,'COLL')
640      end do
641    end if
642    if (ipositron/=0) then
643      write(message, '(a)' ) ' '
644      call wrtout(std_out,message,'COLL')
645      do mu=1,6
646        write(message, '(a,i5,3a,1p,e22.12)' ) &
647 &       ' stress: component',mu,' of ',EPName(abs(ipositron)), &
648 &       ' stress is',electronpositron%stress_ep(mu)
649        call wrtout(std_out,message,'COLL')
650      end do
651    end if
652 
653  end if ! prtvol
654  if (ipositron>=0) then
655    write(message, '(a,a)' )ch10,&
656 &   ' Cartesian components of stress tensor (hartree/bohr^3)'
657    call wrtout(ab_out,message,'COLL')
658    call wrtout(std_out,  message,'COLL')
659    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
660 &   '  sigma(1 1)=',strten(1),'  sigma(3 2)=',strten(4)
661    call wrtout(ab_out,message,'COLL')
662    call wrtout(std_out,  message,'COLL')
663    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
664 &   '  sigma(2 2)=',strten(2),'  sigma(3 1)=',strten(5)
665    call wrtout(ab_out,message,'COLL')
666    call wrtout(std_out,  message,'COLL')
667    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
668 &   '  sigma(3 3)=',strten(3),'  sigma(2 1)=',strten(6)
669    call wrtout(ab_out,message,'COLL')
670    call wrtout(std_out,  message,'COLL')
671    write(message, '(a)' ) ' '
672    call wrtout(ab_out,message,'COLL')
673    call wrtout(std_out,  message,'COLL')
674  end if
675  call timab(37,2,tsec)
676 
677 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=informations 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).

PARENTS

      stress

CHILDREN

      metric,ptabs_fourdp,timab,xmpi_sum

SOURCE

713 subroutine strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd,&
714 &                 rhog2) ! optional argument
715 
716 
717 !This section has been created automatically by the script Abilint (TD).
718 !Do not modify the following lines by hand.
719 #undef ABI_FUNC
720 #define ABI_FUNC 'strhar'
721 !End of the abilint section
722 
723  implicit none
724 
725 !Arguments ------------------------------------
726 !scalars
727  integer,intent(in) :: nfft
728  real(dp),intent(in) :: ehart,gsqcut
729  type(MPI_type),intent(in) :: mpi_enreg
730 !arrays
731  integer,intent(in) :: ngfft(18)
732  real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft)
733  real(dp),intent(in),optional :: rhog2(2,nfft)
734  real(dp),intent(out) :: harstr(6)
735 
736 !Local variables-------------------------------
737 !scalars
738  integer,parameter :: im=2,re=1
739  integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,me_fft,n1,n2,n3,nproc_fft
740  real(dp) :: cutoff,gsquar,rhogsq,tolfix=1.000000001_dp,ucvol
741 !arrays
742  real(dp) :: gcart(3),gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2)
743  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
744  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
745 
746 ! *************************************************************************
747 
748  call timab(568,1,tsec)
749 
750  harstr(:)=zero
751 !ehtest=0.0_dp (used for testing)
752 
753  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
754 
755  irho2=0;if (present(rhog2)) irho2=1
756 
757 !Conduct looping over all fft grid points to find G vecs inside gsqcut
758 !Include G**2 on surface of cutoff sphere as well as inside:
759  cutoff=gsqcut*tolfix
760  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
761  me_fft=ngfft(11)
762  nproc_fft=ngfft(10)
763  id1=n1/2+2
764  id2=n2/2+2
765  id3=n3/2+2
766  ii=0
767 
768  ! Get the distrib associated with this fft_grid
769  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
770 
771  do i3=1,n3
772    ig3=i3-(i3/id3)*n3-1
773    do i2=1,n2
774      ig2=i2-(i2/id2)*n2-1
775      if (fftn2_distrib(i2)==me_fft) then
776        do i1=1,n1
777          ig1=i1-(i1/id1)*n1-1
778 !        ii=ii+1
779          ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
780 !        **     GET RID OF THIS IF STATEMENT LATER for speed if needed
781 !        Avoid G=0:
782 !        if (ii>1) then
783          if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
784 !        Compute cartesian components of G
785          gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+gprimd(1,3)*dble(ig3)
786          gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+gprimd(2,3)*dble(ig3)
787          gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+gprimd(3,3)*dble(ig3)
788 !        Compute |G|^2
789          gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2
790 
791 !        Keep only G**2 inside larger cutoff (not sure this is needed):
792          if (gsquar<=cutoff) then
793 !          take |rho(G)|^2 for complex rhog
794            if (irho2==0) then
795              rhogsq=rhog(re,ii)**2+rhog(im,ii)**2
796            else
797              rhogsq=rhog(re,ii)*rhog2(re,ii)+rhog(im,ii)*rhog2(im,ii)
798            end if
799            harstr(1)=harstr(1)+(rhogsq/gsquar**2)*gcart(1)*gcart(1)
800            harstr(2)=harstr(2)+(rhogsq/gsquar**2)*gcart(2)*gcart(2)
801            harstr(3)=harstr(3)+(rhogsq/gsquar**2)*gcart(3)*gcart(3)
802            harstr(4)=harstr(4)+(rhogsq/gsquar**2)*gcart(3)*gcart(2)
803            harstr(5)=harstr(5)+(rhogsq/gsquar**2)*gcart(3)*gcart(1)
804            harstr(6)=harstr(6)+(rhogsq/gsquar**2)*gcart(2)*gcart(1)
805          end if
806 !        end if
807        end do
808      end if
809    end do
810  end do
811 
812 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel
813 #ifdef FC_IBM
814  write(std_out,*)' strhar : before mpi_comm, harstr=',harstr
815 #endif
816 
817 !Init mpi_comm
818  if(mpi_enreg%nproc_fft>1)then
819    call timab(48,1,tsec)
820    call xmpi_sum(harstr,mpi_enreg%comm_fft ,ierr)
821    call timab(48,2,tsec)
822  end if
823 
824 #ifdef FC_IBM
825 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel
826  write(std_out,*)' strhar : after mpi_comm, harstr=',harstr
827  write(std_out,*)' strhar : ehart,ucvol=',ehart,ucvol
828 #endif
829 
830 !Normalize and add term -ehart/ucvol on diagonal
831  harstr(1)=harstr(1)/pi-ehart/ucvol
832  harstr(2)=harstr(2)/pi-ehart/ucvol
833  harstr(3)=harstr(3)/pi-ehart/ucvol
834  harstr(4)=harstr(4)/pi
835  harstr(5)=harstr(5)/pi
836  harstr(6)=harstr(6)/pi
837 
838  call timab(568,2,tsec)
839 
840 end subroutine strhar