TABLE OF CONTENTS


ABINIT/afterscfloop [ Functions ]

[ Top ] [ Functions ]

NAME

 afterscfloop

FUNCTION

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

INPUTS

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

OUTPUT

  conv_retcode= Non-zero if convergence is not achieved.
  elfr(nfftf,nspden)=electron localization function
  grhor(nfft,nspden,3)= gradient of electron density in electrons/bohr**4, real space
  lrhor(nfft,nspden)= Laplacian of electron density in electrons/bohr**5, real space
  results_gs <type(results_gs_type)>=results (energy and its components,
   forces and its components, the stress tensor) of a ground-state computation
   (should be made a pure output quantity)
  taug(2,nfftf)=Fourier transform of total kinetic energy density
  taur(nfftf,nspden)=total kinetic energy density in real space
  ==== if forces are required ====
   diffor=maximal absolute value of changes in the components of
          force between the input and the output.
   favg(3)=mean of the forces before correction for translational symmetry
   fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
     at input, previous value of forces,
     at output, new value.
     Note : unlike fred, this array has been corrected by enforcing
     the translational symmetry, namely that the sum of force
     on all atoms is zero.
   fred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
   gresid(3,natom)=forces due to the residual of the potential
   grhf(3,natom)=Hellman-Feynman derivatives of the total energy
   grxc(9+3*natom)=d(Exc)/d(xred) if core charges are used
   maxfor=maximal absolute value of the output array force.
   synlgr(3,natom)=symmetrized gradients of energy due to nonlocal contributions
  ==== if stress tensor is required ====
   strten(6)=components of the stress tensor (hartree/bohr^3) for the
    6 unique components of this symmetric 3x3 tensor:
    Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).

SIDE EFFECTS

 computed_forces=1 if forces have been computed, 0 otherwise
 dtefield <type(efield_type)> = variables related to Berry phase
       and electric field calculations (see initberry.f).
       In case dtset%berryopt = 4/6/7/14/16/17, the overlap matrices computed
       in this routine are stored in dtefield%smat in order
       to be used in the electric field calculation.
 dtorbmag <type(orbmag_type)> = variables related to orbital magnetization
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  energies <type(energies_type)>=all part of total energy.
   | entropy(IN)=entropy due to the occupation number smearing (if metal)
   | e_localpsp(IN)=local psp energy (hartree)
   | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree)
   | e_chempot(IN)=energy from spatially varying chemical potential (hartree)
   | e_ewald(IN)=Ewald energy (hartree)
   | e_vdw_dftd(IN)=VdW DFT-D energy
   | e_hartree(IN)=Hartree part of total energy (hartree units)
   | e_corepsp(IN)=psp core-core energy
   | e_hybcomp_E0=energy compensation term for hybrid exchange-correlation energy (hartree) at fixed density
   | e_hybcomp_v0=potential compensation term for hybrid exchange-correlation energy (hartree) at fixed density
   | e_hybcomp_v=potential compensation term for hybrid exchange-correlation energy (hartree) at self-consistent den
   | e_kinetic(IN)=kinetic energy part of total energy.
   | e_nonlocalpsp(IN)=nonlocal pseudopotential part of total energy.
   | e_xc(IN)=exchange-correlation energy (hartree)
   | e_xcdc(IN)=exchange-correlation double-counting energy (hartree)
   | e_paw(IN)=PAW spherical part energy
   | e_pawdc(IN)=PAW spherical part double-counting energy
   | e_elecfield(OUT)=the term of the energy functional that depends explicitely    !!HONG
   |                  on the electric field:
   |                 enefield =  -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j for fixed E/ebar
   |                          =  Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j  for fixed D/d
 etotal=total energy, might be correct by improved polarization computation
 forold(3,natom)=old forces
 xred(3,natom)=reduced dimensionless atomic coordinates
 ===== if dtset%densfor_pred==3 .and. moved_atm_inside==1 =====
   ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases (coarse grid)
   ph1df(2,3*(2*mgfftf+1)*natom)=1-dim structure factor phases (fine PAW grid)
  wvl <type(wvl_data)>=all wavelets data.

NOTES

PARENTS

      scfcv

CHILDREN

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

SOURCE

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

ABINIT/m_afterscfloop [ Modules ]

[ Top ] [ Modules ]

NAME

  m_afterscfloop

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_afterscfloop
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_energies
34  use m_errors
35  use m_abicore
36  use m_efield
37  use m_ab7_mixing
38  use m_hdr
39 
40  use m_time,             only : timab
41  use m_xmpi,             only : xmpi_sum, xmpi_comm_rank,xmpi_comm_size
42  use m_geometry,         only : xred2xcart, metric
43  use m_crystal,          only : prtposcar
44  use m_results_gs ,      only : results_gs_type
45  use m_electronpositron, only : electronpositron_type, electronpositron_calctype, exchange_electronpositron
46  use m_dtset,            only : dtset_copy, dtset_free
47  use m_paw_dmft,         only : paw_dmft_type
48  use m_pawang,           only : pawang_type
49  use m_pawrad,           only : pawrad_type
50  use m_pawtab,           only : pawtab_type
51  use m_pawrhoij,         only : pawrhoij_type
52  use m_paw_an,           only : paw_an_type
53  use m_paw_ij,           only : paw_ij_type
54  use m_pawfgrtab,        only : pawfgrtab_type
55  use m_pawcprj,          only : pawcprj_type,pawcprj_getdim
56  use m_pawfgr,           only : pawfgr_type
57  use m_paw_mkrho,        only : pawmkrho
58  use m_paw_nhat,         only : nhatgrid,wvl_nhatgrid
59  use m_paw_occupancies,  only : pawmkrhoij
60  use m_paw_correlations, only : setnoccmmp
61  use m_orbmag,           only : chern_number,orbmag,orbmag_type
62  use m_fock,             only : fock_type
63  use m_kg,               only : getph
64  use m_spin_current,     only : spin_current
65  use m_mkrho,            only : mkrho, prtrhomxmn
66  use m_elpolariz,        only : elpolariz
67  use m_nonlop_test,      only : nonlop_test
68  use m_common,           only : scprqt
69  use m_xctk,             only : xcden
70  use m_forstr,           only : forstr
71  use m_wvl_rho,          only : wvl_mkrho
72  use m_wvl_psi,          only : wvl_psitohpsi, wvl_tail_corrections
73 
74 #ifdef HAVE_BIGDFT
75  use m_abi2big
76  use BigDFT_API, only : last_orthon, &
77       & kswfn_free_scf_data, denspot_free_history,&
78       & write_energies, total_energies, XC_potential,&
79       & eigensystem_info, applyprojectorsonthefly
80 #endif
81 
82  implicit none
83 
84  private