TABLE OF CONTENTS


ABINIT/etotfor [ Functions ]

[ Top ] [ 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