TABLE OF CONTENTS


ABINIT/afterscfloop [ Functions ]

[ Top ] [ Functions ]

NAME

 afterscfloop

FUNCTION

 Perform all calculations needed after the SCF loop, independent of the
 call to scfcv (with or without atomic displacements), and exclusive
 of print or write purposes, or deallocations.

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  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>
  cpus= cpu time limit in seconds
  deltae=change in energy between the previous and present SCF cycle
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables in this dataset
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs (see NOTES at beginning of scfcv)
   | mkmem=maximum number of k points in core memory
   | mpw = maximum number of plane waves
   | natom=number of atoms in cell
   | nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv)
   | nkpt=number of k points in Brillouin zone
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | nsym=number of symmetries in space group
  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-D2 dispersion (hartree)
  gsqcut=cutoff value on G**2 for (large) sphere inside FFT box.
                       gsqcut=(boxcut**2)*dtset%ecut/(2._dp*(Pi**2)
  hdr <type(hdr_type)>=the header of wf, den and pot files
  indsym(4,nsym,natom)=index showing transformation of atom labels
                       under symmetry operations (computed in symatm)
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  istep=number of the SCF iteration
  kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere
  kxc(nfftf,nkxc)=XC kernel
  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 (see NOTES at beginning of scfcv)
  ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
  moved_atm_inside: if==1, the atoms are allowed to move.
  mpi_enreg=informations about MPI parallelization
  my_natom=number of atoms treated by current processor
  n3xccc=dimension of the xccc3d array (0 or nfft).
  nattyp(dtset%ntypat)=number of atoms of each type
  nfftf= - PAW only - number of FFT grid points for the "fine" grid (see NOTES at beginning of scfcv)
  ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
  ngfftf(18)= - PAW only - contain all needed information about 3D FFT  for the "fine" grid
  nhat(nfftf,nspden*psps%usepaw)= -PAW only- compensation density
  nkxc=dimension of kxc
  npwarr(nkpt)=number of planewaves in basis and on boundary for each k
  nvresid(nfftf,nspden)=array for the residual of the density/potential
  occ(mband*nkpt*nsppol)=occupancies of bands at various k points
  optres=0: the potential residual has been computed in scfcv
         1: the density residual has been computed in scfcv
  paw_an(my_natom*usepaw) <type(paw_an_type)>=paw arrays given on angular mesh
  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
  pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
  pel(3)=reduced coordinates of the electronic polarization (a. u.)
  pel_cg(3) = reduced coordinates of the electronic polarization (a. u.)
             computed in the SCF loop
  ph1df(2,3*(2*mgfftf+1)*natom)= - PAW only - 1-dim structure factor phases for the "fine" grid
      (see NOTES at beginning of scfcv)
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  pion(3)=reduced coordinates of the ionic polarization (a. u.)
  prtfor=1 only if forces have to be printed (0 otherwise)
  prtxml=1 if XML file has to be output
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
   | mpsang= 1+maximum angular momentum
  pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
  pwind_alloc = first dimension of pwind
  pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
  res2=density/potential residual (squared)
  resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins
  residm=maximum value from resid array (except for nbdbuf highest bands)
  rhog(2,nfftf)=Fourier transform of total electron density (including compensation density in PAW)
  rhor(nfftf,nspden)=total electron density (including compensation density in PAW)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  stress_needed=1 if stresses are needed, 0 otherwise
  strsxc(6)=xc correction to stress
  symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
  tollist(12)=list of tolerances
  usecprj=1 if cprj datastructure has been allocated
  vhartr(nfftf)=Hartree potential
  vpsp(nfftf)=array for holding local psp
  vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space
  vxcavg=vxc average
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  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

  conv_retcode= Non-zero if convergence is not achieved.
  elfr(nfftf,nspden)=electron localization function
  grhor(nfft,nspden,3)= gradient of electron density in electrons/bohr**4, real space
  lrhor(nfft,nspden)= Laplacian of electron density in electrons/bohr**5, real space
  results_gs <type(results_gs_type)>=results (energy and its components,
   forces and its components, the stress tensor) of a ground-state computation
   (should be made a pure output quantity)
  taug(2,nfftf)=Fourier transform of total kinetic energy density
  taur(nfftf,nspden)=total kinetic energy density in real space
  ==== if forces are required ====
   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.
   fred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
   gresid(3,natom)=forces due to the residual of the 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 tensor is required ====
   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

 computed_forces=1 if forces have been computed, 0 otherwise
 dtefield <type(efield_type)> = variables related to Berry phase
       and electric field calculations (see initberry.f).
       In case dtset%berryopt = 4/6/7/14/16/17, the overlap matrices computed
       in this routine are stored in dtefield%smat in order
       to be used in the electric field calculation.
 dtorbmag <type(orbmag_type)> = variables related to orbital magnetization
  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=energy compensation term for hybrid exchange-correlation energy (hartree) at fixed density
   | e_hybcomp_v0=potential compensation term for hybrid exchange-correlation energy (hartree) at fixed density
   | e_hybcomp_v=potential compensation term for hybrid exchange-correlation energy (hartree) at self-consistent den
   | 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    !!HONG
   |                  on the electric field:
   |                 enefield =  -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j for fixed E/ebar
   |                          =  Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j  for fixed D/d
 etotal=total energy, might be correct by improved polarization computation
 forold(3,natom)=old forces
 xred(3,natom)=reduced dimensionless atomic coordinates
 ===== if dtset%densfor_pred==3 .and. moved_atm_inside==1 =====
   ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases (coarse grid)
   ph1df(2,3*(2*mgfftf+1)*natom)=1-dim structure factor phases (fine PAW grid)
  wvl <type(wvl_data)>=all wavelets data.

NOTES

PARENTS

      scfcv

CHILDREN

      applyprojectorsonthefly,chern_number,denspot_free_history
      eigensystem_info,elpolariz,energies_copy,exchange_electronpositron
      forstr,getph,hdr_update,kswfn_free_scf_data,last_orthon,metric,mkrho
      nhatgrid,nonlop_test,pawcprj_getdim,pawmkrho,pawmkrhoij,prtposcar
      prtrhomxmn,scprqt,setnoccmmp,spin_current,timab,total_energies
      write_energies,wrtout,wvl_eigen_abi2big,wvl_mkrho,wvl_nhatgrid
      wvl_occ_abi2big,wvl_psitohpsi,wvl_rho_abi2big,wvl_tail_corrections
      wvl_vtrial_abi2big,xcden,xmpi_sum,xred2xcart

SOURCE

 198 #if defined HAVE_CONFIG_H
 199 #include "config.h"
 200 #endif
 201 
 202 #include "abi_common.h"
 203 
 204 subroutine afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,&
 205 & deltae,diffor,dtefield,dtfil,dtorbmag,dtset,eigen,electronpositron,elfr,&
 206 & energies,etotal,favg,fcart,fock,forold,fred,grchempottn,&
 207 & gresid,grewtn,grhf,grhor,grvdw,&
 208 & grxc,gsqcut,hdr,indsym,irrzon,istep,kg,kxc,lrhor,maxfor,mcg,mcprj,mgfftf,&
 209 & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,&
 210 & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,&
 211 & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,prtxml,&
 212 & psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,&
 213 & rhog,rhor,rprimd,stress_needed,strsxc,strten,symrec,synlgr,taug,&
 214 & taur,tollist,usecprj,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,&
 215 & xccc3d,xred,ylm,ylmgr,qvpotzero,conv_retcode)
 216 
 217  use defs_basis
 218  use defs_datatypes
 219  use defs_abitypes
 220  use defs_parameters
 221  use defs_wvltypes
 222  use m_energies
 223  use m_errors
 224  use m_profiling_abi
 225  use m_efield
 226  use m_orbmag
 227  use m_ab7_mixing
 228  use m_hdr
 229 
 230  use m_xmpi,             only : xmpi_sum, xmpi_comm_rank,xmpi_comm_size
 231  use m_geometry,         only : xred2xcart, metric
 232  use m_crystal,          only : prtposcar
 233  use m_results_gs ,      only : results_gs_type
 234  use m_electronpositron, only : electronpositron_type,electronpositron_calctype,exchange_electronpositron
 235  use m_dtset,            only : dtset_copy, dtset_free
 236  use m_paw_dmft,         only : paw_dmft_type
 237  use m_pawang,           only : pawang_type
 238  use m_pawrad,           only : pawrad_type
 239  use m_pawtab,           only : pawtab_type
 240  use m_pawrhoij,         only : pawrhoij_type
 241  use m_paw_an,           only : paw_an_type
 242  use m_paw_ij,           only : paw_ij_type
 243  use m_pawfgrtab,        only : pawfgrtab_type
 244  use m_pawcprj,          only : pawcprj_type,pawcprj_getdim
 245  use m_pawfgr,           only : pawfgr_type
 246  use m_fock,             only : fock_type
 247  use m_kg,               only : getph
 248 
 249 #ifdef HAVE_BIGDFT
 250  use m_abi2big
 251  use BigDFT_API, only : last_orthon, &
 252       & kswfn_free_scf_data, denspot_free_history,&
 253       & write_energies, total_energies, XC_potential,&
 254       & eigensystem_info, applyprojectorsonthefly
 255 #endif
 256 
 257 !This section has been created automatically by the script Abilint (TD).
 258 !Do not modify the following lines by hand.
 259 #undef ABI_FUNC
 260 #define ABI_FUNC 'afterscfloop'
 261  use interfaces_14_hidewrite
 262  use interfaces_18_timing
 263  use interfaces_56_xc
 264  use interfaces_62_wvl_wfs
 265  use interfaces_65_paw
 266  use interfaces_66_nonlocal
 267  use interfaces_67_common
 268  use interfaces_94_scfcv, except_this_one => afterscfloop
 269 !End of the abilint section
 270 
 271  implicit none
 272 
 273 !Arguments ------------------------------------
 274 !scalars
 275  integer,intent(in) :: istep,mcg,mcprj,mgfftf,moved_atm_inside,my_natom,n3xccc,nfftf,ngrvdw,nkxc
 276  integer,intent(in) :: optres,prtfor,prtxml,pwind_alloc,stress_needed,usecprj
 277  integer,intent(inout) :: computed_forces
 278  real(dp),intent(in) :: cpus,deltae,gsqcut,res2,residm
 279  real(dp),intent(in) :: qvpotzero
 280  real(dp),intent(inout) :: diffor,etotal,maxfor,vxcavg
 281  type(MPI_type),intent(inout) :: mpi_enreg
 282  type(datafiles_type),intent(in) :: dtfil
 283  type(dataset_type),intent(inout) :: dtset
 284  type(efield_type),intent(inout) :: dtefield
 285  type(orbmag_type),intent(inout) :: dtorbmag
 286  type(electronpositron_type),pointer :: electronpositron
 287  type(energies_type),intent(inout) :: energies
 288  type(hdr_type),intent(inout) :: hdr
 289  type(pawang_type),intent(in) :: pawang
 290  type(pawfgr_type),intent(in) :: pawfgr
 291  type(pseudopotential_type),intent(in) :: psps
 292  type(results_gs_type),intent(inout) :: results_gs
 293  type(wvl_data),intent(inout) :: wvl
 294  type(fock_type),pointer, intent(inout) :: fock
 295 !arrays
 296  integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom)
 297  integer,intent(in) :: indsym(4,dtset%nsym,dtset%natom)
 298  integer,intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 299  integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%ntypat)
 300  integer,intent(in) :: ngfft(18),ngfftf(18),npwarr(dtset%nkpt)
 301  integer,intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym)
 302  integer,intent(out) :: conv_retcode
 303  real(dp),intent(in) :: grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw)
 304  real(dp),intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 305  real(dp),intent(in) :: pwnsfac(2,pwind_alloc)
 306  real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
 307  real(dp),intent(in) :: rprimd(3,3),tollist(12),vpsp(nfftf)
 308  real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden)
 309  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
 310  real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
 311  real(dp),intent(inout) :: cg(2,mcg)
 312  real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
 313  real(dp),intent(inout) :: forold(3,dtset%natom)
 314  real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw)
 315  real(dp),intent(inout) :: nvresid(nfftf,dtset%nspden),pel(3)
 316  real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),pel_cg(3)
 317  real(dp),intent(inout) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
 318  real(dp),intent(inout) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),pion(3)
 319  real(dp),intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),strsxc(6)
 320  real(dp),intent(inout) :: vhartr(nfftf),vxc(nfftf,dtset%nspden)
 321  real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom)
 322  real(dp),intent(inout) :: favg(3),fcart(3,dtset%natom),fred(3,dtset%natom)
 323  real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom)
 324  real(dp),intent(inout) :: grxc(3,dtset%natom),kxc(nfftf,nkxc),strten(6)
 325  real(dp),intent(inout) :: synlgr(3,dtset%natom)
 326  real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taug(:,:),taur(:,:)
 327  type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj)
 328  type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw)
 329  type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw)
 330  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
 331  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
 332  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw)
 333  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
 334 
 335 !Local variables-------------------------------
 336 !scalars
 337  integer,parameter :: response=0
 338  integer :: bantot,bufsz,choice,cplex,ierr,ifft,igrad,ishift,ispden,nfftotf,ngrad
 339  integer :: optcut,optfor,optgr0,optgr1,optgr2,optrad,quit,shft
 340  integer :: spaceComm_fft,tim_mkrho
 341  logical :: test_gylmgr,test_nfgd,test_rfgd
 342  logical :: wvlbigdft=.false.
 343  real(dp) :: c_fermi,dtaur,dtaurzero
 344  real(dp) :: ucvol
 345  character(len=500) :: message
 346  type(paw_dmft_type) :: paw_dmft
 347 #if defined HAVE_BIGDFT
 348  integer :: ia,ii,mband_cprj
 349  logical :: do_last_ortho,compute_wvl_tail=.false.
 350  real(dp) :: dum,eexctx,eh,ekin,eloc,enl,eproj,esicdc,evxc,exc,ucvol_local
 351 #endif
 352 !arrays
 353  real(dp) :: gmet(3,3),gprimd(3,3),pelev(3),rmet(3,3),tsec(2)
 354  real(dp) :: dmatdum(0,0,0,0)
 355  real(dp),allocatable :: mpibuf(:,:),qphon(:),rhonow(:,:,:),sqnormgrhor(:,:)
 356 #if defined HAVE_BIGDFT
 357  integer,allocatable :: dimcprj_srt(:)
 358  real(dp),allocatable :: hpsi_tmp(:),xcart(:,:)
 359 #endif
 360 
 361 ! *************************************************************************
 362 
 363  DBG_ENTER("COLL")
 364 
 365  call timab(250,1,tsec)
 366  call timab(251,1,tsec)
 367 
 368 !Compute different geometric tensor, as well as ucvol, from rprimd
 369  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 370  nfftotf=product(ngfftf(1:3))
 371 
 372 !MPI FFT communicator
 373  spaceComm_fft=mpi_enreg%comm_fft
 374 
 375 !Recompute structure factor phases if atomic positions have changed
 376  if (moved_atm_inside==1) then
 377    if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
 378      call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
 379    else
 380      ph1d(:,:)=ph1df(:,:)
 381    end if
 382  end if
 383 
 384 !----------------------------------------------------------------------
 385 !Wavelet case: transform psi to KS orbitals
 386 !----------------------------------------------------------------------
 387  if (dtset%usewvl == 1) then
 388 
 389 !  wvlbigdft indicates that the BigDFT workflow will be followed
 390    wvlbigdft=(dtset%wvl_bigdft_comp==1)
 391 
 392 #if defined HAVE_BIGDFT
 393 !  Transform to KS orbitals
 394 
 395 !  Need xcart
 396    ABI_ALLOCATE(xcart,(3, dtset%natom))
 397    call xred2xcart(dtset%natom, rprimd, xcart, xred)
 398    ucvol_local=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp)
 399 
 400 !  do_last_ortho in case of direct minimization, since
 401 !  We never diagonalized the hamiltonian and the eigenvalues are unknown.
 402    if (     wvlbigdft) do_last_ortho=(dtset%iscf==0)
 403    if (.not.wvlbigdft) do_last_ortho=(.false.)
 404    if (do_last_ortho) then
 405      call total_energies(wvl%e%energs, istep, mpi_enreg%me_wvl)
 406      call write_energies(istep,0,wvl%e%energs,zero,zero,"FINAL")
 407      if(.not.wvlbigdft) then
 408 !      If ISCF>10, we exit scfcv at a place where bigdft objects
 409 !      do not contain the KS potential. Hence, we copy vtrial to wvl%den
 410        if(dtset%iscf>=10) call wvl_vtrial_abi2big(1,vtrial,wvl%den)
 411 !      hpsi is lost in hpsitopsi, so we recalculate it (needed for last_orthon).
 412        call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
 413 &       istep,1,-1,mpi_enreg%me_wvl,dtset%natom,&
 414 &       nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
 415 &       dum,.false.,evxc,wvl,wvlbigdft,xcart,strsxc)
 416        if (dtset%iscf==0) then
 417          energies%e_kinetic=ekin    ; energies%e_hartree=eh
 418          energies%e_xc=exc          ; energies%e_localpsp=eloc
 419          energies%e_nonlocalpsp=enl ; energies%e_exactX=eexctx
 420          energies%e_sicdc=esicdc    ; energies%e_xcdc=evxc
 421          energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp &
 422 &         + energies%e_xcdc  + two*energies%e_hartree +energies%e_nonlocalpsp
 423        end if
 424      end if
 425      call last_orthon(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,istep,wvl%wfs%ks,wvl%e%energs%evsum,.true.)
 426      if (mpi_enreg%nproc_wvl == 1) nullify(wvl%wfs%ks%psit)
 427      call eigensystem_info(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,0.d0,&
 428      wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,&
 429      wvl%wfs%ks%orbs,wvl%wfs%ks%psi)
 430 !    Copy eigenvalues from BigDFT object to "eigen"
 431      call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs)
 432 !    Copy occupations from BigDFT objects to ABINIT
 433      call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs)
 434    end if
 435 
 436 !  Tail corrections, pending for wvlbigdft==.false.
 437 !  TODO put it at the end of gstate.
 438 !  WVL - maybe compute the tail corrections to energy
 439    compute_wvl_tail=(dtset%tl_radius>tol12.and.wvlbigdft)
 440    if (compute_wvl_tail) then
 441 !    Use the tails to improve energy precision.
 442      call wvl_tail_corrections(dtset, energies, etotal, mpi_enreg, psps, wvl, xcart)
 443    end if
 444 
 445 !  Clean KSwfn parts only needed in the SCF loop.
 446    call kswfn_free_scf_data(wvl%wfs%ks, (mpi_enreg%nproc_wvl > 1))
 447 !  Clean denspot parts only needed in the SCF loop.
 448    call denspot_free_history(wvl%den%denspot)
 449 
 450 !  If WF have been modified, change the density according to the KS projection.
 451    if ( do_last_ortho ) then
 452 
 453 !    Density from new orthogonalized WFs
 454      call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den)
 455 
 456 !    PAW: has to update cprj, rhoij and compensation charge density
 457      if (psps%usepaw==1) then
 458 !      1-Compute cprj
 459        ABI_ALLOCATE(hpsi_tmp,(size(wvl%wfs%ks%hpsi)))
 460        call applyprojectorsonthefly(mpi_enreg%me_wvl,wvl%wfs%ks%orbs,wvl%descr%atoms,wvl%descr%Glr,&
 461 &       xcart,wvl%descr%h(1),wvl%descr%h(2),wvl%descr%h(3),wvl%wfs%ks%lzd%Glr%wfd,&
 462 &       wvl%projectors%nlpsp,wvl%wfs%ks%psi,hpsi_tmp,eproj,&
 463 &       proj_G=wvl%projectors%G,paw=wvl%descr%paw)
 464        ABI_DEALLOCATE(hpsi_tmp)
 465        do ii=1,mcprj
 466          do ia=1,dtset%natom
 467            !Note that cprj should be allocated (i.e. usepcrj=1 imposed in scfcv)
 468            cprj(ia,ii)%cp(:,:)= wvl%descr%paw%cprj(ia,ii)%cp(:,:)
 469          end do
 470        end do
 471 !      2-Compute rhoij
 472        ABI_ALLOCATE(dimcprj_srt,(dtset%natom))
 473        call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
 474        mband_cprj=mcprj/(dtset%nspinor*dtset%mkmem*dtset%nsppol)
 475        paw_dmft%use_sc_dmft=0 ; paw_dmft%use_dmft=0 ! dmft not used here
 476        call pawmkrhoij(atindx,atindx1,cprj,dimcprj_srt,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,&
 477 &       mcprj,dtset%mkmem,mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspinor,dtset%nsppol,&
 478 &       occ,mpi_enreg%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij,dtfil%unpaw,dtset%usewvl,dtset%wtk)
 479        ABI_DEALLOCATE(dimcprj_srt)
 480 !      3-Symetrize rhoij, compute nhat and add it to rhor
 481        call pawmkrho(dum,1,gprimd,0,indsym,0,mpi_enreg,my_natom,dtset%natom,dtset%nspden,dtset%nsym,&
 482 &       dtset%ntypat,mpi_enreg%paral_kgb,pawang,pawfgr,pawfgrtab,dtset%pawprtvol,pawrhoij,pawrhoij,&
 483 &       pawtab,(/zero,zero,zero/),rhog,rhor,rhor,rprimd,dtset%symafm,symrec,dtset%typat,ucvol_local,&
 484 &       dtset%usewvl,xred,pawnhat=nhat)
 485        call wvl_rho_abi2big(1,rhor,wvl%den)
 486      end if
 487    end if
 488 
 489    ABI_DEALLOCATE(xcart)
 490 
 491 #else
 492    BIGDFT_NOTENABLED_ERROR()
 493 #endif
 494  end if
 495 
 496  call timab(251,2,tsec)
 497  call timab(252,1,tsec)
 498 
 499 !----------------------------------------------------------------------
 500 !Polarization Calculation
 501 !----------------------------------------------------------------------
 502 
 503  if(dtset%berryopt/=0)then
 504    call elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,energies%e_elecfield,gprimd,hdr,&
 505 &   kg,dtset%mband,mcg,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,dtset%nkpt,&
 506 &   npwarr,dtset%nsppol,psps%ntypat,pawrhoij,pawtab,pel,pel_cg,pelev,pion,&
 507 &   psps,pwind,pwind_alloc,pwnsfac,rprimd,ucvol,usecprj,xred)
 508  end if
 509 
 510 !----------------------------------------------------------------------
 511 ! Orbital magnetization calculations
 512 !----------------------------------------------------------------------
 513  if(dtset%orbmag==1) then
 514    call chern_number(atindx1,cg,cprj,dtset,dtorbmag,gmet,gprimd,kg,&
 515 &   mcg,size(cprj,2),mpi_enreg,npwarr,pawang,pawrad,pawtab,pwind,pwind_alloc,&
 516 &   symrec,usecprj,psps%usepaw,xred)
 517  end if
 518 
 519  call timab(252,2,tsec)
 520  call timab(253,1,tsec)
 521 
 522 !----------------------------------------------------------------------
 523 !Gradient and Laplacian of the Density Calculation
 524 !----------------------------------------------------------------------
 525 
 526 !We use routine xcden which get gradient of rhor (grhor), and eventually laplacian of rhor (lrhor).
 527  if(dtset%prtgden/=0 .or. dtset%prtlden/=0)then
 528 
 529 !  Compute gradient of the electron density
 530    ngrad=2
 531    cplex=1
 532    ishift=0
 533    ABI_ALLOCATE(rhonow,(nfftf,dtset%nspden,ngrad*ngrad))
 534    if(dtset%prtlden/=0)then
 535      nullify(lrhor)
 536      ABI_ALLOCATE(lrhor,(nfftf,dtset%nspden))
 537    end if
 538    write(message,'(a,a)') ch10, " Compute gradient of the electron density"
 539    call wrtout(ab_out,message,'COLL')
 540    if(dtset%prtlden/=0) then
 541      write(message,'(a)') " and also Compute Laplacian of the electron density"
 542      call wrtout(ab_out,message,'COLL')
 543    end if
 544    write(message,'(a)') "--------------------------------------------------------------------------------"
 545    call wrtout(ab_out,message,'COLL')
 546 
 547    ABI_ALLOCATE(qphon,(3))
 548    qphon(:)=zero
 549    if(dtset%prtlden/=0) then
 550      call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,dtset%paral_kgb,qphon,rhor,rhonow,lrhonow=lrhor)
 551    else
 552      call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,dtset%paral_kgb,qphon,rhor,rhonow)
 553    end if
 554    ABI_DEALLOCATE(qphon)
 555 
 556 !  Copy gradient which has been output in rhonow to grhor (and free rhonow)
 557    nullify(grhor)
 558    ABI_ALLOCATE(grhor,(nfftf,dtset%nspden,3))
 559    do ispden=1,dtset%nspden
 560      do ifft=1,nfftf
 561        grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4)
 562      end do
 563    end do
 564    ABI_DEALLOCATE(rhonow)
 565 
 566    if(dtset%prtgden/=0) then
 567 !    Print result for grhor
 568      write(message,'(a,a)') ch10, " Result for gradient of the electron density for each direction (1,2,3):"
 569      call wrtout(ab_out,message,'COLL')
 570      write(message,'(a,a,a,a)') ch10," 1rst direction:",ch10,&
 571 &     "--------------------------------------------------------------------------------"
 572      call wrtout(ab_out,message,'COLL')
 573      call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,1),optrhor=2,ucvol=ucvol)
 574      write(message,'(a)') "--------------------------------------------------------------------------------"
 575      call wrtout(ab_out,message,'COLL')
 576      write(message,'(a,a,a,a)') ch10," 2nd direction:",ch10,&
 577 &     "--------------------------------------------------------------------------------"
 578      call wrtout(ab_out,message,'COLL')
 579      call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,2),optrhor=2,ucvol=ucvol)
 580      write(message,'(a)') "--------------------------------------------------------------------------------"
 581      call wrtout(ab_out,message,'COLL')
 582      write(message,'(a,a,a,a)') ch10," 3rd direction:",ch10,&
 583 &     "--------------------------------------------------------------------------------"
 584      call wrtout(ab_out,message,'COLL')
 585      call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,3),optrhor=2,ucvol=ucvol)
 586      write(message,'(a)') "--------------------------------------------------------------------------------"
 587      call wrtout(ab_out,message,'COLL')
 588    end if
 589 
 590    if(dtset%prtlden/=0) then
 591 !    Print result for lrhor
 592      write(message,'(a,a)') ch10, " Result for Laplacian of the electron density :"
 593      call wrtout(ab_out,message,'COLL')
 594      write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------"
 595      call wrtout(ab_out,message,'COLL')
 596      call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,lrhor,optrhor=3,ucvol=ucvol)
 597      write(message,'(a)') "--------------------------------------------------------------------------------"
 598      call wrtout(ab_out,message,'COLL')
 599    end if
 600 
 601    write(message,'(a)') "--------------------------------------------------------------------------------"
 602    call wrtout(ab_out,message,'COLL')
 603  end if
 604 
 605 !----------------------------------------------------------------------
 606 !Kinetic Energy Density Calculation
 607 !----------------------------------------------------------------------
 608 
 609  call timab(253,2,tsec)
 610  call timab(254,1,tsec)
 611 
 612 !We use routine mkrho with option=1 to compute kinetic energy density taur (and taug)
 613  if(dtset%prtkden/=0 .or. dtset%prtelf/=0)then
 614    nullify(taug,taur)
 615 !  tauX are reused in outscfcv for output
 616 !  should be deallocated there
 617    ABI_ALLOCATE(taug,(2,nfftf))
 618    ABI_ALLOCATE(taur,(nfftf,dtset%nspden))
 619    tim_mkrho=5
 620    if(dtset%prtelf/=0) then
 621      write(message,'(a,a)') ch10, " Compute ELF"
 622      call wrtout(ab_out,message,'COLL')
 623      write(message,'(a)') "--------------------------------------------------------------------------------"
 624      call wrtout(ab_out,message,'COLL')
 625    end if
 626    write(message,'(a,a)') ch10, " Compute kinetic energy density"
 627    call wrtout(ab_out,message,'COLL')
 628    paw_dmft%use_sc_dmft=0 ! dmft not used here
 629    paw_dmft%use_dmft=0 ! dmft not used here
 630    call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
 631 &   npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
 632 !  Print result
 633    if(dtset%prtkden/=0) then
 634      write(message,'(a,a)') ch10, "Result for kinetic energy density :"
 635      call wrtout(ab_out,message,'COLL')
 636      write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------"
 637      call wrtout(ab_out,message,'COLL')
 638      call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,taur,optrhor=1,ucvol=ucvol)
 639      write(message,'(a)') "--------------------------------------------------------------------------------"
 640      call wrtout(ab_out,message,'COLL')
 641    end if
 642    ABI_DEALLOCATE(taug)
 643  end if
 644 
 645 !----------------------------------------------------------------------
 646 !Electron Localization Function (ELF) Calculation
 647 !----------------------------------------------------------------------
 648 
 649  call timab(254,2,tsec)
 650  call timab(255,1,tsec)
 651 
 652 !We use routine xcden to compute gradient of electron density (grhor),
 653 !NOTE: If GGA is used, gradient of electron density is already computed
 654 !and it is stored in exchange correlation kernel kxc(:,5:7) (nspden=1) or kxc(:,14:19) (nspden=2).
 655 !In order to save memory and do not have the same quantity twice
 656 !in memory we should use kxc.
 657 !But unfortunately only spin up ans spin down component are stored in kxc.
 658 !So far we thus use grhor which contains all component (nspden=4)
 659 !just like rhonow in xcden instead of kxc.
 660 
 661  if((dtset%prtelf/=0))then
 662    if(dtset%nspden<=2) then
 663 
 664      ngrad=2
 665      cplex=1
 666      if((cplex*dtset%nfft)/=nfftf)then
 667        write(message, '(a,a,a,a)' ) ch10,&
 668 &       ' afterscfloop: ERROR -', ch10, &
 669 &       '   The density is complex, ELF analysis cannot be performed.'
 670        call wrtout(std_out,message,'COLL')
 671 !      MSG_ERROR(message)
 672      end if
 673 
 674      if((dtset%prtgden==0) .and. (dtset%prtlden==0)) then
 675 !      Compute gradient of the electron density
 676        ishift=0
 677        ABI_ALLOCATE(rhonow,(nfftf,dtset%nspden,ngrad*ngrad))
 678        write(message,'(a,a)') ch10, " Compute gradient of the electron density"
 679        call wrtout(ab_out,message,'COLL')
 680        ABI_ALLOCATE(qphon,(3))
 681        qphon(:)=zero
 682        call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,dtset%paral_kgb,qphon,rhor,rhonow)
 683        ABI_DEALLOCATE(qphon)
 684 !      Copy gradient which has been output in rhonow to grhor (and free rhonow)
 685        ABI_ALLOCATE(grhor,(nfftf,dtset%nspden,3))
 686        do ispden=1,dtset%nspden
 687          do ifft=1,nfftf
 688            grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4)
 689          end do
 690        end do
 691        ABI_DEALLOCATE(rhonow)
 692      end if
 693 !    Compute square norm of gradient of the electron density (|grhor|**2)
 694      if(dtset%nspden==1)then
 695        ABI_ALLOCATE(sqnormgrhor,(nfftf,dtset%nspden))
 696        do ifft=1,nfftf
 697          sqnormgrhor(ifft,1) = zero
 698        end do
 699      elseif(dtset%nspden==2)then
 700        ABI_ALLOCATE(sqnormgrhor,(nfftf,dtset%nspden+1))
 701 !      because we not only want (total and up quantities, but also down)
 702 !      Indeed after having token the square norm we can not recover the
 703 !      down quantity by substracting total and up quantities (as we do for usual densities)
 704        do ispden=1,dtset%nspden+1
 705          do ifft=1,nfftf
 706            sqnormgrhor(ifft,ispden) = zero
 707          end do
 708        end do
 709      end if
 710 
 711      do igrad=1,3
 712        do ispden=1,dtset%nspden !total (and up)
 713          do ifft=1,nfftf
 714            sqnormgrhor(ifft,ispden) = sqnormgrhor(ifft,ispden) + grhor(ifft,ispden,igrad)**2
 715          end do
 716        end do
 717        if(dtset%nspden==2)then
 718          do ifft=1,nfftf !down
 719            sqnormgrhor(ifft,3) = sqnormgrhor(ifft,3) + (grhor(ifft,1,igrad)-grhor(ifft,2,igrad))**2
 720          end do
 721        end if
 722      end do
 723 
 724 !    Compute electron localization function (ELF) (here it is elfr)
 725 
 726      nullify(elfr)
 727      if(dtset%nspden==1)then
 728        ABI_ALLOCATE(elfr,(nfftf,dtset%nspden))
 729      elseif(dtset%nspden==2)then
 730        ABI_ALLOCATE(elfr,(nfftf,dtset%nspden+1))
 731 !      1rst is total elf, 2nd is spin-up elf, and 3rd is spin-down elf. (elf_tot /= elf_up + elf_down)
 732      end if
 733      c_fermi = 3.0d0/10.0d0*((3.0d0*pi**2)**(2.0d0/3.0d0))
 734 
 735 !    First compute total elf
 736      ispden=1
 737      do ifft=1,nfftf
 738        dtaurzero = c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0)
 739        dtaur = taur(ifft,ispden)
 740        dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden))
 741 !      Ensure that dtaur is always positive or zero, as it should be.
 742        if(dtaur<0.0d0)dtaur=0.0d0
 743 !      To avoid NaN values we check that dtaurzero is not to small compare to dtaur
 744        if(dtaurzero<(1.0d-20*dtaur)) then
 745          elfr(ifft,ispden) = 0.0d0
 746        else
 747          elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
 748 !        For atomic shell studies we could also add the condition that when (dtaur/dtaurzero)
 749 !        is very close to zero we actually set it exactly to zero (or --> elfr = 1.0d0
 750 !        which is the upper limit of elfr.)
 751        end if
 752      end do
 753 
 754 !    If spin-dependent densities are avalaible, compute spin-dependent elf
 755 !    and offer the possibility to compute total elf in an alternative approach
 756 !    (see doc/theory/ELF)
 757      if(dtset%nspden==2)then
 758 
 759 !      alternative approach to the total elf
 760        if(dtset%prtelf==2)then
 761          ispden=1
 762          do ifft=1,nfftf
 763            dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*( rhor(ifft,ispden+1)**(5.0d0/3.0d0) + &
 764 &           (rhor(ifft,ispden) - rhor(ifft,ispden+1))**(5.0d0/3.0d0) )
 765            dtaur = taur(ifft,ispden)
 766            dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+1)/rhor(ifft,ispden+1)) &
 767 &           - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1)))
 768            if(dtaur<0.0d0)dtaur=0.0d0
 769 !          To avoid NaN values we check that dtaurzero is not to small compare to dtaur
 770            if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then
 771              elfr(ifft,ispden) = 0.0d0
 772            else
 773              elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
 774            end if
 775          end do
 776        end if
 777 
 778 !      elf_up
 779        ispden=2
 780        do ifft=1,nfftf
 781          dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0)
 782          dtaur = taur(ifft,ispden)
 783          dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden))
 784          if(dtaur<0.0d0)dtaur=0.0d0
 785 !        To avoid NaN values we check that dtaurzero is not to small compare to dtaur
 786          if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then
 787            elfr(ifft,ispden) = 0.0d0
 788          else
 789            elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
 790          end if
 791        end do
 792 
 793 !      elf_down
 794        ispden=1
 795        do ifft=1,nfftf
 796          dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*(rhor(ifft,ispden)-rhor(ifft,ispden+1))**(5.0d0/3.0d0)
 797          dtaur = taur(ifft,ispden)-taur(ifft,ispden+1)
 798          dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1)))
 799          if(dtaur<0.0d0)dtaur=0.0d0
 800 !        To avoid NaN values we check that dtaurzero is not to small compare to dtaur
 801          if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then
 802            elfr(ifft,ispden+2) = 0.0d0
 803          else
 804            elfr(ifft,ispden+2) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
 805          end if
 806        end do
 807 
 808      end if !endif dtset%nspden==2
 809 
 810 !    Print result for elfr
 811      call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,elfr,optrhor=4,ucvol=ucvol)
 812 
 813      ABI_DEALLOCATE(grhor)
 814      ABI_DEALLOCATE(sqnormgrhor)
 815 
 816    else
 817      message ='ELF is not yet implemented for non collinear spin cases.'
 818      MSG_WARNING(message)
 819 
 820      ABI_ALLOCATE(elfr,(nfftf,dtset%nspden))
 821      do ispden=1,dtset%nspden
 822        do ifft=1,nfftf
 823          elfr(ifft,ispden) = -2.0d0
 824        end do
 825      end do
 826 !    even if elf is not computed we want to finish the abinit run.
 827 !    To ensure that users won't use the _ELF file which will be produced
 828 !    we set elf to -2.0 (a meaningless value)
 829 
 830    end if ! endif dtset%nspden<=2
 831 
 832    write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------"
 833    call wrtout(ab_out,message,'COLL')
 834    write(message,'(a)') " End of ELF section"
 835    call wrtout(ab_out,message,'COLL')
 836 
 837  end if !endif prtelf/=0
 838 
 839 !######################################################################
 840 !Compute forces (if they were not computed during the elec. iterations)
 841 !and stresses (if requested by user)
 842 !----------------------------------------------------------------------
 843 
 844  call timab(255,2,tsec)
 845  call timab(256,1,tsec)
 846 
 847  optfor=0
 848 
 849  if (computed_forces==0.and.dtset%optforces>0.and.dtset%iscf>=0) then
 850    if (dtset%nstep>0.or.dtfil%ireadwf==1) optfor=1
 851  end if
 852 
 853  if (optfor>0.or.stress_needed>0) then
 854 
 855 !  PAW: eventually, compute g_l(r).Y_lm(r) gradients (if not already done)
 856    if (psps%usepaw==1) then
 857      test_nfgd  =any(pawfgrtab(:)%nfgd==0)
 858      test_rfgd  =any(pawfgrtab(:)%rfgd_allocated==0)
 859      test_gylmgr=any(pawfgrtab(:)%gylmgr_allocated==0)
 860      if (test_nfgd.or.&
 861 &     (test_gylmgr.and.dtset%pawstgylm==1).or.&
 862 &     (test_rfgd.and.stress_needed==1.and.dtset%pawstgylm==1).or.&
 863      (test_rfgd.and.dtset%pawstgylm==0)) then
 864        optcut=0;optgr0=0;optgr1=dtset%pawstgylm;optgr2=0
 865        optrad=1-dtset%pawstgylm;if (stress_needed==1) optrad=1
 866        if (dtset%usewvl==0) then
 867          call nhatgrid(atindx1,gmet,my_natom,dtset%natom,nattyp,ngfftf,dtset%ntypat,&
 868 &         optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
 869 &         comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
 870 &         comm_fft=spaceComm_fft,distribfft=mpi_enreg%distribfft)
 871        else
 872          shft=0
 873 #if defined HAVE_BIGDFT
 874          shft=wvl%descr%Glr%d%n1i*wvl%descr%Glr%d%n2i*wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,4)
 875          call wvl_nhatgrid(atindx1,wvl%descr%atoms%astruct%geocode,&
 876 &         wvl%descr%h,wvl%den%denspot%dpbox%i3s,dtset%natom,dtset%natom,&
 877 &         nattyp,psps%ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,&
 878 &         wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,&
 879 &         wvl%den%denspot%dpbox%n3pi,optcut,optgr0,optgr1,optgr2,optrad,&
 880 &         pawfgrtab,pawtab,psps%gth_params%psppar,rprimd,shft,xred)
 881 #endif
 882        end if
 883      end if
 884    end if
 885 
 886    call forstr(atindx1,cg,cprj,diffor,dtefield,dtset,&
 887 &   eigen,electronpositron,energies,favg,fcart,fock,&
 888 &   forold,fred,grchempottn,gresid,grewtn,&
 889 &   grhf,grvdw,grxc,gsqcut,indsym,&
 890 &   kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,&
 891 &   n3xccc,nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,&
 892 &   npwarr,dtset%ntypat,nvresid,occ,optfor,optres,&
 893 &   paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,&
 894 &   psps,rhog,rhor,rprimd,stress_needed,strsxc,strten,symrec,synlgr,&
 895 &   ucvol,usecprj,vhartr,vpsp,vxc,wvl,xccc3d,xred,ylm,ylmgr,qvpotzero)
 896  end if
 897 
 898  if (optfor==1) computed_forces=1
 899  if (optfor==1) diffor=9.9999999999D99
 900  if (stress_needed==0) strten(:)=9.9999999999D99
 901  if (computed_forces==0) fcart(:,:)=9.9999999999D99
 902  if (dtset%prtstm/=0) strten(:)=zero
 903 
 904  call timab(256,2,tsec)
 905  call timab(257,1,tsec)
 906 
 907 !If SCF convergence was not reached (for dtset%nstep>0),
 908 !print a warning to the output file (non-dummy arguments: dtset%nstep,
 909 !residm, diffor - infos from tollist have been saved inside )
 910  choice=3
 911  call scprqt(choice,cpus,deltae,diffor,dtset,&
 912 & eigen,etotal,favg,fcart,energies%e_fermie,dtfil%fnameabo_app_eig,dtfil%filnam_ds(1),&
 913 & 1,dtset%iscf,istep,dtset%kptns,maxfor,&
 914 & moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,&
 915 & dtset%nstep,occ,optres,prtfor,prtxml,quit,&
 916 & res2,resid,residm,response,tollist,psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,&
 917 & electronpositron=electronpositron, fock=fock)
 918 
 919 !output POSCAR and FORCES files, VASP style, for PHON code and friends.
 920  if (dtset%prtposcar == 1) then
 921    call prtposcar(fcart, dtfil%filnam_ds(4), dtset%natom, dtset%ntypat, rprimd, dtset%typat, ucvol, xred, dtset%znucl)
 922  end if ! prtposcar
 923 
 924  if(allocated(qphon))   then
 925    ABI_DEALLOCATE(qphon)
 926  end if
 927 
 928 !get current operator on wavefunctions
 929  if (dtset%prtspcur == 1) then
 930    call spin_current(cg,dtfil,dtset,gprimd,hdr,kg,mcg,mpi_enreg,psps)
 931  end if
 932 
 933 !Electron-positron stuff: if last calculation was a positron minimization,
 934 !exchange electron and positron data in order to
 935 !get electronic quantities in global variables
 936  if (dtset%positron/=0) then
 937    electronpositron%scf_converged=.false.
 938    if (dtset%positron<0.and.electronpositron_calctype(electronpositron)==1) then
 939      call exchange_electronpositron(cg,cprj,dtset,eigen,electronpositron,energies,fred,mcg,mcprj,&
 940 &     mpi_enreg,my_natom,nfftf,ngfftf,nhat,npwarr,occ,paw_an,pawrhoij,rhog,rhor,strten,usecprj,vhartr)
 941    end if
 942  end if
 943 
 944 !If PAW+U and density mixing, has to update nocc_mmp
 945  if (psps%usepaw==1.and.dtset%usepawu>0.and.(dtset%iscf>0.or.dtset%iscf==-3)) then
 946    call setnoccmmp(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,&
 947 &   dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,dtset%pawprtvol,&
 948 &   pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,&
 949 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 950  end if
 951 
 952 !Update the content of the header (evolving variables)
 953  bantot=hdr%bantot
 954  if (dtset%positron==0) then
 955    call hdr_update(hdr,bantot,etotal,energies%e_fermie,residm,rprimd,occ,&
 956 &   pawrhoij,xred,dtset%amu_orig(:,1),&
 957 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 958  else
 959    call hdr_update(hdr,bantot,electronpositron%e0,energies%e_fermie,residm,rprimd,occ,&
 960 &   pawrhoij,xred,dtset%amu_orig(:,1),&
 961 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 962  end if
 963 
 964 #ifdef HAVE_LOTF
 965  if(dtset%ionmov==23 .and. mpi_enreg%nproc_band>1) then
 966    bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom
 967    ABI_ALLOCATE(mpibuf,(3,bufsz))
 968    mpibuf(:,1:dtset%natom)=fred(:,1:dtset%natom)
 969    mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom)
 970    if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom)
 971    mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/))
 972    call xmpi_sum(mpibuf,mpi_enreg%comm_band,ierr)
 973    fred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_band
 974    fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_band
 975    if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_band
 976    strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_band
 977    ABI_DEALLOCATE(mpibuf)
 978  end if
 979 #endif
 980 
 981 !In case of FFT parallelisation, has to synchronize positions and forces
 982 !to avoid numerical noise
 983  if (mpi_enreg%nproc_fft>1) then
 984    bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom
 985    ABI_ALLOCATE(mpibuf,(3,bufsz))
 986    mpibuf(:,1:dtset%natom)=fred(:,1:dtset%natom)
 987    mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom)
 988    if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom)
 989    mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/))
 990    call xmpi_sum(mpibuf,mpi_enreg%comm_fft,ierr)
 991    fred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_fft
 992    fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_fft
 993    if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_fft
 994    strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_fft
 995    ABI_DEALLOCATE(mpibuf)
 996  end if
 997 
 998 !results_gs%energies   = energies
 999  call energies_copy(energies,results_gs%energies)
