TABLE OF CONTENTS
ABINIT/forces [ Functions ]
NAME
forces
FUNCTION
Assemble gradients of various total energy terms with respect to reduced coordinates, including possible symmetrization, in order to produce forces. fcart(i,iat) = d(Etot)/(d(r(i,iat)))
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, FJ, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx dtefield <type(efield_type)> = variables related to Berry phase dtset <type(dataset_type)>=all input variables in this dataset berryopt = 4/14: electric field is on -> add the contribution of the -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy = 6/16, or 7/17: electric displacement field is on -> add the contribution of the Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy | efield = cartesian coordinates of the electric field in atomic units | dfield = cartesian coordinates of the electric displacement field in atomic units | iatfix(3,natom)=1 for frozen atom along specified direction, 0 for unfrozen | ionmov=governs the movement of atoms (see help file) | densfor_pred=governs the mixed electronic-atomic part of the preconditioner | natom=number of atoms in cell | nconeq=number of atomic constraint equations | nspden=number of spin-density components | nsym=number of symmetries in space group | prtvol=integer controlling volume of printed output | typat(natom)=type integer for each atom in cell | wtatcon(3,natom,nconeq)=weights for atomic constraints 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) grnl(3*natom)=gradients of Etot due to nonlocal contributions 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)=indirect indexing array for atom labels mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization n1xccc=dimension of xccc1d ; 0 if no XC core correction is used n3xccc=dimension of the xccc3d array (0 or nfft). nattyp(ntypat)=number of atoms of each type nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc ntypat=number of types of atoms pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) array psps <type(pseudopotential_type)>=variables related to pseudopotentials rhog(2,nfft)=Fourier transform of charge density (bohr^-3) rhor(nfft,nspden)=array for electron density in electrons/bohr**3 rprimd(3,3)=dimensional primitive translations in real space (bohr) symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates usefock=1 if fock operator is used; 0 otherwise. vresid(nfft,nspden)=potential residual (if non-collinear magn., only trace of it) vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space xred(3,natom)=reduced dimensionless atomic coordinates xred_old(3,natom)=previous reduced dimensionless atomic coordinates
OUTPUT
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
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 d(enl)/d(xred)
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument) fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr) Note : unlike fred, this array has been corrected by enforcing the translational symmetry, namely that the sum of force on all atoms is zero.
NOTES
* Symmetrization of gradients with respect to reduced
coordinates xred is conducted according to the expression
[d(e)/d(t(n,a))]_symmetrized = (1/Nsym) Sum(S) symrec(n,m,S)*
[d(e)/d(t(m,b))]_unsymmetrized
where t(m,b)= (symrel^-1)(m,n)*(t(n,a)-tnons(n)) and tnons
is a possible nonsymmorphic translation. The label "b" here
refers to the atom which gets rotated into "a" under symmetry "S".
symrel is the symmetry matrix in real space, which is the inverse
transpose of symrec. symrec is the symmetry matrix in reciprocal
space. sym_cartesian = R * symrel * R^-1 = G * symrec * G^-1
where the columns of R and G are the dimensional primitive translations
in real and reciprocal space respectively.
* Note the use of "symrec" in the symmetrization expression above.
PARENTS
etotfor,forstr
CHILDREN
atm2fft,constrf,dgemv,fourdp,fred2fcart,fresid,fresidrsp,metric,mkcore mkcore_alt,mkcore_wvl,mklocl,sygrad,timab,xchybrid_ncpp_cc,zerosym
SOURCE
114 #if defined HAVE_CONFIG_H 115 #include "config.h" 116 #endif 117 118 #include "abi_common.h" 119 120 subroutine forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,& 121 & forold,fred,grchempottn,gresid,grewtn,& 122 & grhf,grnl,grvdw,grxc,gsqcut,indsym,& 123 & maxfor,mgfft,mpi_enreg,n1xccc,n3xccc,& 124 & nattyp,nfft,ngfft,ngrvdw,ntypat,& 125 & pawrad,pawtab,ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,usefock,& 126 & vresid,vxc,wvl,wvl_den,xred,& 127 & electronpositron) ! optional argument 128 129 use defs_basis 130 use defs_datatypes 131 use defs_abitypes 132 use defs_wvltypes 133 use m_profiling_abi 134 use m_efield 135 use m_errors 136 137 use m_geometry, only : fred2fcart, metric 138 use m_fock, only : fock_type 139 use m_pawrad, only : pawrad_type 140 use m_pawtab, only : pawtab_type 141 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 142 use libxc_functionals, only : libxc_functionals_is_hybrid 143 use m_fft, only : zerosym 144 145 !This section has been created automatically by the script Abilint (TD). 146 !Do not modify the following lines by hand. 147 #undef ABI_FUNC 148 #define ABI_FUNC 'forces' 149 use interfaces_18_timing 150 use interfaces_53_ffts 151 use interfaces_56_xc 152 use interfaces_64_psp 153 use interfaces_67_common, except_this_one => forces 154 !End of the abilint section 155 156 implicit none 157 158 !Arguments ------------------------------------ 159 !scalars 160 integer,intent(in) :: mgfft,n1xccc,n3xccc,nfft,ngrvdw,ntypat,usefock 161 real(dp),intent(in) :: gsqcut 162 real(dp),intent(out) :: diffor,maxfor 163 type(MPI_type),intent(in) :: mpi_enreg 164 type(efield_type),intent(in) :: dtefield 165 type(dataset_type),intent(in) :: dtset 166 type(electronpositron_type),pointer,optional :: electronpositron 167 type(pseudopotential_type),intent(in) :: psps 168 type(wvl_internal_type), intent(in) :: wvl 169 type(wvl_denspot_type), intent(inout) :: wvl_den 170 type(fock_type),pointer, intent(inout) :: fock 171 !arrays 172 integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom) 173 integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym) 174 real(dp),intent(in) :: grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw),grnl(3*dtset%natom) 175 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom) 176 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rprimd(3,3) 177 real(dp),intent(in) :: vxc(nfft,dtset%nspden) 178 real(dp),intent(inout) :: fcart(3,dtset%natom),forold(3,dtset%natom) 179 real(dp),intent(inout) :: vresid(nfft,dtset%nspden),xred(3,dtset%natom) 180 real(dp),intent(out) :: favg(3),fred(3,dtset%natom),gresid(3,dtset%natom) 181 real(dp),intent(out) :: grhf(3,dtset%natom) 182 real(dp),intent(inout) :: grxc(3,dtset%natom) 183 real(dp),intent(out) :: synlgr(3,dtset%natom) 184 type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw) 185 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 186 187 !Local variables------------------------------- 188 !scalars 189 integer :: coredens_method,fdir,iatom,idir,indx,ipositron,itypat,mu 190 integer :: optatm,optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method 191 real(dp) :: eei_dum1,eei_dum2,ucvol,ucvol_local,vol_element 192 logical :: calc_epaw3_forces, efield_flag 193 logical :: is_hybrid_ncpp 194 !arrays 195 integer :: qprtrb_dum(3) 196 real(dp) :: dummy6(6),ep3(3),fioncart(3),gmet(3,3),gprimd(3,3) 197 real(dp) :: rmet(3,3),strn_dummy6(6),strv_dummy6(6),tsec(2),vprtrb_dum(2) 198 real(dp),allocatable :: atmrho_dum(:),atmvloc_dum(:),dyfrlo_dum(:,:,:) 199 real(dp),allocatable :: dyfrn_dum(:,:,:),dyfrv_dum(:,:,:) 200 real(dp),allocatable :: dyfrx2_dum(:,:,:),eltfrn_dum(:,:),gauss_dum(:,:) 201 real(dp),allocatable :: epawf3red(:,:),fin(:,:),fionred(:,:),grl(:,:),grnl_tmp(:,:),grtn(:,:) 202 real(dp),allocatable :: grtn_indx(:,:),v_dum(:),vxctotg(:,:) 203 real(dp),allocatable :: xccc3d_dum(:) 204 205 ! ************************************************************************* 206 207 call timab(69,1,tsec) 208 209 !Save input value of forces 210 ABI_ALLOCATE(fin,(3,dtset%natom)) 211 fin(:,:)=fcart(:,:) 212 213 !Compute different geometric tensor, as well as ucvol, from rprimd 214 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 215 216 !Check if we're in hybrid norm conserving pseudopotential 217 is_hybrid_ncpp=(psps%usepaw==0 .and. & 218 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) 219 !======================================================================= 220 !========= Local pseudopotential and core charge contributions ========= 221 !======================================================================= 222 223 ABI_ALLOCATE(grl,(3,dtset%natom)) 224 225 !Determine by which method the local ionic potential and/or the pseudo core 226 ! charge density contributions have to be computed 227 !Local ionic potential: 228 ! Method 1: PAW 229 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 230 vloc_method=1;if (psps%usepaw==0) vloc_method=2 231 if (dtset%icoulomb>0) vloc_method=2 232 if (psps%usewvl==1) vloc_method=2 233 !Pseudo core charge density: 234 ! Method 1: PAW, nc_xccc_gspace 235 ! Method 2: Norm-conserving PP, wavelets 236 coredens_method=1;if (psps%usepaw==0) coredens_method=2 237 if (psps%nc_xccc_gspace==1) coredens_method=1 238 if (psps%nc_xccc_gspace==0) coredens_method=2 239 if (psps%usewvl==1) coredens_method=2 240 241 !Local ionic potential and/or pseudo core charge by method 1 242 if (vloc_method==1.or.coredens_method==1) then 243 call timab(550,1,tsec) 244 optv=0;if (vloc_method==1) optv=1 245 optn=0;if (coredens_method==1) optn=n3xccc/nfft 246 optatm=0;optdyfr=0;optgr=1;optstr=0;optn2=1;opteltfr=0 247 if (psps%nc_xccc_gspace==1.and.psps%usepaw==0.and.is_hybrid_ncpp) then 248 MSG_BUG(' Not yet implemented !') 249 end if 250 if (coredens_method==1.and.n3xccc>0) then 251 ABI_ALLOCATE(v_dum,(nfft)) 252 ABI_ALLOCATE(vxctotg,(2,nfft)) 253 v_dum(:)=vxc(:,1);if (dtset%nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2)) 254 call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 255 call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),& 256 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 257 ABI_DEALLOCATE(v_dum) 258 else 259 ABI_ALLOCATE(vxctotg,(0,0)) 260 end if 261 ! Allocate (unused) dummy variables, otherwise some compilers complain 262 ABI_ALLOCATE(gauss_dum,(0,0)) 263 ABI_ALLOCATE(atmrho_dum,(0)) 264 ABI_ALLOCATE(atmvloc_dum,(0)) 265 ABI_ALLOCATE(dyfrn_dum,(0,0,0)) 266 ABI_ALLOCATE(dyfrv_dum,(0,0,0)) 267 ABI_ALLOCATE(eltfrn_dum,(0,0)) 268 call atm2fft(atindx1,atmrho_dum,atmvloc_dum,dyfrn_dum,dyfrv_dum,& 269 & eltfrn_dum,gauss_dum,gmet,gprimd,& 270 & grxc,grl,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,& 271 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb_dum,& 272 & rhog,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,psps%vlspl,& 273 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 274 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 275 ABI_DEALLOCATE(gauss_dum) 276 ABI_DEALLOCATE(atmrho_dum) 277 ABI_DEALLOCATE(atmvloc_dum) 278 ABI_DEALLOCATE(dyfrn_dum) 279 ABI_DEALLOCATE(dyfrv_dum) 280 ABI_DEALLOCATE(eltfrn_dum) 281 ABI_DEALLOCATE(vxctotg) 282 if (n3xccc==0.and.coredens_method==1) grxc=zero 283 call timab(550,2,tsec) 284 end if 285 286 !Local ionic potential by method 2 287 if (vloc_method==2) then 288 option=2 289 ABI_ALLOCATE(dyfrlo_dum,(3,3,dtset%natom)) 290 ABI_ALLOCATE(grtn_indx,(3,dtset%natom)) 291 ABI_ALLOCATE(v_dum,(nfft)) 292 call mklocl(dtset,dyfrlo_dum,eei_dum1,gmet,gprimd,grtn_indx,gsqcut,dummy6,mgfft,& 293 & mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,& 294 & qprtrb_dum,rhog,rhor,rprimd,ucvol,vprtrb_dum,v_dum,wvl,wvl_den,xred) 295 do iatom=1,dtset%natom 296 ! Has to use the indexing array atindx1 297 grl(1:3,atindx1(iatom))=grtn_indx(1:3,iatom) 298 end do 299 ABI_DEALLOCATE(dyfrlo_dum) 300 ABI_DEALLOCATE(grtn_indx) 301 ABI_DEALLOCATE(v_dum) 302 ! If gradients are computed in real space, we need to symmetrize the system before summing. 303 ! Rshaltaf: I changed the following line to include surfaces BC 304 if (dtset%icoulomb == 1 .or. dtset%icoulomb == 2) then 305 ABI_ALLOCATE(grnl_tmp,(3,dtset%natom)) 306 call sygrad(grnl_tmp,dtset%natom,grl,dtset%nsym,symrec,indsym) 307 grl(:, :) = grnl_tmp(:, :) 308 ABI_DEALLOCATE(grnl_tmp) 309 end if 310 end if 311 312 !Pseudo core electron density by method 2 313 if (coredens_method==2) then 314 if (n1xccc/=0) then 315 call timab(53,1,tsec) 316 option=2 317 ABI_ALLOCATE(dyfrx2_dum,(3,3,dtset%natom)) 318 ABI_ALLOCATE(xccc3d_dum,(n3xccc)) 319 if (is_hybrid_ncpp) then 320 call xchybrid_ncpp_cc(dtset,eei_dum1,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,& 321 & dummy6,eei_dum2,xccc3d_dum,grxc=grxc,xcccrc=psps%xcccrc,xccc1d=psps%xccc1d,xred=xred,n1xccc=n1xccc) 322 else 323 if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then 324 call mkcore(dummy6,dyfrx2_dum,grxc,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 325 & ngfft(1),n1xccc, ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,vxc,& 326 & psps%xcccrc,psps%xccc1d,xccc3d_dum,xred) 327 else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then 328 call mkcore_alt(atindx1,dummy6,dyfrx2_dum,grxc,dtset%icoulomb,mpi_enreg,dtset%natom,nfft,& 329 & dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 330 & ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d_dum,xred,pawrad,pawtab,psps%usepaw) 331 else if (psps%usewvl==1.and.psps%usepaw==1) then 332 ucvol_local=ucvol 333 #if defined HAVE_BIGDFT 334 ! ucvol_local=product(wvl_den%denspot%dpbox%hgrids)*real(product(wvl_den%denspot%dpbox%ndims),dp) 335 ! call mkcore_wvl_old(atindx1,dummy6,dyfrx2_dum,wvl%atoms%astruct%geocode,grxc,wvl%h,dtset%natom,& 336 ! & nattyp,nfft,wvl_den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),dtset%nspden,ntypat,& 337 ! & wvl%Glr%d%n1,wvl%Glr%d%n1i,wvl%Glr%d%n2,wvl%Glr%d%n2i,wvl%Glr%d%n3,wvl_den%denspot%dpbox%n3pi,& 338 ! & n3xccc,option,pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d_dum,xred,& 339 ! & mpi_comm_wvl=mpi_enreg%comm_wvl) 340 call mkcore_wvl(atindx1,dummy6,grxc,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,& 341 & n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d_dum,& 342 & psps%xcccrc,xred,wvl_den,wvl,mpi_comm_wvl=mpi_enreg%comm_wvl) 343 #endif 344 end if 345 end if 346 ABI_DEALLOCATE(xccc3d_dum) 347 ABI_DEALLOCATE(dyfrx2_dum) 348 call timab(53,2,tsec) 349 else 350 grxc(:,:)=zero 351 end if 352 end if 353 354 !======================================================================= 355 !===================== Nonlocal contributions ========================== 356 !======================================================================= 357 358 !Only has to apply symmetries 359 ABI_ALLOCATE(grnl_tmp,(3,dtset%natom)) 360 do iatom=1,dtset%natom 361 indx=3*(iatom-1);grnl_tmp(1:3,atindx1(iatom))=grnl(indx+1:indx+3) 362 end do 363 if (dtset%usewvl == 0) then 364 call sygrad(synlgr,dtset%natom,grnl_tmp,dtset%nsym,symrec,indsym) 365 else 366 synlgr = grnl_tmp 367 end if 368 ABI_DEALLOCATE(grnl_tmp) 369 370 !======================================================================= 371 !============ Density/potential residual contributions ================= 372 !======================================================================= 373 374 if (dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1.and.abs(dtset%densfor_pred)<=3) then 375 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,ntypat,1,& 376 & pawtab,rhor,rprimd,ucvol,vresid,xred,xred,psps%znuclpsp) 377 else if (dtset%usewvl==0.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) then 378 call fresidrsp(atindx1,dtset,gmet,gprimd,gresid,gsqcut,mgfft,& 379 & mpi_enreg,psps%mqgrid_vl,nattyp,nfft,ngfft,ntypat,psps,pawtab,ph1d,& 380 & psps%qgrid_vl,ucvol,psps%usepaw,vresid,psps%zionpsp,psps%znuclpsp) 381 else 382 gresid(:,:)=zero 383 end if 384 385 !======================================================================= 386 !======================= Other contributions =========================== 387 !======================================================================= 388 389 !Ewald energy contribution to forces as already been computed in "ewald" 390 391 !Potential residual contribution to forces as already been computed (forstr) 392 393 !Add Berry phase contributions (berryopt == 4/6/7/14/16/17) 394 !(compute the electric field force on the ion cores) 395 efield_flag = (dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. & 396 & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17) 397 calc_epaw3_forces = (efield_flag .and. dtset%optforces /= 0 .and. psps%usepaw == 1) 398 if ( efield_flag ) then 399 ABI_ALLOCATE(fionred,(3,dtset%natom)) 400 fionred(:,:)=zero 401 do iatom=1,dtset%natom 402 itypat=dtset%typat(iatom) 403 ! force on ion due to electric field, cartesian representation 404 fioncart(:)=psps%ziontypat(itypat)*dtset%efield(:) 405 ! form fionred = rprimd^T * fioncart, note that forces transform 406 ! oppositely to coordinates, because they are derivative with respect to 407 ! coordinates 408 call dgemv('T',3,3,one,rprimd,3,fioncart,1,zero,fionred(1:3,iatom),1) 409 ! do mu=1,3 410 ! fionred(mu,iatom)=rprimd(1,mu)*fioncart(1) & 411 !& +rprimd(2,mu)*fioncart(2) & 412 !& +rprimd(3,mu)*fioncart(3) 413 ! end do 414 end do 415 end if 416 417 !(compute additional F3-type force due to projectors for electric field with PAW) 418 if ( efield_flag .and. calc_epaw3_forces ) then 419 ABI_ALLOCATE(epawf3red,(3,dtset%natom)) 420 ! dtefield%epawf3(iatom,idir,fdir) contains 421 epawf3red(:,:)=zero 422 do iatom=1,dtset%natom 423 do fdir = 1, 3 424 do idir = 1, 3 425 ! vol_element is volume/pt for integration of epawf3. volume is BZ volume 426 ! so 1/ucvol, and number of kpts is nstr(idir)*nkstr(idir) 427 vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir)) 428 ep3(idir) = vol_element*dtefield%epawf3(iatom,idir,fdir) 429 end do 430 epawf3red(fdir,iatom) = -ucvol*dot_product(dtset%red_efieldbar(1:3),ep3(1:3)) 431 end do 432 end do ! end loop over iatom 433 end if 434 435 !This was incorrect coding. Bug found by Jiawang Hong 436 !if (dtset%berryopt==4) then 437 !allocate(fionred(3,dtset%natom));fionred(:,:)=zero 438 !iatom = 0 439 !do itypat=1,ntypat 440 !do iattyp=1,nattyp(itypat) 441 !iatom=iatom+1 442 !fioncart(:)=psps%ziontypat(itypat)*dtset%efield(:) 443 !do mu=1,3 444 !fionred(mu,iatom)=rprimd(1,mu)*fioncart(1) & 445 !& +rprimd(2,mu)*fioncart(2) & 446 !& +rprimd(3,mu)*fioncart(3) 447 !end do 448 !end do 449 !end do 450 !end if 451 452 !======================================================================= 453 !======= Assemble the various contributions to the forces ============== 454 !======================================================================= 455 456 !Collect grads of etot wrt reduced coordinates 457 !This gives non-symmetrized Hellman-Feynman reduced gradients 458 ABI_ALLOCATE(grtn,(3,dtset%natom)) 459 grtn(:,:)=grl(:,:)+grchempottn(:,:)+grewtn(:,:)+synlgr(:,:)+grxc(:,:) 460 ! grtn(:,:)=grl(:,:)+grewtn(:,:)+synlgr(:,:)+grxc(:,:) 461 462 if (usefock==1 .and. associated(fock).and.fock%fock_common%optfor) then 463 grtn(:,:)=grtn(:,:)+fock%fock_common%forces(:,:) 464 end if 465 if (ngrvdw==dtset%natom) grtn(:,:)=grtn(:,:)+grvdw(:,:) 466 ! note that fionred is subtracted, because it really is a force and we need to 467 ! turn it back into a gradient. The fred2fcart routine below includes the minus 468 ! sign to convert gradients back to forces 469 if ( efield_flag ) grtn(:,:)=grtn(:,:)-fionred(:,:) 470 ! epawf3red is added, because it actually is a gradient, not a force 471 if ( efield_flag .and. calc_epaw3_forces ) grtn(:,:) = grtn(:,:) + epawf3red(:,:) 472 473 !Symmetrize explicitly for given space group and store in grhf : 474 call sygrad(grhf,dtset%natom,grtn,dtset%nsym,symrec,indsym) 475 476 !If residual forces are too large, there must be a problem: cancel them ! 477 if (dtset%usewvl==0.and.abs(dtset%densfor_pred)>0.and.abs(dtset%densfor_pred)/=5) then 478 do iatom=1,dtset%natom 479 do mu=1,3 480 if (abs(gresid(mu,iatom))>10000._dp*abs(grtn(mu,iatom))) gresid(mu,iatom)=zero 481 end do 482 end do 483 end if 484 485 !Add residual potential correction 486 grtn(:,:)=grtn(:,:)+gresid(:,:) 487 488 !Additional stuff for electron-positron 489 ipositron=0 490 if (present(electronpositron)) then 491 if (associated(electronpositron)) then 492 if (allocated(electronpositron%fred_ep)) ipositron=electronpositron_calctype(electronpositron) 493 end if 494 end if 495 if (abs(ipositron)==1) then 496 grtn(:,:)=grtn(:,:)-grxc(:,:)-grchempottn(:,:)-grewtn(:,:)-gresid(:,:)-two*grl(:,:) 497 ! grtn(:,:)=grtn(:,:)-grxc(:,:)-grewtn(:,:)-gresid(:,:)-two*grl(:,:) 498 grl(:,:)=-grl(:,:);grxc(:,:)=zero;gresid(:,:)=zero 499 if (ngrvdw==dtset%natom) grtn(:,:)=grtn(:,:)-grvdw(:,:) 500 if ( dtset%berryopt== 4 .or. dtset%berryopt== 6 .or. dtset%berryopt== 7 .or. & 501 & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17) then 502 grtn(:,:)=grtn(:,:)+fionred(:,:) 503 fionred(:,:)=zero 504 end if 505 end if 506 if (ipositron>0) grtn(:,:)=grtn(:,:)+electronpositron%fred_ep(:,:) 507 508 !Symmetrize all grads explicitly for given space group: 509 if (dtset%usewvl == 0) then 510 call sygrad(fred,dtset%natom,grtn,dtset%nsym,symrec,indsym) 511 else 512 fred = grtn 513 end if 514 515 !Conversion to cartesian coordinates (bohr) AND 516 !Subtract off average force from each force component 517 !to avoid spurious drifting of atoms across cell. 518 ! notice that fred2fcart multiplies fred by -1 to convert it 519 ! from a gradient (input) to a force (output) 520 521 call fred2fcart(favg,(dtset%jellslab==0 .and. dtset%nzchempot==0),fcart,fred,gprimd,dtset%natom) 522 523 !Compute maximal force and maximal difference 524 maxfor=zero;diffor=zero 525 do iatom=1,dtset%natom 526 do mu=1,3 527 if (dtset%iatfix(mu,iatom) /= 1) then 528 maxfor=max(maxfor,abs(fcart(mu,iatom))) 529 diffor=max(diffor,abs(fcart(mu,iatom)-fin(mu,iatom))) 530 else if (dtset%ionmov==4 .or. dtset%ionmov==5) then 531 ! Make the force vanish on fixed atoms when ionmov=4 or 5 532 ! This is because fixing of atom cannot be imposed at the 533 ! level of a routine similar to brdmin or moldyn for these options. 534 fcart(mu,iatom)=zero 535 end if 536 end do 537 end do 538 539 !Apply any generalized constraints to the forces 540 if (dtset%nconeq>0) call constrf(diffor,fcart,forold,fred,dtset%iatfix,dtset%ionmov,maxfor,& 541 & dtset%natom,dtset%nconeq,dtset%prtvol,rprimd,dtset%wtatcon,xred) 542 543 !======================================================================= 544 !Memory deallocations 545 ABI_DEALLOCATE(grl) 546 ABI_DEALLOCATE(grtn) 547 ABI_DEALLOCATE(fin) 548 if ( efield_flag ) then 549 ABI_DEALLOCATE(fionred) 550 if ( calc_epaw3_forces ) then 551 ABI_DEALLOCATE(epawf3red) 552 end if 553 end if 554 555 call timab(69,2,tsec) 556 557 end subroutine forces