TABLE OF CONTENTS
ABINIT/forstr [ Functions ]
NAME
forstr
FUNCTION
Drives the computation of forces and/or stress tensor
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MB, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx cg(2,mcg)=wavefunctions (may be read from disk instead of input) cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn> dtefield <type(efield_type)> = variables related to Berry phase dtset <type(dataset_type)>=all input variables in this dataset | berryopt = 4: electric field is on -> add the contribution of the | - \Omega E.P term to the total energy | /= 4: electric field is off | from Etot(npw) data (at fixed geometry), used for making | Pulay correction to stress tensor (hartree). Should be <=0. | ecut=cut-off energy for plane wave basis sphere (Ha) | ecutsm=smearing energy for plane wave kinetic energy (Ha) | effmass_free=effective mass for electrons (1. in common case) | efield = cartesian coordinates of the electric field in atomic units | ionmov=governs the movement of atoms (see help file) | densfor_pred=governs the mixed electronic-atomic part of the preconditioner | istwfk(nkpt)=input option parameter that describes the storage of wfs | kptns(3,nkpt)=reduced coordinates of k points in Brillouin zone | mband=maximum number of bands | mgfft=maximum size of 1D FFTs | mkmem=maximum number of k points in core memory | mpw = maximum number of plane waves | natom=number of atoms in cell | nband(nkpt*nsppol)=number of bands to be included in summation at each k point | 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 | nkpt=number of k points in Brillouin zone | nloalg(3)=governs the choice of the algorithm for non-local operator. | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | pawprtvol=control print volume and debugging output for PAW | prtvol=integer controlling volume of printed output | symafm(nsym)=(anti)ferromagnetic part of symmetry operations | tfkinfunc=1 if use of Thomas-Fermi kinetic functional | =2 if use of recursion method | typat(natom)=type integer for each atom in cell | wtk(nkpt)=weights associated with various k points | nsym=number of symmetries in space group energies <type(energies_type)>=all part of total energy. | e_localpsp(IN)=local psp energy (hartree) | e_hartree(IN)=Hartree part of total energy (hartree units) | e_corepsp(IN)=psp core-core energy | e_kinetic(IN)=kinetic energy part of total energy. eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) fock <type(fock_type)>= quantities to calculate Fock exact exchange grchempottn(3,natom)=d(E_chemical potential)/d(xred) (hartree) grewtn(3,natom)=d(Ewald)/d(xred) (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree) gsqcut=cutoff value on G**2 for (large) sphere inside FFT box. gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2) indsym(4,nsym,natom)=index showing transformation of atom labels under symmetry operations (computed in symatm) kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfftf= -PAW ONLY- maximum size of 1D FFTs for the fine grid (mgfftf=mgfft for norm-conserving potential runs) mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor n3xccc=dimension of the xccc3d array (0 or nfftf). nattyp(ntypat)=number of atoms of each type nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for norm-conserving potential runs) ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid (ngfftf=ngfft for norm-conserving potential runs) ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc nhat(nfftf,nspden*psps%usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.f for a description npwarr(nkpt)=number of planewaves in basis and on boundary for each k ntypat=number of types of atoms nvresid(nfftf,nspden)=array for the residual of the density/potential occ(mband*nkpt*nsppol)=occupancies of bands at various k points optfor=1 if computation of forces is required optres=0 if the potential residual has to be used for forces corrections =1 if the density residual has to be used for forces corrections paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid 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 structure factor phases ph1df(2,3*(2*mgfftf+1)*natom)=-PAW only- 1-dim structure factor phases for the fine grid psps <type(pseudopotential_type)>=variables related to pseudopotentials | mpsang= 1+maximum angular momentum rhog(2,nfftf)=Fourier transform of charge density (bohr^-3) rhor(nfftf,nspden)=array for electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) strsxc(6)=xc correction to stress stress_needed=1 if computation of stress tensor is required symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates ucvol=unit cell volume in bohr**3 usecprj=1 if cprj datastructure is stored in memory vhartr(nfftf)=array for holding Hartree potential vpsp(nfftf)=array for holding local psp vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
OUTPUT
==== if (optfor==1) ====
diffor=maximal absolute value of changes in the components of
force between the input and the output.
favg(3)=mean of the forces before correction for translational symmetry
fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
at input, previous value of forces,
at output, new value.
Note : unlike fred, this array has been corrected by enforcing
the translational symmetry, namely that the sum of force
on all atoms is zero.
forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
fred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
gresid(3,natom)=forces due to the residual of the density/potential
grhf(3,natom)=Hellman-Feynman derivatives of the total energy
grxc(9+3*natom)=d(Exc)/d(xred) if core charges are used
maxfor=maximal absolute value of the output array force.
synlgr(3,natom)=symmetrized gradients of energy due to nonlocal contributions
==== if (stress_needed==1) ====
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).
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr) ===== if psps%usepaw==1 pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data (gradients of rhoij for each atom with respect to atomic positions are computed here) wvl <type(wvl_data)>=all wavelets data.
NOTES
Be careful to the meaning of nfft (size of FFT grids):
- In case of norm-conserving calculations the FFT grid is the usual FFT grid.
- In case of PAW calculations:
Two FFT grids are used; one with nfft points (coarse grid) for
the computation of wave functions ; one with nfftf points
(fine grid) for the computation of total density.
PARENTS
afterscfloop,setup_positron
CHILDREN
ctocprj,forces,forstrnps,initylmg,metric,nres2vres,pawcprj_alloc pawcprj_free,pawcprj_getdim,pawgrnl,stress,timab,wvl_nl_gradient xchybrid_ncpp_cc,xred2xcart
SOURCE
168 #if defined HAVE_CONFIG_H 169 #include "config.h" 170 #endif 171 172 #include "abi_common.h" 173 174 subroutine forstr(atindx1,cg,cprj,diffor,dtefield,dtset,eigen,electronpositron,energies,favg,fcart,fock,& 175 & forold,fred,grchempottn,gresid,grewtn,grhf,grvdw,grxc,gsqcut,indsym,& 176 & kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,& 177 & nfftf,ngfftf,ngrvdw,nhat,nkxc,npwarr,& 178 & ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,& 179 & pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,psps,rhog,rhor,rprimd,stress_needed,& 180 & strsxc,strten,symrec,synlgr,ucvol,usecprj,vhartr,vpsp,& 181 & vxc,wvl,xccc3d,xred,ylm,ylmgr,qvpotzero) 182 183 184 use defs_basis 185 use defs_datatypes 186 use defs_abitypes 187 use defs_wvltypes 188 use m_profiling_abi 189 use m_efield 190 use m_errors 191 192 use m_geometry, only : xred2xcart, metric 193 use m_electronpositron, only : electronpositron_type 194 use m_energies, only : energies_type 195 use m_pawang, only : pawang_type 196 use m_pawrad, only : pawrad_type 197 use m_pawtab, only : pawtab_type 198 use m_paw_ij, only : paw_ij_type 199 use m_pawfgrtab, only : pawfgrtab_type 200 use m_pawrhoij, only : pawrhoij_type 201 use m_pawfgr, only : pawfgr_type 202 use m_pawcprj, only : pawcprj_type,pawcprj_free,pawcprj_getdim,pawcprj_alloc 203 use m_fock, only : fock_type 204 use libxc_functionals, only : libxc_functionals_is_hybrid 205 206 !This section has been created automatically by the script Abilint (TD). 207 !Do not modify the following lines by hand. 208 #undef ABI_FUNC 209 #define ABI_FUNC 'forstr' 210 use interfaces_18_timing 211 use interfaces_56_recipspace 212 use interfaces_56_xc 213 use interfaces_62_wvl_wfs 214 use interfaces_65_paw 215 use interfaces_66_nonlocal 216 use interfaces_67_common, except_this_one => forstr 217 !End of the abilint section 218 219 implicit none 220 221 !Arguments ------------------------------------ 222 !scalars 223 integer,intent(in) :: mcg,mcprj,mgfftf,my_natom,n3xccc,nfftf,ngrvdw,nkxc,ntypat,optfor,optres 224 integer,intent(in) :: stress_needed,usecprj 225 real(dp),intent(in) :: gsqcut,qvpotzero,ucvol 226 real(dp),intent(inout) :: diffor,maxfor 227 type(electronpositron_type),pointer :: electronpositron 228 type(MPI_type),intent(inout) :: mpi_enreg 229 type(efield_type),intent(in) :: dtefield 230 type(dataset_type),intent(in) :: dtset 231 type(energies_type),intent(in) :: energies 232 type(pawang_type),intent(in) :: pawang 233 type(pawfgr_type),intent(in) :: pawfgr 234 type(pseudopotential_type),intent(in) :: psps 235 type(wvl_data),intent(inout) :: wvl 236 type(fock_type),pointer, intent(inout) :: fock 237 !arrays 238 integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom) 239 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfftf(18) 240 integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym) 241 real(dp),intent(in) :: cg(2,mcg) 242 real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 243 real(dp),intent(in) :: grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(dtset%nfft,nkxc) 244 real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 245 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 246 real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom) 247 real(dp),intent(in) :: rhog(2,nfftf),rprimd(3,3),strsxc(6),vhartr(nfftf) 248 real(dp),intent(in) :: vpsp(nfftf),vxc(nfftf,dtset%nspden) 249 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 250 real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 251 real(dp),intent(inout) :: forold(3,dtset%natom) 252 real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw),rhor(nfftf,dtset%nspden) 253 real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom) 254 real(dp),intent(inout),target :: nvresid(nfftf,dtset%nspden) 255 real(dp),intent(out) :: favg(3) 256 real(dp),intent(inout) :: fcart(3,dtset%natom),fred(3,dtset%natom) 257 real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom) 258 real(dp),intent(inout) :: grxc(3,dtset%natom),strten(6),synlgr(3,dtset%natom) 259 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj) 260 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 261 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 262 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 263 type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw) 264 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 265 266 !Local variables------------------------------- 267 !scalars 268 269 integer :: comm_grid,ifft,ispden,ncpgr,occopt_,optgr,optgr2,option,optnc,optstr,optstr2,iorder_cprj,ctocprj_choice 270 integer :: idir,iatom,unpaw,mcgbz 271 integer,allocatable :: dimcprj(:) 272 real(dp) ::dum,dum1,ucvol_ 273 logical :: apply_residual 274 !arrays 275 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 276 real(dp) :: kinstr(6),nlstr(6),tsec(2),strdum(6),gmet(3,3),gprimd(3,3),rmet(3,3) 277 real(dp) :: dummy(0) 278 real(dp),allocatable :: grnl(:),vlocal(:,:),vxc_hf(:,:),xcart(:,:),ylmbz(:,:),ylmgrbz(:,:,:) 279 real(dp), ABI_CONTIGUOUS pointer :: resid(:,:) 280 281 282 ! ************************************************************************* 283 284 call timab(910,1,tsec) 285 call timab(911,1,tsec) 286 287 !Do nothing if nothing is required 288 if (optfor==0.and.stress_needed==0) return 289 290 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW) 291 if (dtset%usewvl==0) then 292 if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then 293 MSG_BUG(' wrong values for nfft, nfftf !') 294 end if 295 if ((psps%usepaw==1.and.pawfgr%mgfft/=mgfftf).or.(psps%usepaw==0.and.dtset%mgfft/=mgfftf)) then 296 MSG_BUG('wrong values for mgfft, mgfftf!') 297 end if 298 end if 299 300 !========================================================================== 301 !Here compute terms common to forces and stresses 302 !========================================================================== 303 304 !output only if (optfor==1) but we have to allocate it 305 ABI_ALLOCATE(grnl,(3*dtset%natom*optfor)) 306 grnl(:)=zero 307 308 !Compute nonlocal pseudopotential parts of forces and stress tensor 309 !-involves summations over wavefunctions at all k points 310 if (dtset%tfkinfunc>0.and.stress_needed==1) then 311 kinstr(1:3)=-two/three*energies%e_kinetic/ucvol ; kinstr(4:6)=zero 312 nlstr(1:6)=zero 313 else if (dtset%usewvl==0) then 314 occopt_=0 ! This means that occ are now fixed 315 if(dtset%usefock==1 .and. associated(fock)) then 316 ! if((dtset%optstress/=0).and.(psps%usepaw==1)) then 317 if((psps%usepaw==1).and.((dtset%optstress/=0).or.(dtset%optforces==2))) then 318 if(dtset%optstress==0) then 319 ctocprj_choice=2 320 ncpgr=3 321 end if 322 if(dtset%optstress/=0) then 323 ctocprj_choice=20*optfor+3*dtset%optstress 324 ncpgr=6*dtset%optstress+3*optfor 325 end if 326 if (allocated(fock%fock_BZ%cwaveocc_prj)) then 327 call pawcprj_free(fock%fock_BZ%cwaveocc_prj) 328 ABI_DATATYPE_DEALLOCATE(fock%fock_BZ%cwaveocc_prj) 329 ABI_DATATYPE_ALLOCATE(fock%fock_BZ%cwaveocc_prj,(dtset%natom,fock%fock_BZ%mcprj)) 330 ABI_ALLOCATE(dimcprj,(dtset%natom)) 331 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 332 call pawcprj_alloc(fock%fock_BZ%cwaveocc_prj,ncpgr,dimcprj) 333 ABI_DEALLOCATE(dimcprj) 334 end if 335 iatom=-1;idir=0;iorder_cprj=0;unpaw=26 336 call metric(gmet,gprimd,-1,rmet,rprimd,dum) 337 if (fock%fock_BZ%mkpt/=dtset%mkmem.or.(fock%fock_BZ%mpi_enreg%paral_hf ==1)) then 338 ABI_ALLOCATE(ylmbz,(dtset%mpw*fock%fock_BZ%mkpt,psps%mpsang*psps%mpsang*psps%useylm)) 339 ABI_ALLOCATE(ylmgrbz,(dtset%mpw*fock%fock_BZ%mkpt,3,psps%mpsang*psps%mpsang*psps%useylm)) 340 option=1; mcgbz=dtset%mpw*fock%fock_BZ%mkptband*fock%fock_common%my_nsppol 341 call initylmg(gprimd,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,& 342 & psps%mpsang,dtset%mpw,fock%fock_BZ%nbandocc_bz,fock%fock_BZ%mkpt,& 343 & fock%fock_BZ%npwarr,dtset%nsppol,option,rprimd,ylmbz,ylmgrbz) 344 call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,& 345 & iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcgbz,& 346 & fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,psps%mpsang,& 347 & dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,& 348 & dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,& 349 & dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,& 350 & xred,ylmbz,ylmgrbz) 351 ABI_DEALLOCATE(ylmbz) 352 ABI_DEALLOCATE(ylmgrbz) 353 else 354 call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,& 355 & iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcg,& 356 & fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,mpi_enreg,psps%mpsang,& 357 & dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,& 358 & dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,& 359 & dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,& 360 & xred,ylm,ylmgr) 361 end if 362 end if 363 end if 364 call forstrnps(cg,cprj,dtset%ecut,dtset%ecutsm,dtset%effmass_free,eigen,electronpositron,fock,grnl,& 365 & dtset%istwfk,kg,kinstr,nlstr,dtset%kptns,dtset%mband,mcg,mcprj,dtset%mgfft,dtset%mkmem,& 366 & mpi_enreg,psps%mpsang,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,dtset%ngfft,& 367 & dtset%nkpt,dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,& 368 & dtset%nucdipmom,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,stress_needed,symrec,dtset%typat,& 369 & usecprj,dtset%usefock,dtset%use_gpu_cuda,dtset%wtk,xred,ylm,ylmgr) 370 else if (optfor>0) then !WVL 371 ABI_ALLOCATE(xcart,(3, dtset%natom)) 372 call xred2xcart(dtset%natom, rprimd, xcart, xred) 373 call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart) 374 ABI_DEALLOCATE(xcart) 375 end if 376 377 call timab(911,2,tsec) 378 call timab(912,1,tsec) 379 380 !PAW: add gradients due to Dij derivatives to non-local term 381 if (psps%usepaw==1) then 382 ABI_ALLOCATE(vlocal,(nfftf,dtset%nspden)) 383 384 !$OMP PARALLEL DO COLLAPSE(2) 385 do ispden=1,min(dtset%nspden,2) 386 do ifft=1,nfftf 387 vlocal(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft) 388 end do 389 end do 390 if (dtset%nspden==4) then 391 !$OMP PARALLEL DO COLLAPSE(2) 392 do ispden=3,4 393 do ifft=1,nfftf 394 vlocal(ifft,ispden)=vxc(ifft,ispden) 395 end do 396 end do 397 end if 398 ucvol_=ucvol 399 #if defined HAVE_BIGDFT 400 if (dtset%usewvl==1) ucvol_=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp) 401 #endif 402 optgr=optfor;optgr2=0;optstr=stress_needed;optstr2=0 403 comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl 404 call pawgrnl(atindx1,dtset%nspden,dummy,1,dummy,grnl,gsqcut,mgfftf,my_natom,dtset%natom,& 405 & nattyp,nfftf,ngfftf,nhat,nlstr,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,& 406 & pawang,pawfgrtab,pawrhoij,pawtab,ph1df,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred,& 407 & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=comm_grid) 408 ABI_DEALLOCATE(vlocal) 409 410 end if 411 call timab(912,2,tsec) 412 call timab(913,1,tsec) 413 414 !========================================================================== 415 !Here compute forces (if required) 416 !========================================================================== 417 if (optfor==1) then 418 apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. & 419 & abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5) 420 421 ! If residual is a density residual (and forces from residual asked), 422 ! has to convert it into a potential residual before calling forces routine 423 if (apply_residual) then 424 ABI_ALLOCATE(resid,(nfftf,dtset%nspden)) 425 option=0; if (dtset%densfor_pred<0) option=1 426 optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2 427 call nres2vres(dtset,gsqcut,psps%usepaw,kxc,mpi_enreg,my_natom,nfftf,ngfftf,nhat,& 428 & nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,& 429 & rhor,rprimd,psps%usepaw,resid,xccc3d,xred,vxc) 430 else 431 resid => nvresid 432 end if 433 434 call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,fred,grchempottn,gresid,grewtn,& 435 & grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfftf,& 436 & mpi_enreg,psps%n1xccc,n3xccc,nattyp,& 437 & nfftf,ngfftf,ngrvdw,ntypat,pawrad,pawtab,ph1df,psps,rhog,& 438 & rhor,rprimd,symrec,synlgr,dtset%usefock,resid,vxc,wvl%descr,wvl%den,xred,& 439 & electronpositron=electronpositron) 440 441 if (apply_residual) then 442 ABI_DEALLOCATE(resid) 443 end if 444 end if 445 446 call timab(913,2,tsec) 447 call timab(914,1,tsec) 448 449 !========================================================================== 450 !Here compute stress tensor (if required) 451 !========================================================================== 452 453 if (stress_needed==1.and.dtset%usewvl==0) then 454 ! if (dtset%usefock==1 .and. associated(fock).and.fock%fock_common%optstr.and.psps%usepaw==0) then 455 if (dtset%usefock==1 .and. associated(fock).and.fock%fock_common%optstr) then 456 fock%fock_common%stress(1:3)=fock%fock_common%stress(1:3)-energies%e_fock/ucvol 457 if (n3xccc>0.and.psps%usepaw==0 .and. & 458 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) then 459 ABI_ALLOCATE(vxc_hf,(nfftf,dtset%nspden)) 460 !compute Vxc^GGA(rho_val) 461 call xchybrid_ncpp_cc(dtset,dum,mpi_enreg,nfftf,ngfftf,n3xccc,rhor,rprimd,strdum,dum1,xccc3d,vxc=vxc_hf,optstr=1) 462 end if 463 end if 464 call stress(atindx1,dtset%berryopt,dtefield,energies%e_localpsp,dtset%efield,& 465 & energies%e_hartree,energies%e_corepsp,fock,gsqcut,dtset%ixc,kinstr,mgfftf,& 466 & mpi_enreg,psps%mqgrid_vl,psps%n1xccc,n3xccc,dtset%natom,nattyp,& 467 & nfftf,ngfftf,nlstr,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,psps,pawrad,pawtab,ph1df,& 468 & dtset%prtvol,psps%qgrid_vl,dtset%red_efieldbar,rhog,rprimd,strten,strsxc,symrec,& 469 & dtset%typat,dtset%usefock,psps%usepaw,& 470 & dtset%vdw_tol,dtset%vdw_tol_3bt,dtset%vdw_xc,psps%vlspl,vxc,vxc_hf,psps%xccc1d,xccc3d,psps%xcccrc,xred,& 471 & psps%ziontypat,psps%znucltypat,qvpotzero,& 472 & electronpositron=electronpositron) 473 end if 474 475 !Memory deallocation 476 ABI_DEALLOCATE(grnl) 477 if (allocated(vxc_hf)) then 478 ABI_DEALLOCATE(vxc_hf) 479 end if 480 481 482 call timab(914,2,tsec) 483 call timab(910,2,tsec) 484 485 end subroutine forstr