./94_scfcv/scfcv.F90

ABINIT/scfcv

ABINIT/scfcv [ Functions ]

NAME

 scfcv

FUNCTION

 Self-consistent-field convergence.
 Conducts set of passes or overall iterations of preconditioned
 conjugate gradient algorithm to converge wavefunctions to
 ground state and optionally to compute forces and energy.
 This routine is called to compute forces for given atomic
 positions or else to do non-SCF band structures.

 Copyright (C) 1998-2018 ABINIT group (XG, GMR, AR, MKV, MT, FJ, MB)
 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

  atindx(natom)=index table for atoms (see gstate.f)
  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cpus= cpu time limit in seconds
  dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs for the "coarse" grid (see NOTES below)
   | mkmem =number of k points treated by this node.
   | mpw=maximum dimensioned size of npw.
   | natom=number of atoms in cell.
   | nfft=(effective) number of FFT grid points (for this processor)
   |    for the "coarse" grid (see NOTES below)
   | nkpt=number of k points
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | nsym=number of symmetry elements in space group
  ecore=core psp energy (part of total energy) (hartree)
  fatvshift=factor to multiply dtset%atvshift
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  nattyp(ntypat)= # atoms of each type.
  ndtpawuj=size of dtpawuj
  npwarr(nkpt)=number of planewaves in basis at this k point
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and
     related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
   | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  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
                           (see initberry.f)
  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

  resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins

SIDE EFFECTS

  cg(2,mcg)=updated wavefunctions; if mkmem>=nkpt, these are kept in a disk file.
  dtefield <type(efield_type)> = variables related to Berry phase
  dtorbmag <type(orbmag_type)> = variables related to orbital magnetization
  dtpawuj(ndtpawuj)= data used for the automatic determination of U
     (relevant only for PAW+U) calculations (see initberry.f)
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  electronpositron <type(electronpositron_type)>=quantities for
     the electron-positron annihilation
  hdr <type(hdr_type)>=the header of wf, den and pot files
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  initialized= if 0 the initialization of the gstate run is not yet finished
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  nfftf=(effective) number of FFT grid points (for this processor)
     for the "fine" grid (see NOTES below)
  occ(mband*nkpt*nsppol)=occupation number for each band (often 2) at each k point
  pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  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)
  rhog(2,nfftf)=array for Fourier transform of electron density
  rhor(nfftf,nspden)=array for electron density in el./bohr**3
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  scf_history <type(scf_history_type)>=arrays obtained from previous
     SCF cycles
  symrec(3,3,nsym)=symmetry operations in reciprocal space
  taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic
     energy density
  taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
  wffnew=struct info for wf disk files
  wvl <type(wvl_data)>=all wavelets data
  xred(3,natom)=reduced dimensionless atomic coordinates
  xred_old(3,natom)= at input, previous reduced dimensionless atomic
     coordinates at output, current xred is transferred to xred_old
  conv_retcode=return code, 0 if convergence was achieved

NOTES

 It is worth to explain THE USE OF FFT GRIDS:
 ============================================
 In case of PAW:
 ---------------
    Two FFT grids are used:
    - A "coarse" FFT grid (defined by ecut)
      for the application of the Hamiltonian on the plane waves basis.
      It is defined by nfft, ngfft, mgfft, ...
      Hamiltonian, wave-functions, density related to WFs (rhor here), ...
      are expressed on this grid.
    - A "fine" FFT grid (defined) by ecutdg)
      for the computation of the density inside PAW spheres.
      It is defined by nfftf, ngfftf, mgfftf, ...
      Total density, potentials, ...
      are expressed on this grid.
 In case of norm-conserving:
 ---------------------------
    - Only the usual FFT grid (defined by ecut) is used.
      It is defined by nfft, ngfft, mgfft, ...
      For compatibility reasons, (nfftf,ngfftf,mgfftf)
      are set equal to (nfft,ngfft,mgfft) in that case.

PARENTS

      m_scfcv

CHILDREN

      ab7_mixing_deallocate,ab7_mixing_new,ab7_mixing_use_disk_cache
      afterscfloop,build_vxc,check_kxc,chkpawovlp,cprj_clean,cprj_paw_alloc
      ctocprj,destroy_distribfft,destroy_mpi_enreg,energies_init,energy
      etotfor,extraprho,fftdatar_write_from_hdr,first_rec,fock2ace
      fock_ace_destroy,fock_bz_destroy,fock_common_destroy,fock_destroy
      fock_init,fock_updatecwaveocc,fourdp,fresid,getcut,getmpw,getng,getph
      gshgg_mkncwrite,hdr_update,init_distribfft,init_distribfft_seq
      init_metricrec,initmpi_seq,initylmg,int2char4,kpgio,metric,mkrho,newrho
      newvtr,nhatgrid,odamix,out_geometry_xml,out_resultsgs_xml,outscfcv
      paw2wvl_ij,paw_an_free,paw_an_init,paw_an_nullify,paw_an_reset_flags
      paw_ij_free,paw_ij_init,paw_ij_nullify,paw_ij_reset_flags,pawcprj_alloc
      pawcprj_free,pawcprj_getdim,pawcprj_reorder,pawdenpot,pawdij
      pawfgrtab_free,pawfgrtab_init,pawmknhat,pawmkrho,pawmkrhoij
      pawtab_get_lsize,pawuj_red,prc_mem_free,prtene,psolver_rhohxc,rhotov
      rhotoxc,scf_history_free,scf_history_init,scprqt,setnoccmmp
      setrhoijpbe0,setsym,setup_positron,setvtr,sphereboundary,status,symdij
      symmetrize_xred,timab,transgrid,update_e_field_vars,vtorho,vtorhorec
      vtorhotf,wf_mixing,wrtout,wvl_cprjreorder,wvl_nhatgrid,xcdata_init
      xmpi_isum,xmpi_sum,xmpi_wait

