TABLE OF CONTENTS
ABINIT/etotfor [ Functions ]
NAME
etotfor
FUNCTION
This routine is called to compute the total energy and various parts of it. The routine computes -if requested- the forces.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG, GMR, 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: electric field is on -> add the contribution of the | - \Omega E.P term to the total energy | /= 4: electric field is off | bfield = cartesian coordinates of magnetic field in atomic units | efield = cartesian coordinates of the electric field in atomic units | iatfix(3,natom)=1 for frozen atom along some 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 symmetry elements in space group | occopt=option for occupancies | prtvol=integer controlling volume of printed output | tsmear=smearing energy or temperature (if metal) | typat(natom)=type integer for each atom in cell | wtatcon(3,natom,nconeq)=weights for atomic constraints gmet(3,3)=metric tensor for G vecs (in bohr**-2) fock <type(fock_type)>= quantities to calculate Fock exact exchange grchempottn(3,natom)=grads of spatially-varying chemical potential energy (hartree) grewtn(3,natom)=grads of Ewald energy (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree) gsqcut=cutoff on (k+G)^2 (bohr^-2) indsym(4,nsym,natom)=indirect indexing array for atom labels kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor 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 nhat(nfft,nspden*usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.f for a description ntypat=number of types of atoms in unit cell. nvresid(nfft,nspden)=potential or density residual n1xccc=dimension of xccc1d ; 0 if no XC core correction is used n3xccc=dimension of the xccc3d array (0 or nfft). optene=option for the computation of total energy (-1=no computation; 0=direct scheme; 1=double-counting scheme) optforces=option for the computation of forces optres=0 if residual array (nvresid) contains the potential residual =1 if residual array (nvresid) contains the density residual pawang <type(pawang_type)>=paw angular mesh 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 phase (structure factor) information. psps <type(pseudopotential_type)>=variables related to pseudopotentials rhog(2,nfft)=array for Fourier transform of electron density 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)=symmetry operations in reciprocal space ucvol=unit cell volume usepaw= 0 for non paw calculation; =1 for paw calculation vhartr(nfft)=array for holding Hartree potential vpsp(nfft)=array for holding local psp vxc(nfft,nspden)=array for holding XC potential xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
deltae=change in total energy between the previous and present SCF cycle etotal=total energy (hartree) ===== if optforces==1 diffor=maximum absolute change in component of forces between present and previous SCF cycle. favg(3)=mean of fcart before correction for translational symmetry fcart(3,natom)=cartesian forces from fred (hartree/bohr) fred(3,natom)=symmetrized form of grtn (grads of Etot) (hartree) gresid(3,natom)=forces due to the residual of the density/potential grhf(3,natom)=Hellman-Feynman derivatives of the total energy grxc(3,natom)=d(Exc)/d(xred) derivatives (0 without core charges) maxfor=maximum absolute value of force synlgr(3,natom)=symmetrized form of grads of Enl (hartree)
SIDE EFFECTS
Input/Output: elast=previous value of the energy, needed to compute deltae, then updated. electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation energies <type(energies_type)>=all part of total energy. | entropy(IN)=entropy due to the occupation number smearing (if metal) | e_localpsp(IN)=local psp energy (hartree) | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree) | e_chempot(IN)=energy from spatially varying chemical potential (hartree) | e_ewald(IN)=Ewald energy (hartree) | e_vdw_dftd(IN)=VdW DFT-D energy | e_hartree(IN)=Hartree part of total energy (hartree units) | e_corepsp(IN)=psp core-core energy | e_hybcomp_E0(IN)=energy compensation energy for the hybrid functionals at frozen density | e_hybcomp_v0(IN)=potential compensation energy for the hybrid functionals at frozen density | e_hybcomp_v (IN)=potential compensation energy for the hybrid functionals at self-consistent density | e_kinetic(IN)=kinetic energy part of total energy. | e_nonlocalpsp(IN)=nonlocal pseudopotential part of total energy. | e_xc(IN)=exchange-correlation energy (hartree) | e_xcdc(IN)=exchange-correlation double-counting energy (hartree) | e_paw(IN)=PAW spherical part energy | e_pawdc(IN)=PAW spherical part double-counting energy | e_elecfield(OUT)=the term of the energy functional that depends explicitely | on the electric field: enefield = -ucvol*E*P | e_magfield(OUT)=the term of the energy functional that depends explicitely | on the magnetic field: e_magfield = -ucvol*E*P | e_entropy(OUT)=entropy energy due to the occupation number smearing (if metal) | this value is %entropy * dtset%tsmear (hartree). ===== if optforces==1 forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr) grnl(3*natom)=gradients of Etot due to nonlocal contributions Input for norm-conserving psps, output for PAW ===== 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)
NOTES
In case of PAW calculations: All computations are done on the fine FFT grid. All variables (nfft,ngfft,mgfft) refer to this fine FFT grid. All arrays (densities/potentials...) are computed on this fine FFT grid. ! Developpers have to be careful when introducing others arrays: they have to be stored on the fine FFT grid. In case of norm-conserving calculations the FFT grid is the usual FFT grid.
PARENTS
scfcv
CHILDREN
forces,nres2vres,pawgrnl,timab
SOURCE
150 #if defined HAVE_CONFIG_H 151 #include "config.h" 152 #endif 153 154 #include "abi_common.h" 155 156 subroutine etotfor(atindx1,deltae,diffor,dtefield,dtset,& 157 & elast,electronpositron,energies,& 158 & etotal,favg,fcart,fock,forold,fred,gmet,grchempottn,gresid,grewtn,grhf,grnl,grvdw,& 159 & grxc,gsqcut,indsym,kxc,maxfor,mgfft,mpi_enreg,my_natom,nattyp,& 160 & nfft,ngfft,ngrvdw,nhat,nkxc,ntypat,nvresid,n1xccc,n3xccc,optene,optforces,optres,& 161 & pawang,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,red_ptot,psps,rhog,rhor,rmet,rprimd,& 162 & symrec,synlgr,ucvol,usepaw,vhartr,vpsp,vxc,wvl,wvl_den,xccc3d,xred) 163 164 use defs_basis 165 use defs_datatypes 166 use defs_abitypes 167 use defs_wvltypes 168 use m_efield 169 use m_profiling_abi 170 use m_fock, only : fock_type 171 use m_pawang, only : pawang_type 172 use m_pawrad, only : pawrad_type 173 use m_pawtab, only : pawtab_type 174 use m_pawfgrtab, only : pawfgrtab_type 175 use m_pawrhoij, only : pawrhoij_type 176 use m_energies, only : energies_type 177 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 178 179 !This section has been created automatically by the script Abilint (TD). 180 !Do not modify the following lines by hand. 181 #undef ABI_FUNC 182 #define ABI_FUNC 'etotfor' 183 use interfaces_18_timing 184 use interfaces_65_paw 185 use interfaces_67_common, except_this_one => etotfor 186 !End of the abilint section 187 188 implicit none 189 190 !Arguments ------------------------------------ 191 !scalars 192 integer,intent(in) :: my_natom,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nkxc,ntypat,optene,optforces 193 integer,intent(in) :: optres,usepaw 194 real(dp),intent(in) :: gsqcut 195 real(dp),intent(inout) :: elast,ucvol 196 real(dp),intent(out) :: deltae,diffor,etotal,maxfor 197 type(MPI_type),intent(in) :: mpi_enreg 198 type(efield_type),intent(in) :: dtefield 199 type(dataset_type),intent(in) :: dtset 200 type(electronpositron_type),pointer :: electronpositron 201 type(energies_type),intent(inout) :: energies 202 type(pawang_type),intent(in) :: pawang 203 type(pseudopotential_type),intent(in) :: psps 204 type(wvl_internal_type), intent(in) :: wvl 205 type(wvl_denspot_type), intent(inout) :: wvl_den 206 type(fock_type),pointer, intent(inout) :: fock 207 !arrays 208 integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom) 209 integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym) 210 real(dp),intent(in) :: gmet(3,3),grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc) 211 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),red_ptot(3) 212 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rmet(3,3),rprimd(3,3) 213 real(dp),intent(in) :: vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden) 214 real(dp),intent(in) :: xccc3d(n3xccc) 215 real(dp),intent(inout) :: forold(3,dtset%natom),grnl(3*dtset%natom) 216 real(dp),intent(inout) :: nhat(nfft,dtset%nspden*psps%usepaw) 217 real(dp),intent(inout),target :: nvresid(nfft,dtset%nspden) 218 real(dp),intent(inout) :: xred(3,dtset%natom) 219 real(dp),intent(out) :: favg(3),fred(3,dtset%natom) 220 real(dp),intent(inout) :: fcart(3,dtset%natom) 221 real(dp),intent(out) :: gresid(3,dtset%natom),grhf(3,dtset%natom) 222 real(dp),intent(inout) :: grxc(3,dtset%natom) 223 real(dp),intent(out) :: synlgr(3,dtset%natom) 224 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 225 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 226 type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw) 227 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 228 229 !Local variables------------------------------- 230 !scalars 231 integer :: comm_grid,dimnhat,ifft,ipositron,ispden,optgr,optgr2,option,optnc,optstr,optstr2,iir,jjr,kkr 232 logical :: apply_residual 233 real(dp) :: eenth,ucvol_ 234 !arrays 235 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 236 real(dp) :: tsec(2),A(3,3),A1(3,3),A_new(3,3),efield_new(3) 237 real(dp) :: dummy(0),nhat_dum(0,0) 238 real(dp),allocatable :: vlocal(:,:) 239 real(dp), ABI_CONTIGUOUS pointer :: resid(:,:) 240 241 ! ********************************************************************* 242 243 call timab(80,1,tsec) 244 245 ipositron=electronpositron_calctype(electronpositron) 246 247 if (optene>-1) then 248 249 ! When the finite-temperature VG broadening scheme is used, 250 ! the total entropy contribution "tsmear*entropy" has a meaning, 251 ! and gather the two last terms of Eq.8 of VG paper 252 ! Warning : might have to be changed for fixed moment calculations 253 if(dtset%occopt>=3 .and. dtset%occopt<=8) then 254 if (abs(dtset%tphysel) < tol10) then 255 energies%e_entropy = - dtset%tsmear * energies%entropy 256 else 257 energies%e_entropy = - dtset%tphysel * energies%entropy 258 end if 259 else 260 energies%e_entropy = zero 261 end if 262 263 ! Turn it into an electric enthalpy, refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009), 264 ! the missing volume is added here 265 energies%e_elecfield=zero 266 if ((dtset%berryopt==4.or.dtset%berryopt==14).and.ipositron/=1) then 267 energies%e_elecfield=-dot_product(dtset%red_efieldbar,red_ptot) !!ebar_i p_i 268 eenth=zero 269 do iir=1,3 270 do jjr=1,3 271 eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !!g^{-1})_ij ebar_i ebar_j 272 end do 273 end do 274 energies%e_elecfield=energies%e_elecfield-eenth*ucvol/(8._dp*pi) 275 end if 276 277 ! Turn it into an internal energy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009), 278 ! but a little different: U=E_ks + (vol/8*pi) * g^{-1})_ij ebar_i ebar_j 279 if ((dtset%berryopt==6.or.dtset%berryopt==16).and.ipositron/=1) then 280 energies%e_elecfield=zero 281 eenth=zero 282 do iir=1,3 283 do jjr=1,3 284 eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! g^{-1})_ij ebar_i ebar_j 285 end do 286 end do 287 energies%e_elecfield=energies%e_elecfield+eenth*ucvol/(8._dp*pi) 288 end if 289 290 ! Calculate internal energy and electric enthalpy for mixed BC case. 291 if (dtset%berryopt==17.and.ipositron/=1) then 292 energies%e_elecfield=zero 293 A(:,:)=(four_pi/ucvol)*rmet(:,:) 294 A1(:,:)=A(:,:) ; A_new(:,:)=A(:,:) 295 efield_new(:)=dtset%red_efield(:) 296 eenth=zero 297 do kkr=1,3 298 if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction 299 ! step 1 add -ebar*p 300 eenth=eenth-dtset%red_efieldbar(kkr)*red_ptot(kkr) 301 ! step 2 chang to e_new (change e to ebar) 302 efield_new(kkr)=dtset%red_efieldbar(kkr) 303 ! step 3 chang matrix A to A1 304 do iir=1,3 305 do jjr=1,3 306 if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr) 307 if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) & 308 & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr) 309 if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr) 310 end do 311 end do 312 A(:,:)=A1(:,:) ; A_new(:,:)=A1(:,:) 313 end if 314 end do ! end fo kkr 315 do iir=1,3 316 do jjr=1,3 317 eenth= eenth+half*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr) 318 end do 319 end do 320 energies%e_elecfield=energies%e_elecfield+eenth 321 end if ! berryopt==17 322 323 ! Turn it into a magnetic enthalpy, by adding orbital electronic contribution 324 energies%e_magfield = zero 325 ! if (dtset%berryopt == 5 .and. ipositron/=1) then 326 ! emag = dot_product(mag_cart,dtset%bfield) 327 ! energies%e_magfield = emag 328 ! end if 329 330 ! Compute total (free)- energy by direct scheme 331 if (optene==0) then 332 etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc + & 333 & energies%e_localpsp + energies%e_corepsp + energies%e_fock+& 334 & energies%e_entropy + energies%e_elecfield + energies%e_magfield+& 335 & energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_hybcomp_v 336 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 337 if (usepaw==0) etotal = etotal + energies%e_nonlocalpsp 338 if (usepaw/=0) etotal = etotal + energies%e_paw 339 end if 340 341 ! Compute total (free) energy by double-counting scheme 342 if (optene==1) then 343 etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc & 344 & - energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc+ energies%e_fock- energies%e_fockdc & 345 & + energies%e_entropy + energies%e_elecfield + energies%e_magfield & 346 & + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 347 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 348 if (usepaw/=0) etotal = etotal + energies%e_pawdc 349 end if 350 351 ! Additional stuff for electron-positron 352 if (dtset%positron/=0) then 353 if (ipositron==0) then 354 energies%e_electronpositron =zero 355 energies%edc_electronpositron=zero 356 else 357 energies%e_electronpositron =electronpositron%e_hartree+electronpositron%e_xc 358 energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc 359 if (usepaw==1) then 360 energies%e_electronpositron =energies%e_electronpositron +electronpositron%e_paw 361 energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc 362 end if 363 end if 364 if (optene==0) electronpositron%e0=etotal 365 if (optene==1) electronpositron%e0=etotal-energies%edc_electronpositron 366 etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron 367 end if 368 369 ! Compute energy residual 370 deltae=etotal-elast 371 elast=etotal 372 end if !optene/=-1 373 374 call timab(80,2,tsec) 375 376 !------Compute forces----------------------------------------------------- 377 378 if (optforces==1) then 379 380 ! PAW: add gradients due to Dij derivatives to non-local term 381 if (usepaw==1) then 382 ABI_ALLOCATE(vlocal,(nfft,dtset%nspden)) 383 do ispden=1,min(dtset%nspden,2) 384 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vlocal,vpsp,vxc) 385 do ifft=1,nfft 386 vlocal(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden) 387 end do 388 end do 389 390 if(dtset%nspden==4)then 391 do ispden=3,4 392 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vlocal,vxc) 393 do ifft=1,nfft 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 dimnhat=0;optgr=1;optgr2=0;optstr=0;optstr2=0 403 comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl 404 call pawgrnl(atindx1,dimnhat,dummy,1,dummy,grnl,gsqcut,mgfft,my_natom, & 405 & dtset%natom, nattyp,nfft,ngfft,nhat_dum,dummy,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,& 406 & pawang,pawfgrtab,pawrhoij,pawtab,ph1d,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=mpi_enreg%comm_fft,& 408 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,paral_kgb=mpi_enreg%paral_kgb) 409 ABI_DEALLOCATE(vlocal) 410 end if 411 412 apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. & 413 & abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5) 414 415 ! If residual is a density residual (and forces from residual asked), 416 ! has to convert it into a potential residual before calling forces routine 417 if (apply_residual) then 418 ABI_ALLOCATE(resid,(nfft,dtset%nspden)) 419 option=0; if (dtset%densfor_pred<0) option=1 420 optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2 421 call nres2vres(dtset,gsqcut,usepaw,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,& 422 & nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,& 423 & rhor,rprimd,usepaw,resid,xccc3d,xred,vxc) 424 else 425 resid => nvresid 426 end if 427 call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,fred,grchempottn,gresid,grewtn,& 428 & grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfft,mpi_enreg,& 429 & n1xccc,n3xccc,nattyp,nfft,ngfft,ngrvdw,ntypat,& 430 & pawrad,pawtab,ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,dtset%usefock,resid,vxc,wvl,wvl_den,xred,& 431 & electronpositron=electronpositron) 432 if (apply_residual) then 433 ABI_DEALLOCATE(resid) 434 end if 435 436 ! Returned fred are full symmetrized gradients of Etotal 437 ! wrt reduced coordinates xred, d(Etotal)/d(xred) 438 ! Forces are contained in array fcart 439 440 else ! if optforces==0 441 fcart=zero 442 fred=zero 443 favg=zero 444 diffor=zero 445 gresid=zero 446 grhf=zero 447 maxfor=zero 448 synlgr=zero 449 end if 450 451 call timab(80,2,tsec) 452 453 end subroutine etotfor