Branch data 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 : 6329 : subroutine stress(atindx1,berryopt,dtefield,eei,efield,ehart,eii,fock,gsqcut,extfpmd,&
180 : 6329 : & ixc,kinstr,mggastr,mgfft,mpi_enreg,mqgrid,n1xccc,n3xccc,natom,nattyp,&
181 : 6329 : & nfft,ngfft,nlstr,nspden,nsym,ntypat,psps,pawrad,pawtab,ph1d,&
182 : 6329 : & prtvol,qgrid,red_efieldbar,rhog,rprimd,strten,strscondft,strsxc,symrec,&
183 : 6329 : & typat,usefock,usekden,usepaw,usevxctau,vdw_tol,vdw_tol_3bt,vdw_xc,&
184 : 6329 : & 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 : 6329 : real(dp),allocatable :: dummy(:),dyfr_dum(:,:,:),gr_dum(:,:),rhog_ep(:,:),v_dum(:)
230 : 6329 : real(dp),allocatable :: vxctotg(:,:)
231 : : character(len=10) :: EPName(1:2)=(/"Electronic","Positronic"/)
232 : : ! *************************************************************************
233 : :
234 : 6329 : 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 : 6329 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
239 : :
240 [ + + ]: 6329 : opt_hybr=0;if (allocated(vxc_hf)) opt_hybr=1
241 : 6329 : 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 [ + + ]: 6329 : vloc_method=1;if (usepaw==0) vloc_method=2
253 [ - + ]: 6329 : 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 [ + + ]: 6329 : if (psps%nc_xccc_gspace==1) coredens_method=1
259 [ + + ]: 6329 : if (psps%nc_xccc_gspace==0) coredens_method=2
260 [ - + ]: 6329 : if (psps%usewvl==1) coredens_method=2
261 : 6329 : coretau_method=0
262 [ + + ]: 6329 : 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 [ + + - + ]: 6329 : if (vloc_method==1.or.coredens_method==1.or.coretau_method==1) then
268 : 1417 : call timab(551,1,tsec)
269 : : ! Compute Vxc in reciprocal space
270 [ + - + + ]: 1417 : if (coredens_method==1.and.n3xccc>0) then
271 [ - + + - : 2880 : ABI_MALLOC(v_dum,(nfft))
+ - + - -
+ - + ]
272 [ - + + - : 2880 : ABI_MALLOC(vxctotg,(2,nfft))
+ - + - -
+ - + ]
273 [ + + + + : 35056098 : v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
+ + ]
274 : 960 : 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 : 960 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
277 [ - + ]: 960 : 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 [ + + ]: 1417 : optv=0;if (vloc_method==1) optv=1
283 [ + - ]: 1417 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
284 : 1417 : optatm=0;optdyfr=0;opteltfr=0;optgr=0;optstr=1;optn2=1
285 [ + - ]: 1417 : 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 : 1417 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
292 : : end if
293 [ + + + - ]: 1417 : if (n3xccc==0.and.coredens_method==1) corstr=zero
294 [ - + ]: 1417 : ABI_FREE(vxctotg)
295 [ + + + - : 1417 : 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 : 1417 : call timab(551,2,tsec)
315 : : end if
316 : :
317 : : !Local ionic potential by method 2
318 [ + + ]: 6329 : if (vloc_method==2) then
319 : 4915 : option=3
320 [ - + + - : 14745 : ABI_MALLOC(dyfr_dum,(3,3,natom))
+ - + - -
+ - + ]
321 [ - + + - : 14745 : ABI_MALLOC(gr_dum,(3,natom))
+ - + - -
+ - + ]
322 [ - + + - : 14745 : 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 : 4915 : & qprtrb_dum,dum_rcut,rhog,rprimd,ucvol,vcutgeo,vlspl,vprtrb_dum,v_dum)
326 [ - + ]: 4915 : ABI_FREE(dyfr_dum)
327 [ - + ]: 4915 : ABI_FREE(gr_dum)
328 [ - + ]: 4915 : ABI_FREE(v_dum)
329 : : end if
330 : :
331 : : !Pseudo core electron density by method 2
332 [ + + ]: 6329 : if (coredens_method==2.or.coretau_method==2) then
333 [ + + ]: 4912 : if (n1xccc/=0) then
334 : 1817 : call timab(55,1,tsec)
335 : 1817 : option=3
336 [ - + + - : 5451 : ABI_MALLOC(dyfr_dum,(3,3,natom))
+ - + - -
+ - + ]
337 [ - + + - : 5451 : ABI_MALLOC(gr_dum,(3,natom))
+ - + - -
+ - + ]
338 [ - + + - : 5451 : ABI_MALLOC(v_dum,(nfft))
+ - + - -
+ - + ]
339 [ + - ]: 1817 : if (coredens_method==2) then
340 [ + - + - ]: 1817 : if (psps%usewvl==0.and.usepaw==0.and.icoulomb==0) then
341 [ + + ]: 1817 : 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 : 1787 : & 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 [ + + + - : 1817 : 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 [ - + ]: 1817 : ABI_FREE(dyfr_dum)
363 [ - + ]: 1817 : ABI_FREE(gr_dum)
364 [ - + ]: 1817 : ABI_FREE(v_dum)
365 : 1817 : 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 : 6329 : call strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd)
376 : :
377 : : !=======================================================================
378 : : !======================= Ewald contribution ============================
379 : : !=======================================================================
380 : :
381 : 6329 : call timab(38,1,tsec)
382 : 6329 : call ewald2(gmet,natom,ntypat,rmet,rprimd,ewestr,typat,ucvol,xred,zion)
383 : :
384 : : !=======================================================================
385 : : !================== VdW DFT-D contribution ============================
386 : : !=======================================================================
387 : :
388 [ + + ]: 6329 : 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 [ + + ]: 6319 : 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 : 6329 : 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 : 6329 : strsii=-eii/ucvol
453 : : !qvpotzero is non zero, only when usepotzero=1
454 : 6329 : 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 [ + + ]: 6329 : & berryopt==14 .or. berryopt==16 .or. berryopt==17)
461 [ + + ]: 6370 : 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 [ + + ]: 44303 : strten(:)=kinstr(:)+ewestr(:)+corstr(:)+strscondft(:)+strsxc(:)+harstr(:)+lpsstr(:)+nlstr(:)+mggastr(:)
540 : :
541 [ + + + - ]: 6329 : 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 [ + + ]: 6329 : if ( efield_flag ) then
549 [ + + ]: 175 : strten(:)=strten(:)+Maxstr(:)
550 [ + + + + ]: 79 : if ( calc_epaw3_stress ) strten(:) = strten(:) + epaws3red(:)
551 : : end if
552 [ + + + + ]: 6431 : if (vdw_xc>=5.and.vdw_xc<=7) strten(:)=strten(:)+vdwstr(:)
553 : :
554 : : !Additional stuff for electron-positron
555 : 6329 : ipositron=0
556 [ + - ]: 6329 : if (present(electronpositron)) then
557 [ + + ]: 6329 : 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 [ - + ]: 6329 : 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 [ + + + + ]: 6335 : if (ipositron>0) strten(:)=strten(:)+electronpositron%stress_ep(:)
580 : :
581 : : !Symmetrize resulting tensor if nsym>1
582 [ + + ]: 6329 : if (nsym>1) then
583 : 4854 : call stresssym(gprimd,nsym,strten,symrec)
584 : : end if
585 : :
586 : : !Set to zero very small values of stress
587 [ + + ]: 44303 : do mu=1,6
588 [ + + ]: 44303 : if (abs(strten(mu))<tol) strten(mu)=zero
589 : : end do
590 : :
591 : : !Include diagonal terms, save uncorrected stress for output
592 [ + + ]: 25316 : do mu=1,3
593 : 18987 : uncorr(mu)=strten(mu)+strsii
594 : 25316 : strten(mu)=uncorr(mu)
595 : : end do
596 : :
597 : : !Adding the extfpmd continous contribution to stress tensor
598 [ + + ]: 6329 : 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 [ + + + - ]: 6329 : 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 [ + - ]: 6085 : if (ipositron>=0) then
717 : 6329 : write(message, '(a,a)' )ch10,&
718 : 12658 : & ' Cartesian components of stress tensor (hartree/bohr^3)'
719 : 6329 : call wrtout(ab_out,message,'COLL')
720 : 6329 : call wrtout(std_out, message,'COLL')
721 : : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
722 : 6329 : & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4)
723 : 6329 : call wrtout(ab_out,message,'COLL')
724 : 6329 : call wrtout(std_out, message,'COLL')
725 : : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
726 : 6329 : & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5)
727 : 6329 : call wrtout(ab_out,message,'COLL')
728 : 6329 : call wrtout(std_out, message,'COLL')
729 : : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
730 : 6329 : & ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6)
731 : 6329 : call wrtout(ab_out,message,'COLL')
732 : 6329 : call wrtout(std_out, message,'COLL')
733 : 6329 : write(message, '(a)' ) ' '
734 : 6329 : call wrtout(ab_out,message,'COLL')
735 : 6329 : call wrtout(std_out, message,'COLL')
736 : : end if
737 : : ABI_NVTX_END_RANGE()
738 : 6329 : call timab(37,2,tsec)
739 : :
740 : 6329 : 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 : 6329 : 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 : 6329 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
792 : 6329 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
793 : : ! *************************************************************************
794 : :
795 : 6329 : call timab(568,1,tsec)
796 : :
797 : 6329 : harstr(:)=zero
798 : : !ehtest=0.0_dp (used for testing)
799 : :
800 : 6329 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
801 : :
802 [ - + ]: 6329 : 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 : 6329 : cutoff=gsqcut*tolfix
807 : 6329 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
808 : 6329 : me_fft=ngfft(11)
809 : 6329 : nproc_fft=ngfft(10)
810 : 6329 : id1=n1/2+2
811 : 6329 : id2=n2/2+2
812 : 6329 : id3=n3/2+2
813 : 6329 : ii=0
814 : :
815 : : ! Get the distrib associated with this fft_grid
816 : 6329 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
817 : :
818 [ + + ]: 148336 : do i3=1,n3
819 : 142007 : ig3=i3-(i3/id3)*n3-1
820 [ + + ]: 4011998 : do i2=1,n2
821 : 3863662 : ig2=i2-(i2/id2)*n2-1
822 [ + + ]: 4005669 : if (fftn2_distrib(i2)==me_fft) then
823 [ + + ]: 141193999 : do i1=1,n1
824 : 137582233 : ig1=i1-(i1/id1)*n1-1
825 : : ! ii=ii+1
826 : 137582233 : 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 [ + + ]: 137582233 : if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
831 : : ! Compute cartesian components of G
832 : 137576126 : gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+gprimd(1,3)*dble(ig3)
833 : 137576126 : gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+gprimd(2,3)*dble(ig3)
834 : 137576126 : gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+gprimd(3,3)*dble(ig3)
835 : : ! Compute |G|^2
836 : 137576126 : 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 [ + + ]: 141187892 : if (gsquar<=cutoff) then
840 : : ! take |rho(G)|^2 for complex rhog
841 [ + - ]: 55585523 : if (irho2==0) then
842 : 55585523 : 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 : 55585523 : harstr(1)=harstr(1)+(rhogsq/gsquar**2)*gcart(1)*gcart(1)
847 : 55585523 : harstr(2)=harstr(2)+(rhogsq/gsquar**2)*gcart(2)*gcart(2)
848 : 55585523 : harstr(3)=harstr(3)+(rhogsq/gsquar**2)*gcart(3)*gcart(3)
849 : 55585523 : harstr(4)=harstr(4)+(rhogsq/gsquar**2)*gcart(3)*gcart(2)
850 : 55585523 : harstr(5)=harstr(5)+(rhogsq/gsquar**2)*gcart(3)*gcart(1)
851 : 55585523 : 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 [ + + ]: 6329 : 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 : 6329 : harstr(1)=harstr(1)/pi-ehart/ucvol
879 : 6329 : harstr(2)=harstr(2)/pi-ehart/ucvol
880 : 6329 : harstr(3)=harstr(3)/pi-ehart/ucvol
881 : 6329 : harstr(4)=harstr(4)/pi
882 : 6329 : harstr(5)=harstr(5)/pi
883 : 6329 : harstr(6)=harstr(6)/pi
884 : :
885 : 6329 : call timab(568,2,tsec)
886 : :
887 : 6329 : end subroutine strhar
888 : : !!***
889 : :
890 : : end module m_stress
891 : : !!***
|