SOURCE

 156 #if defined HAVE_CONFIG_H
 157 #include "config.h"
 158 #endif
 159 
 160 #include "abi_common.h"
 161 
 162 subroutine scfcv(atindx,atindx1,cg,cpus,dmatpawu,dtefield,dtfil,dtorbmag,dtpawuj,&
 163 &  dtset,ecore,eigen,electronpositron,fatvshift,hdr,indsym,&
 164 &  initialized,irrzon,kg,mcg,mpi_enreg,my_natom,nattyp,ndtpawuj,nfftf,npwarr,occ,&
 165 &  paw_dmft,pawang,pawfgr,pawrad,pawrhoij,pawtab,phnons,psps,pwind,&
 166 &  pwind_alloc,pwnsfac,rec_set,resid,results_gs,rhog,rhor,rprimd,&
 167 &  scf_history,symrec,taug,taur,wffnew,wvl,xred,xred_old,ylm,ylmgr,conv_retcode)
 168 
 169 
 170  use defs_basis
 171  use defs_datatypes
 172  use defs_abitypes
 173  use defs_wvltypes
 174  use defs_parameters
 175  use defs_rectypes
 176  use m_xmpi
 177  use m_profiling_abi
 178  use m_wffile
 179  use m_rec
 180  use m_ab7_mixing
 181  use m_errors
 182  use m_efield
 183  use m_orbmag
 184  use mod_prc_memory
 185  use m_nctk
 186  use m_hdr
 187  use m_xcdata
 188 
 189  use m_fstrings,         only : int2char4, sjoin
 190  use m_geometry,         only : metric
 191  use m_time,             only : abi_wtime, sec2str
 192  use m_exit,             only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in
 193  use m_abi_etsf,         only : abi_etsf_init
 194  use m_mpinfo,           only : destroy_mpi_enreg, iwrite_fftdatar
 195  use m_ioarr,            only : ioarr, fftdatar_write_from_hdr
 196  use m_results_gs ,      only : results_gs_type
 197  use m_scf_history,      only : scf_history_type, scf_history_init, scf_history_free
 198  use m_energies,         only : energies_type, energies_init, energies_copy
 199  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
 200  use m_pawang,           only : pawang_type
 201  use m_pawrad,           only : pawrad_type
 202  use m_pawtab,           only : pawtab_type,pawtab_get_lsize
 203  use m_paw_an,           only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify,&
 204 &                               paw_an_reset_flags
 205  use m_pawfgrtab,        only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
 206  use m_pawrhoij,         only : pawrhoij_type
 207  use m_pawcprj,          only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_reorder, pawcprj_getdim
 208  use m_pawdij,           only : pawdij, symdij,pawdijhat
 209  use m_pawfgr,           only : pawfgr_type
 210  use m_paw_ij,           only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify,&
 211 &                               paw_ij_reset_flags
 212  use m_paw_dmft,         only : paw_dmft_type
 213  use m_fock,             only : fock_type,fock_init,fock_destroy,fock_ACE_destroy,fock_common_destroy,&
 214 &                               fock_BZ_destroy,fock_update_exc,fock_updatecwaveocc
 215  use gwls_hamiltonian,   only : build_vxc
 216 #if defined HAVE_BIGDFT
 217  use BigDFT_API,         only : cprj_clean,cprj_paw_alloc
 218 #endif
 219  use m_io_kss,           only : gshgg_mkncwrite
 220  use m_outxml,           only : out_resultsgs_XML, out_geometry_XML
 221  use m_kg,               only : getcut, getmpw, kpgio, getph
 222 
 223 !This section has been created automatically by the script Abilint (TD).
 224 !Do not modify the following lines by hand.
 225 #undef ABI_FUNC
 226 #define ABI_FUNC 'scfcv'
 227  use interfaces_14_hidewrite
 228  use interfaces_18_timing
 229  use interfaces_32_util
 230  use interfaces_41_geometry
 231  use interfaces_41_xc_lowlevel
 232  use interfaces_43_wvl_wrappers
 233  use interfaces_51_manage_mpi
 234  use interfaces_52_fft_mpi_noabirule
 235  use interfaces_53_ffts
 236  use interfaces_56_recipspace
 237  use interfaces_56_xc
 238  use interfaces_62_poisson
 239  use interfaces_65_paw
 240  use interfaces_66_nonlocal
 241  use interfaces_66_wfs
 242  use interfaces_67_common
 243  use interfaces_68_recursion
 244  use interfaces_68_rsprc
 245  use interfaces_79_seqpar_mpi
 246  use interfaces_94_scfcv, except_this_one => scfcv
 247 !End of the abilint section
 248 
 249  implicit none
 250 
 251 !Arguments ------------------------------------
 252 !scalars
 253  integer,intent(in) :: mcg,my_natom,ndtpawuj,pwind_alloc
 254  integer,intent(inout) :: initialized,nfftf
 255  integer,intent(out) :: conv_retcode
 256  real(dp),intent(in) :: cpus,ecore,fatvshift
 257  type(MPI_type),intent(inout) :: mpi_enreg
 258  type(datafiles_type),intent(in) :: dtfil
 259  type(dataset_type),intent(inout) :: dtset
 260  type(efield_type),intent(inout) :: dtefield
 261  type(orbmag_type),intent(inout) :: dtorbmag
 262  type(electronpositron_type),pointer:: electronpositron
 263  type(hdr_type),intent(inout) :: hdr
 264  type(pawang_type),intent(in) :: pawang
 265  type(pawfgr_type),intent(inout) :: pawfgr
 266  type(pseudopotential_type),intent(in) :: psps
 267  type(recursion_type),intent(inout) :: rec_set
 268  type(results_gs_type),intent(inout) :: results_gs
 269  type(scf_history_type),intent(inout) :: scf_history
 270  type(wffile_type),intent(inout) :: wffnew
 271  type(wvl_data),intent(inout) :: wvl
 272 !arrays
 273  integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom)
 274  integer,intent(in) :: indsym(4,dtset%nsym,dtset%natom)
 275 !no_abirules
 276  integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 277   !(nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise)
 278  integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem)
 279  integer, intent(in) :: nattyp(psps%ntypat),npwarr(dtset%nkpt),pwind(pwind_alloc,2,3)
 280  integer, intent(in) :: symrec(3,3,dtset%nsym)
 281  real(dp), intent(inout) :: cg(2,mcg),dmatpawu(:,:,:,:)
 282  real(dp), intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
 283  real(dp), intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 284  real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 285   !(nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise)
 286  real(dp), intent(in) :: pwnsfac(2,pwind_alloc)
 287  real(dp), intent(in) :: rprimd(3,3)
 288  real(dp), pointer :: rhog(:,:),rhor(:,:)
 289  real(dp), pointer :: taug(:,:),taur(:,:)
 290  real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
 291  real(dp), intent(inout) :: xred(3,dtset%natom)
 292  real(dp), intent(inout) :: xred_old(3,dtset%natom)
 293  real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
 294  real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
 295  type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtpawuj)
 296  type(pawrhoij_type), intent(inout) :: pawrhoij(my_natom*psps%usepaw)
 297  type(pawrad_type), intent(in) :: pawrad(psps%ntypat*psps%usepaw)
 298  type(pawtab_type), intent(in) :: pawtab(psps%ntypat*psps%usepaw)
 299  type(paw_dmft_type), intent(inout) :: paw_dmft
 300 
 301 !Local variables -------------------------
 302 !scalars
 303  integer,parameter :: level=110,response=0,cplex1=1
 304  integer :: afford,bantot,choice
 305  integer :: computed_forces,cplex,cplex_hf,ctocprj_choice,dbl_nnsclo,dielop,dielstrt,dimdmat
 306  integer :: forces_needed,errid,has_dijhat,has_dijnd,has_vhartree,has_dijfock,history_size,usefock
 307 !integer :: dtset_iprcel
 308  integer :: iatom,ider,idir,ierr,iexit,ii,ikpt,impose_dmat,denpot
 309  integer :: initialized0,iorder_cprj,ipert,ipositron,isave_den,isave_kden,iscf10,ispden
 310  integer :: ispmix,istep,istep_fock_outer,istep_mix,istep_updatedfock,itypat,izero,lmax_diel,lpawumax,mband_cprj,mcprj,me,me_wvl
 311  integer :: mgfftdiel,mgfftf,moved_atm_inside,moved_rhor,my_nspinor,n1xccc
 312  integer :: n3xccc,ncpgr,nfftdiel,nfftmix,nfftmix_per_nfft,nfftotf,ngrvdw,nhatgrdim,nk3xc,nkxc
 313  integer :: npawmix,npwdiel,nstep,nzlmopt,optcut,optcut_hf,optene,optgr0,optgr0_hf
 314  integer :: optgr1,optgr2,optgr1_hf,optgr2_hf,option,optrad,optrad_hf,optres,optxc,prtfor,prtxml,quit
 315  integer :: quit_sum,req_cplex_dij,rdwrpaw,shft,spaceComm,spaceComm_fft,spaceComm_wvl,spaceComm_grid
 316  integer :: spare_mem
 317  integer :: stress_needed,sz1,sz2,tim_mkrho,unit_out
 318  integer :: usecprj,usexcnhat,use_hybcomp
 319  integer :: my_quit,quitsum_request,timelimit_exit,usecg,wfmixalg
 320  integer ABI_ASYNC :: quitsum_async
 321  real(dp) :: boxcut,compch_fft,compch_sph,deltae,diecut,diffor,ecut
 322  real(dp) :: ecutf,ecutsus,edum,elast,etotal,evxc,fermie,gsqcut,hyb_mixing,hyb_mixing_sr
 323  real(dp) :: maxfor,res2,residm,ucvol,ucvol_local,val_max
 324  real(dp) :: val_min,vxcavg,vxcavg_dum
 325  real(dp) :: zion,wtime_step,now,prev
 326  character(len=10) :: tag
 327  character(len=1500) :: message
 328  !character(len=500) :: dilatmx_errmsg
 329  character(len=fnlen) :: fildata
 330  type(MPI_type) :: mpi_enreg_diel
 331  type(xcdata_type) :: xcdata
 332  type(energies_type) :: energies
 333  type(ab7_mixing_object) :: mix
 334  logical,parameter :: VERBOSE=.FALSE.
 335  logical :: dummy_nhatgr
 336  logical :: finite_efield_flag=.false.
 337  logical :: recompute_cprj=.false.,reset_mixing=.false.
 338  logical,save :: tfw_activated=.false.
 339  logical :: wvlbigdft=.false.
 340  logical :: non_magnetic_xc
 341 !type(energies_type),pointer :: energies_wvl  ! TO BE ACTIVATED LATER
 342 !arrays
 343  integer :: ngfft(18),ngfftdiel(18),ngfftf(18),ngfftmix(18),npwarr_diel(1)
 344  integer :: npwtot_diel(1)
 345  integer, save :: scfcv_jdtset = 0 ! To simulate iapp behavior
 346  integer, save :: scfcv_itime = 1 ! To simulate iapp behavior
 347  integer,allocatable :: dimcprj(:),dimcprj_srt(:)
 348  integer,allocatable :: gbound_diel(:,:),irrzondiel(:,:,:),kg_diel(:,:)
 349  integer,allocatable :: l_size_atm(:)
 350  integer,allocatable :: indsym_dum(:,:,:),symrec_dum(:,:,:)
 351  logical,pointer :: lmselect_ep(:,:)
 352  real(dp) :: dielar(7),dphase(3),dummy2(6),favg(3),gmet(3,3),gprimd(3,3)
 353  real(dp) :: kpt_diel(3),pel(3),pel_cg(3),pelev(3),pion(3),ptot(3),qpt(3),red_ptot(3) !!REC
 354  real(dp) :: rhodum(1),rmet(3,3),strsxc(6),strten(6),tollist(12)
 355  real(dp) :: tsec(2),vnew_mean(dtset%nspden),vres_mean(dtset%nspden)
 356  real(dp) :: efield_old_cart(3), ptot_cart(3)
 357  real(dp) :: red_efield2(3),red_efield2_old(3)
 358  real(dp) :: vpotzero(2)
 359 ! red_efield1(3),red_efield2(3) is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009)
 360 ! red_efield1(3) for fixed ebar calculation, red_efield2(3) for fixed reduced d calculation, in mixed BC
 361 ! red_efieldbar_lc(3) is local reduced electric field, defined by Eq.(28) of Nat. Phys. suppl. (2009)
 362 ! pbar(3) and dbar(3) are reduced polarization and displacement field, defined by Eq.(27) and (29) Nat. Phys. suppl. (2009)
 363  real(dp),parameter :: k0(3)=(/zero,zero,zero/)
 364  real(dp),allocatable :: dielinv(:,:,:,:,:),dtn_pc(:,:)
 365  real(dp),allocatable :: fcart(:,:),forold(:,:),fred(:,:),gresid(:,:)
 366  real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grhf(:,:),grnl(:),grvdw(:,:),grxc(:,:)
 367  real(dp),allocatable :: kxc(:,:),nhat(:,:),nhatgr(:,:,:),nvresid(:,:)
 368  real(dp),allocatable :: ph1d(:,:),ph1ddiel(:,:),ph1df(:,:)
 369  real(dp),allocatable :: phnonsdiel(:,:,:),rhowfg(:,:),rhowfr(:,:),shiftvector(:)
 370  real(dp),allocatable :: susmat(:,:,:,:,:),synlgr(:,:)
 371  real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:)
 372  real(dp),allocatable :: vxc(:,:),vxc_hybcomp(:,:),vxctau(:,:,:),workr(:,:),xccc3d(:),ylmdiel(:,:)
 373  real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:)
 374  type(scf_history_type) :: scf_history_wf
 375  type(pawcprj_type),allocatable :: cprj(:,:)
 376  type(paw_an_type),allocatable :: paw_an(:)
 377  type(paw_ij_type),allocatable :: paw_ij(:)
 378  type(pawfgrtab_type),allocatable,save :: pawfgrtab(:)
 379  type(pawrhoij_type),pointer :: pawrhoij_ep(:)
 380  type(fock_type),pointer :: fock
 381 
 382 ! *********************************************************************
 383 
 384  _IBM6("Hello, I'm running on IBM6")
 385 
 386  DBG_ENTER("COLL")
 387 
 388  call timab(238,1,tsec)
 389  call timab(54,1,tsec)
 390 
 391  call status(0,dtfil%filstat,iexit,level,'enter')
 392 
 393  ! enable time limit handler if not done in callers.
 394  if (enable_timelimit_in(ABI_FUNC) == ABI_FUNC) then
 395    write(std_out,*)"Enabling timelimit check in function: ",trim(ABI_FUNC)," with timelimit: ",trim(sec2str(get_timelimit()))
 396  end if
 397 
 398 ! Initialise non_magnetic_xc for rhohxc
 399  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
 400 
 401 !######################################################################
 402 !Initializations - Memory allocations
 403 !----------------------------------------------------------------------
 404  lmax_diel = 0
 405 
 406 !MPI communicators
 407  if ((xmpi_paral==1).and.(mpi_enreg%paral_hf==1)) then
 408    spaceComm=mpi_enreg%comm_kpt
 409  else
 410    spaceComm=mpi_enreg%comm_cell
 411  end if
 412  me=xmpi_comm_rank(spaceComm)
 413  spaceComm_fft=mpi_enreg%comm_fft
 414  spaceComm_wvl=mpi_enreg%comm_wvl
 415  me_wvl=mpi_enreg%me_wvl
 416  spaceComm_grid=mpi_enreg%comm_fft
 417  if(dtset%usewvl==1) spaceComm_grid=mpi_enreg%comm_wvl
 418  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 419 
 420 !Save some variables from dataset definition
 421  nstep=dtset%nstep
 422  ecut=dtset%ecut
 423  ecutf=ecut; if (psps%usepaw==1) ecutf=dtset%pawecutdg
 424  if (psps%usepaw==1.and.pawfgr%usefinegrid==1) ecutf=dtset%pawecutdg
 425  iscf10=mod(dtset%iscf,10)
 426  tollist(1)=dtset%tolmxf;tollist(2)=dtset%tolwfr
 427  tollist(3)=dtset%toldff;tollist(4)=dtset%toldfe
 428  tollist(6)=dtset%tolvrs;tollist(7)=dtset%tolrff
 429  tollist(8)=dtset%vdw_df_threshold
 430  dielstrt=0
 431  finite_efield_flag=(dtset%berryopt == 4  .or. &
 432 & dtset%berryopt == 6  .or. &
 433 & dtset%berryopt == 7  .or. &
 434 & dtset%berryopt == 14 .or. &
 435 & dtset%berryopt == 16 .or. &
 436 & dtset%berryopt == 17)
 437 
 438 !Get FFT grid(s) sizes (be careful !)
 439 !See NOTES in the comments at the beginning of this file.
 440  ngfft(:)=dtset%ngfft(:)
 441  if (psps%usepaw==1) then
 442    mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:)
 443  else
 444    mgfftf=dtset%mgfft;ngfftf(:)=ngfft(:)
 445  end if
 446 
 447 !Calculate zion: the total positive charge acting on the valence electrons
 448  zion=zero
 449  do iatom=1,dtset%natom
 450    zion=zion+psps%ziontypat(dtset%typat(iatom))
 451  end do
 452 
 453 !Compute different geometric tensor, as well as ucvol, from rprimd
 454  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 455 
 456 !Fock: be sure that the pointer is initialized to Null.
 457  usefock=dtset%usefock
 458  nullify(fock)
 459 
 460 !Special care in case of WVL
 461 !wvlbigdft indicates that the BigDFT workflow will be followed
 462  wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
 463 !if (wvlbigdft) then   ! TO BE ACTIVATED LATER
 464 !  ABI_ALLOCATE(energies_wvl,)
 465 !end if
 466  ucvol_local = ucvol
 467 #if defined HAVE_BIGDFT
 468  if (dtset%usewvl == 1) then
 469 !  We need to tune the volume when wavelets are used because, not
 470 !  all FFT points are used.
 471 !  ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
 472    ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp)
 473  end if
 474 #endif
 475 
 476 !Some variables need to be initialized/nullified at start
 477  nullify(grhor,lrhor,elfr)
 478  quit=0 ; dbl_nnsclo=0 ; conv_retcode=0
 479  dielop=0 ; strsxc=zero
 480  deltae=zero ; elast=zero ;
 481  vpotzero(:)=zero
 482  ! JWZ April 12 2018: Intel 18 compiler seems to require maxfor initialized,
 483  ! else it dies in scprqt in some scenarios
 484  maxfor=zero
 485  !
 486  results_gs%residm=zero;results_gs%res2=zero
 487  results_gs%deltae=zero;results_gs%diffor=zero
 488  call energies_init(energies)
 489  if (dtset%positron/=0.and.initialized/=0) then
 490    energies%e0_electronpositron =results_gs%energies%e0_electronpositron
 491    energies%e_electronpositron  =results_gs%energies%e_electronpositron
 492    energies%edc_electronpositron=results_gs%energies%edc_electronpositron
 493    maxfor=zero
 494  end if
 495 
 496  ! Initialize fermi level.
 497  !if (dtset%nstep==0 .or. dtset%iscf < 0) then
 498  if ((dtset%nstep==0 .or. dtset%iscf < 0) .and. dtset%plowan_compute==0) then
 499    energies%e_fermie = results_gs%energies%e_fermie
 500    results_gs%fermie = results_gs%energies%e_fermie
 501    write(std_out,*)"in scfcv: results_gs%fermie: ",results_gs%fermie
 502  end if
 503 
 504  select case(dtset%usepotzero)
 505  case(0,1)
 506    energies%e_corepsp   = ecore / ucvol
 507    energies%e_corepspdc = zero
 508  case(2)
 509      ! No need to include the PspCore energy since it is already included in the
 510      ! local pseudopotential  (vpsp)
 511    energies%e_corepsp   = zero
 512    energies%e_corepspdc = zero
 513  end select
 514  if(wvlbigdft) energies%e_corepsp = zero
 515 
 516  fermie=energies%e_fermie
 517  isave_den=0; isave_kden=0 !initial index of density protection file
 518  optres=merge(0,1,dtset%iscf<10)
 519  usexcnhat=0;mcprj=0
 520  initialized0=initialized
 521  if (dtset%tfkinfunc==12) tfw_activated=.true.
 522  ipert=0;idir=0;cplex=1
 523  istep_mix=1
 524  istep_fock_outer=1
 525  ipositron=electronpositron_calctype(electronpositron)
 526  unit_out=0;if (dtset%prtvol >= 10) unit_out=ab_out
 527  nfftotf=product(ngfftf(1:3))
 528 
 529  usecprj=0 ; iorder_cprj=0
 530  if (psps%usepaw==1) then
 531    if (associated(electronpositron)) then
 532      if (dtset%positron/=0.and.electronpositron%dimcprj>0) usecprj=1
 533    end if
 534    if (dtset%prtnabla>0) usecprj=1
 535    if (dtset%extrapwf>0) usecprj=1
 536    if (dtset%usewvl==1)  usecprj=1
 537    if (dtset%pawfatbnd>0)usecprj=1
 538    if (dtset%prtdos==3)  usecprj=1
 539    if (nstep==0) usecprj=0
 540    if (usefock==1)  usecprj=1
 541  end if
 542 
 543 !Stresses and forces flags
 544  forces_needed=0;prtfor=0
 545  if ((dtset%optforces==1.or.dtset%ionmov==4.or.abs(tollist(3))>tiny(0._dp))) then
 546    if (dtset%iscf>0.and.nstep>0) forces_needed=1
 547    if (nstep==0) forces_needed=2
 548    prtfor=1
 549  else if (dtset%iscf>0.and.dtset%optforces==2) then
 550    forces_needed=2
 551  end if
 552 
 553  stress_needed=0
 554  if (dtset%optstress>0.and.dtset%iscf>0.and.dtset%prtstm==0.and. (nstep>0.or.dtfil%ireadwf==1)) stress_needed=1
 555  if (dtset%optstress>0.and.dtset%iscf>0.and.psps%usepaw==1 &
 556 & .and.finite_efield_flag.and.(nstep>0.or.dtfil%ireadwf==1)) stress_needed=1
 557 
 558 !This is only needed for the tddft routine, and does not
 559 !correspond to the intented use of results_gs (should be only
 560 !for output of scfcv
 561  etotal = results_gs%etotal
 562 
 563 !Entering a scfcv loop, printing data to XML file if required.
 564  prtxml=0;if (me==0.and.dtset%prtxml==1) prtxml=1
 565  if (prtxml == 1) then
 566 !  scfcv() will handle a scf loop, so we output the scfcv markup.
 567    write(ab_xml_out, "(A)") '    <scfcvLoop>'
 568    write(ab_xml_out, "(A)") '      <initialConditions>'
 569 !  We output the geometry of the dataset given in argument.
 570 !  xred and rprimd are given independently since dtset only
 571 !  stores original and final values.
 572    call out_geometry_XML(dtset, 4, dtset%natom, rprimd, xred)
 573    write(ab_xml_out, "(A)") '      </initialConditions>'
 574  end if
 575 
 576 !Examine tolerance criteria, and eventually  print a line to the output
 577 !file (with choice=1, the only non-dummy arguments of scprqt are
 578 !nstep, tollist and iscf - still, diffor and res2 are here initialized to 0)
 579  choice=1 ; diffor=zero ; res2=zero
 580  ABI_ALLOCATE(fcart,(3,dtset%natom))
 581  ABI_ALLOCATE(fred,(3,dtset%natom))
 582  fred(:,:)=zero
 583  fcart(:,:)=results_gs%fcart(:,:) ! This is a side effect ...
 584 !results_gs should not be used as input of scfcv
 585 !HERE IS PRINTED THE FIRST LINE OF SCFCV
 586 
 587  call scprqt(choice,cpus,deltae,diffor,dtset,&
 588 & eigen,etotal,favg,fcart,energies%e_fermie,dtfil%fnameabo_app_eig,&
 589 & dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,dtset%kptns,&
 590 & maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,&
 591 & occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,&
 592 & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode)
 593 
 594 !Various allocations (potentials, gradients, ...)
 595  ABI_ALLOCATE(forold,(3,dtset%natom))
 596  ABI_ALLOCATE(grchempottn,(3,dtset%natom))
 597  ABI_ALLOCATE(gresid,(3,dtset%natom))
 598  ABI_ALLOCATE(grewtn,(3,dtset%natom))
 599  ABI_ALLOCATE(grnl,(3*dtset%natom))
 600  ABI_ALLOCATE(grxc,(3,dtset%natom))
 601  ABI_ALLOCATE(synlgr,(3,dtset%natom))
 602  ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
 603  ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
 604  ABI_ALLOCATE(vhartr,(nfftf))
 605  ABI_ALLOCATE(vtrial,(nfftf,dtset%nspden))
 606  ABI_ALLOCATE(vpsp,(nfftf))
 607  ABI_ALLOCATE(vxc,(nfftf,dtset%nspden))
 608  ABI_ALLOCATE(vxctau,(nfftf,dtset%nspden,4*dtset%usekden))
 609 
 610  wfmixalg=dtset%fockoptmix/100
 611  use_hybcomp=0
 612  if(mod(dtset%fockoptmix,100)==11)use_hybcomp=1
 613  ABI_ALLOCATE(vxc_hybcomp,(nfftf,dtset%nspden*use_hybcomp))
 614 
 615  ngrvdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) ngrvdw=dtset%natom
 616  ABI_ALLOCATE(grvdw,(3,ngrvdw))
 617  grchempottn(:,:)=zero
 618  forold(:,:)=zero ; gresid(:,:)=zero ; pel(:)=zero
 619  vtrial(:,:)=zero; vxc(:,:)=zero
 620  n1xccc=0;if (psps%n1xccc/=0) n1xccc=psps%n1xccc
 621  n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
 622  ABI_ALLOCATE(xccc3d,(n3xccc))
 623 
 624 !Allocations/initializations for PAW only
 625  lpawumax=-1
 626 
 627  !Allocate fake cprj -> valgrind is happier and so am I
 628  ABI_DATATYPE_ALLOCATE(cprj,(1,1))
 629  ABI_ALLOCATE(dimcprj,(1))
 630  dimcprj(1) = 1
 631  call pawcprj_alloc(cprj,0,dimcprj)
 632  ABI_DEALLOCATE(dimcprj)
 633  if(psps%usepaw==1) then
 634 !  Variables/arrays related to the fine FFT grid
 635    ABI_ALLOCATE(nhat,(nfftf,dtset%nspden*psps%usepaw))
 636    if (nstep==0) nhat=zero
 637    ABI_DATATYPE_ALLOCATE(pawfgrtab,(my_natom))
 638    if (my_natom>0) then
 639      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,&
 640 &     mpi_atmtab=mpi_enreg%my_atmtab)
 641      call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat,&
 642 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 643      ABI_DEALLOCATE(l_size_atm)
 644    end if
 645    compch_fft=-1.d5
 646    usexcnhat=maxval(pawtab(:)%usexcnhat)
 647    if (usexcnhat==0.and.dtset%ionmov==4.and.dtset%iscf<10) then
 648      message = 'You cannot simultaneously use ionmov=4 and such a PAW psp file !'
 649      MSG_ERROR(message)
 650    end if
 651 
 652 !  Variables/arrays related to the PAW spheres
 653    ABI_DATATYPE_ALLOCATE(paw_ij,(my_natom))
 654    ABI_DATATYPE_ALLOCATE(paw_an,(my_natom))
 655    call paw_an_nullify(paw_an)
 656    call paw_ij_nullify(paw_ij)
 657    has_dijhat=0;if (dtset%iscf==22) has_dijhat=1
 658    has_vhartree=0; if (dtset%prtvha > 0 .or. dtset%prtvclmb > 0) has_vhartree=1
 659    has_dijfock=0; if (usefock==1) has_dijfock=1
 660    has_dijnd=0; req_cplex_dij=1
 661    if(any(abs(dtset%nucdipmom)>tol8)) then
 662      has_dijnd=1; req_cplex_dij=2
 663    end if
 664    call paw_an_init(paw_an,dtset%natom,dtset%ntypat,0,dtset%nspden,&
 665 &   cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,has_vxc_ex=1,has_vhartree=has_vhartree,&
 666 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 667    call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,&
 668 &   dtset%pawspnorb,dtset%natom,dtset%ntypat,dtset%typat,pawtab,&
 669 &   has_dij=1,has_dijfock=has_dijfock,has_dijhartree=1,has_dijnd=has_dijnd,has_dijso=1,has_dijhat=has_dijhat,&
 670 &   has_pawu_occ=1,has_exexch_pot=1,req_cplex_dij=req_cplex_dij,&
 671 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 672    if(dtset%usewvl==1) then
 673      call paw2wvl_ij(1,paw_ij,wvl%descr)
 674    end if
 675    compch_sph=-1.d5
 676    ABI_ALLOCATE(dimcprj,(dtset%natom))
 677    ABI_ALLOCATE(dimcprj_srt,(dtset%natom))
 678    call pawcprj_getdim(dimcprj    ,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'R')
 679    call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
 680    do itypat=1,dtset%ntypat
 681      if (pawtab(itypat)%usepawu>0) lpawumax=max(pawtab(itypat)%lpawu,lpawumax)
 682    end do
 683    if (dtset%usedmatpu/=0.and.lpawumax>0) then
 684      if (2*lpawumax+1/=size(dmatpawu,1).or.2*lpawumax+1/=size(dmatpawu,2)) then
 685        message = 'Incorrect size for dmatpawu!'
 686        MSG_BUG(message)
 687      end if
 688    end if
 689 
 690 !  Allocation of projected WF (optional)
 691    if (usecprj==1) then
 692      iorder_cprj=0
 693      mband_cprj=dtset%mband;if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band
 694      mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
 695      !Was allocated above for valgrind sake so should always be true (safety)
 696      if (allocated(cprj)) then
 697        call pawcprj_free(cprj)
 698        ABI_DATATYPE_DEALLOCATE(cprj)
 699      end if
 700      ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj))
 701      ncpgr=0
 702      if (usefock==1) then
 703        ctocprj_choice = 1
 704        if (dtset%optforces == 1) then
 705          ncpgr = 3 ; ctocprj_choice = 2
 706        end if
 707 !       if (dtset%optstress /= 0) then
 708 !         ncpgr = 6 ; ctocprj_choice = 3
 709 !       end if
 710      end if
 711      call pawcprj_alloc(cprj,ncpgr,dimcprj_srt)
 712 #if defined HAVE_BIGDFT
 713      if (dtset%usewvl==1) then
 714        ABI_DATATYPE_ALLOCATE(wvl%descr%paw%cprj,(dtset%natom,mcprj))
 715        call cprj_paw_alloc(wvl%descr%paw%cprj,0,dimcprj_srt)
 716      end if
 717 #endif
 718    end if
 719 
 720 !  Other variables for PAW
 721    nullify(pawrhoij_ep);if(associated(electronpositron))pawrhoij_ep=>electronpositron%pawrhoij_ep
 722    nullify(lmselect_ep);if(associated(electronpositron))lmselect_ep=>electronpositron%lmselect_ep
 723  else
 724    ABI_ALLOCATE(dimcprj,(0))
 725    ABI_ALLOCATE(dimcprj_srt,(0))
 726    ABI_ALLOCATE(nhat,(0,0))
 727    ABI_DATATYPE_ALLOCATE(paw_ij,(0))
 728    ABI_DATATYPE_ALLOCATE(paw_an,(0))
 729    ABI_DATATYPE_ALLOCATE(pawfgrtab,(0))
 730  end if ! PAW
 731 
 732 !Several parameters and arrays for the SCF mixing:
 733 !These arrays are needed only in the self-consistent case
 734  if (dtset%iscf>=0) then
 735    dielar(1)=dtset%diecut;dielar(2)=dtset%dielng
 736    dielar(3)=dtset%diemac;dielar(4)=dtset%diemix
 737    dielar(5)=dtset%diegap;dielar(6)=dtset%dielam
 738    dielar(7)=dtset%diemix;if (dtset%iscf>=10) dielar(7)=dtset%diemixmag
 739    ABI_ALLOCATE(nvresid,(nfftf,dtset%nspden))
 740    if (nstep==0) nvresid=zero
 741    ABI_ALLOCATE(dtn_pc,(3,dtset%natom))
 742 !  The next arrays are needed if iscf==5 and ionmov==4,
 743 !  but for the time being, they are always allocated
 744    ABI_ALLOCATE(grhf,(3,dtset%natom))
 745 !  Additional allocation for mixing within PAW
 746    npawmix=0
 747    if(psps%usepaw==1) then
 748      do iatom=1,my_natom
 749        itypat=pawrhoij(iatom)%itypat
 750        pawrhoij(iatom)%use_rhoijres=1
 751        sz1=pawrhoij(iatom)%cplex*pawtab(itypat)%lmn2_size
 752        sz2=pawrhoij(iatom)%nspden
 753        ABI_ALLOCATE(pawrhoij(iatom)%rhoijres,(sz1,sz2))
 754        do ispden=1,pawrhoij(iatom)%nspden
 755          pawrhoij(iatom)%rhoijres(:,ispden)=zero
 756        end do
 757        ABI_ALLOCATE(pawrhoij(iatom)%kpawmix,(pawtab(itypat)%lmnmix_sz))
 758        pawrhoij(iatom)%lmnmix_sz=pawtab(itypat)%lmnmix_sz
 759        pawrhoij(iatom)%kpawmix=pawtab(itypat)%kmix
 760        npawmix=npawmix+pawrhoij(iatom)%nspden*pawtab(itypat)%lmnmix_sz*pawrhoij(iatom)%cplex
 761      end do
 762    end if
 763    if (dtset%iscf > 0) then
 764      denpot = AB7_MIXING_POTENTIAL
 765      if (dtset%iscf > 10) denpot = AB7_MIXING_DENSITY
 766      if (psps%usepaw==1.and.dtset%pawmixdg==0 .and. dtset%usewvl==0) then
 767        ispmix=AB7_MIXING_FOURRIER_SPACE;nfftmix=dtset%nfft;ngfftmix(:)=ngfft(:)
 768      else
 769        ispmix=AB7_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:)
 770      end if
 771      !TRangel: added to avoid segfaults with Wavelets
 772      nfftmix_per_nfft=0;if(nfftf>0) nfftmix_per_nfft=(1-nfftmix/nfftf)
 773      call ab7_mixing_new(mix, iscf10, denpot, ispmix, nfftmix, dtset%nspden, npawmix, errid, message, dtset%npulayit)
 774      if (errid /= AB7_NO_ERROR) then
 775        MSG_ERROR(message)
 776      end if
 777      if (dtset%mffmem == 0) then
 778        call ab7_mixing_use_disk_cache(mix, dtfil%fnametmp_fft)
 779      end if
 780 !   else if (dtset%iscf==0.and.dtset%usewvl==1) then
 781 !     ispmix=AB7_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:)
 782    end if
 783  else
 784    ABI_ALLOCATE(nvresid,(0,0))
 785    ABI_ALLOCATE(dtn_pc,(0,0))
 786    ABI_ALLOCATE(grhf,(0,0))
 787  end if ! iscf>0
 788 
 789 !Here, allocate arrays for computation of susceptibility and dielectric matrix or for TDDFT
 790 
 791  if( (nstep>0 .and. dtset%iscf>=0) .or. dtset%iscf==-1 ) then !MF
 792 
 793 !  Here, for TDDFT, artificially set iprcel . Also set a variable to reduce the memory needs.
 794    afford=1
 795    if(dtset%iscf==-1) then
 796 !    dtset%iprcel=21
 797      afford=0
 798    end if
 799 
 800 !  First compute dimensions
 801    if(dtset%iprcel>=21 .or. dtset%iscf==-1)then
 802 !    With dielop=1, the matrices will be computed when istep=dielstrt
 803 !    With dielop=2, the matrices will be computed when istep=dielstrt and 1
 804      dielop=1
 805      if(dtset%iprcel>=41)dielop=2
 806      if((dtset%iprcel >= 71).and.(dtset%iprcel<=79)) dielop=0 !RSkerker preconditioner do not need the susceptibility matrix
 807 !    Immediate computation of dielectric matrix
 808      dielstrt=1
 809 !    Or delayed computation
 810      if(modulo(dtset%iprcel,100)>21 .and. modulo(dtset%iprcel,100)<=29)dielstrt=modulo(dtset%iprcel,100)-20
 811      if(modulo(dtset%iprcel,100)>31 .and. modulo(dtset%iprcel,100)<=39)dielstrt=modulo(dtset%iprcel,100)-30
 812      if(modulo(dtset%iprcel,100)>41 .and. modulo(dtset%iprcel,100)<=49)dielstrt=modulo(dtset%iprcel,100)-40
 813      if(modulo(dtset%iprcel,100)>51 .and. modulo(dtset%iprcel,100)<=59)dielstrt=modulo(dtset%iprcel,100)-50
 814      if(modulo(dtset%iprcel,100)>61 .and. modulo(dtset%iprcel,100)<=69)dielstrt=modulo(dtset%iprcel,100)-60
 815 !    Get diecut, and the fft grid to be used for the susceptibility computation
 816      diecut=abs(dtset%diecut)
 817      if( dtset%diecut<0.0_dp )then
 818        ecutsus=ecut
 819      else
 820        ecutsus= ( sqrt(ecut) *0.5_dp + sqrt(diecut) *0.25_dp )**2
 821      end if
 822 !    Impose sequential calculation
 823      ngfftdiel(1:3)=0 ; ngfftdiel(7)=100 ; ngfftdiel(9)=0; ngfftdiel(8)=dtset%ngfft(8);ngfftdiel(10:18)=0
 824      if(dtset%iscf==-1)ngfftdiel(7)=102
 825 
 826 !    The dielectric stuff is performed in sequential mode; set mpi_enreg_diel accordingly
 827      call initmpi_seq(mpi_enreg_diel)
 828      call getng(dtset%boxcutmin,ecutsus,gmet,k0,mpi_enreg_diel%me_fft,mgfftdiel,nfftdiel,ngfftdiel,&
 829 &     mpi_enreg_diel%nproc_fft,dtset%nsym,mpi_enreg_diel%paral_kgb,dtset%symrel,&
 830 &     use_gpu_cuda=dtset%use_gpu_cuda)
 831 !    Update the fft distribution
 832      call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all')
 833 
 834 !    Compute the size of the dielectric matrix
 835      kpt_diel(1:3)=(/ 0.0_dp, 0.0_dp, 0.0_dp /)
 836      call getmpw(diecut,dtset%exchn2n3d,gmet,(/1/),kpt_diel,mpi_enreg_diel,npwdiel,1)
 837      lmax_diel=0
 838      if (psps%usepaw==1) then
 839        do ii=1,dtset%ntypat
 840          lmax_diel=max(lmax_diel,pawtab(ii)%lcut_size)
 841        end do
 842      end if
 843    else
 844      npwdiel=1
 845      mgfftdiel=1
 846      nfftdiel=1
 847      lmax_diel=0
 848      afford=0
 849    end if
 850 
 851 !  Now, performs allocation
 852    ABI_ALLOCATE(dielinv,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden))
 853    ABI_ALLOCATE(susmat,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden))
 854    ABI_ALLOCATE(kg_diel,(3,npwdiel))
 855    ABI_ALLOCATE(gbound_diel,(2*mgfftdiel+8,2))
 856    ABI_ALLOCATE(irrzondiel,(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 857    ABI_ALLOCATE(phnonsdiel,(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 858    ABI_ALLOCATE(ph1ddiel,(2,3*(2*mgfftdiel+1)*dtset%natom*psps%usepaw))
 859    ABI_ALLOCATE(ylmdiel,(npwdiel,lmax_diel**2))
 860 !  Then, compute the values of different arrays
 861    if(dielop>=1)then
 862 !    Note : npwarr_diel is dummy, npwtot_diel is dummy
 863 !    This kpgio call for going from the suscep FFT grid to the diel sphere
 864      npwarr_diel(1)=npwdiel
 865 
 866      call kpgio(diecut,dtset%exchn2n3d,gmet,(/1/),kg_diel,&
 867 &     kpt_diel,1,(/1/),1,'COLL',mpi_enreg_diel,npwdiel,&
 868 &     npwarr_diel,npwtot_diel,dtset%nsppol)
 869      call sphereboundary(gbound_diel,1,kg_diel,mgfftdiel,npwdiel)
 870 
 871      if (dtset%nsym>1 .and. dtset%iscf>=0 ) then
 872 !      Should replace this initialization of irrzondiel and phnonsdiel through setsym by a direct call to irrzg
 873        ABI_ALLOCATE(indsym_dum,(4,dtset%nsym,dtset%natom))
 874        ABI_ALLOCATE(symrec_dum,(3,3,dtset%nsym))
 875        call setsym(indsym_dum,irrzondiel,dtset%iscf,dtset%natom,&
 876 &       nfftdiel,ngfftdiel,dtset%nspden,dtset%nsppol,dtset%nsym,phnonsdiel,&
 877 &       dtset%symafm,symrec_dum,dtset%symrel,dtset%tnons,dtset%typat,xred)
 878        ABI_DEALLOCATE(indsym_dum)
 879        ABI_DEALLOCATE(symrec_dum)
 880      end if
 881      if (psps%usepaw==1) then
 882        call getph(atindx,dtset%natom,ngfftdiel(1),ngfftdiel(2),&
 883 &       ngfftdiel(3),ph1ddiel,xred)
 884        call initylmg(gprimd,kg_diel,kpt_diel,1,mpi_enreg_diel,&
 885 &       lmax_diel,npwdiel,dtset%nband,1,npwarr_diel,dtset%nsppol,0,&
 886 &       rprimd,ylmdiel,rhodum)
 887      end if
 888    end if
 889 
 890    if(dtset%iprcel>=21 .or. dtset%iscf==-1)then
 891      call destroy_mpi_enreg(mpi_enreg_diel)
 892    end if
 893 
 894  else
 895    npwdiel=1
 896    mgfftdiel=1
 897    nfftdiel=1
 898    afford = 0
 899    ABI_ALLOCATE(susmat,(0,0,0,0,0))
 900    ABI_ALLOCATE(kg_diel,(0,0))
 901    ABI_ALLOCATE(gbound_diel,(0,0))
 902    ABI_ALLOCATE(irrzondiel,(0,0,0))
 903    ABI_ALLOCATE(phnonsdiel,(0,0,0))
 904    ABI_ALLOCATE(ph1ddiel,(0,0))
 905    ABI_ALLOCATE(ylmdiel,(0,0))
 906  end if
 907 
 908  nkxc=0
 909 !TDDFT - For a first coding
 910  if (dtset%iscf==-1 .and. dtset%nspden==1) nkxc=2
 911  if (dtset%iscf==-1 .and. dtset%nspden==2) nkxc=3
 912 !Eventually need kxc-LDA when susceptibility matrix has to be computed
 913  if (dtset%iscf>0.and.modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79)) nkxc=2*min(dtset%nspden,2)-1
 914 !Eventually need kxc-LDA for residual forces (when density mixing is selected)
 915  if (dtset%iscf>=10.and.dtset%usewvl==0.and.forces_needed>0 .and. &
 916 & abs(dtset%densfor_pred)>=1.and.abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5) then
 917    if (dtset%xclevel==1.or.dtset%densfor_pred>=0) nkxc=2*min(dtset%nspden,2)-1
 918    if (dtset%xclevel==2.and.dtset%nspden==1.and.dtset%densfor_pred<0) nkxc=7
 919    if (dtset%xclevel==2.and.dtset%nspden==2.and.dtset%densfor_pred<0) nkxc=19
 920  end if
 921  if (nkxc>0) then
 922    call check_kxc(dtset%ixc,dtset%optdriver)
 923  end if
 924  ABI_ALLOCATE(kxc,(nfftf,nkxc))
 925 
 926 !This flag will be set to 1 just before an eventual change of atomic
 927 !positions inside the iteration, and set to zero when the consequences
 928 !of this change are taken into account.
 929  moved_atm_inside=0
 930 !This flag will be set to 1 if the forces are computed inside the iteration.
 931  computed_forces=0
 932 
 933  if(dtset%wfoptalg==2)then
 934    ABI_ALLOCATE(shiftvector,((dtset%mband+2)*dtset%nkpt))
 935    val_min=-1.0_dp
 936    val_max=zero
 937  else
 938    ABI_ALLOCATE(shiftvector,(1))
 939  end if
 940 
 941  call status(0,dtfil%filstat,iexit,level,'berryphase    ')
 942 
 943 !!PAW+DMFT: allocate structured datatype paw_dmft if dtset%usedmft=1
 944 !call init_sc_dmft(dtset%dmftbandi,dtset%dmftbandf,dtset%mband,dtset%nkpt,&
 945 !&  dtset%nsppol,dtset%usedmft,paw_dmft,dtset%usedmft)
 946 !call print_sc_dmft(paw_dmft)
 947 
 948 !!Electric field initializations: initialize pel_cg(:) and p_ion(:)
 949  call update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,&
 950 & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,&
 951 & dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,ngfft,dtset%nkpt,npwarr,&
 952 & dtset%ntypat,pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,&
 953 & pwind,pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,&
 954 & 0,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr)
 955 
 956  if (dtset%iscf==22) energies%h0=zero
 957 
 958  call timab(54,2,tsec)
 959 
 960 !##################################################################
 961 !PERFORM ELECTRONIC ITERATIONS
 962 !##################################################################
 963 
 964 !Offer option of computing total energy with existing
 965 !wavefunctions when nstep<=0, else do nstep iterations
 966 !Note that for non-self-consistent calculations, this loop will be exited
 967 !after the first call to vtorho
 968 !Pass through the first routines even when nstep==0
 969 
 970  quitsum_request = xmpi_request_null; timelimit_exit = 0
 971  istep_updatedfock=0
 972 
 973 ! start SCF loop
 974  do istep=1,max(1,nstep)
 975    call status(istep,dtfil%filstat,iexit,level,'loop istep    ')
 976 
 977    ! Handle time limit condition.
 978    if (istep == 1) prev = abi_wtime()
 979    if (istep  > 1) then
 980      now = abi_wtime()
 981      wtime_step = now - prev
 982      prev = now
 983      call wrtout(std_out, sjoin(" scfcv: previous iteration took ",sec2str(wtime_step)))
 984 
 985      if (have_timelimit_in(ABI_FUNC)) then
 986        if (istep > 2) then
 987          call xmpi_wait(quitsum_request,ierr)
 988          if (quitsum_async > 0) then
 989            write(message,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())),". Will exit istep loop in scfcv."
 990            MSG_COMMENT(message)
 991            call wrtout(ab_out, message, "COLL")
 992            timelimit_exit = 1
 993            exit
 994          end if
 995        end if
 996 
 997        my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1
 998        call xmpi_isum(my_quit,quitsum_async,spacecomm,quitsum_request,ierr)
 999      end if
1000    end if
1001 
1002    call timab(240,1,tsec)
1003    if (moved_atm_inside==1 .or. istep==1) then
1004 !    ##############################################################
1005 !    The following steps are done once for a given set of atomic
1006 !    coordinates or for the nstep=1 case
1007 !    --------------------------------------------------------------
1008 
1009 !    Eventually symmetrize atomic coordinates over space group elements:
1010      call symmetrize_xred(indsym,dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,xred)
1011 
1012      if (dtset%usewvl == 0) then
1013 !      Get cut-off for g-vectors
1014        if (psps%usepaw==1) then
1015          call wrtout(std_out,' FFT (fine) grid used in SCF cycle:','COLL')
1016        end if
1017        call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf)
1018 
1019 !      Compute structure factor phases and large sphere cut-off (gsqcut):
1020        call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
1021 
1022        if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
1023          call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
1024        else
1025          ph1df(:,:)=ph1d(:,:)
1026        end if
1027      end if
1028 
1029 !    Initialization of atomic data for PAW
1030      if (psps%usepaw==1) then
1031 !      Check for non-overlapping spheres
1032        call chkpawovlp(dtset%natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred)
1033 
1034 !      Identify parts of the rectangular grid where the density has to be calculated
1035        optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
1036        if (forces_needed==1.or.(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0)) then
1037          optgr1=dtset%pawstgylm;if (stress_needed==1) optrad=1; if (dtset%pawprtwf==1) optrad=1
1038        end if
1039 
1040        if(dtset%usewvl==0) then
1041          call nhatgrid(atindx1,gmet,my_natom,dtset%natom,&
1042 &         nattyp,ngfftf,psps%ntypat,optcut,optgr0,optgr1,optgr2,optrad,&
1043 &         pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
1044 &         comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1045 &         comm_fft=spaceComm_fft,distribfft=mpi_enreg%distribfft)
1046        else
1047          shft=0
1048 #if defined HAVE_BIGDFT
1049          shft=wvl%descr%Glr%d%n1i*wvl%descr%Glr%d%n2i*wvl%den%denspot%dpbox%nscatterarr(me_wvl,4)
1050          call wvl_nhatgrid(atindx1,wvl%descr%atoms%astruct%geocode,&
1051 &         wvl%descr%h,wvl%den%denspot%dpbox%i3s,dtset%natom,dtset%natom,&
1052 &         nattyp,psps%ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,&
1053 &         wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,&
1054 &         wvl%den%denspot%dpbox%n3pi,optcut,optgr0,optgr1,optgr2,optrad,&
1055 &         pawfgrtab,pawtab,psps%gth_params%psppar,rprimd,shft,xred)
1056 #endif
1057        end if
1058      end if
1059 
1060 !    If we are inside SCF cycle or inside dynamics over ions,
1061 !    we have to translate the density of previous iteration
1062      moved_rhor=0
1063 
1064      if (initialized/=0.and.dtset%usewvl == 0.and.ipositron/=1.and. &
1065 &     (abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)) then
1066        moved_rhor=1
1067        if (abs(dtset%densfor_pred)==2) then
1068          option=2
1069          ABI_ALLOCATE(workr,(nfftf,dtset%nspden))
1070          call fresid(dtset,gresid,mpi_enreg,nfftf,ngfftf,&
1071 &         psps%ntypat,option,pawtab,rhor,rprimd,&
1072 &         ucvol,workr,xred,xred_old,psps%znuclpsp)
1073          rhor=workr
1074          ABI_DEALLOCATE(workr)
1075        else if (abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6) then
1076          scf_history%icall=scf_history%icall+1
1077          call extraprho(atindx,atindx1,cg,dtset,gmet,gprimd,gsqcut,&
1078 &         scf_history%icall,kg,mcg,mgfftf,mpi_enreg,psps%mqgrid_vl,&
1079 &         my_natom,nattyp,nfftf,ngfftf,npwarr,psps%ntypat,pawrhoij,pawtab,&
1080 &         ph1df,psps,psps%qgrid_vl,rhor,rprimd,scf_history,ucvol,&
1081 &         psps%usepaw,xred,xred_old,ylm,psps%ziontypat,psps%znuclpsp)
1082        end if
1083        call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1084      end if
1085 
1086    end if ! moved_atm_inside==1 .or. istep==1
1087 
1088    !Initialize/Update data in the case of an Exact-exchange (Hartree-Fock) or hybrid XC calculation
1089    hyb_mixing=zero;hyb_mixing_sr=zero
1090    if (usefock==1) then
1091      if (istep==1) then
1092        ! Initialize data_type fock for the calculation
1093        cplex_hf=cplex
1094        if (psps%usepaw==1) cplex_hf=dtset%pawcpxocc
1095        call fock_init(atindx,cplex_hf,dtset,fock,gsqcut,kg,mpi_enreg,nattyp,npwarr,pawang,pawfgr,pawtab,rprimd)
1096        if (fock%fock_common%usepaw==1) then
1097          optcut_hf = 0 ! use rpaw to construct local_pawfgrtab
1098          optgr0_hf = 0; optgr1_hf = 0; optgr2_hf = 0 ! dont need gY terms locally
1099          optrad_hf = 1 ! do store r-R
1100          call nhatgrid(atindx1,gmet,dtset%natom,dtset%natom,nattyp,ngfftf,psps%ntypat,&
1101 &         optcut_hf,optgr0_hf,optgr1_hf,optgr2_hf,optrad_hf,fock%fock_common%pawfgrtab,pawtab,&
1102 &         rprimd,dtset%typat,ucvol,xred,typord=1)
1103          iatom=-1;idir=0
1104          call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,iatom,idir,&
1105 &         iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,&
1106 &         dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft, dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,&
1107 &         dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,&
1108 &         xred,ylm,ylmgr)
1109        end if
1110        if(wfmixalg/=0)then
1111          spare_mem=0
1112          if(spare_mem==1)history_size=wfmixalg ! Not yet coded
1113          if(spare_mem==0)history_size=2*(wfmixalg-1)+1
1114 !        Specific case of simple mixing : always history_size=1
1115          if(wfmixalg==2)history_size=1
1116          scf_history_wf%history_size=history_size
1117          usecg=2
1118          call scf_history_init(dtset,mpi_enreg,usecg,scf_history_wf)
1119        end if
1120      end if
1121 
1122      !Fock energy
1123      energies%e_exactX=zero
1124      if (fock%fock_common%optfor) then
1125        fock%fock_common%forces=zero
1126      end if
1127 
1128      if (istep==1 .or. istep_updatedfock==fock%fock_common%nnsclo_hf) then
1129 
1130        istep_updatedfock=1
1131 
1132        !Possibly mix the wavefunctions from different steps before computing the Fock operator
1133        if(wfmixalg/=0 .and. .not. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) )then
1134          call wf_mixing(atindx1,cg,cprj,dtset,istep_fock_outer,mcg,mcprj,mpi_enreg,&
1135 &         nattyp,npwarr,pawtab,scf_history_wf)
1136          istep_fock_outer=istep_fock_outer+1
1137 
1138 !DEBUG
1139          if(.false.)then
1140          !Update the density, from the newly mixed cg and cprj.
1141          !Be careful: in PAW, rho does not include the compensation density (added later) !
1142            tim_mkrho=6
1143            if (psps%usepaw==1) then
1144              ABI_ALLOCATE(rhowfg,(2,dtset%nfft))
1145              ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden))
1146 !          write(std_out,*) "mkrhogstate"
1147            !From this call, rho does not include the compensation density
1148              call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
1149 &             mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1150              call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
1151 
1152 !          2-Compute rhoij
1153              call pawmkrhoij(atindx,atindx1,cprj,dimcprj,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,&
1154 &             mcprj,dtset%mkmem,mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspinor,dtset%nsppol,&
1155 &             occ,dtset%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij,dtfil%unpaw,&
1156 &             dtset%usewvl,dtset%wtk)
1157 
1158 !          3-Symetrize rhoij, compute nhat and add it to rhor
1159 !          Note pawrhoij_unsym and pawrhoij are the same, which means that pawrhoij cannot be distributed over different atomic sites.
1160              cplex=1;ipert=0;idir=0;qpt(:)=zero
1161              call pawmkrho(compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
1162 &             my_natom,dtset%natom,dtset%nspden,dtset%nsym,dtset%ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
1163 &             dtset%pawprtvol,pawrhoij,pawrhoij,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
1164 &             symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat,rhog=rhog)
1165              ABI_DEALLOCATE(rhowfg)
1166              ABI_DEALLOCATE(rhowfr)
1167 
1168            else
1169 !DEBUG
1170              write(std_out,*)' scfcv : recompute the density after the wf mixing '
1171 !ENDDEBUG
1172              call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
1173 &             mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1174 !DEBUG
1175 !          write(std_out,*)' scfcv : for debugging purposes, set rhor to zero '
1176 !          rhor=zero
1177 !ENDDEBUG
1178              if(dtset%usekden==1)then
1179                call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
1180 &               mpi_enreg,npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1181              end if
1182            end if
1183          end if
1184        end if
1185 !ENDDEBUG
1186 
1187        ! Update data relative to the occupied states in fock
1188        call fock_updatecwaveocc(cg,cprj,dtset,fock,indsym,mcg,mcprj,mpi_enreg,nattyp,npwarr,occ,ucvol)
1189        ! Possibly (re)compute the ACE operator
1190        if(fock%fock_common%use_ACE/=0) then
1191          call fock2ACE(cg,cprj,fock,dtset%istwfk,kg,dtset%kptns,dtset%mband,mcg,mcprj,dtset%mgfft,&
1192 &         dtset%mkmem,mpi_enreg,psps%mpsang,&
1193 &         dtset%mpw,dtset%natom,dtset%natom,dtset%nband,dtset%nfft,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspden,&
1194 &         dtset%nspinor,dtset%nsppol,dtset%ntypat,occ,dtset%optforces,paw_ij,pawtab,ph1d,psps,rprimd,&
1195 &         dtset%typat,usecprj,dtset%use_gpu_cuda,dtset%wtk,xred,ylm)
1196        end if
1197 
1198        !Should place a test on whether there should be the final exit of the istep loop.
1199        !This test should use focktoldfe.
1200        !This should update the info in fock%fock_common%fock_converged.
1201        !For the time being, fock%fock_common%fock_converged=.false. , so the loop end with the maximal value of nstep always,
1202        !except when nnsclo_hf==1 (so the Fock operator is always updated), in which case, the usual exit tests (toldfe, tolvrs, etc)
1203        !work fine.
1204        !if(fock%fock_common%nnsclo_hf==1 .and. fock%fock_common%use_ACE==0)then
1205        if(fock%fock_common%nnsclo_hf==1)then
1206          fock%fock_common%fock_converged=.TRUE.
1207        end if
1208 
1209 !DEBUG
1210 !      endif
1211 !ENDDEBUG
1212 
1213        !Depending on fockoptmix, possibly restart the mixing procedure for the potential
1214        if(mod(dtset%fockoptmix,10)==1)then
1215          istep_mix=1
1216        end if
1217      else
1218        istep_updatedfock=istep_updatedfock+1
1219      end if
1220 
1221      !Used locally
1222      hyb_mixing=fock%fock_common%hyb_mixing ; hyb_mixing_sr=fock%fock_common%hyb_mixing_sr
1223 
1224    end if ! usefock
1225 
1226 !  Initialize/update data in the electron-positron case
1227    if (dtset%positron<0.or.(dtset%positron>0.and.istep==1)) then
1228      call setup_positron(atindx,atindx1,cg,cprj,dtefield,dtfil,dtset,ecore,eigen,&
1229 &     etotal,electronpositron,energies,fock,forces_needed,fred,gmet,gprimd,&
1230 &     grchempottn,grewtn,grvdw,gsqcut,hdr,initialized0,indsym,istep,istep_mix,kg,&
1231 &     kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfftf,ngrvdw,nhat,&
1232 &     nkxc,npwarr,nvresid,occ,optres,paw_ij,pawang,pawfgr,pawfgrtab,&
1233 &     pawrad,pawrhoij,pawtab,ph1df,ph1d,psps,rhog,rhor,rprimd,&
1234 &     stress_needed,strsxc,symrec,ucvol,usecprj,vhartr,vpsp,vxc,&
1235 &     xccc3d,xred,ylm,ylmgr)
1236      ipositron=electronpositron_calctype(electronpositron)
1237    end if
1238 
1239    if ((moved_atm_inside==1 .or. istep==1).or.&
1240 &   (dtset%positron<0.and.istep_mix==1).or.&
1241 &   (mod(dtset%fockoptmix,100)==11 .and. istep_updatedfock==1)) then
1242 !    PAW only: we sometimes have to compute compensation density
1243 !    and eventually add it to density from WFs
1244      nhatgrdim=0
1245      dummy_nhatgr = .False.
1246 !    This is executed only in the positron case.
1247      if (psps%usepaw==1.and.(dtset%positron>=0.or.ipositron/=1) &
1248 &     .and.((usexcnhat==0) &
1249 &     .or.(dtset%xclevel==2.and.(dtfil%ireadwf/=0.or.dtfil%ireadden/=0.or.initialized/=0)) &
1250 &     .or.(dtfil%ireadwf/=0.and.dtfil%ireadden==0.and.initialized==0))) then
1251        call timab(558,1,tsec)
1252        nhatgrdim=0;if (dtset%xclevel==2) nhatgrdim=usexcnhat*dtset%pawnhatxc
1253        ider=2*nhatgrdim;izero=0
1254        if (nhatgrdim>0)   then
1255          ABI_ALLOCATE(nhatgr,(cplex*nfftf,dtset%nspden,3*nhatgrdim))
1256        else
1257          ABI_ALLOCATE(nhatgr,(0,0,0))
1258          dummy_nhatgr = .True.
1259        end if
1260        call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
1261 &       nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,&
1262 &       nhatgr,nhat,pawrhoij,pawrhoij,pawtab,k0,rprimd,ucvol_local,dtset%usewvl,xred,&
1263 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1264 &       comm_fft=spaceComm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
1265 &       distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
1266        if (dtfil%ireadwf/=0.and.dtfil%ireadden==0.and.initialized==0) then
1267          rhor(:,:)=rhor(:,:)+nhat(:,:)
1268          if(dtset%usewvl==0) then
1269            call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1270          end if
1271        end if
1272        call timab(558,2,tsec)
1273      end if
1274 
1275 !    The following steps have been gathered in the setvtr routine:
1276 !    - get Ewald energy and Ewald forces
1277 !    - compute local ionic pseudopotential vpsp
1278 !    - eventually compute 3D core electron density xccc3d
1279 !    - eventually compute vxc and vhartr
1280 !    - set up vtrial
1281 
1282      optene = 4 * optres
1283      if(dtset%iscf==-3) optene=4
1284      if (wvlbigdft) optene = 1 ! VH needed for the WF mixing
1285 
1286      if (.not.allocated(nhatgr))  then
1287        ABI_ALLOCATE(nhatgr,(nfftf,dtset%nspden,3*nhatgrdim))
1288        dummy_nhatgr = .True.
1289      end if
1290 
1291      call setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,&
1292 &     istep,kxc,mgfftf,moved_atm_inside,moved_rhor,mpi_enreg,&
1293 &     nattyp,nfftf,ngfftf,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,psps%ntypat,&
1294 &     n1xccc,n3xccc,optene,pawrad,pawtab,ph1df,psps,rhog,rhor,rmet,rprimd,&
1295 &     strsxc,ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,&
1296 &     xccc3d,xred,electronpositron=electronpositron,&
1297 &     taug=taug,taur=taur,vxc_hybcomp=vxc_hybcomp,vxctau=vxctau,add_tfw=tfw_activated)
1298 
1299      ! set the zero of the potentials here
1300      if(dtset%usepotzero==2) then
1301        vpsp(:) = vpsp(:) + ecore / ( zion * ucvol )
1302      end if
1303 
1304      if(dtset%optdriver==RUNL_GWLS) then
1305        call build_vxc(vxc,nfftf,dtset%nspden)
1306      end if
1307 
1308      if ((nhatgrdim>0.and.nstep>0).or.dummy_nhatgr) then
1309        ABI_DEALLOCATE(nhatgr)
1310      end if
1311 
1312 !    Recursion Initialisation
1313      if(dtset%userec==1 .and. istep==1)  then
1314        rec_set%quitrec = 0
1315 !      --At any step calculate the metric
1316        call Init_MetricRec(rec_set%inf,rec_set%nl%nlpsp,rmet,ucvol,rprimd,xred,dtset%ngfft(1:3),dtset%natom,rec_set%debug)
1317        call destroy_distribfft(rec_set%mpi%distribfft)
1318        call init_distribfft(rec_set%mpi%distribfft,'c',rec_set%mpi%nproc_fft,rec_set%ngfftrec(2),rec_set%ngfftrec(3))
1319        call init_distribfft(rec_set%mpi%distribfft,'f',rec_set%mpi%nproc_fft,dtset%ngfft(2),dtset%ngfft(3))
1320        if(initialized==0) then
1321          call first_rec(dtset,psps,rec_set)
1322        end if
1323      end if
1324 
1325 !    End the condition of atomic position change or istep==1
1326    end if
1327    call timab(240,2,tsec)
1328    call timab(241,1,tsec)
1329 
1330 !  ######################################################################
1331 !  The following steps are done at every iteration
1332 !  ----------------------------------------------------------------------
1333 !  PAW: Compute energies and potentials in the augmentation regions (spheres)
1334 !  Compute pseudopotential strengths (Dij quantities)
1335    if (psps%usepaw==1)then
1336 
1337 !    Local exact exch.: impose occ. matrix if required
1338      if (dtset%useexexch>0) then
1339        call setrhoijpbe0(dtset,initialized0,istep,istep_mix,&
1340 &       spaceComm,my_natom,dtset%natom,dtset%ntypat,pawrhoij,pawtab,dtset%typat,&
1341 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1342      end if
1343 
1344 !    Computation of on-site densities/potentials/energies
1345      nzlmopt=0;if (istep_mix==2.and.dtset%pawnzlm>0) nzlmopt=-1
1346      if (istep_mix>2) nzlmopt=dtset%pawnzlm
1347      call paw_an_reset_flags(paw_an) ! Force the recomputation of on-site potentials
1348      call paw_ij_reset_flags(paw_ij,self_consistent=.true.) ! Force the recomputation of Dij
1349      option=0;if (dtset%iscf>0.and.dtset%iscf<10.and.nstep>0) option=1
1350      call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,ipert,dtset%ixc,my_natom,dtset%natom,&
1351 &     dtset%nspden,psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,pawrad,&
1352 &     pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,&
1353 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1354 &     hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,&
1355 &     electronpositron=electronpositron,vpotzero=vpotzero)
1356 
1357 !    Correct the average potential with the calculated constant vpotzero
1358 !    Correct the total energies accordingly
1359 !    vpotzero(1) = -beta/ucvol
1360 !    vpotzero(2) = -1/ucvol sum_ij rho_ij gamma_ij
1361      write(message,'(a,f14.6,2x,f14.6)') &
1362 &     ' average electrostatic smooth potential [Ha] , [eV]',SUM(vpotzero(:)),SUM(vpotzero(:))*Ha_eV
1363      call wrtout(std_out,message,'COLL')
1364      vtrial(:,:)=vtrial(:,:)+SUM(vpotzero(:))
1365      if(option/=1)then
1366 !      Fix the direct total energy (non-zero only for charged systems)
1367        energies%e_paw=energies%e_paw-SUM(vpotzero(:))*dtset%charge
1368 !      Fix the double counting total energy accordingly (for both charged AND
1369 !      neutral systems)
1370        energies%e_pawdc=energies%e_pawdc-SUM(vpotzero(:))*zion+vpotzero(2)*dtset%charge
1371      end if
1372 
1373 !    PAW+U: impose density matrix if required
1374      if (dtset%usepawu>0.and.(ipositron/=1)) then
1375        impose_dmat=0
1376        if ((istep<=abs(dtset%usedmatpu)).and.(dtset%usedmatpu<0.or.initialized0==0)) impose_dmat=1
1377        if (impose_dmat==1.or.dtset%dmatudiag/=0) then
1378          dimdmat=0;if (impose_dmat==1) dimdmat=2*lpawumax+1
1379          call setnoccmmp(0,dimdmat,&
1380 &         dmatpawu(1:dimdmat,1:dimdmat,1:dtset%nsppol*dtset%nspinor,1:dtset%natpawu*impose_dmat),&
1381 &         dtset%dmatudiag,impose_dmat,indsym,my_natom,dtset%natom,dtset%natpawu,&
1382 &         dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,dtset%pawprtvol,&
1383 &         pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,&
1384 &         comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1385 !        Reinitalize mixing if PAW+U and occupation matrix now allowed to change
1386 !        For experimental purpose...
1387          if ((dtset%userib==1234).and.(istep==abs(dtset%usedmatpu)).and. &
1388 &         (dtset%usedmatpu<0.or.initialized0==0)) reset_mixing=.true.
1389        end if
1390      end if
1391 
1392 !    Dij computation
1393      call timab(561,1,tsec)
1394 
1395      call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,dtset%natom,nfftf,nfftotf,&
1396 &     dtset%nspden,psps%ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,&
1397 &     pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,dtset%spnorbscl,&
1398 &     ucvol_local,dtset%charge,vtrial,vxc,xred,&
1399 &     natvshift=dtset%natvshift,atvshift=dtset%atvshift,fatvshift=fatvshift,&
1400 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1401 &     mpi_comm_grid=spaceComm_grid,&
1402 &     hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,&
1403 &     electronpositron_calctype=ipositron,&
1404 &     electronpositron_pawrhoij=pawrhoij_ep,&
1405 &     electronpositron_lmselect=lmselect_ep,&
1406 &     nucdipmom=dtset%nucdipmom)
1407 
1408 !    Symetrize Dij
1409      call symdij(gprimd,indsym,ipert,my_natom,dtset%natom,dtset%nsym,&
1410 &     psps%ntypat,0,paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
1411 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1412      if (has_dijhat==1) then
1413        call symdij(gprimd,indsym,ipert,my_natom,dtset%natom,dtset%nsym,&
1414 &       psps%ntypat,1,paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
1415 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1416      end if
1417      if(dtset%usewvl==1) then
1418        call paw2wvl_ij(3,paw_ij,wvl%descr)
1419      end if
1420 
1421      call timab(561,2,tsec)
1422    end if
1423 
1424 !  Write out occupancies to dtpawuj-dataset
1425    if (dtset%usepawu>0.and.dtset%macro_uj>0.and.istep>1.and.ipositron/=1) then
1426      call pawuj_red(dtset,dtpawuj,fatvshift,my_natom,dtset%natom,dtset%ntypat,&
1427      paw_ij,pawrad,pawtab,ndtpawuj,&
1428 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1429    end if
1430 
1431    call timab(241,2,tsec)
1432 
1433 
1434 !  No need to continue and call vtorho, when nstep==0
1435    if(nstep==0)exit
1436 
1437 !  ######################################################################
1438 !  The following steps are done only when nstep>0
1439 !  ----------------------------------------------------------------------
1440    call timab(56,1,tsec)
1441 
1442    if(dtset%iscf>=0)then
1443      write(message, '(a,a,i4)' )ch10,' ITER STEP NUMBER  ',istep
1444      call wrtout(std_out,message,'COLL')
1445    end if
1446 
1447    if (dtset%useria == -4242) then
1448      call gshgg_mkncwrite(istep, dtset, dtfil, psps, hdr, pawtab, pawfgr, paw_ij, mpi_enreg, &
1449      rprimd, xred, eigen, npwarr, kg, ylm, ngfft, dtset%nfft, ngfftf, nfftf, vtrial) !electronpositron) ! Optional arguments
1450    end if
1451 
1452 !  The next flag says whether the xred have to be changed in the current iteration
1453    moved_atm_inside=0
1454    ! /< Hack to remove iapp from scfcv
1455    ! for ionmov 4|5 ncycle=1
1456    ! Hence iapp = itime
1457    if ( dtset%jdtset /= scfcv_jdtset ) then
1458      ! new dtset -> reinitialize
1459      scfcv_jdtset = dtset%jdtset
1460      scfcv_itime = 0
1461    end if
1462    if ( istep .eq. 1 ) scfcv_itime = scfcv_itime + 1
1463    if(dtset%ionmov==4 .and. mod(scfcv_itime,2)/=1 .and. dtset%iscf>=0 ) moved_atm_inside=1
1464    if(dtset%ionmov==5 .and. scfcv_itime/=1 .and. istep==1 .and. dtset%iscf>=0) moved_atm_inside=1
1465    ! /< Hack to remove iapp from scfcv
1466 
1467 !  Thomas-Fermi scheme might use a different toldfe criterion
1468    if (dtset%tfkinfunc>0.and.dtset%tfkinfunc/=2) then
1469      tollist(4)=dtset%toldfe;if (.not.tfw_activated) tollist(4)=dtset%tfw_toldfe
1470    end if
1471 
1472 !  The next flag says whether the forces have to be computed in the current iteration
1473    computed_forces=0
1474    if ((dtset%optforces==1 .and. dtset%usewvl == 0).or.(moved_atm_inside==1)) computed_forces=1
1475    if (abs(tollist(3))>tiny(0._dp)) computed_forces=1
1476    if (dtset%iscf<0) computed_forces=0
1477    if ((istep==1).and.(dtset%optforces/=1)) then
1478      if (moved_atm_inside==1) then
1479        write(message,'(5a)')&
1480 &       'Although the computation of forces during electronic iterations',ch10,&
1481 &       'was not required by user, it is done (required by the',ch10,&
1482 &       'choice of ionmov input parameter).'
1483        MSG_WARNING(message)
1484      end if
1485      if (abs(tollist(3))+abs(tollist(7))>tiny(0._dp)) then
1486        write(message,'(5a)')&
1487 &       'Although the computation of forces during electronic iterations',ch10,&
1488 &       'was not required by user, it is done (required by the',ch10,&
1489 &       '"toldff" or "tolrff" tolerance criteria).'
1490        MSG_WARNING(message)
1491      end if
1492    end if
1493    if ((istep==1).and.(dtset%optforces==1).and. dtset%usewvl == 1) then
1494      write(message,'(5a)')&
1495 &     'Although the computation of forces during electronic iterations',ch10,&
1496 &     'was required by user, it has been disable since the tolerence',ch10,&
1497 &     'is not on forces (force computation is expensive in wavelets).'
1498      MSG_WARNING(message)
1499    end if
1500 
1501    call timab(56,2,tsec)
1502 
1503 !  ######################################################################
1504 !  Compute the density rho from the trial potential
1505 !  ----------------------------------------------------------------------
1506 
1507    call timab(242,1,tsec)
1508 !  Compute the density from the trial potential
1509    if (dtset%tfkinfunc==0) then
1510 
1511      if(VERBOSE)then
1512        call wrtout(std_out,'*. Compute the density from the trial potential (vtorho)',"COLL")
1513      end if
1514 
1515      call vtorho(afford,atindx,atindx1,cg,compch_fft,cprj,cpus,dbl_nnsclo,&
1516 &     dielop,dielstrt,dmatpawu,dphase,dtefield,dtfil,dtset,&
1517 &     eigen,electronpositron,energies,etotal,gbound_diel,&
1518 &     gmet,gprimd,grnl,gsqcut,hdr,indsym,irrzon,irrzondiel,&
1519 &     istep,istep_mix,kg,kg_diel,kxc,lmax_diel,mcg,mcprj,mgfftdiel,mpi_enreg,&
1520 &     my_natom,dtset%natom,nattyp,nfftf,nfftdiel,ngfftdiel,nhat,nkxc,&
1521 &     npwarr,npwdiel,res2,psps%ntypat,nvresid,occ,&
1522 &     computed_forces,optres,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,&
1523 &     pawrhoij,pawtab,phnons,phnonsdiel,ph1d,ph1ddiel,psps,fock,&
1524 &     pwind,pwind_alloc,pwnsfac,resid,residm,rhog,rhor,rmet,rprimd,&
1525 &     susmat,symrec,taug,taur,ucvol_local,usecprj,wffnew,vtrial,vxctau,wvl,xred,&
1526 &     ylm,ylmgr,ylmdiel)
1527 
1528    else if (dtset%tfkinfunc==1.or.dtset%tfkinfunc==11.or.dtset%tfkinfunc==12) then
1529      MSG_WARNING('THOMAS FERMI')
1530      call vtorhotf(dtfil,dtset,energies%e_kinetic,energies%e_nonlocalpsp,&
1531 &     energies%entropy,energies%e_fermie,gprimd,grnl,irrzon,mpi_enreg,&
1532 &     dtset%natom,nfftf,dtset%nspden,dtset%nsppol,dtset%nsym,phnons,&
1533 &     rhog,rhor,rprimd,ucvol,vtrial)
1534      residm=zero
1535      energies%e_eigenvalues=zero
1536    end if
1537 
1538 !  Recursion method
1539    if(dtset%userec==1)then
1540      call vtorhorec(dtset,&
1541 &     energies%e_kinetic,energies%e_nonlocalpsp,energies%entropy,energies%e_eigenvalues,&
1542 &     energies%e_fermie,grnl,initialized,irrzon,nfftf,phnons,&
1543 &     rhog,rhor,vtrial,rec_set,istep-nstep,rprimd,gprimd)
1544      residm=zero
1545    end if
1546 
1547    ! Update Fermi level in energies
1548    results_gs%fermie = energies%e_fermie
1549 
1550    if(dtset%wfoptalg==2)then
1551      do ikpt=1,dtset%nkpt
1552        shiftvector(1+(ikpt-1)*(dtset%mband+2))=val_min
1553        shiftvector(2+(ikpt-1)*(dtset%mband+2):ikpt*(dtset%mband+2)-1)=&
1554 &       eigen((ikpt-1)*dtset%mband+1:ikpt*dtset%mband)
1555        shiftvector(ikpt*(dtset%mband+2))=val_max
1556      end do
1557    end if
1558 
1559    call timab(242,2,tsec)
1560 
1561     !if (dtset%useria == -4242) then
1562     !  call gshgg_mkncwrite(istep, dtset, dtfil, psps, hdr, pawtab, pawfgr, paw_ij, mpi_enreg, &
1563     !     rprimd, xred, eigen, npwarr, kg, ylm, ngfft, dtset%nfft, ngfftf, nfftf, vtrial) !electronpositron) ! Optional arguments
1564     !end if
1565 
1566 !  ######################################################################
1567 !  Skip out of step loop if non-SCF (completed)
1568 !  ----------------------------------------------------------------------
1569 
1570 !  Indeed, nstep loops have been done inside vtorho
1571    if (dtset%iscf<0) exit
1572 
1573 !  ######################################################################
1574 !  In case of density mixing or wavelet handling, compute the total energy
1575 !  ----------------------------------------------------------------------
1576    call timab(60,1,tsec)
1577    if (dtset%iscf>=10 .or. wvlbigdft) then
1578      optene = 1  ! use double counting scheme (default)
1579      if (wvlbigdft.and.dtset%iscf==0) optene = 0 ! use direct scheme
1580      if (dtset%iscf==22) optene = -1
1581 
1582 !    Add the Fock contribution to E_xc and E_xcdc if required
1583      if (usefock==1) then
1584        energies%e_fockdc=two*energies%e_fock
1585      end if
1586 
1587 !    if the mixing is the ODA mixing, compute energy and new density here
1588      if (dtset%iscf==22) then
1589        call odamix(deltae,dtset,&
1590 &       elast,energies,etotal,gprimd,gsqcut,kxc,mpi_enreg,&
1591 &       my_natom,nfftf,ngfftf,nhat,nkxc,psps%ntypat,nvresid,n3xccc,optres,&
1592 &       paw_ij,paw_an,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,&
1593 &       red_ptot,psps,rhog,rhor,rprimd,strsxc,ucvol,psps%usepaw,&
1594 &       usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,xccc3d,xred,&
1595 &       taug=taug,taur=taur,vxctau=vxctau,add_tfw=tfw_activated)
1596      end if
1597 !    If the density mixing is required, compute the total energy here
1598      call etotfor(atindx1,deltae,diffor,dtefield,dtset,&
1599 &     elast,electronpositron,energies,&
1600 &     etotal,favg,fcart,fock,forold,fred,gmet,grchempottn,gresid,grewtn,grhf,grnl,grvdw,&
1601 &     grxc,gsqcut,indsym,kxc,maxfor,mgfftf,mpi_enreg,my_natom,&
1602 &     nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,psps%ntypat,nvresid,n1xccc,n3xccc,&
1603 &     optene,computed_forces,optres,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,&
1604 &     ph1df,red_ptot,psps,rhog,rhor,rmet,rprimd,symrec,synlgr,ucvol,&
1605 &     psps%usepaw,vhartr,vpsp,vxc,wvl%descr,wvl%den,xccc3d,xred)
1606     !if (wvlbigdft) energies_copy(energies,energies_wvl) ! TO BE ACTIVATED LATER
1607    end if
1608    call timab(60,2,tsec)
1609 
1610 !  ######################################################################
1611 !  In case of density mixing, check the exit criterion
1612 !  ----------------------------------------------------------------------
1613    if (dtset%iscf>=10.or.(wvlbigdft.and.dtset%iscf>0)) then
1614 !    Check exit criteria
1615      call timab(52,1,tsec)
1616      choice=2
1617      if(paw_dmft%use_dmft==1) then
1618        call prtene(dtset,energies,std_out,psps%usepaw)
1619      end if
1620      call scprqt(choice,cpus,deltae,diffor,dtset,&
1621 &     eigen,etotal,favg,fcart,energies%e_fermie,dtfil%fnameabo_app_eig,&
1622 &     dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,dtset%kptns,&
1623 &     maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,&
1624 &     occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,&
1625 &     psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,&
1626 &     electronpositron=electronpositron,fock=fock)
1627      call timab(52,2,tsec)
1628 
1629 !    Check if we need to exit the loop
1630      call timab(244,1,tsec)
1631      if (dtset%tfkinfunc>10.and.(.not.tfw_activated).and.quit==1) then
1632        quit=0;tfw_activated=.true.;reset_mixing=.true.
1633      end if
1634      if(dtset%userec==1.and.rec_set%quitrec==2)quit=1
1635      if (istep==nstep) quit=1
1636      quit_sum=quit
1637      call xmpi_sum(quit_sum,spaceComm,ierr)
1638      if (quit_sum>0) quit=1
1639      call timab(244,2,tsec)
1640 
1641 !    If criteria in scprqt say to quit, then exit the loop over istep.
1642      if (quit==1) exit
1643    end if
1644 
1645 !  ######################################################################
1646 !  Mix the total density (if required)
1647 !  ----------------------------------------------------------------------
1648    call timab(68,1,tsec)
1649 
1650    if (dtset%iscf>=10 .and.dtset%iscf/=22.and. .not. wvlbigdft ) then
1651 
1652 !    If LDA dielectric matrix is used for preconditionning, has to update here Kxc
1653      if (nkxc>0.and.modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79) &
1654 &     .and.((istep==1.or.istep==dielstrt).or.(dtset%iprcel>=100))) then
1655        optxc=10
1656        call xcdata_init(xcdata,dtset=dtset)
1657 !      to be adjusted for the call to rhotoxc
1658        nk3xc=1
1659        if(dtset%icoulomb==0 .and. dtset%usewvl==0) then
1660          call rhotoxc(edum,kxc,mpi_enreg,nfftf,&
1661 &         ngfftf,nhat,psps%usepaw,nhatgr,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
1662 &         optxc,dtset%paral_kgb,rhor,rprimd,dummy2,0,vxc,vxcavg_dum,xccc3d,xcdata,&
1663 &         add_tfw=tfw_activated,taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau)
1664        else if(.not. wvlbigdft) then
1665 !        WVL case:
1666          call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
1667 &         dtset%icoulomb, dtset%ixc, &
1668 &         mpi_enreg, nfftf, ngfftf,&
1669 &         nhat,psps%usepaw,&
1670 &         dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, &
1671 &         usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, vxcavg,&
1672 &         wvl%descr,wvl%den,&
1673 &         wvl%e,xccc3d,dtset%xclevel,dtset%xc_denpos)
1674        end if
1675      end if
1676 
1677      call newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,&
1678 &     dtset,etotal,fcart,pawfgr%fintocoa,&
1679 &     gmet,grhf,gsqcut,initialized,ispmix,istep_mix,kg_diel,kxc,&
1680 &     mgfftf,mix,pawfgr%coatofin,moved_atm_inside,mpi_enreg,my_natom,nattyp,nfftf,&
1681 &     nfftmix,nfftmix_per_nfft,ngfftf,ngfftmix,nkxc,npawmix,npwdiel,nvresid,psps%ntypat,&
1682 &     n1xccc,pawrhoij,pawtab,ph1df,psps,rhog,rhor,&
1683 &     rprimd,susmat,psps%usepaw,vtrial,wvl%descr,wvl%den,xred)
1684    end if   ! iscf>=10
1685 
1686    call timab(68,2,tsec)
1687 
1688 !  ######################################################################
1689 !  Additional computation in case of an electric field or electric displacement field
1690 !  ----------------------------------------------------------------------
1691 
1692    call timab(239,1,tsec)
1693 
1694    call update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,&
1695 &   efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,&
1696 &   dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,ngfft,dtset%nkpt,npwarr,&
1697 &   dtset%ntypat,pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,&
1698 &   pwind,pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,&
1699 &   1,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr)
1700 
1701    call timab(239,2,tsec)
1702 
1703 !  ######################################################################
1704 !  Compute the new potential from the trial density
1705 !  ----------------------------------------------------------------------
1706 
1707    call timab(243,1,tsec)
1708    if(VERBOSE)then
1709      call wrtout(std_out,'*. Compute the new potential from the trial density',"COLL")
1710    end if
1711 
1712 !  Set XC computation flag
1713    optxc=1
1714    if (nkxc>0) then
1715 ! MJV 2017 May 25: you should not be able to get here with iscf < 0
1716      if (dtset%iscf<0) optxc=2
1717      if (modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79).and. &
1718 &     dtset%iscf<10.and. &
1719 &     (dtset%iprcel>=100.or.istep==1.or.istep==dielstrt)) optxc=2
1720      if (dtset%iscf>=10.and.dtset%densfor_pred/=0.and.abs(dtset%densfor_pred)/=5) optxc=2
1721      if (optxc==2.and.dtset%xclevel==2.and.nkxc==2*min(dtset%nspden,2)-1) optxc=12
1722    end if
1723 
1724    if (dtset%iscf/=22) then
1725 !    PAW: eventually recompute compensation density (and gradients)
1726      nhatgrdim=0
1727      if ( allocated(nhatgr) ) then
1728        ABI_DEALLOCATE(nhatgr)
1729      end if
1730      if (psps%usepaw==1) then
1731        ider=-1;if (dtset%iscf>=10.and.((dtset%xclevel==2.and.dtset%pawnhatxc>0).or.usexcnhat==0)) ider=0
1732        if (dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) ider=ider+2
1733        if (ipositron==1) ider=-1
1734        if (ider>0) then
1735          nhatgrdim=1
1736          ABI_ALLOCATE(nhatgr,(nfftf,dtset%nspden,3))
1737        else
1738          ABI_ALLOCATE(nhatgr,(0,0,0))
1739        end if
1740        if (ider>=0) then
1741          call timab(558,1,tsec)
1742          izero=0
1743          call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,nfftf,ngfftf,&
1744 &         nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,nhatgr,nhat,&
1745 &         pawrhoij,pawrhoij,pawtab,k0,rprimd,ucvol_local,dtset%usewvl,xred,&
1746 &         comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1747 &         comm_fft=spaceComm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
1748 &         distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
1749          call timab(558,2,tsec)
1750        end if
1751      else
1752        ABI_ALLOCATE(nhatgr,(0,0,0))
1753      end if
1754 
1755 !    Compute new potential from the trial density
1756 
1757      optene=2*optres;if(psps%usepaw==1) optene=2
1758 
1759      call rhotov(dtset,energies,gprimd,gsqcut,istep,kxc,mpi_enreg,nfftf,ngfftf, &
1760 &     nhat,nhatgr,nhatgrdim,nkxc,nvresid,n3xccc,optene,optres,optxc,&
1761 &     rhog,rhor,rprimd,strsxc,ucvol_local,psps%usepaw,usexcnhat,&
1762 &     vhartr,vnew_mean,vpsp,vres_mean,res2,vtrial,vxcavg,vxc,wvl,xccc3d,xred,&
1763 &     electronpositron=electronpositron,taug=taug,taur=taur,vxctau=vxctau,&
1764 &     vxc_hybcomp=vxc_hybcomp,add_tfw=tfw_activated)
1765 
1766    end if
1767 
1768    call timab(243,2,tsec)
1769    call timab(60,1,tsec)
1770 
1771 !  This is inside the loop, its not equivalent to the line 1821
1772    if(moved_atm_inside==1) xred_old(:,:)=xred(:,:)
1773 
1774    if (dtset%iscf<10) then
1775 
1776      if(VERBOSE)then
1777        call wrtout(std_out,'Check exit criteria in case of potential mixing',"COLL")
1778      end if
1779 
1780 !    If the potential mixing is required, compute the total energy here
1781 !    PAW: has to compute here spherical terms
1782      if (psps%usepaw==1) then
1783        nzlmopt=0;if (istep_mix==1.and.dtset%pawnzlm>0) nzlmopt=-1
1784        if (istep_mix>1) nzlmopt=dtset%pawnzlm
1785        call paw_an_reset_flags(paw_an) ! Force the recomputation of on-site potentials
1786        option=2
1787        call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,ipert,&
1788 &       dtset%ixc,my_natom,dtset%natom,dtset%nspden,&
1789 &       psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,&
1790 &       paw_ij,pawang,dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,&
1791 &       pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,&
1792 &       hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1793 &       electronpositron=electronpositron)
1794 
1795      end if
1796 
1797 !    Add the Fock contribution to E_xc and E_xcdc if required
1798      if (usefock==1) then
1799        energies%e_fockdc=two*energies%e_fock
1800      end if
1801 
1802      if (.not.wvlbigdft) then
1803        call etotfor(atindx1,deltae,diffor,dtefield,dtset,&
1804 &       elast,electronpositron,energies,&
1805 &       etotal,favg,fcart,fock,forold,fred,gmet,grchempottn,gresid,grewtn,grhf,grnl,grvdw,&
1806 &       grxc,gsqcut,indsym,kxc,maxfor,mgfftf,mpi_enreg,my_natom,&
1807 &       nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,dtset%ntypat,nvresid,n1xccc, &
1808 &       n3xccc,0,computed_forces,optres,pawang,pawfgrtab,pawrad,pawrhoij,&
1809 &       pawtab,ph1df,red_ptot,psps,rhog,rhor,rmet,rprimd,symrec,synlgr,ucvol,&
1810 &       psps%usepaw,vhartr,vpsp,vxc,wvl%descr,wvl%den,xccc3d,xred)
1811      end if
1812 
1813    end if
1814    call timab(60,2,tsec)
1815 
1816 !  ######################################################################
1817 !  Check exit criteria in case of potential mixing or direct minimization
1818 !  ----------------------------------------------------------------------
1819    if ((dtset%iscf<10.and.(.not.wvlbigdft)) .or. dtset%iscf == 0) then
1820 !    Check exit criteria
1821      call timab(52,1,tsec)
1822      choice=2
1823      call scprqt(choice,cpus,deltae,diffor,dtset,&
1824 &     eigen,etotal,favg,fcart,energies%e_fermie,dtfil%fnameabo_app_eig,&
1825 &     dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,dtset%kptns,&
1826 &     maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,&
1827 &     occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,&
1828 &     psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,&
1829 &     electronpositron=electronpositron,fock=fock)
1830      call timab(52,2,tsec)
1831 
1832 !    Check if we need to exit the loop
1833      call timab(244,1,tsec)
1834      if (dtset%tfkinfunc>10.and.(.not.tfw_activated).and.quit==1) then
1835        quit=0;tfw_activated=.true.;reset_mixing=.true.
1836      end if
1837      if (istep==nstep.and.psps%usepaw==1) quit=1
1838      if(dtset%userec==1 .and. rec_set%quitrec==2) quit=1
1839      quit_sum=quit
1840      call xmpi_sum(quit_sum,spaceComm,ierr)
1841      if (quit_sum > 0) quit=1
1842 
1843 !    If criteria in scprqt say to quit, then exit the loop over istep.
1844      if (quit==1) then
1845        do ispden=1,dtset%nspden
1846          vtrial(:,ispden)=vtrial(:,ispden)+nvresid(:,ispden)+vres_mean(ispden)
1847        end do
1848        call timab(244,2,tsec) ! Due to the exit instruction, two timab calls are needed
1849        exit ! exit the loop over istep
1850      end if
1851      call timab(244,2,tsec) ! Due to the exit instruction, two timab calls are needed
1852    end if
1853 
1854 !  ######################################################################
1855 !  Mix the potential (if required) - Check exit criteria
1856 !  ----------------------------------------------------------------------
1857 
1858    call timab(245,1,tsec)
1859    if (dtset%iscf<10 .and. dtset%iscf>0 .and. .not. wvlbigdft) then
1860 
1861      if(VERBOSE)then
1862        call wrtout(std_out,'*. Mix the potential (if required) - Check exit criteria',"COLL")
1863      end if
1864 
1865 !    Precondition the residual and forces, then determine the new vtrial
1866 !    (Warning: the (H)xc potential may have been subtracted from vtrial)
1867 
1868      call newvtr(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,&
1869 &     dtn_pc,dtset,etotal,fcart,pawfgr%fintocoa,&
1870 &     gmet,grhf,gsqcut,initialized,ispmix,&
1871 &     istep_mix,kg_diel,kxc,mgfftf,mix,pawfgr%coatofin,&
1872 &     moved_atm_inside,mpi_enreg,my_natom,nattyp,nfftf,nfftmix,&
1873 &     ngfftf,ngfftmix,nkxc,npawmix,npwdiel,&
1874 &     nstep,psps%ntypat,n1xccc,&
1875 &     pawrhoij,ph1df,psps,rhor,rprimd,susmat,psps%usepaw,&
1876 &     vhartr,vnew_mean,vpsp,nvresid,vtrial,vxc,xred,&
1877 &     nfftf,&
1878 &     pawtab,rhog,wvl)
1879 
1880    end if   ! iscf<10
1881 
1882 !  ######################################################################
1883 !  END MINIMIZATION ITERATIONS
1884 !  ######################################################################
1885 
1886    if(VERBOSE)then
1887      call wrtout(std_out,'*. END MINIMIZATION ITERATIONS',"COLL")
1888    end if
1889 
1890 !  The initialisation of the gstate run should be done when this point is reached
1891    initialized=1
1892 
1893    !if (dtset%useria == -4242) then
1894    !  call gshgg_mkncwrite(istep, dtset, dtfil, psps, hdr, pawtab, pawfgr, paw_ij, mpi_enreg, &
1895    !     rprimd, xred, eigen, npwarr, kg, ylm, ngfft, dtset%nfft, ngfftf, nfftf, vtrial) !electronpositron) ! Optional arguments
1896    !end if
1897 
1898 !  This is to save the density for restart.
1899    if (iwrite_fftdatar(mpi_enreg)) then
1900 
1901      if(dtset%prtden<0.or.dtset%prtkden<0) then
1902 !      Update the content of the header (evolving variables)
1903 !      Don't use parallelism over atoms because only me=0 accesses here
1904        bantot=hdr%bantot
1905        if (dtset%positron==0) then
1906          call hdr_update(hdr,bantot,etotal,energies%e_fermie,residm,&
1907 &         rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1))
1908        else
1909          call hdr_update(hdr,bantot,electronpositron%e0,energies%e_fermie,residm,&
1910 &         rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1))
1911        end if
1912      end if
1913 
1914      if (dtset%prtden<0) then
1915        if (mod(istep-1,abs(dtset%prtden))==0) then
1916          isave_den=isave_den+1
1917          rdwrpaw=0
1918          call int2char4(mod(isave_den,2),tag)
1919          ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
1920          fildata=trim(dtfil%fnametmp_app_den)//'_'//trim(tag)
1921          if (dtset%iomode == IO_MODE_ETSF) fildata = nctk_ncify(fildata)
1922          call fftdatar_write_from_hdr("density",fildata,dtset%iomode,hdr,ngfftf,cplex1,nfftf,&
1923          dtset%nspden,rhor,mpi_enreg,eigen=eigen)
1924        end if
1925      end if
1926 
1927      if (dtset%prtkden<0) then
1928        if (mod(istep-1,abs(dtset%prtkden))==0) then
1929          isave_kden=isave_kden+1
1930          rdwrpaw=0
1931          call int2char4(mod(isave_kden,2),tag)
1932          ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
1933          fildata=trim(dtfil%fnametmp_app_kden)//'_'//trim(tag)
1934          if (dtset%iomode == IO_MODE_ETSF) fildata = nctk_ncify(fildata)
1935          ! output the Laplacian of density
1936          call fftdatar_write_from_hdr("kinedr",fildata,dtset%iomode,hdr,ngfftf,cplex1,nfftf,&
1937          dtset%nspden,taur,mpi_enreg,eigen=eigen)
1938        end if
1939      end if
1940 
1941    end if
1942 
1943    ABI_DEALLOCATE(nhatgr)
1944 
1945    istep_mix=istep_mix+1
1946    if (reset_mixing) then
1947      istep_mix=1;reset_mixing=.false.
1948    end if
1949    if (ipositron/=0) electronpositron%istep_scf=electronpositron%istep_scf+1
1950 
1951    call timab(245,2,tsec)
1952 
1953  end do ! istep
1954 
1955  ! Avoid pending requests if itime == ntime.
1956  call xmpi_wait(quitsum_request,ierr)
1957  if (timelimit_exit == 1) istep = istep - 1
1958 
1959  call timab(246,1,tsec)
1960 
1961  if (dtset%iscf > 0) then
1962    call ab7_mixing_deallocate(mix)
1963  end if
1964 
1965  if (usefock==1)then
1966    if(wfmixalg/=0)then
1967      call scf_history_free(scf_history_wf)
1968    end if
1969  end if
1970 
1971 
1972  if (quit==1.and.nstep==1) initialized=1
1973 
1974 !######################################################################
1975 !Case nstep==0: compute energy based on incoming wf
1976 !----------------------------------------------------------------------
1977 
1978  if(nstep==0) then
1979    optene=2*psps%usepaw+optres
1980    energies%entropy=results_gs%energies%entropy  !MT20070219: entropy is not recomputed in routine energy
1981    if (.not.allocated(nhatgr) ) then
1982      ABI_ALLOCATE(nhatgr,(0,0,0))
1983    end if
1984    call energy(cg,compch_fft,dtset,electronpositron,&
1985 &   energies,eigen,etotal,gsqcut,indsym,irrzon,kg,mcg,mpi_enreg,my_natom,&
1986 &   nfftf,ngfftf,nhat,nhatgr,nhatgrdim,npwarr,n3xccc,&
1987 &   occ,optene,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrhoij,pawtab,&
1988 &   phnons,ph1d,psps,resid,rhog,rhor,rprimd,strsxc,symrec,taug,taur,usexcnhat,&
1989 &   vhartr,vtrial,vpsp,vxc,vxctau,wvl%wfs,wvl%descr,wvl%den,wvl%e,xccc3d,xred,ylm,&
1990 &   add_tfw=tfw_activated)
1991    if (nhatgrdim>0)  then
1992      ABI_DEALLOCATE(nhatgr)
1993    end if
1994 
1995  end if ! nstep==0
1996 
1997 !######################################################################
1998 !Additional steps after SC iterations, including force, stress, polarization calculation
1999 !----------------------------------------------------------------------
2000 
2001  if (dtset%userec==1) then
2002    call prtene(dtset,energies,ab_out,psps%usepaw)
2003    call prtene(dtset,energies,std_out,psps%usepaw)
2004  end if
2005 
2006 !if (wvlbigdft) call energies_copy(energies_wvl,energies) ! TO BE ACTIVATED LATER
2007 
2008 !PAW: if cprj=<p_lmn|Cnk> are in memory,
2009 !need to reorder them (from atom-sorted to unsorted)
2010  if (psps%usepaw==1.and.usecprj==1) then
2011    iorder_cprj=1
2012    call pawcprj_reorder(cprj,atindx1)
2013    if (dtset%positron/=0) then
2014      if (electronpositron%dimcprj>0) then
2015        call pawcprj_reorder(electronpositron%cprj_ep,atindx1)
2016      end if
2017    end if
2018    if (dtset%usewvl==1) then
2019      call wvl_cprjreorder(wvl%descr,atindx1)
2020    end if
2021  end if
2022 
2023 !PAW: if cprj=<p_lmn|Cnk> are not in memory,need to compute them in some cases
2024  recompute_cprj = psps%usepaw ==1 .and. usecprj==0 .and. &
2025 & (dtset%prtwant  ==2 .or. &
2026 & dtset%prtwant  ==3  .or. &
2027 & dtset%prtnabla > 0  .or. &
2028 & dtset%prtdos   ==3  .or. &
2029 & dtset%berryopt /=0  .or. &
2030 & dtset%orbmag /=0    .or. &
2031 & dtset%kssform  ==3  .or. &
2032 & dtset%pawfatbnd> 0  .or. &
2033 & dtset%pawprtwf > 0  )
2034 
2035  if (recompute_cprj) then
2036    usecprj=1
2037    mband_cprj=dtset%mband/mpi_enreg%nproc_band
2038    mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
2039    call pawcprj_free(cprj)
2040    ABI_DATATYPE_DEALLOCATE(cprj) ! Was previously allocated (event if size = 0,0)
2041    ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj))
2042    ncpgr = 0 ; ctocprj_choice = 1
2043    if (finite_efield_flag) then
2044      if (forces_needed /= 0 .and. stress_needed == 0) then
2045        ncpgr = 3 ; ctocprj_choice = 2
2046      else if (forces_needed /= 0 .and. stress_needed /= 0) then
2047        ncpgr = 9 ; ctocprj_choice = 23
2048      else if (forces_needed == 0 .and. stress_needed /= 0) then
2049        ncpgr = 6 ; ctocprj_choice = 3
2050      end if
2051    end if
2052    call pawcprj_alloc(cprj,ncpgr,dimcprj)
2053    iatom=0 ; iorder_cprj=1 ! cprj are not ordered
2054    call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,&
2055 &   iatom,idir,iorder_cprj,dtset%istwfk,kg,dtset%kptns,&
2056 &   mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,&
2057 &   dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,&
2058 &   dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,&
2059 &   dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,&
2060 &   ucvol,dtfil%unpaw,xred,ylm,ylmgr)
2061  end if
2062 
2063  call timab(246,2,tsec)
2064  call timab(247,1,tsec)
2065 
2066 !SHOULD CLEAN THE ARGS OF THIS ROUTINE
2067  call afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,&
2068 & deltae,diffor,dtefield,dtfil,dtorbmag,dtset,eigen,electronpositron,elfr,&
2069 & energies,etotal,favg,fcart,fock,forold,fred,grchempottn,&
2070 & gresid,grewtn,grhf,grhor,grvdw,&
2071 & grxc,gsqcut,hdr,indsym,irrzon,istep,kg,kxc,lrhor,maxfor,mcg,mcprj,mgfftf,&
2072 & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,&
2073 & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,&
2074 & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,&
2075 & prtxml,psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,&
2076 & rhog,rhor,rprimd,stress_needed,strsxc,strten,symrec,synlgr,taug,&
2077 & taur,tollist,usecprj,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,&
2078 & xccc3d,xred,ylm,ylmgr,dtset%charge*SUM(vpotzero(:)),conv_retcode)
2079 
2080 !Before leaving the present routine, save the current value of xred.
2081  xred_old(:,:)=xred(:,:)
2082 
2083  call timab(247,2,tsec)
2084 
2085 !######################################################################
2086 !All calculations in scfcv are finished. Printing section
2087 !----------------------------------------------------------------------
2088 
2089  call timab(248,1,tsec)
2090 
2091  call outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,&
2092 & dtset,ecut,eigen,electronpositron,elfr,etotal,&
2093 & gmet,gprimd,grhor,hdr,kg,lrhor,dtset%mband,mcg,mcprj,dtset%mgfft,&
2094 & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,my_natom,dtset%natom,nattyp,&
2095 & nfftf,ngfftf,nhat,dtset%nkpt,npwarr,dtset%nspden,&
2096 & dtset%nsppol,dtset%nsym,psps%ntypat,n3xccc,occ,paw_dmft,pawang,pawfgr,pawfgrtab,&
2097 & pawrad,pawrhoij,pawtab,paw_an,paw_ij,dtset%prtvol,psps,results_gs,&
2098 & rhor,rprimd,taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl%den,xccc3d,xred)
2099 
2100  if(associated(elfr))then
2101    ABI_DEALLOCATE(elfr)
2102    nullify(elfr)
2103  end if
2104 
2105  if(associated(grhor))then
2106    ABI_DEALLOCATE(grhor)
2107    nullify(grhor)
2108  end if
2109 
2110  if(associated(lrhor))then
2111    ABI_DEALLOCATE(lrhor)
2112    nullify(lrhor)
2113  end if
2114 
2115  if(dtset%prtkden/=0 .or. dtset%prtelf/=0)then
2116    ABI_DEALLOCATE(taur)
2117  end if
2118 
2119  call timab(248,2,tsec)
2120  call timab(249,1,tsec)
2121 
2122 !Transfer eigenvalues and occupation computed by BigDFT in afterscfloop to eigen.
2123 #if defined HAVE_BIGDFT
2124  if (dtset%usewvl == 1) then
2125    if (dtset%nsppol == 1) then
2126      eigen = wvl%wfs%ks%orbs%eval
2127      occ = wvl%wfs%ks%orbs%occup
2128    else
2129      eigen(1:wvl%wfs%ks%orbs%norbu) = wvl%wfs%ks%orbs%eval(1:wvl%wfs%ks%orbs%norbu)
2130      eigen(dtset%mband + 1:dtset%mband + wvl%wfs%ks%orbs%norbd) = &
2131 &     wvl%wfs%ks%orbs%eval(wvl%wfs%ks%orbs%norbu + 1:wvl%wfs%ks%orbs%norb)
2132      occ(1:wvl%wfs%ks%orbs%norbu) = wvl%wfs%ks%orbs%occup(1:wvl%wfs%ks%orbs%norbu)
2133      occ(dtset%mband + 1:dtset%mband + wvl%wfs%ks%orbs%norbd) = &
2134 &     wvl%wfs%ks%orbs%occup(wvl%wfs%ks%orbs%norbu + 1:wvl%wfs%ks%orbs%norb)
2135    end if
2136  end if
2137 #endif
2138 
2139 !######################################################################
2140 !Deallocate memory and save results
2141 !----------------------------------------------------------------------
2142 
2143  call prc_mem_free()
2144 
2145  ABI_DEALLOCATE(fcart)
2146  ABI_DEALLOCATE(fred)
2147  ABI_DEALLOCATE(forold)
2148  ABI_DEALLOCATE(grchempottn)
2149  ABI_DEALLOCATE(gresid)
2150  ABI_DEALLOCATE(grewtn)
2151  ABI_DEALLOCATE(grnl)
2152  ABI_DEALLOCATE(grvdw)
2153  ABI_DEALLOCATE(grxc)
2154  ABI_DEALLOCATE(synlgr)
2155  ABI_DEALLOCATE(ph1d)
2156  ABI_DEALLOCATE(ph1df)
2157  ABI_DEALLOCATE(vhartr)
2158  ABI_DEALLOCATE(vtrial)
2159  ABI_DEALLOCATE(vpsp)
2160  ABI_DEALLOCATE(vxc)
2161  ABI_DEALLOCATE(vxc_hybcomp)
2162  ABI_DEALLOCATE(vxctau)
2163  ABI_DEALLOCATE(xccc3d)
2164  ABI_DEALLOCATE(kxc)
2165  ABI_DEALLOCATE(shiftvector)
2166  ABI_DEALLOCATE(dtn_pc)
2167  ABI_DEALLOCATE(grhf)
2168  ABI_DEALLOCATE(nvresid)
2169 
2170  if((nstep>0.and.dtset%iscf>0).or.dtset%iscf==-1) then
2171    ABI_DEALLOCATE(dielinv)
2172  end if
2173  ABI_DEALLOCATE(gbound_diel)
2174  ABI_DEALLOCATE(irrzondiel)
2175  ABI_DEALLOCATE(kg_diel)
2176  ABI_DEALLOCATE(phnonsdiel)
2177  ABI_DEALLOCATE(susmat)
2178  ABI_DEALLOCATE(ph1ddiel)
2179  ABI_DEALLOCATE(ylmdiel)
2180  !end if
2181 
2182  if (psps%usepaw==1) then
2183    if (dtset%iscf>0) then
2184      do iatom=1,my_natom
2185        pawrhoij(iatom)%lmnmix_sz=0
2186        pawrhoij(iatom)%use_rhoijres=0
2187        ABI_DEALLOCATE(pawrhoij(iatom)%kpawmix)
2188        ABI_DEALLOCATE(pawrhoij(iatom)%rhoijres)
2189      end do
2190    end if
2191    if (recompute_cprj.or.usecprj==1) then
2192      usecprj=0;mcprj=0
2193      if (recompute_cprj.or.dtset%mkmem/=0) then
2194        call pawcprj_free(cprj)
2195      end if
2196    end if
2197    call paw_an_free(paw_an)
2198    call paw_ij_free(paw_ij)
2199    call pawfgrtab_free(pawfgrtab)
2200    if(dtset%usewvl==1) then
2201 #if defined HAVE_BIGDFT
2202      call cprj_clean(wvl%descr%paw%cprj)
2203      ABI_DATATYPE_DEALLOCATE(wvl%descr%paw%cprj)
2204 #endif
2205      call paw2wvl_ij(2,paw_ij,wvl%descr)
2206    end if
2207  end if
2208  ABI_DATATYPE_DEALLOCATE(pawfgrtab)
2209  ABI_DATATYPE_DEALLOCATE(paw_an)
2210  ABI_DATATYPE_DEALLOCATE(paw_ij)
2211  ABI_DEALLOCATE(nhat)
2212  ABI_DEALLOCATE(dimcprj_srt)
2213  ABI_DEALLOCATE(dimcprj)
2214  call pawcprj_free(cprj)
2215  ABI_DATATYPE_DEALLOCATE(cprj)
2216 
2217 ! Deallocate exact exchange data at the end of the calculation
2218  if (usefock==1) then
2219    if (fock%fock_common%use_ACE/=0) then
2220      call fock_ACE_destroy(fock%fockACE)
2221    end if
2222    call fock_common_destroy(fock%fock_common)
2223    call fock_BZ_destroy(fock%fock_BZ)
2224    call fock_destroy(fock)
2225    nullify(fock)
2226  end if
2227 
2228  if (prtxml == 1) then
2229 !  We output the final result given in results_gs
2230    write(ab_xml_out, "(A)") '      <finalConditions>'
2231    call out_resultsgs_XML(dtset, 4, results_gs, psps%usepaw)
2232    write(ab_xml_out, "(A)") '      </finalConditions>'
2233    write(ab_xml_out, "(A)") '    </scfcvLoop>'
2234  end if
2235 
2236  call status(0,dtfil%filstat,iexit,level,'exit')
2237 
2238  call timab(249,2,tsec)
2239  call timab(238,2,tsec)
2240 
2241  DBG_EXIT("COLL")
2242 
2243 end subroutine scfcv

TABLE OF CONTENTS

ABINIT/scfcv [ Functions ]