1000  results_gs%etotal     =etotal
1001  results_gs%deltae     =deltae
1002  results_gs%diffor     =diffor
1003  results_gs%residm     =residm
1004  results_gs%res2       =res2
1005  results_gs%fcart(:,:) =fcart(:,:)
1006  results_gs%fred(:,:)  =fred(:,:)
1007  results_gs%grchempottn(:,:)=grchempottn(:,:)
1008  results_gs%gresid(:,:)=gresid(:,:)
1009  results_gs%grewtn(:,:)=grewtn(:,:)
1010  results_gs%grxc(:,:)  =grxc(:,:)
1011  results_gs%pel(1:3)   =pel(1:3)
1012  results_gs%strten(1:6)=strten(1:6)
1013  results_gs%synlgr(:,:)=synlgr(:,:)
1014  results_gs%vxcavg     =vxcavg
1015  if (ngrvdw>0) results_gs%grvdw(1:3,1:ngrvdw)=grvdw(1:3,1:ngrvdw)
1016  if (dtset%nstep == 0 .and. dtset%occopt>=3.and.dtset%occopt<=8) then
1017    results_gs%etotal = results_gs%etotal - dtset%tsmear * results_gs%entropy
1018  end if
1019 
1020 !This call is only for testing purpose:
1021 !test of the nonlop routine (DFPT vs Finite Differences)
1022  if (dtset%useria==112233) then
1023    call nonlop_test(cg,eigen,dtset%istwfk,kg,dtset%kptns,dtset%mband,mcg,dtset%mgfft,dtset%mkmem,&
1024 &   mpi_enreg,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,dtset%ngfft,dtset%nkpt,&
1025 &   dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,paw_ij,pawtab,&
1026 &   ph1d,psps,rprimd,dtset%typat,xred)
1027  end if
1028 
1029  call timab(257,2,tsec)
1030  call timab(250,2,tsec)
1031 
1032  DBG_EXIT("COLL")
1033 
1034 #if !defined HAVE_BIGDFT
1035  if (.false.) write(std_out,*) vtrial(1,1)
1036 #endif
1037 
1038 end subroutine afterscfloop