Line data Source code
1 : !!****m* ABINIT/m_stress
2 : !! NAME
3 : !! m_stress
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, FJ, MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : module m_stress
26 :
27 : use defs_basis
28 : use m_efield
29 : use m_abicore
30 : use m_errors
31 : use m_xmpi
32 : use m_extfpmd
33 :
34 : use defs_abitypes, only : MPI_type
35 : use m_time, only : timab
36 : use m_geometry, only : metric, stresssym
37 : use m_fock, only : fock_type
38 : use m_ewald, only : ewald2
39 : use defs_datatypes, only : pseudopotential_type
40 : use m_pawrad, only : pawrad_type
41 : use m_pawtab, only : pawtab_type
42 : use m_electronpositron, only : electronpositron_type,electronpositron_calctype
43 : use m_fft, only : zerosym, fourdp
44 : use m_mpinfo, only : ptabs_fourdp
45 : use m_vdw_dftd2, only : vdw_dftd2
46 : use m_vdw_dftd3, only : vdw_dftd3
47 : use m_atm2fft, only : atm2fft
48 : use m_mklocl, only : mklocl_recipspace
49 : use m_mkcore, only : mkcore, mkcore_alt
50 :
51 : #if defined(HAVE_GPU_MARKERS)
52 : use m_nvtx_data
53 : #endif
54 :
55 : implicit none
56 :
57 : private
58 : !!***
59 :
60 : public :: stress
61 : !!***
62 :
63 : contains
64 : !!***
65 :
66 : !!****f* ABINIT/stress
67 : !!
68 : !! NAME
69 : !! stress
70 : !!
71 : !! FUNCTION
72 : !! Compute the stress tensor
73 : !! strten(i,j) = (1/ucvol)*d(Etot)/(d(eps(i,j)))
74 : !! where Etot is energy per unit cell, ucvol is the unstrained unit cell
75 : !! volume, r(i,iat) is the ith position of atom iat,
76 : !! and eps(i,j) is an infinitesimal strain which maps each
77 : !! point r to r(i) -> r(i) + Sum(j) [eps(i,j)*r(j)].
78 : !!
79 : !! INPUTS
80 : !! atindx1(natom)=index table for atoms, inverse of atindx
81 : !! berryopt = 4/14: electric field is on -> add the contribution of the
82 : !! -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy
83 : !! = 6/16, or 7/17: electric displacement field is on -> add the contribution of the
84 : !! Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy
85 : !! from Etot(npw) data (at fixed geometry), used for making
86 : !! Pulay correction to stress tensor (hartree). Should be <=0.
87 : !! dtefield <type(efield_type)> = variables related to Berry phase
88 : !! eei=local pseudopotential part of Etot (hartree)
89 : !! efield = cartesian coordinates of the electric field in atomic units
90 : !! ehart=Hartree energy (hartree)
91 : !! eii=pseudoion core correction energy part of Etot (hartree)
92 : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
93 : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
94 : !! gsqcut=cutoff value on G**2 for (large) sphere inside FFT box.
95 : !! gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2)
96 : !! ixc = choice of exchange-correlation functional
97 : !! kinstr(6)=kinetic energy part of stress tensor
98 : !! mggastr(6)=meta-GGA part of stress tensor (hartree/bohr^3)
99 : !! Only non-local contribution from Div(V_tau.Grad(Psi))
100 : !! mgfft=maximum size of 1D FFTs
101 : !! mpi_enreg=information about MPI parallelization
102 : !! mqgrid=dimensioned number of q grid points for local psp spline
103 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
104 : !! n3xccc=dimension of the xccc3d array (0 or nfft).
105 : !! natom=number of atoms in cell
106 : !! nattyp(ntypat)=number of atoms of each type
107 : !! nfft=(effective) number of FFT grid points (for this processor)
108 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
109 : !! nlstr(6)=nonlocal part of stress tensor
110 : !! nspden=number of spin-density components
111 : !! nsym=number of symmetries in space group
112 : !! ntypat=number of types of atoms
113 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
114 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
115 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
116 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) array
117 : !! prtvol=integer controlling volume of printed output
118 : !! qgrid(mqgrid)=q point array for local psp spline fits
119 : !! red_efieldbar(3) = efield in reduced units relative to reciprocal lattice
120 : !! rhog(2,nfft)=Fourier transform of charge density (bohr^-3)
121 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
122 : !! strscondft(6)=cDFT correction to stress
123 : !! strsxc(6)=xc correction to stress
124 : !! symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
125 : !! typat(natom)=type integer for each atom in cell
126 : !! usefock=1 if fock operator is used; 0 otherwise.
127 : !! usekden=1 is kinetic energy density has to be taken into account, 0 otherwise
128 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
129 : !! usevxctau=1 if if XC functional depends on kinetic energy density
130 : !! vdw_tol= Van der Waals tolerance
131 : !! vdw_tol_3bt= Van der Waals tolerance on the 3-body term (only effective
132 : !! vdw_xc=6)
133 : !! vdw_xc= Van der Waals correction flag
134 : !! vlspl(mqgrid,2,ntypat)=local psp spline
135 : !! vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space
136 : !! vxctau(nfft,nspden,4*usevxctau)=(only for meta-GGA) derivative of XC energy density
137 : !! wrt kinetic energy density (depsxcdtau)
138 : !! vxc_hf(nfft,nspden)=exchange-correlation potential (hartree) in real space for Hartree-Fock corrections
139 : !! xccc1d(n1xccc*(1-usepaw),6,ntypat)=1D core charge function and five derivatives,
140 : !! for each type of atom, from psp (used in Norm-conserving only)
141 : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
142 : !! xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction
143 : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
144 : !! xred(3,natom)=reduced dimensionless atomic coordinates
145 : !! zion(ntypat)=valence charge of each type of atom
146 : !! znucl(ntypat)=atomic number of atom type
147 : !!
148 : !! OUTPUT
149 : !! strten(6)=components of the stress tensor (hartree/bohr^3) for the
150 : !! 6 unique components of this symmetric 3x3 tensor:
151 : !! Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
152 : !! The diagonal components of the returned stress tensor are
153 : !! CORRECTED for the Pulay stress.
154 : !!
155 : !! SIDE EFFECTS
156 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
157 : !!
158 : !! NOTES
159 : !! * Concerning the stress tensor:
160 : !! See O. H. Nielsen and R. M. Martin, PRB 32, 3792 (1985) [[cite:Nielsen1985a]].
161 : !! Note that first term in equation (2) should have minus sign
162 : !! (for kinetic energy contribution to stress tensor).
163 : !! Normalizations in this code differ somewhat from those employed
164 : !! by Nielsen and Martin.
165 : !! For the stress tensor contribution from the nonlocal Kleinman-Bylander
166 : !! separable pseudopotential, see D. M. Bylander, L. Kleinman, and
167 : !! S. Lee, PRB 42, 1394 (1990) [[cite:Bylander1990]].
168 : !! Again normalization conventions differ somewhat.
169 : !! See Doug Allan s notes starting page 795 (13 Jan 1992).
170 : !! * This subroutine calls different subroutines to compute the stress
171 : !! tensor contributions from the following parts of the total energy:
172 : !! (1) kinetic energy, (2) exchange-correlation energy,
173 : !! (3) Hartree energy, (4) local pseudopotential energy,
174 : !! (5) pseudoion core correction energy, (6) nonlocal pseudopotential energy,
175 : !! (7) Ewald energy.
176 : !!
177 : !! SOURCE
178 :
179 6302 : subroutine stress(atindx1,berryopt,dtefield,eei,efield,ehart,eii,fock,gsqcut,extfpmd,&
180 6302 : & ixc,kinstr,mggastr,mgfft,mpi_enreg,mqgrid,n1xccc,n3xccc,natom,nattyp,&
181 6302 : & nfft,ngfft,nlstr,nspden,nsym,ntypat,psps,pawrad,pawtab,ph1d,&
182 6302 : & prtvol,qgrid,red_efieldbar,rhog,rprimd,strten,strscondft,strsxc,symrec,&
183 6302 : & typat,usefock,usekden,usepaw,usevxctau,vdw_tol,vdw_tol_3bt,vdw_xc,&
184 6302 : & vlspl,vxc,vxctau,vxc_hf,xccc1d,xccc3d,xcctau3d,xcccrc,xred,zion,znucl,qvpotzero,&
185 : & electronpositron) ! optional argument
186 :
187 : !Arguments ------------------------------------
188 : !scalars
189 : integer,intent(in) :: berryopt,ixc,mgfft,mqgrid,n1xccc,n3xccc,natom,nfft,nspden
190 : integer,intent(in) :: nsym,ntypat,prtvol,usefock,usekden,usepaw,usevxctau,vdw_xc
191 : real(dp),intent(in) :: eei,ehart,eii,gsqcut,vdw_tol,vdw_tol_3bt,qvpotzero
192 : type(efield_type),intent(in) :: dtefield
193 : type(extfpmd_type),pointer,intent(inout) :: extfpmd
194 : type(pseudopotential_type),intent(in) :: psps
195 : type(electronpositron_type),pointer,optional :: electronpositron
196 : type(MPI_type),intent(in) :: mpi_enreg
197 : type(fock_type),pointer, intent(inout) :: fock
198 : !arrays
199 : integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),symrec(3,3,nsym)
200 : integer,intent(in) :: typat(natom)
201 : real(dp),intent(in) :: efield(3),kinstr(6),mggastr(6),nlstr(6)
202 : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid)
203 : real(dp),intent(in) :: red_efieldbar(3),rhog(2,nfft),strscondft(6),strsxc(6)
204 : real(dp),intent(in) :: vlspl(mqgrid,2,ntypat),vxc(nfft,nspden),vxctau(nfft,nspden,4*usevxctau)
205 : real(dp),allocatable,intent(in) :: vxc_hf(:,:)
206 : real(dp),intent(in) :: xccc1d(n1xccc*(1-usepaw),6,ntypat),xcccrc(ntypat)
207 : real(dp),intent(in) :: xred(3,natom),zion(ntypat),znucl(ntypat)
208 : real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*usekden),rprimd(3,3)
209 : real(dp),intent(out) :: strten(6)
210 : type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw)
211 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
212 :
213 : !Local variables-------------------------------
214 : !scalars
215 : integer :: coredens_method,coretau_method,iatom,icoulomb,idir,ii,ipositron,mu,nkpt=1
216 : integer :: optatm,optdyfr,opteltfr,opt_hybr,optgr,option,optn,optn2,optstr,optv,sdir,vloc_method
217 : real(dp),parameter :: tol=1.0d-15
218 : real(dp) :: e_dum,dum_rcut=zero,strsii,ucvol,vol_element
219 : character(len=500) :: message
220 : logical :: calc_epaw3_stress, efield_flag
221 : !arrays
222 : integer :: qprtrb_dum(3),icutcoul=3
223 : real(dp) :: corstr(6),dumstr(6),ep3(3),epaws3red(6),ewestr(6),gmet(3,3),vcutgeo(3)
224 : real(dp) :: gprimd(3,3),harstr(6),lpsstr(6),rmet(3,3),taustr(6),tsec(2),uncorr(3)
225 : real(dp) :: vdwstr(6),vprtrb_dum(2)
226 : real(dp) :: Maxstr(6),ModE !Maxwell-stress constribution, and magnitude of efield
227 : real(dp) :: dummy_in(0)
228 : real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0)
229 6302 : real(dp),allocatable :: dummy(:),dyfr_dum(:,:,:),gr_dum(:,:),rhog_ep(:,:),v_dum(:)
230 6302 : real(dp),allocatable :: vxctotg(:,:)
231 : character(len=10) :: EPName(1:2)=(/"Electronic","Positronic"/)
232 : ! *************************************************************************
233 :
234 6302 : call timab(37,1,tsec)
235 : ABI_NVTX_START_RANGE(NVTX_STRESS)
236 :
237 : !Compute different geometric tensor, as well as ucvol, from rprimd
238 6302 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
239 :
240 6302 : opt_hybr=0;if (allocated(vxc_hf)) opt_hybr=1
241 6302 : icoulomb=0 ! not yet compatible with icoulomb
242 :
243 : !=======================================================================
244 : !========= Local pseudopotential and core charge contributions =========
245 : !=======================================================================
246 :
247 : !Determine by which method the local ionic potential and/or the pseudo core charge density
248 : ! contributions have to be computed
249 : !Local ionic potential:
250 : ! Method 1: PAW
251 : ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
252 6302 : vloc_method=1;if (usepaw==0) vloc_method=2
253 6302 : if (psps%usewvl==1) vloc_method=2
254 : !Pseudo core charge density:
255 : ! Method 1: PAW, nc_xccc_gspace
256 : ! Method 2: Norm-conserving PP, wavelets
257 : coredens_method=1;if (usepaw==0) coredens_method=2
258 6302 : if (psps%nc_xccc_gspace==1) coredens_method=1
259 6302 : if (psps%nc_xccc_gspace==0) coredens_method=2
260 6302 : if (psps%usewvl==1) coredens_method=2
261 6302 : coretau_method=0
262 6302 : if (usekden==1) then
263 60 : coretau_method=1;if (psps%nc_xccc_gspace==0) coretau_method=2
264 : end if
265 :
266 : !Local ionic potential and/or pseudo core charge by method 1
267 6302 : if (vloc_method==1.or.coredens_method==1.or.coretau_method==1) then
268 1398 : call timab(551,1,tsec)
269 : ! Compute Vxc in reciprocal space
270 1398 : if (coredens_method==1.and.n3xccc>0) then
271 2823 : ABI_MALLOC(v_dum,(nfft))
272 2823 : ABI_MALLOC(vxctotg,(2,nfft))
273 34572887 : v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
274 941 : call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0)
275 : call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
276 941 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
277 941 : ABI_FREE(v_dum)
278 : else
279 457 : ABI_MALLOC(vxctotg,(0,0))
280 : end if
281 : ! Compute contribution to stresses from Vloc and/or pseudo core density
282 1398 : optv=0;if (vloc_method==1) optv=1
283 1398 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
284 1398 : optatm=0;optdyfr=0;opteltfr=0;optgr=0;optstr=1;optn2=1
285 1398 : if (vloc_method==1.or.coredens_method==1) then
286 : call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,&
287 : & dummy_out5,dummy_in,gmet,gprimd,dummy_out6,dummy_out7,gsqcut,&
288 : & mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
289 : & psps,pawtab,ph1d,qgrid,qprtrb_dum,dum_rcut,rhog,rprimd,corstr,lpsstr,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,vlspl,&
290 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
291 1398 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
292 : end if
293 1398 : if (n3xccc==0.and.coredens_method==1) corstr=zero
294 1398 : ABI_FREE(vxctotg)
295 1398 : if (usekden==1.and.usevxctau==1.and.coretau_method==1..and.n3xccc>0) then
296 : ! Compute contribution to stresses from pseudo kinetic energy core density
297 29 : optv=0;optn=1;optn2=4
298 87 : ABI_MALLOC(v_dum,(nfft))
299 87 : ABI_MALLOC(vxctotg,(2,nfft))
300 1353053 : v_dum(:)=vxctau(:,1,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxctau(:,2,1))
301 29 : call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0)
302 : call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
303 29 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
304 29 : ABI_FREE(v_dum)
305 : call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,&
306 : & dummy_out5,dummy_in,gmet,gprimd,dummy_out6,dummy_out7,gsqcut,&
307 : & mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
308 : & psps,pawtab,ph1d,qgrid,qprtrb_dum,dum_rcut,rhog,rprimd,taustr,dumstr,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,vlspl,&
309 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
310 29 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
311 203 : corstr(1:6)=corstr(1:6)+taustr(1:6)
312 29 : ABI_FREE(vxctotg)
313 : end if
314 1398 : call timab(551,2,tsec)
315 : end if
316 :
317 : !Local ionic potential by method 2
318 6302 : if (vloc_method==2) then
319 4907 : option=3
320 14721 : ABI_MALLOC(dyfr_dum,(3,3,natom))
321 14721 : ABI_MALLOC(gr_dum,(3,natom))
322 14721 : ABI_MALLOC(v_dum,(nfft))
323 : call mklocl_recipspace(dyfr_dum,eei,gmet,gprimd,gr_dum,gsqcut,icutcoul,lpsstr,mgfft,&
324 : & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,nkpt,ntypat,option,ph1d,qgrid,&
325 4907 : & qprtrb_dum,dum_rcut,rhog,rprimd,ucvol,vcutgeo,vlspl,vprtrb_dum,v_dum)
326 4907 : ABI_FREE(dyfr_dum)
327 4907 : ABI_FREE(gr_dum)
328 4907 : ABI_FREE(v_dum)
329 : end if
330 :
331 : !Pseudo core electron density by method 2
332 6302 : if (coredens_method==2.or.coretau_method==2) then
333 4904 : if (n1xccc/=0) then
334 1809 : call timab(55,1,tsec)
335 1809 : option=3
336 5427 : ABI_MALLOC(dyfr_dum,(3,3,natom))
337 5427 : ABI_MALLOC(gr_dum,(3,natom))
338 5427 : ABI_MALLOC(v_dum,(nfft))
339 1809 : if (coredens_method==2) then
340 1809 : if (psps%usewvl==0.and.usepaw==0.and.icoulomb==0) then
341 1809 : if(opt_hybr==0) then
342 : call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),&
343 : & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc,&
344 1779 : & xcccrc,xccc1d,xccc3d,xred)
345 : else
346 : call mkcore(corstr,dyfr_dum,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,ngfft(1),&
347 : & n1xccc,ngfft(2),ngfft(3),option,rprimd,typat,ucvol,vxc_hf,&
348 30 : & xcccrc,xccc1d,xccc3d,xred)
349 : end if
350 0 : else if (psps%usewvl==0.and.(usepaw==1.or.icoulomb==1)) then
351 : call mkcore_alt(atindx1,corstr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,&
352 : & nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
353 0 : & ucvol,vxc,xcccrc,xccc1d,xccc3d,xred,pawrad,pawtab,usepaw)
354 : end if
355 : end if
356 1809 : if (usekden==1.and.usevxctau==1.and.coretau_method==2) then
357 : call mkcore_alt(atindx1,taustr,dyfr_dum,gr_dum,icoulomb,mpi_enreg,natom,nfft,&
358 : & nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
359 : & ucvol,vxctau(:,:,1),xcccrc,xccc1d,xcctau3d,xred,pawrad,pawtab,usepaw,&
360 3 : & usekden=.true.)
361 : end if
362 1809 : ABI_FREE(dyfr_dum)
363 1809 : ABI_FREE(gr_dum)
364 1809 : ABI_FREE(v_dum)
365 1809 : call timab(55,2,tsec)
366 : else
367 3095 : corstr(:)=zero
368 : end if
369 : end if
370 :
371 : !=======================================================================
372 : !======================= Hartree energy contribution ===================
373 : !=======================================================================
374 :
375 6302 : call strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd)
376 :
377 : !=======================================================================
378 : !======================= Ewald contribution ============================
379 : !=======================================================================
380 :
381 6302 : call timab(38,1,tsec)
382 6302 : call ewald2(gmet,natom,ntypat,rmet,rprimd,ewestr,typat,ucvol,xred,zion)
383 :
384 : !=======================================================================
385 : !================== VdW DFT-D contribution ============================
386 : !=======================================================================
387 :
388 6302 : if (vdw_xc==5) then
389 : call vdw_dftd2(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_tol,&
390 10 : & xred,znucl,str_vdw_dftd2=vdwstr)
391 6292 : elseif (vdw_xc==6.or.vdw_xc==7) then
392 : call vdw_dftd3(e_dum,ixc,natom,ntypat,0,typat,rprimd,vdw_xc,&
393 7 : & vdw_tol,vdw_tol_3bt,xred,znucl,str_vdw_dftd3=vdwstr)
394 : end if
395 :
396 6302 : call timab(38,2,tsec)
397 :
398 : !HONG no Berry phase contribution if using reduced ebar or d according to
399 : !HONG PRL 89, 117602 (2002) [[cite:Souza2002]]
400 : !HONG Nature Physics: M. Stengel et.al. (2009)) [[cite:Stengel1999]]
401 : !=======================================================================
402 : !=================== Berry phase contribution ==========================
403 : !=======================================================================
404 :
405 : !if (berryopt==4) then
406 : !berrystr_tmp(:,:) = zero
407 : !Diagonal:
408 : !do mu = 1, 3
409 : !do ii = 1, 3
410 : !berrystr_tmp(mu,mu) = berrystr_tmp(mu,mu) - &
411 : !& efield(mu)*rprimd(mu,ii)*(pel(ii) + pion(ii))/ucvol
412 : !end do
413 : !end do
414 : !Off-diagonal (symmetrized before adding it to strten):
415 : !do ii = 1, 3
416 : !berrystr_tmp(3,2) = berrystr_tmp(3,2) &
417 : !& - efield(3)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol
418 : !berrystr_tmp(2,3) = berrystr_tmp(2,3) &
419 : !& - efield(2)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol
420 : !berrystr_tmp(3,1) = berrystr_tmp(3,1) &
421 : !& - efield(3)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol
422 : !berrystr_tmp(1,3) = berrystr_tmp(1,3) &
423 : !& - efield(1)*rprimd(3,ii)*(pel(ii) + pion(ii))/ucvol
424 : !berrystr_tmp(2,1) = berrystr_tmp(2,1) &
425 : !& - efield(2)*rprimd(1,ii)*(pel(ii) + pion(ii))/ucvol
426 : !berrystr_tmp(1,2) = berrystr_tmp(1,2) &
427 : !& - efield(1)*rprimd(2,ii)*(pel(ii) + pion(ii))/ucvol
428 : !end do
429 : !berrystr(1) = berrystr_tmp(1,1)
430 : !berrystr(2) = berrystr_tmp(2,2)
431 : !berrystr(3) = berrystr_tmp(3,3)
432 : !berrystr(4) = (berrystr_tmp(3,2) + berrystr_tmp(2,3))/two
433 : !berrystr(5) = (berrystr_tmp(3,1) + berrystr_tmp(1,3))/two
434 : !berrystr(6) = (berrystr_tmp(2,1) + berrystr_tmp(1,2))/two
435 : !end if
436 :
437 : !=======================================================================
438 : !================= Other (trivial) contributions =======================
439 : !=======================================================================
440 :
441 : !Nonlocal part of stress has already been computed
442 : !(in forstrnps(norm-conserving) or pawgrnl(PAW))
443 :
444 : !Kinetic part of stress has already been computed
445 : !(in forstrnps)
446 :
447 : !cDFT part of stress tensor has already been computed in "constrained_residual"
448 :
449 : !XC part of stress tensor has already been computed in "strsxc"
450 :
451 : !ii part of stress (diagonal) is trivial!
452 6302 : strsii=-eii/ucvol
453 : !qvpotzero is non zero, only when usepotzero=1
454 6302 : strsii=strsii+qvpotzero/ucvol
455 :
456 : !======================================================================
457 : !HONG Maxwell stress when electric/displacement field is non-zero=====
458 : !======================================================================
459 : efield_flag = (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. &
460 6302 : & berryopt==14 .or. berryopt==16 .or. berryopt==17)
461 6343 : calc_epaw3_stress = (efield_flag .and. usepaw == 1)
462 : if ( efield_flag ) then
463 100 : ModE=dot_product(efield,efield)
464 100 : do ii=1,3
465 100 : Maxstr(ii)=two*efield(ii)*efield(ii)-ModE
466 : end do
467 25 : Maxstr(4)=two*efield(3)*efield(2)
468 25 : Maxstr(5)=two*efield(3)*efield(1)
469 25 : Maxstr(6)=two*efield(2)*efield(1)
470 : ! Converting to units of Ha/Bohr^3
471 : ! Maxstr(:)=Maxstr(:)*e_Cb*Bohr_Ang*1.0d-10/(Ha_J*8.0d0*pi)
472 :
473 175 : Maxstr(:)=Maxstr(:)*eps0*Ha_J*Bohr_Ang*1.0d-10/(8.0d0*pi*e_Cb**2)
474 :
475 25 : write(message, '(a,a)' )ch10,&
476 50 : & ' Cartesian components of Maxwell stress tensor (hartree/bohr^3)'
477 25 : call wrtout(ab_out,message,'COLL')
478 25 : call wrtout(std_out, message,'COLL')
479 : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
480 25 : & ' Maxstr(1 1)=',Maxstr(1),' Maxstr(3 2)=',Maxstr(4)
481 25 : call wrtout(ab_out,message,'COLL')
482 25 : call wrtout(std_out, message,'COLL')
483 : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
484 25 : & ' Maxstr(2 2)=',Maxstr(2),' Maxstr(3 1)=',Maxstr(5)
485 25 : call wrtout(ab_out,message,'COLL')
486 25 : call wrtout(std_out, message,'COLL')
487 : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
488 25 : & ' Maxstr(3 3)=',Maxstr(3),' Maxstr(2 1)=',Maxstr(6)
489 25 : call wrtout(ab_out,message,'COLL')
490 25 : call wrtout(std_out, message,'COLL')
491 25 : write(message, '(a)' ) ' '
492 25 : call wrtout(ab_out,message,'COLL')
493 25 : call wrtout(std_out, message,'COLL')
494 :
495 : end if
496 :
497 : ! compute additional F3-type stress due to projectors for electric field with PAW
498 25 : if ( efield_flag .and. calc_epaw3_stress ) then
499 63 : do sdir = 1, 6
500 54 : ep3(:) = zero
501 216 : do idir = 1, 3
502 162 : vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir))
503 540 : do iatom = 1, natom
504 486 : ep3(idir) = ep3(idir) + vol_element*dtefield%epaws3(iatom,idir,sdir)
505 : end do ! end loop over atoms
506 : end do ! end loop over idir (components of P)
507 : ! note no appearance of ucvol here unlike in forces, stress definition includes
508 : ! division by ucvol, which cancels the factor in -ucvol e . p
509 225 : epaws3red(sdir) = -dot_product(red_efieldbar(1:3),ep3(1:3))
510 : end do
511 :
512 : ! write(message, '(a,a)' )ch10,&
513 : !& ' Cartesian components of PAW sigma_3 stress tensor (hartree/bohr^3)'
514 : ! call wrtout(ab_out,message,'COLL')
515 : ! call wrtout(std_out, message,'COLL')
516 : ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
517 : !& ' epaws3red(1 1)=',epaws3red(1),' epaws3red(3 2)=',epaws3red(4)
518 : ! call wrtout(ab_out,message,'COLL')
519 : ! call wrtout(std_out, message,'COLL')
520 : ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
521 : !& ' epaws3red(2 2)=',epaws3red(2),' epaws3red(3 1)=',epaws3red(5)
522 : ! call wrtout(ab_out,message,'COLL')
523 : ! call wrtout(std_out, message,'COLL')
524 : ! write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
525 : !& ' epaws3red(3 3)=',epaws3red(3),' epaws3red(2 1)=',epaws3red(6)
526 : ! call wrtout(ab_out,message,'COLL')
527 : ! call wrtout(std_out, message,'COLL')
528 : ! write(message, '(a)' ) ' '
529 : ! call wrtout(ab_out,message,'COLL')
530 : ! call wrtout(std_out, message,'COLL')
531 :
532 : end if
533 :
534 : !=======================================================================
535 : !===== Assemble the various contributions to the stress tensor =========
536 : !=======================================================================
537 : !In cartesian coordinates (symmetric storage)
538 :
539 44114 : strten(:)=kinstr(:)+ewestr(:)+corstr(:)+strscondft(:)+strsxc(:)+harstr(:)+lpsstr(:)+nlstr(:)+mggastr(:)
540 :
541 6302 : if (usefock==1 .and. associated(fock)) then
542 121 : if (fock%fock_common%optstr) then
543 847 : strten(:)=strten(:)+fock%fock_common%stress(:)
544 : end if
545 : end if
546 :
547 : !Add contributions for constant E or D calculation.
548 6302 : if ( efield_flag ) then
549 175 : strten(:)=strten(:)+Maxstr(:)
550 79 : if ( calc_epaw3_stress ) strten(:) = strten(:) + epaws3red(:)
551 : end if
552 6404 : if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)+vdwstr(:)
553 :
554 : !Additional stuff for electron-positron
555 6302 : ipositron=0
556 6302 : if (present(electronpositron)) then
557 6302 : if (associated(electronpositron)) then
558 1 : if (allocated(electronpositron%stress_ep)) ipositron=electronpositron_calctype(electronpositron)
559 : end if
560 : end if
561 1 : if (abs(ipositron)==1) then
562 7 : strten(:)=strten(:)-harstr(:)-ewestr(:)-corstr(:)-lpsstr(:)
563 1 : harstr(:)=zero;ewestr(:)=zero;corstr(:)=zero;strsii=zero
564 10 : lpsstr(:)=-lpsstr(:);lpsstr(1:3)=lpsstr(1:3)-two*eei/ucvol
565 7 : strten(:)=strten(:)+lpsstr(:)
566 1 : if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)-vdwstr(:)
567 1 : if (vdw_xc>=5.and.vdw_xc<=7) vdwstr(:)=zero
568 : end if
569 6302 : if (abs(ipositron)==2) then
570 0 : ABI_MALLOC(rhog_ep,(2,nfft))
571 0 : ABI_MALLOC(dummy,(6))
572 0 : call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,1,ngfft,0)
573 0 : rhog_ep=-rhog_ep
574 0 : call strhar(electronpositron%e_hartree,gsqcut,dummy,mpi_enreg,nfft,ngfft,rhog_ep,rprimd)
575 0 : strten(:)=strten(:)+dummy(:);harstr(:)=harstr(:)+dummy(:)
576 0 : ABI_FREE(rhog_ep)
577 0 : ABI_FREE(dummy)
578 : end if
579 6308 : if (ipositron>0) strten(:)=strten(:)+electronpositron%stress_ep(:)
580 :
581 : !Symmetrize resulting tensor if nsym>1
582 6302 : if (nsym>1) then
583 4850 : call stresssym(gprimd,nsym,strten,symrec)
584 : end if
585 :
586 : !Set to zero very small values of stress
587 44114 : do mu=1,6
588 44114 : if (abs(strten(mu))<tol) strten(mu)=zero
589 : end do
590 :
591 : !Include diagonal terms, save uncorrected stress for output
592 25208 : do mu=1,3
593 18906 : uncorr(mu)=strten(mu)+strsii
594 25208 : strten(mu)=uncorr(mu)
595 : end do
596 :
597 : !Adding the extfpmd continous contribution to stress tensor
598 6302 : if(associated(extfpmd)) then
599 28 : strten(1:3)=strten(1:3)-(2./3.)*extfpmd%e_kinetic/extfpmd%ucvol
600 : end if
601 :
602 : !=======================================================================
603 : !================ Print out info about stress tensor ===================
604 : !=======================================================================
605 6302 : if (prtvol>=10.and.ipositron>=0) then
606 244 : write(message, '(a)' ) ' '
607 244 : call wrtout(std_out,message,'COLL')
608 1708 : do mu=1,6
609 : write(message, '(a,i5,a,1p,e22.12)' )&
610 1464 : & ' stress: component',mu,' of hartree stress is',harstr(mu)
611 1708 : call wrtout(std_out,message,'COLL')
612 : end do
613 244 : write(message, '(a)' ) ' '
614 244 : call wrtout(std_out,message,'COLL')
615 1708 : do mu=1,6
616 : write(message, '(a,i5,a,1p,e22.12)' )&
617 1464 : & ' stress: component',mu,' of loc psp stress is',lpsstr(mu)
618 1708 : call wrtout(std_out,message,'COLL')
619 : end do
620 244 : write(message, '(a)' ) ' '
621 244 : call wrtout(std_out,message,'COLL')
622 1708 : do mu=1,6
623 : write(message, '(a,i5,a,1p,e22.12)' )&
624 1464 : & ' stress: component',mu,&
625 2928 : & ' of kinetic stress is',kinstr(mu)
626 1708 : call wrtout(std_out,message,'COLL')
627 : end do
628 244 : if (usekden>0) then
629 0 : write(message, '(a)' ) ' '
630 0 : call wrtout(std_out,message,'COLL')
631 0 : do mu=1,6
632 : write(message, '(a,i5,a,1p,e22.12)' )&
633 0 : & ' stress: component',mu,&
634 0 : & ' of metaGGA stress is',mggastr(mu)
635 0 : call wrtout(std_out,message,'COLL')
636 : end do
637 : end if
638 244 : write(message, '(a)' ) ' '
639 244 : call wrtout(std_out,message,'COLL')
640 1708 : do mu=1,6
641 : write(message, '(a,i5,a,1p,e22.12)' )&
642 1464 : & ' stress: component',mu,' of nonlocal ps stress is',nlstr(mu)
643 1708 : call wrtout(std_out,message,'COLL')
644 : end do
645 244 : write(message, '(a)' ) ' '
646 244 : call wrtout(std_out,message,'COLL')
647 1708 : do mu=1,6
648 : write(message, '(a,i5,a,1p,e22.12)' )&
649 1464 : & ' stress: component',mu,' of core xc stress is',corstr(mu)
650 1708 : call wrtout(std_out,message,'COLL')
651 : end do
652 244 : write(message, '(a)' ) ' '
653 244 : call wrtout(std_out,message,'COLL')
654 1708 : do mu=1,6
655 : write(message, '(a,i5,a,1p,e22.12)' )&
656 1464 : & ' stress: component',mu,&
657 2928 : & ' of Ewald energ stress is',ewestr(mu)
658 1708 : call wrtout(std_out,message,'COLL')
659 : end do
660 244 : write(message, '(a)' ) ' '
661 244 : call wrtout(std_out,message,'COLL')
662 1708 : do mu=1,6
663 : write(message, '(a,i5,a,1p,e22.12)' ) &
664 1464 : & ' stress: component',mu,' of xc stress is',strsxc(mu)
665 1708 : call wrtout(std_out,message,'COLL')
666 : end do
667 :
668 1708 : if( any( abs(strscondft(:))>tol8 ) )then
669 0 : write(message, '(a)' ) ' '
670 0 : call wrtout(std_out,message,'COLL')
671 0 : do mu=1,6
672 : write(message, '(a,i5,a,1p,e22.12)' ) &
673 0 : & ' stress: component',mu,' of cDFT stress is',strscondft(mu)
674 0 : call wrtout(std_out,message,'COLL')
675 : end do
676 : endif
677 :
678 244 : if (vdw_xc>=5.and.vdw_xc<=7) then
679 0 : write(message, '(a)' ) ' '
680 0 : call wrtout(std_out,message,'COLL')
681 0 : do mu=1,6
682 : write(message, '(a,i5,a,1p,e22.12)' )&
683 0 : & ' stress: component',mu,&
684 0 : & ' of VdW DFT-D stress is',vdwstr(mu)
685 0 : call wrtout(std_out,message,'COLL')
686 : end do
687 : end if
688 244 : write(message, '(a)' ) ' '
689 244 : call wrtout(std_out,message,'COLL')
690 : write(message, '(a,1p,e22.12)' ) &
691 244 : & ' stress: ii (diagonal) part is',strsii
692 244 : call wrtout(std_out,message,'COLL')
693 : if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. &
694 : & berryopt==14 .or. berryopt==16 .or. berryopt==17) then !!HONG
695 0 : write(message, '(a)' ) ' '
696 0 : call wrtout(std_out,message,'COLL')
697 0 : do mu = 1, 6
698 : write(message, '(a,i2,a,1p,e22.12)' )&
699 0 : & ' stress: component',mu,' of Maxwell stress is',&
700 0 : & Maxstr(mu)
701 0 : call wrtout(std_out,message,'COLL')
702 : end do
703 : end if
704 244 : if (ipositron/=0) then
705 0 : write(message, '(a)' ) ' '
706 0 : call wrtout(std_out,message,'COLL')
707 0 : do mu=1,6
708 : write(message, '(a,i5,3a,1p,e22.12)' ) &
709 0 : & ' stress: component',mu,' of ',EPName(abs(ipositron)), &
710 0 : & ' stress is',electronpositron%stress_ep(mu)
711 0 : call wrtout(std_out,message,'COLL')
712 : end do
713 : end if
714 :
715 : end if ! prtvol
716 6058 : if (ipositron>=0) then
717 6302 : write(message, '(a,a)' )ch10,&
718 12604 : & ' Cartesian components of stress tensor (hartree/bohr^3)'
719 6302 : call wrtout(ab_out,message,'COLL')
720 6302 : call wrtout(std_out, message,'COLL')
721 : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
722 6302 : & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4)
723 6302 : call wrtout(ab_out,message,'COLL')
724 6302 : call wrtout(std_out, message,'COLL')
725 : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
726 6302 : & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5)
727 6302 : call wrtout(ab_out,message,'COLL')
728 6302 : call wrtout(std_out, message,'COLL')
729 : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
730 6302 : & ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6)
731 6302 : call wrtout(ab_out,message,'COLL')
732 6302 : call wrtout(std_out, message,'COLL')
733 6302 : write(message, '(a)' ) ' '
734 6302 : call wrtout(ab_out,message,'COLL')
735 6302 : call wrtout(std_out, message,'COLL')
736 : end if
737 : ABI_NVTX_END_RANGE()
738 6302 : call timab(37,2,tsec)
739 :
740 6302 : end subroutine stress
741 : !!***
742 :
743 : !!****f* ABINIT/strhar
744 : !!
745 : !! NAME
746 : !! strhar
747 : !!
748 : !! FUNCTION
749 : !! Compute Hartree energy contribution to stress tensor (Cartesian coordinates).
750 : !!
751 : !! INPUTS
752 : !! ehart=Hartree energy (hartree)
753 : !! gsqcut=cutoff value on $G^2$ for (large) sphere inside fft box.
754 : !! $gsqcut=(boxcut^2)*ecut/(2._dp*(\pi^2))$
755 : !! mpi_enreg=information about MPI parallelization
756 : !! nfft=(effective) number of FFT grid points (for this processor)
757 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
758 : !! rhog(2,nfft)=Fourier transform of charge density (bohr^-3)
759 : !! rhog(2,nfft)= optional argument: Fourier transform of a second charge density (bohr^-3)
760 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
761 : !!
762 : !! OUTPUT
763 : !! harstr(6)=components of Hartree part of stress tensor
764 : !! (Cartesian coordinates, symmetric tensor) in hartree/bohr^3
765 : !! Definition of symmetric tensor storage: store 6 unique components
766 : !! in the order 11, 22, 33, 32, 31, 21 (suggested by Xavier Gonze).
767 : !!
768 : !! SOURCE
769 :
770 6302 : subroutine strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd,&
771 : & rhog2) ! optional argument
772 :
773 : !Arguments ------------------------------------
774 : !scalars
775 : integer,intent(in) :: nfft
776 : real(dp),intent(in) :: ehart,gsqcut
777 : type(MPI_type),intent(in) :: mpi_enreg
778 : !arrays
779 : integer,intent(in) :: ngfft(18)
780 : real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft)
781 : real(dp),intent(in),optional :: rhog2(2,nfft)
782 : real(dp),intent(out) :: harstr(6)
783 :
784 : !Local variables-------------------------------
785 : !scalars
786 : integer,parameter :: im=2,re=1
787 : integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,me_fft,n1,n2,n3,nproc_fft
788 : real(dp) :: cutoff,gsquar,rhogsq,tolfix=1.000000001_dp,ucvol
789 : !arrays
790 : real(dp) :: gcart(3),gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2)
791 6302 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
792 6302 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
793 : ! *************************************************************************
794 :
795 6302 : call timab(568,1,tsec)
796 :
797 6302 : harstr(:)=zero
798 : !ehtest=0.0_dp (used for testing)
799 :
800 6302 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
801 :
802 6302 : irho2=0;if (present(rhog2)) irho2=1
803 :
804 : !Conduct looping over all fft grid points to find G vecs inside gsqcut
805 : !Include G**2 on surface of cutoff sphere as well as inside:
806 6302 : cutoff=gsqcut*tolfix
807 6302 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
808 6302 : me_fft=ngfft(11)
809 6302 : nproc_fft=ngfft(10)
810 6302 : id1=n1/2+2
811 6302 : id2=n2/2+2
812 6302 : id3=n3/2+2
813 6302 : ii=0
814 :
815 : ! Get the distrib associated with this fft_grid
816 6302 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
817 :
818 147626 : do i3=1,n3
819 141324 : ig3=i3-(i3/id3)*n3-1
820 3993293 : do i2=1,n2
821 3845667 : ig2=i2-(i2/id2)*n2-1
822 3986991 : if (fftn2_distrib(i2)==me_fft) then
823 140673931 : do i1=1,n1
824 137080160 : ig1=i1-(i1/id1)*n1-1
825 : ! ii=ii+1
826 137080160 : ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
827 : ! ** GET RID OF THIS IF STATEMENT LATER for speed if needed
828 : ! Avoid G=0:
829 : ! if (ii>1) then
830 137080160 : if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
831 : ! Compute cartesian components of G
832 137074080 : gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+gprimd(1,3)*dble(ig3)
833 137074080 : gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+gprimd(2,3)*dble(ig3)
834 137074080 : gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+gprimd(3,3)*dble(ig3)
835 : ! Compute |G|^2
836 137074080 : gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2
837 :
838 : ! Keep only G**2 inside larger cutoff (not sure this is needed):
839 140667851 : if (gsquar<=cutoff) then
840 : ! take |rho(G)|^2 for complex rhog
841 55483159 : if (irho2==0) then
842 55483159 : rhogsq=rhog(re,ii)**2+rhog(im,ii)**2
843 : else
844 0 : rhogsq=rhog(re,ii)*rhog2(re,ii)+rhog(im,ii)*rhog2(im,ii)
845 : end if
846 55483159 : harstr(1)=harstr(1)+(rhogsq/gsquar**2)*gcart(1)*gcart(1)
847 55483159 : harstr(2)=harstr(2)+(rhogsq/gsquar**2)*gcart(2)*gcart(2)
848 55483159 : harstr(3)=harstr(3)+(rhogsq/gsquar**2)*gcart(3)*gcart(3)
849 55483159 : harstr(4)=harstr(4)+(rhogsq/gsquar**2)*gcart(3)*gcart(2)
850 55483159 : harstr(5)=harstr(5)+(rhogsq/gsquar**2)*gcart(3)*gcart(1)
851 55483159 : harstr(6)=harstr(6)+(rhogsq/gsquar**2)*gcart(2)*gcart(1)
852 : end if
853 : ! end if
854 : end do
855 : end if
856 : end do
857 : end do
858 :
859 : !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel
860 : #ifdef FC_IBM
861 : write(std_out,*)' strhar : before mpi_comm, harstr=',harstr
862 : #endif
863 :
864 : !Init mpi_comm
865 6302 : if(mpi_enreg%nproc_fft>1)then
866 436 : call timab(48,1,tsec)
867 436 : call xmpi_sum(harstr,mpi_enreg%comm_fft ,ierr)
868 436 : call timab(48,2,tsec)
869 : end if
870 :
871 : #ifdef FC_IBM
872 : !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel
873 : write(std_out,*)' strhar : after mpi_comm, harstr=',harstr
874 : write(std_out,*)' strhar : ehart,ucvol=',ehart,ucvol
875 : #endif
876 :
877 : !Normalize and add term -ehart/ucvol on diagonal
878 6302 : harstr(1)=harstr(1)/pi-ehart/ucvol
879 6302 : harstr(2)=harstr(2)/pi-ehart/ucvol
880 6302 : harstr(3)=harstr(3)/pi-ehart/ucvol
881 6302 : harstr(4)=harstr(4)/pi
882 6302 : harstr(5)=harstr(5)/pi
883 6302 : harstr(6)=harstr(6)/pi
884 :
885 6302 : call timab(568,2,tsec)
886 :
887 6302 : end subroutine strhar
888 : !!***
889 :
890 : end module m_stress
891 : !!***
|
