TABLE OF CONTENTS


ABINIT/e_eigen [ Functions ]

[ Top ] [ Functions ]

NAME

  e_eigen

FUNCTION

  Computes eigenvalues energy from eigen, occ, kpt, wtk

COPYRIGHT

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

INPUTS

  eigen(nkpt*nsppol)=eigenvalues
  mband= maximum number of bands
  nband(nkpt*nsppol)= number of bands for each k-point and spin
  nkpt= number of k-points
  nsppol= number of spin polarization
  occ(mband*nkpt*nsppol)=occupations
  wtk(nkpt)= k-point weights

OUTPUT

  e_eigenvalues= eigenvalues energy

SIDE EFFECTS

NOTES

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2053 subroutine e_eigen(eigen,e_eigenvalues,mband,nband,nkpt,nsppol,occ,wtk)
2054 
2055  use defs_basis
2056  use m_errors
2057 
2058 !This section has been created automatically by the script Abilint (TD).
2059 !Do not modify the following lines by hand.
2060 #undef ABI_FUNC
2061 #define ABI_FUNC 'e_eigen'
2062 !End of the abilint section
2063 
2064  implicit none
2065 
2066 !Arguments ------------------------------------
2067  integer , intent(in)  :: mband,nkpt,nsppol
2068  integer , intent(in)  :: nband(nkpt*nsppol)
2069  real(dp) , intent(in)  :: eigen(mband*nkpt*nsppol)
2070  real(dp) , intent(in)  :: occ(mband*nkpt*nsppol)
2071  real(dp) , intent(in)  :: wtk(nkpt)
2072  real(dp) , intent(out) :: e_eigenvalues
2073 
2074 !Local variables-------------------------------
2075  integer :: ib,iband,ii,ikpt,isppol,nband_k
2076  real(dp) :: wtk_k
2077 
2078 ! *************************************************************************
2079 
2080    DBG_ENTER("COLL")
2081    ii=0;ib=0
2082    do isppol=1,nsppol
2083      do ikpt=1,nkpt
2084        ii=ii+1
2085        nband_k=nband(ii) ;  wtk_k=wtk(ii)
2086        do iband=1,nband_k
2087          ib=ib+1
2088          if(abs(occ(ib)) > tol8) then
2089            e_eigenvalues = e_eigenvalues + wtk_k*occ(ib)*eigen(ib)
2090          end if
2091        end do
2092      end do
2093    end do
2094 
2095    DBG_EXIT("COLL")
2096 
2097  end subroutine e_eigen

ABINIT/vtorho [ Functions ]

[ Top ] [ Functions ]

NAME

 vtorho

FUNCTION

 This routine compute the new density from a fixed potential (vtrial)
 but might also simply compute eigenvectors and eigenvalues.
 The main part of it is a wf update over all k points.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MF, AR, MM, MT, FJ, MB, MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  afford=used to dimension susmat
  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
  dbl_nnsclo=if 1, will double the value of dtset%nnsclo
  dielop= if positive, the dielectric matrix must be computed.
  dielstrt=number of the step at which the dielectric preconditioning begins.
  dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only)
  dtefield <type(efield_type)> = variables related to Berry phase
      calculations (see initberry.f)
  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
   | mkmem =number of k points treated by this node.
   | mpw=maximum dimensioned size of npw
   | nfft=(effective) number of FFT grid points (for this processor)
   | 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
   | typat= array of types of the natoms
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  etotal=total energy (Ha) - only needed for tddft
  fock <type(fock_type)>= quantities to calculate Fock exact exchange
  gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for the dielectric matrix
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
   (3x3 tensor) and grads wrt atomic coordinates (3*natom)
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  hdr <type(hdr_type)>=the header of wf, den and pot files
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for diel matrix
                                     nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
  istep=index of the number of steps in the routine scfcv
  istep_mix=index of the number of steps for the SCF mixing (can be <istep)
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
  kxc(nfftf,nkxc)=exchange-correlation kernel, needed only if nkxc/=0 .
  lmax_diel=1+max. value of l angular momentum used for dielectric matrix
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
  mpi_enreg=informations about MPI parallelization
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nattyp(ntypat)= # atoms of each type.
  nfftf= -PAW ONLY- number of FFT grid points for the fine grid
         (nfftf=nfft for norm-conserving potential runs)
  nfftdiel=number of fft grid points for the computation of the diel matrix
  ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix,
                see ~abinit/doc/variables/vargs.htm#ngfft
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  npwarr(nkpt)=number of planewaves in basis at this k point
  npwdiel=size of the susmat array.
  ntypat=number of types of atoms in unit cell.
  optforces=option for the computation of forces (0: no force;1: forces)
  optres=0: the new value of the density is computed in place of the input value
         1: only the density residual is computed ; the input density is kept
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawang <type(pawang)>=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
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
                                    nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
  phnonsdiel(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases,
   for diel matr
                                     nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  ph1ddiel(2,3*(2*mgfftdiel+1)*natom*usepaw)=one-dimensional structure factor information
                                             for the dielectric matrix
  psps <type(pseudopotential_type)>=variables related to 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)
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional primitive vectors
  symrec(3,3,nsym)=symmetry operations in reciprocal space
  ucvol=unit cell volume in bohr**3.
  usecprj=1 if cprj datastructure is stored in memory
  wffnew,unit numbers for wf disk files.
  vtrial(nfftf,nspden)=INPUT potential Vtrial(r).
  vxctau=(only for meta-GGA): derivative of XC energy density with respect to
    kinetic energy density (depsxcdtau). The arrays vxctau(nfft,nspden,4) contains also
    the gradient of vxctau (gvxctau) in vxctau(:,:,2:4)
  xred(3,natom)=reduced dimensionless atomic coordinates
  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
  ylmdiel(npwdiel,lmax_diel**2)= real spherical harmonics for each G and k point
                                 for the dielectric matrix

OUTPUT

  compch_fft=-PAW only- compensation charge inside spheres computed over fine fft grid
  dphase(3) : dphase(idir) = accumulated change in the string-averaged
     Zak phase along the idir-th direction caused by the update of all
     the occupied Bloch states at all the k-points (only if finite electric field)
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  resid(mband*nkpt*nsppol)=residuals for each band over all k points.
  residm=maximum value from resid array (except for nbdbuf highest bands)
  susmat(2,npwdiel*afford,nspden,npwdiel,nspden)=
   the susceptibility (or density-density response) matrix in reciprocal space
  === if optforces>0 ===
    grnl(3*natom)=stores grads of nonlocal energy wrt length scales
  ==== if optres==1
    nres2=square of the norm of the residual
    nvresid(nfftf,nspden)=density residual
  ==== if psps%usepaw==1
    cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
                               cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
    nhat(nfftf,nspden*psps%usepaw)=compensation charge density on rectangular grid in real space

SIDE EFFECTS

  cg(2,mpw*dtset%nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions.
   At output contains updated wavefunctions coefficients;
    if nkpt>1, these are kept in a disk file.
  energies <type(energies_type)>=storage for energies computed here :
   | e_eigenvalues=Sum of the eigenvalues - Band energy (Hartree)
   | e_kinetic=kinetic energy part of total energy
   | e_nonlocalpsp=nonlocal pseudopotential part of total energy
   | e_fermie=fermi energy (Hartree)
  occ(mband*nkpt*nsppol)=occupation number for each band for each k.
      (input if insulator - occopt<3 - ; output if metallic)
  pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  rhog(2,nfftf)=Fourier transform of total electron density
  rhor(nfftf,nspden)=total electron density (el/bohr**3)
  taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
  taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
  wvl <type(wvl_data)>=wavelets structures in case of wavelets basis.
  ==== if (usepaw==1) ====
    cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
                               cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.

PARENTS

      scfcv

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

NOTES

  Be careful to the meaning of nfft (size of FFT grids):
   - In case of norm-conserving calculations the FFT grid is the usual FFT grid.
   - In case of PAW calculations:
     Two FFT grids are used; one with nfft points (coarse grid) for
     the computation of wave functions ; one with nfftf points
     (fine grid) for the computation of total density.

  The total electronic density (rhor,rhog) is divided into two terms:
   - The density related to WFs =Sum[Psi**2]
   - The compensation density (nhat) - only in PAW

  The parallelisation needed for the electric field should be
  made an independent subroutine, so that this routine could be put
  back in the 95_drive directory.

SOURCE

 180 #if defined HAVE_CONFIG_H
 181 #include "config.h"
 182 #endif
 183 
 184 #include "abi_common.h"
 185 
 186 subroutine vtorho(afford,atindx,atindx1,cg,compch_fft,cprj,cpus,dbl_nnsclo,&
 187 &           dielop,dielstrt,dmatpawu,dphase,dtefield,dtfil,dtset,&
 188 &           eigen,electronpositron,energies,etotal,gbound_diel,&
 189 &           gmet,gprimd,grnl,gsqcut,hdr,indsym,irrzon,irrzondiel,&
 190 &           istep,istep_mix,kg,kg_diel,kxc,lmax_diel,mcg,mcprj,mgfftdiel,mpi_enreg,&
 191 &           my_natom,natom,nattyp,nfftf,nfftdiel,ngfftdiel,nhat,nkxc,&
 192 &           npwarr,npwdiel,nres2,ntypat,nvresid,occ,optforces,&
 193 &           optres,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrhoij,pawtab,&
 194 &           phnons,phnonsdiel,ph1d,ph1ddiel,psps,fock,&
 195 &           pwind,pwind_alloc,pwnsfac,resid,residm,rhog,rhor,&
 196 &           rmet,rprimd,susmat,symrec,taug,taur,&
 197 &           ucvol,usecprj,wffnew,vtrial,vxctau,wvl,xred,ylm,ylmgr,ylmdiel)
 198 
 199  use defs_basis
 200  use defs_datatypes
 201  use defs_abitypes
 202  use defs_wvltypes
 203  use m_profiling_abi
 204  use m_xmpi
 205  use m_ab7_mixing
 206  use m_errors
 207  use m_wffile
 208  use m_efield
 209  use m_cgtools
 210  use m_gemm_nonlop
 211 
 212  use m_geometry,           only : xred2xcart
 213  use m_occ,                only : newocc
 214  use m_dtset,              only : testsusmat
 215  use m_pawang,             only : pawang_type
 216  use m_pawtab,             only : pawtab_type
 217  use m_paw_ij,             only : paw_ij_type
 218  use m_pawfgrtab,          only : pawfgrtab_type
 219  use m_pawrhoij,           only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_io, pawrhoij_get_nspden
 220  use m_pawcprj,            only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_getdim
 221  use m_pawfgr,             only : pawfgr_type
 222  use m_energies,           only : energies_type
 223  use m_hamiltonian,        only : init_hamiltonian,destroy_hamiltonian,&
 224 &                                 load_spin_hamiltonian,load_k_hamiltonian,gs_hamiltonian_type
 225  use m_bandfft_kpt,        only : bandfft_kpt,bandfft_kpt_type,bandfft_kpt_set_ikpt,&
 226 &                                 bandfft_kpt_savetabs,bandfft_kpt_restoretabs
 227  use m_electronpositron,   only : electronpositron_type,electronpositron_calctype
 228  use m_paw_dmft,           only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft,saveocc_dmft
 229  use m_crystal,            only : crystal_init, crystal_free, crystal_t
 230  use m_oper,               only : oper_type,init_oper,destroy_oper
 231  use m_io_tools,           only : flush_unit
 232  use m_abi2big,            only : wvl_occ_abi2big,wvl_rho_abi2big,wvl_occopt_abi2big
 233  use m_fock,               only : fock_type,fock_ACE_type,fock_updateikpt,fock_calc_ene
 234  use m_invovl,             only : make_invovl
 235  use m_tddft,              only : tddft
 236  use m_kg,                 only : mkkin
 237  use m_suscep_stat,        only : suscep_stat
 238 
 239 #if defined HAVE_BIGDFT
 240  use BigDFT_API,           only : last_orthon,evaltoocc,write_energies
 241 #endif
 242 
 243 !This section has been created automatically by the script Abilint (TD).
 244 !Do not modify the following lines by hand.
 245 #undef ABI_FUNC
 246 #define ABI_FUNC 'vtorho'
 247  use interfaces_14_hidewrite
 248  use interfaces_18_timing
 249  use interfaces_32_util
 250  use interfaces_53_ffts
 251  use interfaces_56_recipspace
 252  use interfaces_62_wvl_wfs
 253  use interfaces_65_paw
 254  use interfaces_66_nonlocal
 255  use interfaces_66_wfs
 256  use interfaces_67_common
 257  use interfaces_68_dmft
 258  use interfaces_79_seqpar_mpi, except_this_one => vtorho
 259 !End of the abilint section
 260 
 261  implicit none
 262 
 263 !Arguments -------------------------------
 264  integer, intent(in) :: afford,dbl_nnsclo,dielop,dielstrt,istep,istep_mix,lmax_diel,mcg,mcprj,mgfftdiel
 265  integer, intent(in) :: my_natom,natom,nfftf,nfftdiel,nkxc,npwdiel
 266  integer, intent(in) :: ntypat,optforces,optres,pwind_alloc,usecprj
 267  real(dp), intent(in) :: cpus,etotal,gsqcut,ucvol
 268  real(dp), intent(out) :: compch_fft,nres2,residm
 269  type(MPI_type), intent(inout) :: mpi_enreg
 270  type(datafiles_type), intent(in) :: dtfil
 271  type(dataset_type), intent(inout) :: dtset
 272  type(efield_type), intent(inout) :: dtefield
 273  type(electronpositron_type),pointer :: electronpositron
 274  type(energies_type), intent(inout) :: energies
 275  type(hdr_type), intent(inout) :: hdr
 276  type(paw_dmft_type), intent(inout)  :: paw_dmft
 277  type(pawang_type), intent(in) :: pawang
 278  type(pawfgr_type), intent(in) :: pawfgr
 279  type(pseudopotential_type), intent(in) :: psps
 280  type(fock_type),pointer, intent(inout) :: fock
 281  type(wffile_type), intent(inout) :: wffnew
 282  type(wvl_data), intent(inout) :: wvl
 283  integer, intent(in) :: atindx(natom),atindx1(natom),gbound_diel(2*mgfftdiel+8,2)
 284  integer, intent(in) :: indsym(4,dtset%nsym,natom)
 285  integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 286  integer, intent(in) :: irrzondiel(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 287  integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),kg_diel(3,npwdiel),nattyp(ntypat),ngfftdiel(18),npwarr(dtset%nkpt)
 288  integer, intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym)
 289  real(dp), intent(in) :: dmatpawu(:,:,:,:),gmet(3,3),gprimd(3,3),ph1d(2,3*(2*dtset%mgfft+1)*natom)
 290  real(dp), intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*psps%usepaw)
 291  real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 292  real(dp), intent(in) :: phnonsdiel(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 293  real(dp), intent(in) :: pwnsfac(2,pwind_alloc),rmet(3,3),rprimd(3,3)
 294  real(dp), intent(inout) :: vtrial(nfftf,dtset%nspden)
 295  real(dp), intent(inout) :: xred(3,natom)
 296  real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
 297  real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
 298  real(dp), intent(in) :: ylmdiel(npwdiel,lmax_diel**2)
 299  real(dp), intent(out) :: dphase(3),grnl(3*natom)
 300  real(dp), intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
 301  real(dp), intent(out) :: nhat(nfftf,dtset%nspden*psps%usepaw)
 302  real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
 303  real(dp), intent(out) :: nvresid(nfftf,dtset%nspden)
 304  real(dp), intent(out) :: susmat(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden)
 305  real(dp), intent(inout) :: cg(2,mcg)
 306  real(dp), intent(inout) :: kxc(nfftf,nkxc),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 307  real(dp), intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden)
 308  real(dp), intent(inout) :: taug(2,nfftf*dtset%usekden),taur(nfftf,dtset%nspden*dtset%usekden)
 309  real(dp), intent(inout),optional :: vxctau(nfftf,dtset%nspden,4*dtset%usekden)
 310  type(pawcprj_type),allocatable,intent(inout) :: cprj(:,:)
 311  type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw)
 312  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
 313  type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw)
 314  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
 315 
 316 !Local variables-------------------------------
 317 !scalars
 318  integer,parameter :: level=111,tim_mkrho=2
 319  integer,save :: nwarning=0
 320  integer :: bdtot_index,counter,cplex,dimffnl,enunit,iband,iband1,ibdkpt
 321  integer :: ibg,icg,ider,idir,ierr,ifft,ifor,ifor1,ii,ikg,ikpt
 322  integer :: ikpt_loc,ikpt1,my_ikpt,ikxc,ilm,imagn,index1,iorder_cprj,ipert,iplex
 323  integer :: iscf,ispden,isppol,istwf_k,mband_cprj,mbdkpsp,mb2dkpsp
 324  integer :: mcgq,mcprj_local,mcprj_tmp,me_distrb,mkgq,mpi_comm_sphgrid
 325  integer :: mwarning,my_nspinor,n1,n2,n3,n4,n5,n6,nband_eff
 326  integer :: nband_k,nband_cprj_k,nbuf,neglect_pawhat,nfftot,nkpg,nkpt1,nnn,nnsclo_now
 327  integer :: nproc_distrb,npw_k,nspden_rhoij,option,prtvol
 328  integer :: spaceComm_distrb,usecprj_local,usefock_ACE,usetimerev
 329  logical :: berryflag,computesusmat,fixed_occ
 330  logical :: locc_test,paral_atom,remove_inv,usefock,with_vxctau
 331  logical :: do_last_ortho,wvlbigdft=.false.
 332  real(dp) :: dmft_ldaocc
 333  real(dp) :: edmft,ebandlda,ebanddmft,ebandldatot,ekindmft,ekindmft2,ekinlda
 334  real(dp) :: min_occ,vxcavg_dum,strsxc(6)
 335  character(len=500) :: message
 336  type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
 337  type(gs_hamiltonian_type) :: gs_hamk
 338 !arrays
 339  integer,allocatable :: kg_k(:,:)
 340  real(dp) :: dielar(7),dphase_k(3),kpoint(3),qpt(3),rhodum(1),tsec(2),ylmgr_dum(0,0,0)
 341  real(dp),allocatable :: EigMin(:,:),buffer1(:),buffer2(:),cgq(:,:)
 342  real(dp),allocatable :: cgrkxc(:,:),cgrvtrial(:,:),doccde(:)
 343  real(dp),allocatable :: dphasek(:,:),eig_k(:),ek_k(:),ek_k_nd(:,:,:),eknk(:),eknk_nd(:,:,:,:,:)
 344  real(dp),allocatable :: enl_k(:),enlnk(:),focknk(:),fockfornk(:,:,:),ffnl(:,:,:,:),grnl_k(:,:), xcart(:,:)
 345  real(dp),allocatable :: grnlnk(:,:),kinpw(:),kpg_k(:,:),occ_k(:),ph3d(:,:,:)
 346  real(dp),allocatable :: pwnsfacq(:,:),resid_k(:),rhoaug(:,:,:,:),rhowfg(:,:),rhowfr(:,:)
 347  real(dp),allocatable :: vlocal(:,:,:,:),vlocal_tmp(:,:,:),vxctaulocal(:,:,:,:,:),ylm_k(:,:),zshift(:)
 348  complex(dpc),allocatable :: nucdipmom_k(:)
 349  type(pawcprj_type),allocatable :: cprj_tmp(:,:)
 350  type(oper_type) :: lda_occup
 351  type(pawrhoij_type),pointer :: pawrhoij_unsym(:)
 352  type(crystal_t) :: cryst_struc
 353  integer :: idum1(0),idum3(0,0,0)
 354  real(dp) :: rdum2(0,0),rdum4(0,0,0,0)
 355 
 356 !Variables for BigDFT
 357 #if defined HAVE_BIGDFT
 358  integer :: occopt_bigdft
 359 #endif
 360 
 361 ! *********************************************************************
 362 
 363  DBG_ENTER("COLL")
 364 
 365 !Keep track of total time spent in vtorho
 366  call timab(980,1,tsec)
 367  call timab(981,1,tsec)
 368 
 369 !Structured debugging if prtvol==-level
 370  prtvol=dtset%prtvol
 371 
 372 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
 373  wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
 374 
 375 !Several inits
 376  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
 377  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
 378  usecprj_local=0;if (psps%usepaw==1) usecprj_local=1
 379  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 380  paral_atom=(my_natom/=natom)
 381  compch_fft=-1.d5
 382 
 383 !Check that usekden is not 0 if want to use vxctau
 384  with_vxctau = (dtset%usekden/=0)
 385 
 386 !Check that fock is present if want to use fock option
 387  usefock = (dtset%usefock==1 .and. associated(fock))
 388  usefock_ACE=0
 389  if (usefock) usefock_ACE=fock%fock_common%use_ACE
 390 
 391 
 392 !Init MPI
 393  spaceComm_distrb=mpi_enreg%comm_cell
 394  if (mpi_enreg%paral_kgb==1) spaceComm_distrb=mpi_enreg%comm_kpt
 395  if (mpi_enreg%paral_hf ==1) spaceComm_distrb=mpi_enreg%comm_kpt
 396  nproc_distrb=xmpi_comm_size(spaceComm_distrb)
 397  me_distrb=xmpi_comm_rank(spaceComm_distrb)
 398  mpi_comm_sphgrid=mpi_enreg%comm_fft
 399  if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
 400  if (mpi_enreg%me_img/=0) nwarning=nwarning+1
 401 
 402 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
 403  if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
 404    MSG_BUG('wrong values for nfft, nfftf!')
 405  end if
 406 
 407 !Test optforces (to prevent memory overflow)
 408  if (optforces/=0.and.optforces/=1) then
 409    write(message,'(a,i0)')' wrong value for optforces = ',optforces
 410    MSG_BUG(message)
 411  end if
 412 
 413  iscf=dtset%iscf
 414  fixed_occ=(dtset%occopt<3.or.electronpositron_calctype(electronpositron)==1)
 415  if(.not. wvlbigdft) then
 416    energies%e_eigenvalues = zero
 417    energies%e_kinetic     = zero
 418    energies%e_nonlocalpsp = zero
 419    if (usefock) then
 420      energies%e_fock=zero
 421      energies%e_fockdc=zero
 422    end if
 423    grnl(:)=zero
 424    resid(:) = zero ! JWZ 13 May 2010. resid and eigen need to be fully zeroed each time before use
 425    eigen(:) = zero
 426    bdtot_index=0
 427    ibg=0;icg=0
 428    mbdkpsp=dtset%mband*dtset%nkpt*dtset%nsppol
 429    if(paw_dmft%use_dmft==1) mb2dkpsp=2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol
 430  end if
 431 
 432  if(dtset%usewvl==0) then
 433    ABI_ALLOCATE(eknk,(mbdkpsp))
 434    ABI_ALLOCATE(enlnk,(mbdkpsp))
 435    ABI_ALLOCATE(eknk_nd,(dtset%nsppol,dtset%nkpt,2,dtset%mband,dtset%mband*paw_dmft%use_dmft))
 436    ABI_ALLOCATE(EigMin,(2,dtset%mband))
 437    ABI_ALLOCATE(grnlnk,(3*natom,mbdkpsp*optforces))
 438    if (usefock) then
 439      ABI_ALLOCATE(focknk,(mbdkpsp))
 440      focknk=zero
 441      if (optforces>0)then
 442        ABI_ALLOCATE(fockfornk,(3,natom,mbdkpsp))
 443        fockfornk=zero
 444      end if
 445    end if
 446    eknk(:)=zero;enlnk(:)=zero
 447    if (optforces>0) grnlnk(:,:)=zero
 448    if(paw_dmft%use_dmft==1) eknk_nd=zero
 449  end if !usewvl==0
 450 
 451 !Initialize rhor if needed; store old rhor
 452  if(iscf>=0 .or. iscf==-3) then
 453    if (optres==1) nvresid=rhor
 454 !  NC and plane waves
 455    if (psps%usepaw==0 .and. dtset%usewvl==0) then
 456      rhor=zero
 457 !    PAW
 458    elseif(psps%usepaw==1) then
 459      ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden))
 460      ABI_ALLOCATE(rhowfg,(2,dtset%nfft))
 461      rhowfr(:,:)=zero
 462    end if
 463  end if
 464 
 465 !Set max number of non-self-consistent loops nnsclo_now for use in vtowfk
 466  if(iscf<0)then
 467    ! ===== Non self-consistent =====
 468    nnsclo_now=dtset%nstep
 469  else
 470    ! ===== Self-consistent =====
 471    if(dtset%nnsclo>0) then
 472    ! ===== Self-consistent + imposed =====
 473      nnsclo_now=dtset%nnsclo
 474    else
 475    ! ===== Self-consistent + default =====
 476      nnsclo_now=1
 477      if (dtset%usewvl==0) then
 478      ! ----- Plane waves -----
 479        if (istep<=2.and.iscf/=0) nnsclo_now=2
 480      else
 481      ! ----- Wavelets -----
 482        if (iscf==0) then
 483          nnsclo_now=0
 484        else if (istep<=2) then
 485          nnsclo_now=3
 486        else if (istep<=4) then
 487          nnsclo_now=2
 488        end if
 489      end if
 490    end if
 491    ! ===== Double is required =====
 492    if(dbl_nnsclo==1)then
 493      nnsclo_now=nnsclo_now*2
 494    end if
 495  end if
 496  if(dtset%wfoptalg==2)nnsclo_now=40  ! UNDER DEVELOPMENT
 497 
 498  write(message, '(a,i0,a,3(i0,1x))' ) ' vtorho : nnsclo_now=',nnsclo_now,&
 499 & ', note that nnsclo,dbl_nnsclo,istep=',dtset%nnsclo,dbl_nnsclo,istep
 500  call wrtout(std_out,message,'COLL')
 501 
 502 !==== Initialize most of the Hamiltonian ====
 503 !Allocate all arrays and initialize quantities that do not depend on k and spin.
 504  call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,&
 505 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
 506 & paw_ij=paw_ij,ph1d=ph1d,usecprj=usecprj_local,electronpositron=electronpositron,fock=fock,&
 507 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
 508 & nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
 509 
 510 !Initializations for PAW (projected wave functions)
 511  mcprj_local=0 ; mband_cprj=0
 512  if (psps%usepaw==1) then
 513    mband_cprj=dtset%mband
 514    if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band
 515    iorder_cprj=0 ; mcprj_local=mcprj
 516    if (usecprj==0) then
 517      mcprj_local=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
 518      !This is a check but should always be true since scfcv allocated cprj anyway
 519      if (allocated(cprj)) then
 520        !Was allocated in scfcv so we just destroy and reconstruct it as desired
 521        call pawcprj_free(cprj)
 522        ABI_DATATYPE_DEALLOCATE(cprj)
 523      end if
 524      ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj_local))
 525      call pawcprj_alloc(cprj,0,gs_hamk%dimcprj)
 526    end if
 527  end if
 528 
 529  call timab(981,2,tsec)
 530 
 531 !===================================================================
 532 ! WAVELETS - Branching with a separate VTORHO procedure
 533 !===================================================================
 534 
 535  if (dtset%usewvl == 1) then
 536 #ifndef HAVE_BIGDFT
 537    BIGDFT_NOTENABLED_ERROR()
 538 #else
 539 
 540 !  do_last_ortho in case of diagonalization scheme
 541    if (     wvlbigdft) do_last_ortho=(dtset%iscf/=0)
 542    if (.not.wvlbigdft) do_last_ortho=(.true.)
 543 
 544    ABI_ALLOCATE(xcart,(3, dtset%natom))
 545    call xred2xcart(dtset%natom, rprimd, xcart, xred)
 546 
 547    if(wvlbigdft) then
 548 !    NSCF loop for wvlbigdt:
 549      call wvl_nscf_loop_bigdft()
 550    else
 551 !    NSCF loop for WVL: (not wvlbigdft)
 552      call wvl_nscf_loop()
 553    end if
 554 
 555 !  Eventually orthogonalize WFs now
 556    if (do_last_ortho) then
 557      call write_energies(ii,0,wvl%e%energs,0.d0,0.d0,"final")
 558      call last_orthon(me_distrb, nproc_distrb, ii, wvl%wfs%ks, wvl%e%energs%evsum, .true.)
 559      if(wvlbigdft) energies%e_xcdc = wvl%e%energs%evxc
 560 !    If occupation numbers are not changed...
 561      if (fixed_occ .or. (iscf<0 .and. iscf/=-3)) then
 562        call wvl_comm_eigen()
 563      end if
 564 !    ... or update occupations:
 565      if( ( .not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then
 566        if(wvlbigdft) then
 567          call wvl_occ_bigdft()
 568        else
 569 !        Communicate eigenvalues:
 570          call wvl_comm_eigen()
 571 !        Update occ and Fermi level
 572          call wvl_occ()
 573        end if
 574      end if
 575 !    This might accelerate convergence:
 576      wvl%wfs%ks%diis%energy_min=one
 577      wvl%wfs%ks%diis%alpha=two
 578    end if !do_last_ortho
 579 
 580 !  Compute eigenvalues energy
 581    if(.not. wvlbigdft .and. nnsclo_now>0) then
 582      call e_eigen(eigen,energies%e_eigenvalues,dtset%mband,dtset%nband,dtset%nkpt,&
 583 &     dtset%nsppol,occ,dtset%wtk)
 584    else
 585      energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp &
 586 &     + energies%e_xcdc  + two*energies%e_hartree +energies%e_nonlocalpsp
 587    end if
 588 
 589    if (optforces == 1) then ! not compatible with iscf=0 and wvlbigdftcomp=1 + PAW
 590      call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart)
 591    end if
 592 
 593 !  For iscf<0 we do not update the density
 594    if (dtset%iscf>=0) then !(dtset%iscf>=0 .and. .not. wvlbigdft ) then
 595      call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den)
 596    end if
 597    ABI_DEALLOCATE(xcart)
 598 
 599 !  Note in WVL+NC: the rest will be skipped.
 600 !  For PAW: we will compute Rho_ij at the end.
 601    !if(wvlbigdft) return
 602 #endif
 603  else
 604 
 605 !===================================================================
 606 ! PLANE WAVES - Standard VTORHO procedure
 607 !===================================================================
 608 
 609 !  Electric fields: set flag to turn on various behaviors
 610    berryflag = (dtset%berryopt == 4 .or. dtset%berryopt == 14 .or. &
 611 &   dtset%berryopt == 6 .or. dtset%berryopt == 16 .or. &
 612 &   dtset%berryopt == 7 .or. dtset%berryopt == 17)
 613 
 614 !  Electric field: allocate dphasek
 615    nkpt1 = dtset%nkpt
 616    if ( berryflag ) then
 617      ABI_ALLOCATE(dphasek,(3,dtset%nkpt*dtset%nsppol))
 618      dphasek(:,:) = zero
 619      nkpt1 = dtefield%mkmem_max
 620    end if
 621 
 622    ABI_ALLOCATE(rhoaug,(n4,n5,n6,gs_hamk%nvloc))
 623    ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamk%nvloc))
 624    if(with_vxctau) then
 625      ABI_ALLOCATE(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4))
 626    end if
 627 
 628 !  LOOP OVER SPINS
 629    do isppol=1,dtset%nsppol
 630      call timab(982,1,tsec)
 631 
 632      ikpt_loc = 0
 633      ikg=0
 634 
 635 !    Set up local potential vlocal with proper dimensioning, from vtrial
 636 !    Also take into account the spin.
 637      if(dtset%nspden/=4)then
 638        if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then
 639          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal,2)
 640          if(with_vxctau) then
 641            do ii=1,4
 642              call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,&
 643 &             vxctau(:,:,ii),vxctaulocal(:,:,:,:,ii),2)
 644            end do
 645          end if
 646        else
 647          ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden))
 648          call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
 649          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal,2)
 650          ABI_DEALLOCATE(cgrvtrial)
 651        end if
 652      else
 653        ABI_ALLOCATE(vlocal_tmp,(n4,n5,n6))
 654        if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then
 655          do ispden=1,dtset%nspden
 656            call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal_tmp,2)
 657            vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
 658          end do
 659        else
 660          ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden))
 661          call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
 662          do ispden=1,dtset%nspden
 663            call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal_tmp,2)
 664            vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
 665          end do
 666          ABI_DEALLOCATE(cgrvtrial)
 667        end if
 668        ABI_DEALLOCATE(vlocal_tmp)
 669      end if ! nspden
 670 
 671      rhoaug(:,:,:,:)=zero
 672 
 673 !    Continue to initialize the Hamiltonian
 674      call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.)
 675      if (with_vxctau) then
 676        call load_spin_hamiltonian(gs_hamk,isppol,vxctaulocal=vxctaulocal)
 677      end if
 678 
 679      call timab(982,2,tsec)
 680 
 681 !    BIG FAT k POINT LOOP
 682 !    MVeithen: I had to modify the structure of this loop in order to implement MPI // of the electric field
 683 !    Note that the loop here differs from the similar one in berryphase_new.F90.
 684 !    here, ikpt_loc numbers the kpts treated by the current processor.
 685 !    in berryphase_new.F90, ikpt_loc ALSO includes info about value of isppol.
 686 
 687      ikpt = 0
 688      do while (ikpt_loc < nkpt1)
 689 
 690        call timab(997,1,tsec)
 691 
 692        if ( .not.berryflag ) then
 693          ikpt_loc = ikpt_loc + 1
 694          ikpt = ikpt_loc
 695          my_ikpt = mpi_enreg%my_kpttab(ikpt)
 696        else
 697          if (ikpt_loc < dtset%mkmem) ikpt = ikpt + 1
 698          if ((ikpt > dtset%nkpt).and.(ikpt_loc < dtset%mkmem)) exit
 699          my_ikpt=ikpt
 700        end if
 701 
 702        dphase_k(:) = zero
 703        counter=100*ikpt+isppol
 704        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 705        nband_cprj_k=nband_k/mpi_enreg%nproc_band
 706        istwf_k=dtset%istwfk(ikpt)
 707        npw_k=npwarr(ikpt)
 708 
 709        mcgq=1 ; mkgq=1
 710        if (.not.berryflag) then
 711          if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
 712            eigen(1+bdtot_index : nband_k+bdtot_index) = zero
 713            resid(1+bdtot_index : nband_k+bdtot_index) = zero
 714            bdtot_index=bdtot_index+nband_k
 715            cycle
 716          end if
 717        else
 718          if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) &
 719 &         .and.(ikpt_loc <= dtset%mkmem)) then
 720            eigen(1+bdtot_index : nband_k+bdtot_index) = zero
 721            resid(1+bdtot_index : nband_k+bdtot_index) = zero
 722            bdtot_index = bdtot_index + nband_k
 723            cycle
 724          end if
 725          ikpt_loc = ikpt_loc + 1
 726          mcgq = dtset%mpw*my_nspinor*nband_k*dtefield%nneigh(ikpt)
 727          ikg = dtefield%kgindex(ikpt)
 728          mkgq = 6*dtset%mpw
 729        end if ! berryflag
 730 
 731        call timab(997,2,tsec)
 732 
 733 !      In case of MPI // of a finite field calculation
 734 !      build the cgq array that stores the wavefunctions for the
 735 !      neighbours of ikpt, and the pwnsfacq array that stores the
 736 !      corresponding phase factors (in case of tnons)
 737        ABI_ALLOCATE(cgq,(2,mcgq))
 738        ABI_ALLOCATE(pwnsfacq,(2,mkgq))
 739        if ( berryflag ) then
 740          call cgq_builder(berryflag,cg,cgq,dtefield,dtset,ikpt,ikpt_loc,isppol,mcg,mcgq,&
 741 &         me_distrb,mkgq,mpi_enreg,my_nspinor,nband_k,nproc_distrb,&
 742 &         npwarr,pwnsfac,pwnsfacq,pwind_alloc,spaceComm_distrb)
 743          if (ikpt_loc > dtset%mkmem) then
 744            ABI_DEALLOCATE(cgq)
 745            ABI_DEALLOCATE(pwnsfacq)
 746            cycle
 747          end if
 748        end if !berryopt
 749 
 750        call timab(984,1,tsec)
 751 
 752        if (mpi_enreg%paral_kgb==1) my_bandfft_kpt => bandfft_kpt(my_ikpt)
 753        call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
 754 !      my_ikpt = ikpt
 755 !       if (mpi_enreg%paral_kgb==1) then
 756 !        my_ikpt = mpi_enreg%my_kpttab(ikpt)
 757 !         my_bandfft_kpt => bandfft_kpt(my_ikpt)
 758 !         call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
 759 !       end if
 760 
 761        ABI_ALLOCATE(eig_k,(nband_k))
 762        ABI_ALLOCATE(ek_k,(nband_k))
 763        ABI_ALLOCATE(enl_k,(nband_k))
 764        ABI_ALLOCATE(ek_k_nd,(2,nband_k,nband_k*paw_dmft%use_dmft))
 765        ABI_ALLOCATE(occ_k,(nband_k))
 766        ABI_ALLOCATE(resid_k,(nband_k))
 767        ABI_ALLOCATE(zshift,(nband_k))
 768        ABI_ALLOCATE(grnl_k,(3*natom,nband_k*optforces))
 769 
 770        eig_k(:)=zero
 771        ek_k(:)=zero
 772        enl_k(:)=zero
 773        if(paw_dmft%use_dmft==1) ek_k_nd(:,:,:)=zero
 774        if (optforces>0) grnl_k(:,:)=zero
 775        kpoint(:)=dtset%kptns(:,ikpt)
 776        occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
 777        resid_k(:)=zero
 778        zshift(:)=dtset%eshift
 779 
 780        ABI_ALLOCATE(kg_k,(3,npw_k))
 781        ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
 782        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
 783        if (psps%useylm==1) then
 784          do ilm=1,psps%mpsang*psps%mpsang
 785            ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
 786          end do
 787        end if
 788 
 789 !      Set up remaining of the Hamiltonian
 790 
 791 !      Compute (1/2) (2 Pi)**2 (k+G)**2:
 792        ABI_ALLOCATE(kinpw,(npw_k))
 793 !       call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k)
 794        call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
 795 
 796 !      Compute (k+G) vectors (only if useylm=1)
 797        nkpg=3*optforces*dtset%nloalg(3)
 798        ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
 799        if ((mpi_enreg%paral_kgb/=1.or.istep<=1).and.nkpg>0) then
 800          call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
 801        end if
 802 
 803 !      Compute nonlocal form factors ffnl at all (k+G):
 804        ider=0;idir=0;dimffnl=1
 805        ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
 806        if (mpi_enreg%paral_kgb/=1.or.istep<=1) then
 807          call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
 808 &         gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
 809 &         psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
 810 &         npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
 811 &         psps%usepaw,psps%useylm,ylm_k,ylmgr)
 812        end if
 813 
 814 !     compute and load nuclear dipole Hamiltonian at current k point
 815        if(any(abs(gs_hamk%nucdipmom)>0.0)) then
 816          if(allocated(nucdipmom_k)) then
 817            ABI_DEALLOCATE(nucdipmom_k)
 818          end if
 819          ABI_ALLOCATE(nucdipmom_k,(npw_k*(npw_k+1)/2))
 820          call mknucdipmom_k(gmet,kg_k,kpoint,natom,gs_hamk%nucdipmom,nucdipmom_k,npw_k,rprimd,ucvol,xred)
 821          call load_k_hamiltonian(gs_hamk,nucdipmom_k=nucdipmom_k)
 822        end if
 823 
 824 
 825 !      Load k-dependent part in the Hamiltonian datastructure
 826 !       - Compute 3D phase factors
 827 !       - Prepare various tabs in case of band-FFT parallelism
 828 !       - Load k-dependent quantities in the Hamiltonian
 829        ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
 830 
 831        if(usefock_ACE/=0) then
 832          call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,&
 833 &         kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,fockACE_k=fock%fockACE(ikpt,isppol),ph3d_k=ph3d,&
 834 &         compute_ph3d=(mpi_enreg%paral_kgb/=1.or.istep<=1),&
 835 &         compute_gbound=(mpi_enreg%paral_kgb/=1))
 836        else
 837          call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,&
 838 &         kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,&
 839 &         compute_ph3d=(mpi_enreg%paral_kgb/=1.or.istep<=1),&
 840 &         compute_gbound=(mpi_enreg%paral_kgb/=1))
 841        end if
 842 
 843 !      Load band-FFT tabs (transposed k-dependent arrays)
 844        if (mpi_enreg%paral_kgb==1) then
 845          if (istep<=1) then
 846            call prep_bandfft_tabs(gs_hamk,ikpt,dtset%mkmem,mpi_enreg)
 847          end if
 848          call load_k_hamiltonian(gs_hamk,npw_fft_k=my_bandfft_kpt%ndatarecv, &
 849 &         gbound_k =my_bandfft_kpt%gbound, &
 850 &         kinpw_k  =my_bandfft_kpt%kinpw_gather, &
 851 &         kg_k     =my_bandfft_kpt%kg_k_gather, &
 852 &         kpg_k    =my_bandfft_kpt%kpg_k_gather, &
 853          ffnl_k   =my_bandfft_kpt%ffnl_gather, &
 854          ph3d_k   =my_bandfft_kpt%ph3d_gather)
 855        end if
 856 
 857 !      Build inverse of overlap matrix for chebfi
 858        if(psps%usepaw == 1 .and. dtset%wfoptalg == 1 .and. istep <= 1) then
 859          call make_invovl(gs_hamk, dimffnl, ffnl, ph3d, mpi_enreg)
 860        end if
 861 
 862        ! Setup gemm_nonlop
 863        if (gemm_nonlop_use_gemm) then
 864          !set the global variable indicating to gemm_nonlop where to get its data from
 865          gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt
 866          if (istep <= 1) then
 867            !Init the arrays
 868            call make_gemm_nonlop(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%lmnmax, &
 869 &           gs_hamk%ntypat, gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%istwf_k, gs_hamk%ucvol, gs_hamk%ffnl_k,&
 870 &           gs_hamk%ph3d_k)
 871          end if
 872        end if
 873 
 874 #if defined HAVE_GPU_CUDA
 875        if (dtset%use_gpu_cuda==1) then
 876          call gpu_update_ffnl_ph3d(ph3d,size(ph3d),ffnl,size(ffnl))
 877        end if
 878 #endif
 879 
 880        call timab(984,2,tsec)
 881 
 882 !      Update the value of ikpt,isppol in fock_exchange and allocate the memory space to perform HF calculation.
 883        if (usefock) then
 884          call fock_updateikpt(fock%fock_common,ikpt,isppol)
 885        end if
 886        if ((psps%usepaw==1).and.(usefock)) then
 887          if ((fock%fock_common%optfor).and.(usefock_ACE==0)) then
 888            fock%fock_common%forces_ikpt=zero
 889          end if
 890        end if
 891 
 892 !      Compute the eigenvalues, wavefunction, residuals,
 893 !      contributions to kinetic energy, nonlocal energy, forces,
 894 !      and update of rhor to this k-point and this spin polarization.
 895        call vtowfk(cg,cgq,cprj,cpus,dphase_k,dtefield,dtfil,&
 896 &       dtset,eig_k,ek_k,ek_k_nd,enl_k,fixed_occ,grnl_k,gs_hamk,&
 897 &       ibg,icg,ikpt,iscf,isppol,kg_k,kinpw,mband_cprj,mcg,mcgq,mcprj_local,mkgq,&
 898 &       mpi_enreg,dtset%mpw,natom,nband_k,dtset%nkpt,nnsclo_now,npw_k,npwarr,&
 899 &       occ_k,optforces,prtvol,pwind,pwind_alloc,pwnsfac,pwnsfacq,resid_k,&
 900 &       rhoaug,paw_dmft,dtset%wtk(ikpt),zshift)
 901 
 902        call timab(985,1,tsec)
 903 
 904 #if defined HAVE_GPU_CUDA
 905        if(dtset%use_gpu_cuda==1) then
 906          call gpu_finalize_ffnl_ph3d()
 907        end if
 908 #endif
 909        ABI_DEALLOCATE(ffnl)
 910        ABI_DEALLOCATE(kinpw)
 911        ABI_DEALLOCATE(kg_k)
 912        ABI_DEALLOCATE(kpg_k)
 913        if(allocated(nucdipmom_k)) then
 914          ABI_DEALLOCATE(nucdipmom_k)
 915        end if
 916        ABI_DEALLOCATE(ylm_k)
 917        ABI_DEALLOCATE(ph3d)
 918        ABI_DEALLOCATE(cgq)
 919        ABI_DEALLOCATE(pwnsfacq)
 920 
 921 !      electric field
 922        if (berryflag) then
 923          dphasek(:,ikpt + (isppol - 1)*dtset%nkpt) = dphase_k(:)
 924 
 925 !        The overlap matrices for all first neighbours of ikpt are no more up to date
 926          do idir = 1, 3
 927            do ifor = 1, 2
 928              ikpt1 = dtefield%ikpt_dk(dtefield%i2fbz(ikpt),ifor,idir)
 929              ikpt1 = dtefield%indkk_f2ibz(ikpt1,1)
 930              ifor1 = -1*ifor + 3   ! ifor = 1 -> ifor1 = 2 & ifor = 2 -> ifor1 = 1
 931              dtefield%sflag(:,ikpt1+(isppol-1)*dtset%nkpt,ifor1,idir) = 0
 932            end do
 933          end do
 934        end if  ! berryflag
 935 
 936 !      Save eigenvalues (hartree), residuals (hartree**2)
 937        eigen(1+bdtot_index : nband_k+bdtot_index) = eig_k(:)
 938        eknk (1+bdtot_index : nband_k+bdtot_index) = ek_k (:)
 939        if(usefock) then
 940          focknk (1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%eigen_ikpt (:)
 941          if (optforces>0) fockfornk(:,:,1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%forces_ikpt(:,:,:)
 942        end if
 943        if(paw_dmft%use_dmft==1) eknk_nd(isppol,ikpt,:,:,:) = ek_k_nd(:,:,:)
 944        resid(1+bdtot_index : nband_k+bdtot_index) = resid_k(:)
 945        if (optforces>0) grnlnk(:,1+bdtot_index : nband_k+bdtot_index) = grnl_k(:,:)
 946        enlnk(1+bdtot_index : nband_k+bdtot_index) = enl_k(:)
 947 
 948        if(iscf>0 .or. iscf==-3)then
 949 !        Accumulate sum over k points for band, nonlocal and kinetic energies,
 950 !        also accumulate gradients of Enonlocal:
 951          do iband=1,nband_k
 952            if (abs(occ_k(iband))>tol8) then
 953              energies%e_kinetic     = energies%e_kinetic     + dtset%wtk(ikpt)*occ_k(iband)*ek_k(iband)
 954              energies%e_eigenvalues = energies%e_eigenvalues + dtset%wtk(ikpt)*occ_k(iband)*eig_k(iband)
 955              energies%e_nonlocalpsp = energies%e_nonlocalpsp + dtset%wtk(ikpt)*occ_k(iband)*enl_k(iband)
 956              if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ_k(iband)*grnl_k(:,iband)
 957              if (usefock) energies%e_fock=energies%e_fock + half*fock%fock_common%eigen_ikpt(iband)*occ_k(iband)*dtset%wtk(ikpt)
 958            end if
 959          end do
 960 
 961 !        Calculate Fock contribution to the total energy if required
 962          if ((psps%usepaw==1).and.(usefock)) then
 963            if ((fock%fock_common%optfor).and.(usefock_ACE==0)) then
 964              call fock_calc_ene(dtset,fock%fock_common,energies%e_exactX,ikpt,nband_k,occ_k)
 965            end if
 966          end if
 967        end if
 968 
 969        call timab(985,2,tsec)
 970 
 971        ABI_DEALLOCATE(eig_k)
 972        ABI_DEALLOCATE(ek_k)
 973        ABI_DEALLOCATE(ek_k_nd)
 974        ABI_DEALLOCATE(grnl_k)
 975        ABI_DEALLOCATE(occ_k)
 976        ABI_DEALLOCATE(resid_k)
 977        ABI_DEALLOCATE(zshift)
 978        ABI_DEALLOCATE(enl_k)
 979 
 980 !      Keep track of total number of bands (all k points so far, even for k points not treated by me)
 981        bdtot_index=bdtot_index+nband_k
 982 
 983 !      Also shift array memory if dtset%mkmem/=0
 984        if (dtset%mkmem/=0) then
 985          ibg=ibg+my_nspinor*nband_cprj_k
 986          icg=icg+npw_k*my_nspinor*nband_k
 987          ikg=ikg+npw_k
 988        end if
 989 
 990      end do ! End big k point loop
 991 
 992      call timab(986,1,tsec)
 993 
 994      if (fixed_occ .and. mpi_enreg%paral_kgb==1) then
 995        call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr) !Sum the contributions over bands/FFT/spinors
 996      end if
 997 
 998 !    Transfer density on augmented fft grid to normal fft grid in real space
 999 !    Also take into account the spin.
1000      if(iscf>0.or.iscf==-3)then
1001        if (psps%usepaw==0) then
1002          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,1),1)
1003          if(dtset%nspden==4)then
1004            do imagn=2,4
1005              call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,imagn),1)
1006            end do
1007          end if
1008        else
1009          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,1),1)
1010          if(dtset%nspden==4)then
1011            do imagn=2,4
1012              call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,imagn),1)
1013            end do
1014          end if
1015        end if
1016      end if
1017 
1018      call timab(986,2,tsec)
1019 
1020    end do ! End loop over spins
1021 
1022    call timab(988,1,tsec)
1023    if (usefock) then
1024      if(fock%fock_common%optfor) then
1025        call xmpi_sum(fock%fock_common%forces,mpi_enreg%comm_kpt,ierr)
1026      end if
1027    end if
1028 !  Electric field: compute string-averaged change in Zak phase
1029 !  along each direction, store it in dphase(idir)
1030 !  ji: it is not convenient to do this anymore. Remove. Set dphase(idir)=0.0_dp.
1031 !  eventually, dphase(idir) will have to go...
1032    if (berryflag) then
1033      dphase(:) = zero
1034 !    In case of MPI // of a finite field calculation, send dphasek to all cpus
1035      call xmpi_sum(dphasek,spaceComm_distrb,ierr)
1036      ABI_DEALLOCATE(dphasek)
1037    end if ! berryflag
1038    ABI_DEALLOCATE(rhoaug)
1039    ABI_DEALLOCATE(vlocal)
1040    if(with_vxctau) then
1041      ABI_DEALLOCATE(vxctaulocal)
1042    end if
1043 
1044    call timab(988,2,tsec)
1045 
1046    ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
1047    doccde(:)=zero
1048 
1049 !  Treat now varying occupation numbers, in the self-consistent case
1050    if((.not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then
1051 
1052 !    Parallel case
1053      if (mpi_enreg%nproc_kpt>1) then
1054 
1055        call timab(989,1,tsec)
1056 
1057 !      If needed, exchange the values of eigen,resid,eknk,enlnk,grnlnk
1058        ABI_ALLOCATE(buffer1,((4+3*natom*optforces+dtset%usefock+3*natom*dtset%usefock*optforces)*mbdkpsp))
1059        if(paw_dmft%use_dmft==1) then
1060          ABI_ALLOCATE(buffer2,(mb2dkpsp*paw_dmft%use_dmft))
1061        end if
1062 !      Pack eigen,resid,eknk,enlnk,grnlnk in buffer1
1063        buffer1(1          :  mbdkpsp)=eigen(:)
1064        buffer1(1+  mbdkpsp:2*mbdkpsp)=resid(:)
1065        buffer1(1+2*mbdkpsp:3*mbdkpsp)=eknk(:)
1066        buffer1(1+3*mbdkpsp:4*mbdkpsp)=enlnk(:)
1067        index1=4*mbdkpsp
1068        if (optforces>0) then
1069          buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(grnlnk,(/(3*natom)*mbdkpsp/) )
1070          index1=index1+3*natom*mbdkpsp
1071        end if
1072        if (usefock) then
1073          buffer1(1+index1:index1+mbdkpsp)=focknk(:)
1074          if (optforces>0) then
1075            index1=index1+mbdkpsp
1076            buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(fockfornk,(/(3*natom)*mbdkpsp/) )
1077          end if
1078        end if
1079        if(paw_dmft%use_dmft==1) then
1080          nnn=0
1081          do ikpt=1,dtset%nkpt
1082            do isppol=1,dtset%nsppol
1083              do iband=1,dtset%mband
1084                do iband1=1,dtset%mband
1085                  do iplex=1,2
1086                    nnn=nnn+1
1087                    buffer2(nnn)=eknk_nd(isppol,ikpt,iplex,iband,iband1)
1088                  end do
1089                end do
1090              end do
1091            end do
1092          end do
1093          if(nnn.ne.mb2dkpsp)  then
1094            write(message,*)' BUG in vtorho2, buffer2',nnn,mb2dkpsp
1095            MSG_BUG(message)
1096          end if
1097        end if
1098 !      Build sum of everything
1099        call timab(48,1,tsec)
1100        call xmpi_sum(buffer1,mpi_enreg%comm_kpt,ierr)
1101        if(mpi_enreg%paral_kgb/=1.and.paw_dmft%use_dmft==1) then
1102          call xmpi_sum(buffer2,mpi_enreg%comm_kpt,ierr)
1103        end if
1104        call timab(48,2,tsec)
1105 
1106 !      Unpack eigen,resid,eknk,enlnk,grnlnk in buffer1
1107        eigen(:) =buffer1(1          :  mbdkpsp)
1108        resid(:) =buffer1(1+  mbdkpsp:2*mbdkpsp)
1109        eknk(:)  =buffer1(1+2*mbdkpsp:3*mbdkpsp)
1110        enlnk(:) =buffer1(1+3*mbdkpsp:4*mbdkpsp)
1111        index1=4*mbdkpsp
1112        if (optforces>0) then
1113          grnlnk(:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3*natom,mbdkpsp/) )
1114        end if
1115        if (usefock) then
1116          focknk(:)=buffer1(1+index1:index1+mbdkpsp)
1117          if (optforces>0) then
1118            index1=index1+mbdkpsp
1119            fockfornk(:,:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3,natom,mbdkpsp/) )
1120          end if
1121        end if
1122        if(paw_dmft%use_dmft==1) then
1123          nnn=0
1124          do ikpt=1,dtset%nkpt
1125            do isppol=1,dtset%nsppol
1126              do iband=1,dtset%mband
1127                do iband1=1,dtset%mband
1128                  do iplex=1,2
1129                    nnn=nnn+1
1130                    eknk_nd(isppol,ikpt,iplex,iband,iband1)=buffer2(nnn)
1131                  end do
1132                end do
1133              end do
1134            end do
1135          end do
1136        end if
1137        if(allocated(buffer2))  then
1138          ABI_DEALLOCATE(buffer2)
1139        end if
1140        ABI_DEALLOCATE(buffer1)
1141        call timab(989,2,tsec)
1142 
1143      end if ! nproc_kpt>1
1144 
1145 !    Compute the new occupation numbers from eigen
1146      call timab(990,1,tsec)
1147      call newocc(doccde,eigen,energies%entropy,energies%e_fermie,dtset%spinmagntarget,&
1148 &     dtset%mband,dtset%nband,dtset%nelect,dtset%nkpt,dtset%nspinor,&
1149 &     dtset%nsppol,occ,dtset%occopt,prtvol,dtset%stmbias,dtset%tphysel,dtset%tsmear,dtset%wtk)
1150      call timab(990,2,tsec)
1151 
1152 !    !=========  DMFT call begin ============================================
1153      dmft_ldaocc=0
1154      if(paw_dmft%use_dmft==1.and.psps%usepaw==1.and.dtset%nbandkss==0) then
1155        call timab(991,1,tsec)
1156 
1157        if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(pawtab(:)%upawu)>=tol8.or.  &
1158 &       sum(pawtab(:)%jpawu)>tol8).and.dtset%dmft_entropy==0) energies%entropy=zero
1159 
1160 !      ==  0 to a dmft calculation and do not use lda occupations
1161 !      ==  1 to a lda calculation with the dmft loop
1162        if(dtset%dmftcheck==-1) dmft_ldaocc=1
1163 
1164 !      ==  initialise occnd
1165        paw_dmft%occnd=zero
1166 
1167        if(dmft_ldaocc==0) then
1168          if(dtset%occopt/=3) then
1169            message = ' occopt should be equal to 3 in dmft'
1170            MSG_ERROR(message)
1171          end if
1172 !        ==  initialise edmft
1173          if(paw_dmft%use_dmft>=1) edmft = zero
1174 
1175          !  Compute residm to check the value
1176          ibdkpt=1
1177          residm=zero
1178          do isppol=1,dtset%nsppol
1179            do ikpt=1,dtset%nkpt
1180              nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1181              nband_eff=max(1,nband_k-dtset%nbdbuf)
1182              residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
1183              ibdkpt=ibdkpt+nband_k
1184            end do
1185          end do
1186 
1187          !  Test residm
1188          if (paw_dmft%use_dmft>0 .and. residm>tol4 .and. dtset%dmftcheck>=0) then
1189            if(dtset%dmft_entropy>0)  then
1190              write(message,'(a,e12.3)')&
1191 &             ' WARNING: Wavefunctions not converged : DFT+DMFT calculation cannot be carried out safely ',residm
1192              call wrtout(std_out,message,'COLL')
1193            else
1194              write(message,'(a,e12.3)')&
1195 &             ' ERROR: Wavefunctions not converged : DFT+DMFT calculation cannot be carried out safely ',residm
1196              call wrtout(std_out,message,'COLL')
1197              write(message,'(a,i0)')'  Action: increase nline and nnsclo',dtset%nstep
1198              MSG_ERROR(message)
1199            end if
1200 
1201          else if (paw_dmft%use_dmft>0 .and. residm>tol10.and. dtset%dmftcheck>=0) then
1202            write(message,'(3a)')ch10,&
1203 &           '  Wavefunctions not converged : DFT+DMFT calculation might not be carried out safely ',ch10
1204            MSG_WARNING(message)
1205          end if
1206 
1207 !        ==  allocate paw_dmft%psichi and paw_dmft%eigen_lda
1208          call init_dmft(dmatpawu,dtset,energies%e_fermie,dtfil%fnameabo_app,dtset%nspinor,paw_dmft,pawtab,psps,dtset%typat)
1209          call print_dmft(paw_dmft,dtset%pawprtvol)
1210 
1211 !        ==  gather crystal structure date into data "cryst_struc"
1212          remove_inv=.false.
1213          if(dtset%nspden==4) remove_inv=.true.
1214          call crystal_init(dtset%amu_orig(:,1),cryst_struc,dtset%spgroup,natom,dtset%npsp,ntypat, &
1215 &         dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
1216 &         dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,&
1217 &         dtset%symrel,dtset%tnons,dtset%symafm)
1218 
1219 !        ==  compute psichi
1220          call xmpi_barrier(spaceComm_distrb)
1221          call init_oper(paw_dmft,lda_occup)
1222          call flush_unit(std_out)
1223          call timab(620,1,tsec)
1224          call datafordmft(cryst_struc,cprj,gs_hamk%dimcprj,dtset,eigen,energies%e_fermie,&
1225 &         lda_occup,dtset%mband,mband_cprj,dtset%mkmem,mpi_enreg,&
1226 &         dtset%nkpt,my_nspinor,dtset%nsppol,occ,&
1227 &         paw_dmft,paw_ij,pawang,pawtab,psps,usecprj_local,dtfil%unpaw)
1228          call timab(620,2,tsec)
1229          call flush_unit(std_out)
1230 
1231 !        ==  solve dmft loop
1232          call xmpi_barrier(spaceComm_distrb)
1233 
1234          call dmft_solve(cryst_struc,istep,lda_occup,paw_dmft,pawang,pawtab,dtset%pawprtvol)
1235          edmft=paw_dmft%edmft
1236          energies%e_paw=energies%e_paw+edmft
1237          energies%e_pawdc=energies%e_pawdc+edmft
1238          call flush_unit(std_out)
1239 !        paw_dmft%occnd(:,:,:,:,:)=0.5_dp
1240 
1241 !        call print_dmft(paw_dmft,dtset%pawprtvol)
1242          if(dtset%paral_kgb==1) then
1243            write(message,'(5a)')ch10,&
1244 &           ' Parallelization over bands is not yet compatible with self-consistency in DMFT ',ch10,&
1245 &           ' Calculation of density does not taken into account non diagonal occupations',ch10
1246            call wrtout(std_out,message,'COLL')
1247            call wrtout(ab_out,message,'COLL')
1248 !          MSG_ERROR(message)
1249            if(dtset%nstep>1) then
1250              write(message,'(a,i0)')'  Action: use nstep=1 instead of nstep=',dtset%nstep
1251              MSG_ERROR(message)
1252            end if
1253            residm=zero
1254          end if
1255          if(dtset%nspinor==2) then
1256 !          call flush_unit(ab_out)
1257 !          write(message,'(3a)')&
1258 !          &         ' Self consistent DFT+DMFT with nspinor==2 is not possible yet ',ch10,&
1259 !          &         ' Calculation are restricted to nstep =1'
1260 !          !         MSG_ERROR(message)
1261 !          if(dtset%nstep>1) then
1262 !          write(message,'(a,i0)')' Action: use nstep=1 instead of nstep=',dtset%nstep
1263 !          !           MSG_ERROR(message)
1264 !          endif
1265          end if
1266 
1267          if(me_distrb==0) then
1268            call saveocc_dmft(paw_dmft)
1269          end if
1270          call destroy_dmft(paw_dmft)
1271 
1272 !        ==  destroy crystal_t cryst_struc
1273          call crystal_free(cryst_struc)
1274          call destroy_oper(lda_occup)
1275        end if ! dmft_ldaocc
1276        call timab(991,2,tsec)
1277 
1278      end if ! usedmft
1279 
1280      if(dtset%nbandkss/=0) then
1281        write(message,'(a,i3,2a,i3,4a)') &
1282 &       " dtset%nbandkss ==",dtset%nbandkss,ch10,&
1283 &       " and dtset%usedmft ==",dtset%usedmft,ch10,&
1284 &       " a DFT loop is carried out without DMFT.",ch10,&
1285 &       " Only psichi's will be written at convergence of the DFT loop."
1286        call wrtout(std_out,message,'COLL')
1287      end if
1288 !    !=========  DMFT call end   ============================================
1289 
1290      call timab(992,1,tsec)
1291 
1292 !    Compute eeig, ek,enl and grnl from the new occ, and the shared eknk,enlnk,grnlnk
1293      energies%e_eigenvalues = zero
1294      energies%e_kinetic     = zero
1295      energies%e_nonlocalpsp = zero
1296      if (usefock) then
1297        energies%e_fock     = zero
1298        if (optforces>0) fock%fock_common%forces=zero
1299      end if
1300      if (optforces>0) grnl(:)=zero
1301      if(paw_dmft%use_dmft>=1) then
1302        ebandlda               = zero
1303        ebanddmft              = zero
1304        ebandldatot            = zero
1305        ekindmft               = zero
1306        ekindmft2              = zero
1307        ekinlda                = zero
1308      end if
1309 
1310 !    Compute new energy terms due to non diagonal occupations and DMFT.
1311      bdtot_index=1
1312      do isppol=1,dtset%nsppol
1313        do ikpt=1,dtset%nkpt
1314          nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1315          do iband=1,nband_k
1316 
1317            locc_test = abs(occ(bdtot_index))>tol8
1318 !          dmft
1319            if(paw_dmft%use_dmft>=1.and.dtset%nbandkss==0) then
1320              if(paw_dmft%band_in(iband)) then
1321                if( paw_dmft%use_dmft == 1 .and. dmft_ldaocc == 1 ) then ! test of the code
1322                  paw_dmft%occnd(1,iband,iband,ikpt,isppol)=occ(bdtot_index)
1323                end if
1324                locc_test = abs(paw_dmft%occnd(1,iband,iband,ikpt,isppol))+&
1325 &               abs(paw_dmft%occnd(2,iband,iband,ikpt,isppol))>tol8
1326              end if
1327            end if
1328 
1329            if (locc_test) then
1330 !            dmft
1331              if(paw_dmft%use_dmft==1.and.dtset%nbandkss==0) then
1332                ebandldatot=ebandldatot+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1333                if(paw_dmft%band_in(iband)) then
1334                  ebandlda=ebandlda+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1335                  ekinlda=ekinlda+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1336                  occ(bdtot_index)=paw_dmft%occnd(1,iband,iband,ikpt,isppol)
1337                  ebanddmft=ebanddmft+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1338                  ekindmft=ekindmft+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1339                end if
1340              end if
1341 
1342              energies%e_eigenvalues = energies%e_eigenvalues + &
1343 &             dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1344              energies%e_kinetic = energies%e_kinetic + &
1345 &             dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1346              energies%e_nonlocalpsp = energies%e_nonlocalpsp + &
1347 &             dtset%wtk(ikpt)*occ(bdtot_index)*enlnk(bdtot_index)
1348              if (usefock) then
1349                energies%e_fock=energies%e_fock + half*focknk(bdtot_index)*occ(bdtot_index)*dtset%wtk(ikpt)
1350                if (optforces>0) fock%fock_common%forces(:,:)=fock%fock_common%forces(:,:)+&
1351 &               dtset%wtk(ikpt)*occ(bdtot_index)*fockfornk(:,:,bdtot_index)
1352              end if
1353              if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ(bdtot_index)*grnlnk(:,bdtot_index)
1354            end if
1355            bdtot_index=bdtot_index+1
1356            if(paw_dmft%use_dmft==1.and.dtset%nbandkss==0) then
1357              do iband1=1,nband_k
1358                if(paw_dmft%band_in(iband).and.paw_dmft%band_in(iband1)) then
1359 !                write(std_out,*) "II+", isppol,ikpt,iband,iband1
1360                  ekindmft2=ekindmft2  +  dtset%wtk(ikpt)*paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*&
1361 &                 eknk_nd(isppol,ikpt,1,iband,iband1)
1362                  if(dtset%nspinor==2) then
1363                    ekindmft2=ekindmft2  -  dtset%wtk(ikpt)*paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*&
1364 &                   eknk_nd(isppol,ikpt,2,iband,iband1)
1365                  end if
1366 !                write(std_out,*) "II", occnd(1,iband,iband1,ikpt,isppol),eknk_nd(isppol,ikpt,iband,iband1)
1367                end if
1368              end do
1369            end if
1370          end do
1371        end do
1372      end do
1373 
1374      if(paw_dmft%use_dmft==1) then
1375        energies%e_kinetic = energies%e_kinetic -ekindmft+ekindmft2
1376        if(abs(dtset%pawprtvol)>=2) then
1377          write(message,'(4a,7(2x,a,2x,e14.7,a),a)') &
1378 &         "-----------------------------------------------",ch10,&
1379 &         "--- Energy for DMFT and tests (in Ha)  ",ch10,&
1380 &         "--- Ebandldatot    (Ha.) = ",ebandldatot,ch10,&
1381 &         "--- Ebandlda       (Ha.) = ",ebandlda,ch10,&
1382 &         "--- Ebanddmft      (Ha.) = ",ebanddmft,ch10,&
1383 &         "--- Ekinlda        (Ha.) = ",ekinlda,ch10, &
1384 &         "--- Ekindmftdiag   (Ha.) = ",ekindmft,ch10,&
1385 &         "--- Ekindmftnondiag(Ha.) = ",ekindmft2,ch10,&
1386 &         "--- Edmft=         (Ha.) = ",edmft,ch10,&
1387 &         "-----------------------------------------------"
1388          call wrtout(std_out,message,'COLL')
1389        end if
1390        if(paw_dmft%use_dmft==1.and.mpi_enreg%paral_kgb==1) paw_dmft%use_dmft=0
1391      end if
1392 
1393      if (psps%usepaw==0) then
1394        call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1395 &       rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1396      else
1397        call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1398 &       rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1399      end if
1400      call timab(992,2,tsec)
1401 
1402 !    Treat fixed occupation numbers or non-self-consistent case
1403    else
1404 
1405      if (mpi_enreg%nproc_kpt>1) then
1406 
1407        call timab(989,1,tsec)
1408 
1409        nbuf=2*mbdkpsp+dtset%nfft*dtset%nspden+3+3*natom*optforces
1410 !      * If Hartree-Fock calculation, the exact exchange energy is k-dependent.
1411        if(dtset%usefock==1) then
1412          nbuf=nbuf+1
1413          if (optforces>0) nbuf=nbuf+3*natom
1414        end if
1415        if(iscf==-1 .or. iscf==-2)nbuf=2*mbdkpsp
1416        ABI_ALLOCATE(buffer1,(nbuf))
1417 !      Pack eigen,resid,rho[wf]r,grnl,enl,ek
1418        buffer1(1:mbdkpsp)=eigen(:)
1419        buffer1(1+mbdkpsp:2*mbdkpsp)=resid(:)
1420        index1=2*mbdkpsp
1421        if(iscf/=-1 .and. iscf/=-2)then
1422          if (psps%usepaw==0) then
1423            buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhor,(/dtset%nfft*dtset%nspden/))
1424          else
1425            buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhowfr,(/dtset%nfft*dtset%nspden/))
1426          end if
1427          index1=index1+dtset%nfft*dtset%nspden
1428          buffer1(index1+1) = energies%e_kinetic
1429          buffer1(index1+2) = energies%e_eigenvalues
1430          buffer1(index1+3) = energies%e_nonlocalpsp
1431          index1=index1+3
1432 !        * If Hartree-Fock calculation, save e_fock in buffer1
1433          if (dtset%usefock==1) then
1434            buffer1(index1+1) = energies%e_fock
1435            index1=index1+1
1436            if (optforces>0)then
1437              buffer1(index1+1:index1+3*natom)=reshape(fock%fock_common%forces,(/3*natom/))
1438              index1=index1+3*natom
1439            end if
1440          end if
1441          if (optforces>0) buffer1(index1+1:index1+3*natom)=grnl(1:3*natom)
1442        end if
1443 
1444 !      Build sum of everything
1445        call timab(48,1,tsec)
1446        call xmpi_sum(buffer1,nbuf,mpi_enreg%comm_kpt ,ierr)
1447        call timab(48,2,tsec)
1448 
1449 !      Unpack the final result
1450        eigen(:)=buffer1(1:mbdkpsp)
1451        resid(:)=buffer1(1+mbdkpsp:2*mbdkpsp)
1452        index1=2*mbdkpsp
1453        if(iscf/=-1 .and. iscf/=-2)then
1454          if (psps%usepaw==0) then
1455            ii=1
1456            do ispden=1,dtset%nspden
1457              do ifft=1,dtset%nfft
1458                rhor(ifft,ispden)=buffer1(index1+ii)
1459                ii=ii+1
1460              end do
1461            end do
1462          else
1463            ii=1
1464            do ispden=1,dtset%nspden
1465              do ifft=1,dtset%nfft
1466                rhowfr(ifft,ispden)=buffer1(index1+ii)
1467                ii=ii+1
1468              end do
1469            end do
1470          end if
1471          index1=index1+dtset%nfft*dtset%nspden
1472          energies%e_kinetic = buffer1(index1+1)
1473          energies%e_eigenvalues = buffer1(index1+2)
1474          energies%e_nonlocalpsp = buffer1(index1+3)
1475          index1=index1+3
1476 !        * If Hartree-Fock calculation, save e_fock in buffer1
1477          if (dtset%usefock==1) then
1478            energies%e_fock = buffer1(index1+1)
1479            index1=index1+1
1480            if (optforces>0) then
1481              fock%fock_common%forces(:,:)=reshape(buffer1(index1+1:index1+3*natom),(/3,natom/))
1482              index1=index1+3*natom
1483            end if
1484          end if
1485          if (optforces>0) grnl(1:3*natom)=buffer1(index1+1:index1+3*natom)
1486        end if
1487        ABI_DEALLOCATE(buffer1)
1488        call timab(989,2,tsec)
1489 
1490      end if ! nproc_kpt>1
1491 
1492 
1493 !    Compute the highest occupied eigenenergy
1494      if(iscf/=-1 .and. iscf/=-2)then
1495        call timab(993,1,tsec)
1496        energies%e_fermie = -huge(one)
1497        bdtot_index=1
1498        do isppol=1,dtset%nsppol
1499          do ikpt=1,dtset%nkpt
1500            nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1501            do iband=1,nband_k
1502              if(abs(occ(bdtot_index))>tol8 .and. eigen(bdtot_index)>energies%e_fermie+tol10) then
1503                energies%e_fermie=eigen(bdtot_index)
1504              end if
1505              bdtot_index=bdtot_index+1
1506            end do
1507          end do
1508        end do
1509        call xmpi_max(energies%e_fermie,spaceComm_distrb,ierr)
1510        call timab(993,2,tsec)
1511      end if
1512 
1513 !    If needed, compute rhog, and symmetrizes the density
1514      if (iscf > 0 .or. iscf==-3 ) then
1515 !      energies%e_fermie=zero  ! Actually, should determine the maximum of the valence band XG20020802
1516        nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)
1517 
1518        call timab(994,1,tsec)
1519        if (psps%usepaw==0) then
1520          call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,&
1521 &         dtset%nsppol,dtset%nsym,dtset%paral_kgb,phnons,rhog  ,rhor  ,rprimd,dtset%symafm,dtset%symrel)
1522        else
1523          call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,&
1524 &         dtset%nsppol,dtset%nsym,dtset%paral_kgb,phnons,rhowfg,rhowfr,rprimd,dtset%symafm,dtset%symrel)
1525        end if
1526        call timab(994,2,tsec)
1527 !      We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2
1528 !      we also have the spin-up density, symmetrized, in rhor(:,2).
1529      end if
1530 
1531    end if !  End of test on varying or fixed occupation numbers
1532 
1533    call timab(994,1,tsec)
1534 
1535 
1536 !  Compute the kinetic energy density
1537    if(dtset%usekden==1 .and. (iscf > 0 .or. iscf==-3 ) )then
1538      call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1539 &     taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1540    end if
1541 
1542    ABI_DEALLOCATE(eknk)
1543    if (usefock) then
1544      ABI_DEALLOCATE(focknk)
1545      if (optforces>0)then
1546        ABI_DEALLOCATE(fockfornk)
1547      end if
1548    end if
1549    ABI_DEALLOCATE(eknk_nd)
1550    ABI_DEALLOCATE(grnlnk)
1551    ABI_DEALLOCATE(enlnk)
1552 
1553 !  In the non-self-consistent case, print eigenvalues and residuals
1554    if(iscf<=0)then
1555      option=2 ; enunit=1 ; vxcavg_dum=zero
1556      call prteigrs(eigen,enunit,energies%e_fermie,dtfil%fnameabo_app_eig,&
1557 &     ab_out,iscf,dtset%kptns,dtset%kptopt,dtset%mband,dtset%nband,&
1558 &     dtset%nkpt,nnsclo_now,dtset%nsppol,occ,dtset%occopt,option,&
1559 &     dtset%prteig,prtvol,resid,dtset%tolwfr,vxcavg_dum,dtset%wtk)
1560    end if
1561 
1562 !  Find largest residual over bands, k points, and spins, except for nbdbuf highest bands
1563    ibdkpt=1
1564    residm=zero
1565    do isppol=1,dtset%nsppol
1566      do ikpt=1,dtset%nkpt
1567        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1568        nband_eff=max(1,nband_k-dtset%nbdbuf)
1569        residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
1570        ibdkpt=ibdkpt+nband_k
1571      end do
1572    end do
1573 
1574  end if !usewvl==0
1575 
1576 !===================================================================
1577 ! End of PLANE WAVES section
1578 !===================================================================
1579 
1580 !In the self-consistent case, diagnose lack of unoccupied state (for each spin and k-point).
1581 !Print a warning if the number of such messages already written does not exceed mwarning.
1582  mwarning=5
1583  if(nwarning<mwarning .and. iscf>=0)then
1584    nwarning=nwarning+1
1585    bdtot_index=1
1586    do isppol=1,dtset%nsppol
1587      do ikpt=1,dtset%nkpt
1588        min_occ=two
1589        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1590        do iband=1,nband_k
1591          if(occ(bdtot_index)<min_occ)min_occ=occ(bdtot_index)
1592          bdtot_index=bdtot_index+1
1593        end do
1594        if(min_occ>0.01_dp)then
1595          if(dtset%nsppol==1)then
1596            write(message, '(a,i0,3a,f7.3,5a)' )&
1597 &           'For k-point number ',ikpt,',',ch10,&
1598 &           'The minimal occupation factor is',min_occ,'.',ch10,&
1599 &           'An adequate monitoring of convergence requires it to be  at most 0.01_dp.',ch10,&
1600 &           'Action: increase slightly the number of bands.'
1601          else
1602            write(message, '(a,i0,3a,i0,a,f7.3,5a)' )&
1603 &           'For k-point number ',ikpt,', and',ch10,&
1604 &           'for spin polarization',isppol,'the minimal occupation factor is',min_occ,'.',ch10,&
1605 &           'An adequate monitoring of convergence requires it to be at most 0.01_dp.',ch10,&
1606 &           'Action: increase slightly the number of bands.'
1607          end if
1608          MSG_WARNING(message)
1609          exit ! It is enough if one lack of adequate occupation is identified, so exit.
1610        end if
1611      end do
1612    end do
1613  end if
1614 
1615  if (iscf>0.or.iscf==-3 .or. (dtset%usewvl==1 .and. iscf==0)) then
1616 
1617 !  PAW: Build new rhoij quantities from new occ then symetrize them
1618 !  Compute and add the compensation density to rhowfr to get the total density
1619    if (psps%usepaw==1) then
1620      call timab(555,1,tsec)
1621      if (paral_atom) then
1622        ABI_DATATYPE_ALLOCATE(pawrhoij_unsym,(natom))
1623        nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
1624        call pawrhoij_alloc(pawrhoij_unsym,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,&
1625 &       dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0)
1626      else
1627        pawrhoij_unsym => pawrhoij
1628      end if
1629      if (usecprj_local==1) then
1630        call pawmkrhoij(atindx,atindx1,cprj,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,&
1631 &       mcprj_local,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,dtset%nspinor,dtset%nsppol,&
1632 &       occ,dtset%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij_unsym,dtfil%unpaw,&
1633 &       dtset%usewvl,dtset%wtk)
1634      else
1635        mcprj_tmp=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
1636        ABI_DATATYPE_ALLOCATE(cprj_tmp,(natom,mcprj_tmp))
1637        call pawcprj_alloc(cprj_tmp,0,gs_hamk%dimcprj)
1638        call ctocprj(atindx,cg,1,cprj_tmp,gmet,gprimd,0,0,0,dtset%istwfk,kg,dtset%kptns,&
1639 &       mcg,mcprj_tmp,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,&
1640 &       dtset%natom,nattyp,dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg,&
1641 &       npwarr,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,&
1642 &       ucvol,dtfil%unpaw,xred,ylm,ylmgr_dum)
1643        call pawmkrhoij(atindx,atindx1,cprj_tmp,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,&
1644 &       dtset%mband,mband_cprj,mcprj_tmp,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,&
1645 &       dtset%nspinor,dtset%nsppol,occ,dtset%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij_unsym,&
1646 &       dtfil%unpaw,dtset%usewvl,dtset%wtk)
1647        call pawcprj_free(cprj_tmp)
1648        ABI_DATATYPE_DEALLOCATE(cprj_tmp)
1649      end if
1650      call timab(555,2,tsec)
1651 !    Build symetrized packed rhoij and compensated pseudo density
1652      cplex=1;ipert=0;idir=0;qpt(:)=zero
1653      if(dtset%usewvl==0) then
1654        call pawmkrho(compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
1655 &       my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
1656 &       dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
1657 &       symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat,rhog=rhog)
1658      else
1659 !      here do not pass rhog, we do not use it
1660        call pawmkrho(compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
1661 &       my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
1662 &       dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
1663 &       symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat)
1664 !      In WVL: copy density to BigDFT object:
1665        call wvl_rho_abi2big(1,rhor,wvl%den)
1666      end if
1667      if (paral_atom) then
1668        call pawrhoij_free(pawrhoij_unsym)
1669        ABI_DATATYPE_DEALLOCATE(pawrhoij_unsym)
1670      end if
1671    end if
1672 
1673    if(paw_dmft%use_dmft==1) then
1674 !    == check noccmmp
1675      call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,dtset%nsppol,0,ntypat,&
1676 &     paw_ij,pawang,dtset%pawprtvol,pawrhoij,pawtab,rdum2,idum1,dtset%typat,0,dtset%usepawu,&
1677 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1678    end if
1679 
1680 !  Find and print minimum and maximum total electron density and locations
1681 !  Compute density residual (if required) and its squared norm
1682    if (iscf>=0) then
1683      if (psps%usepaw==0) then
1684        call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,ucvol=ucvol)
1685      else
1686        call prtrhomxmn(std_out,mpi_enreg,nfftf,pawfgr%ngfft,dtset%nspden,1,rhor,ucvol=ucvol)
1687      end if
1688      if (optres==1) then
1689        nvresid=rhor-nvresid
1690        call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid,mpi_comm_sphgrid=mpi_comm_sphgrid)
1691      end if
1692    end if
1693 
1694  end if ! iscf>0 or iscf=-3
1695 
1696  if(psps%usepaw==1.and.(iscf>=0.or.iscf==-3))  then
1697    ABI_DEALLOCATE(rhowfr)
1698    ABI_DEALLOCATE(rhowfg)
1699  end if
1700 
1701  call timab(994,2,tsec)
1702 
1703  if (iscf==-1) then
1704 !  Eventually compute the excited states within tddft
1705    call timab(995,1,tsec)
1706    if (psps%usepaw==1) then
1707 !    In case of PAW calculation, have to transfer kxc from the fine to the coarse grid:
1708      ABI_ALLOCATE(cgrkxc,(dtset%nfft,nkxc))
1709      do ikxc=1,nkxc
1710        call transgrid(1,mpi_enreg,1,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrkxc(:,ikxc),kxc(:,ikxc))
1711      end do
1712      call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,&
1713 &     kg,cgrkxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,&
1714 &     ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew)
1715      ABI_DEALLOCATE(cgrkxc)
1716    else
1717      call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,&
1718 &     kg,kxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,&
1719 &     ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew)
1720    end if
1721    call timab(995,2,tsec)
1722 
1723  else
1724 !  Eventually compute the susceptibility matrix and the
1725 !  dielectric matrix when istep_mix is equal to 1 or dielstrt
1726    call timab(996,1,tsec)
1727    call testsusmat(computesusmat,dielop,dielstrt,dtset,istep_mix) !test if the matrix is to be computed
1728    if(computesusmat) then
1729      dielar(1)=dtset%diecut;dielar(2)=dtset%dielng
1730      dielar(3)=dtset%diemac;dielar(4)=dtset%diemix
1731      dielar(5)=dtset%diegap;dielar(6)=dtset%dielam
1732      dielar(7)=dtset%diemix;if (iscf>=10) dielar(7)=dtset%diemixmag
1733      usetimerev=1
1734      if (psps%usepaw==1.and.dtset%pawspnorb>0.and.dtset%kptopt/=1.and.dtset%kptopt/=2) usetimerev=0
1735      neglect_pawhat=1-dtset%pawsushat
1736      call suscep_stat(atindx,atindx1,cg,cprj,dielar,&
1737 &     gs_hamk%dimcprj,doccde,eigen,gbound_diel,gprimd,&
1738 &     irrzondiel,dtset%istwfk,kg,kg_diel,lmax_diel,&
1739 &     dtset%mband,mcg,mcprj_local,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,natom,dtset%nband,&
1740 &     neglect_pawhat,nfftdiel,ngfftdiel,&
1741 &     dtset%nkpt,npwarr,npwdiel,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,&
1742 &     occ,dtset%occopt,pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,&
1743 &     susmat,dtset%symafm,dtset%symrel,dtset%tnons,dtset%typat,ucvol,&
1744 &     dtfil%unpaw,usecprj_local,psps%usepaw,usetimerev,dtset%wtk,ylmdiel)
1745    end if
1746    call timab(996,2,tsec)
1747 
1748  end if ! end condition on iscf
1749 
1750  call destroy_hamiltonian(gs_hamk)
1751 
1752  if(dtset%usewvl==0) then
1753    ABI_DEALLOCATE(EigMin)
1754    ABI_DEALLOCATE(doccde)
1755 #if defined HAVE_GPU_CUDA
1756    if(dtset%use_gpu_cuda==1) then
1757      call gpu_finalize_ham_data()
1758    end if
1759 #endif
1760  end if
1761 
1762  call timab(980,2,tsec)
1763 
1764  DBG_EXIT("COLL")
1765 
1766  contains

ABINIT/wvl_comm_eigen [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_comm_eigen

FUNCTION

  Computes occupations for the wavelet case
  Using BigDFT routines

COPYRIGHT

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

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

 for the wvlbigdft case, see the routine 'wvl_occ_bigdft'

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2269 subroutine wvl_comm_eigen()
2270 
2271  use defs_basis
2272  use m_errors
2273 
2274  use m_abi2big, only : wvl_eigen_abi2big
2275 #if defined HAVE_BIGDFT
2276  use BigDFT_API, only: eigensystem_info
2277 #endif
2278 
2279 !This section has been created automatically by the script Abilint (TD).
2280 !Do not modify the following lines by hand.
2281 #undef ABI_FUNC
2282 #define ABI_FUNC 'wvl_comm_eigen'
2283 !End of the abilint section
2284 
2285  implicit none
2286 
2287 !Arguments ------------------------------------
2288 
2289 !Local variables-------------------------------
2290 #if defined HAVE_BIGDFT
2291  integer:: ikpt,norb,shift
2292 #endif
2293 
2294 ! *************************************************************************
2295 
2296    DBG_ENTER("COLL")
2297 
2298 #if defined HAVE_BIGDFT
2299    if(wvlbigdft) then
2300 !  Communicates eigenvalues to all procs.
2301 !  This will print out the eigenvalues and Fermi level.
2302      call eigensystem_info(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl,0.d0,&
2303 &     wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,&
2304 &     wvl%wfs%ks%orbs,wvl%wfs%ks%psi)
2305    else
2306 !  Send all eigenvalues to all procs.
2307 !  I simply communicate eigenvalues: I do not print them into screen, nor calculate Fermi-level.
2308      norb=wvl%wfs%ks%orbs%norb
2309      if (mpi_enreg%nproc_wvl > 1) then
2310        shift=1
2311        do ikpt = 1, wvl%wfs%ks%orbs%nkpts
2312          call xmpi_bcast(wvl%wfs%ks%orbs%eval(shift:shift+norb-1),wvl%wfs%ks%orbs%ikptproc(ikpt),mpi_enreg%comm_wvl,ierr)
2313          shift=shift+norb
2314        end do
2315      end if
2316    end if
2317 
2318 !Copy eigenvalues from BigDFT object to "eigen"
2319    call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs)
2320 
2321 #else
2322    BIGDFT_NOTENABLED_ERROR()
2323 #endif
2324 
2325    DBG_EXIT("COLL")
2326 
2327  end subroutine wvl_comm_eigen
2328 
2329 end subroutine vtorho

ABINIT/wvl_nscf_loop [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_nscf_loop

FUNCTION

  Non-self-consistent field cycle in Wavelets
  See also "wvl_nscf_loop_bigdft"

COPYRIGHT

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

INPUTS

  nnsclo= number of non-self consistent field iterations

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

1800 subroutine wvl_nscf_loop()
1801 
1802  use defs_basis
1803  use defs_abitypes
1804  use defs_wvltypes
1805 
1806  use m_errors
1807  use m_energies, only : energies_type
1808  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
1809 
1810 !This section has been created automatically by the script Abilint (TD).
1811 !Do not modify the following lines by hand.
1812 #undef ABI_FUNC
1813 #define ABI_FUNC 'wvl_nscf_loop'
1814  use interfaces_62_wvl_wfs
1815 !End of the abilint section
1816 
1817  implicit none
1818 
1819 !Arguments ------------------------------------
1820 ! integer, intent(in)                    :: istep,mcprj,nfft,nnsclo
1821 ! real(dp), intent(inout)                :: residm
1822 ! type(dataset_type), intent(in)         :: dtset
1823 ! type(MPI_type), intent(in)             :: mpi_enreg
1824 ! type(energies_type), intent(inout)     :: energies
1825 ! type(wvl_data), intent(inout)          :: wvl
1826 ! !arrays
1827 ! real(dp), intent(inout)                :: xcart(3, dtset%natom)
1828 ! real(dp), dimension(6), intent(out)    :: strsxc
1829 ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj
1830 
1831 !Local variables-------------------------------
1832  integer :: inonsc,ii
1833  integer,parameter :: iscf_=-1       !do not do a SCF cycle
1834  logical,parameter :: do_scf=.false. !do not do a SCF cycle
1835  logical,parameter :: wvlbigdft=.false.
1836  real(dp) :: dum,eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc
1837 
1838 ! *************************************************************************
1839 
1840    DBG_ENTER("COLL")
1841 
1842    if(nnsclo_now>0) then
1843      do inonsc=1,nnsclo_now
1844        call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
1845 &       istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,&
1846 &       nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
1847 &       dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
1848        call wvl_hpsitopsi(cprj,dtset,energies,inonsc,mcprj_local,mpi_enreg, &
1849 &       residm,wvl,xcart)
1850        if(residm<dtset%tolwfr) exit !Exit loop if converged
1851      end do
1852 
1853    else
1854      do ii=1, dtset%nline
1855 !      Direct minimization technique: no diagonalization
1856        call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
1857 &       istep,ii,iscf_,mpi_enreg%me_wvl,dtset%natom,&
1858 &       nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
1859 &       dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
1860        call wvl_hpsitopsi(cprj,dtset,energies,ii,mcprj_local,mpi_enreg, &
1861 &       residm,wvl,xcart)
1862        if(residm<dtset%tolwfr) exit !Exit loop if converged
1863      end do
1864    end if
1865 
1866 !  Update energies depending on new WF
1867    energies%e_kinetic=ekin
1868    energies%e_nonlocalpsp=enl
1869    energies%e_exactX=eexctx
1870    energies%e_sicdc=esicdc
1871 
1872 !  Eventually update energies depending on density
1873    if (dtset%iscf<10) then
1874      energies%e_localpsp=eloc
1875      energies%e_hartree=eh
1876      energies%e_xc=exc ; energies%e_xcdc=evxc
1877    else if (nnsclo_now==0) then
1878      energies%e_localpsp=eloc
1879    end if
1880 
1881 !  End of nscf iterations
1882    if (do_last_ortho) then
1883 !    !Don't update energies (nscf cycle has been done); just recompute potential
1884      inonsc=nnsclo_now;if (nnsclo_now==0) inonsc=dtset%nline
1885      call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
1886 &     istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,&
1887 &     nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
1888 &     dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
1889    end if
1890 
1891    DBG_EXIT("COLL")
1892 
1893  end subroutine wvl_nscf_loop

ABINIT/wvl_nscf_loop_bigdft [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_nscf_loop_bigdft

FUNCTION

  Non-self-consistent field cycle in Wavelets
  It follows the BigDFT scheme.
  See also "wvl_nscf_loop"

COPYRIGHT

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

INPUTS

  nnsclo= number of non-self consistent field iterations

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

1929 #if defined HAVE_CONFIG_H
1930 #include "config.h"
1931 #endif
1932 
1933 #include "abi_common.h"
1934 
1935 subroutine wvl_nscf_loop_bigdft()
1936 
1937  use defs_basis
1938  use defs_abitypes
1939  use defs_wvltypes
1940  use m_errors
1941 
1942  use m_energies, only : energies_type
1943  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
1944 
1945 !This section has been created automatically by the script Abilint (TD).
1946 !Do not modify the following lines by hand.
1947 #undef ABI_FUNC
1948 #define ABI_FUNC 'wvl_nscf_loop_bigdft'
1949  use interfaces_62_wvl_wfs
1950 !End of the abilint section
1951 
1952  implicit none
1953 
1954 !Arguments ------------------------------------
1955 ! integer, intent(in)                    :: istep,mcprj,nfft,nnsclo
1956 ! real(dp), intent(inout)                :: residm
1957 ! real(dp), intent(out)                  :: nres2
1958 ! type(dataset_type), intent(in)         :: dtset
1959 ! type(MPI_type), intent(in)             :: mpi_enreg
1960 ! type(energies_type), intent(inout)     :: energies
1961 ! type(wvl_data), intent(inout)          :: wvl
1962  !arrays
1963 ! real(dp), intent(inout)                :: xcart(3, dtset%natom)
1964 ! real(dp), dimension(6), intent(out)    :: strsxc
1965 ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj
1966 
1967 !Local variables-------------------------------
1968  integer :: inonsc
1969  integer,parameter :: iscf_=-1       !do not do a SCF cycle
1970  logical,parameter :: do_scf=.false. !do not do a SCF cycle
1971  logical,parameter :: wvlbigdft=.true.
1972  real(dp) :: eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc
1973 
1974 ! *************************************************************************
1975 
1976    DBG_ENTER("COLL")
1977 
1978    call wvl_hpsitopsi(cprj,dtset, energies, istep, mcprj_local,mpi_enreg, &
1979 &   residm, wvl,xcart)
1980 
1981    if (nnsclo_now>2) then
1982      do inonsc = 2, nnsclo_now-1
1983        call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, &
1984 &       energies%e_hartree, energies%e_kinetic, energies%e_localpsp, &
1985 &       energies%e_nonlocalpsp, energies%e_sicdc, istep, inonsc, iscf_, &
1986 &       mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,&
1987 &       dtset%nspden, nres2, do_scf,energies%e_xcdc, &
1988 &       wvl, wvlbigdft, xcart, strsxc)
1989        call wvl_hpsitopsi(cprj,dtset, energies, inonsc, mcprj_local,mpi_enreg, &
1990 &       residm, wvl,xcart)
1991      end do
1992    end if
1993 
1994 !  End of nscf iterations
1995    if (do_last_ortho.and.nnsclo_now<=1) then
1996 !    !Don't update energies (nscf cycle has been done); just recompute potential
1997      call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc, &
1998 &     istep, 1, iscf_, mpi_enreg%me_wvl, dtset%natom, nfftf, &
1999 &     mpi_enreg%nproc_wvl,dtset%nspden, nres2, do_scf,evxc, &
2000 &     wvl, wvlbigdft, xcart, strsxc)
2001    else if (do_last_ortho.and.nnsclo_now>1) then
2002 !    !Update energies and potential (nscf cycles are not finished)
2003      call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, &
2004 &     energies%e_hartree,energies%e_kinetic, energies%e_localpsp, &
2005 &     energies%e_nonlocalpsp, energies%e_sicdc, istep, nnsclo_now, iscf_, &
2006 &     mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,&
2007 &     dtset%nspden, nres2, do_scf,energies%e_xcdc, &
2008 &     wvl, wvlbigdft, xcart, strsxc)
2009    end if
2010 
2011    DBG_EXIT("COLL")
2012 
2013  end subroutine wvl_nscf_loop_bigdft

ABINIT/wvl_occ [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_occ

FUNCTION

  Computes occupations for the wavelet case

COPYRIGHT

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

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

 for the wvlbigdft case, see the routine 'wvl_occ_bigdft'

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2130 subroutine wvl_occ()
2131 
2132  use defs_basis
2133  use m_errors
2134  use defs_wvltypes
2135 
2136 !This section has been created automatically by the script Abilint (TD).
2137 !Do not modify the following lines by hand.
2138 #undef ABI_FUNC
2139 #define ABI_FUNC 'wvl_occ'
2140 !End of the abilint section
2141 
2142  implicit none
2143 
2144 !Local variables-------------------------------
2145  real(dp):: doccde_(dtset%mband*dtset%nkpt*dtset%nsppol)
2146 ! *************************************************************************
2147 
2148    DBG_ENTER("COLL")
2149 
2150 !  Compute the new occupation numbers from eigen
2151    call newocc(doccde_,eigen,energies%entropy,energies%e_fermie,dtset%spinmagntarget,&
2152 &   dtset%mband,dtset%nband,dtset%nelect,dtset%nkpt,dtset%nspinor,&
2153 &   dtset%nsppol,occ,dtset%occopt,prtvol,&
2154 &   dtset%stmbias,dtset%tphysel,dtset%tsmear,dtset%wtk)
2155 
2156 ! Copy occupations and efermi to BigDFT variables
2157    call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs)
2158 
2159 #if defined HAVE_BIGDFT
2160 !  Copy Fermi level to BigDFT variable:
2161    wvl%wfs%ks%orbs%efermi=energies%e_fermie
2162 #endif
2163 
2164    DBG_EXIT("COLL")
2165 
2166  end subroutine wvl_occ

ABINIT/wvl_occ_bigdft [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_occ_bigdft

FUNCTION

  Computes occupations for the wavelet case
  Using BigDFT routines

COPYRIGHT

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

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

 for the wvlbigdft case, see the routine 'wvl_occ_bigdft'

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2200 subroutine wvl_occ_bigdft()
2201 
2202  use defs_basis
2203  use m_errors
2204 
2205 !This section has been created automatically by the script Abilint (TD).
2206 !Do not modify the following lines by hand.
2207 #undef ABI_FUNC
2208 #define ABI_FUNC 'wvl_occ_bigdft'
2209 !End of the abilint section
2210 
2211  implicit none
2212 
2213 ! *************************************************************************
2214 
2215    DBG_ENTER("COLL")
2216 
2217 ! Transfer occopt from ABINIT to BigDFT
2218 #if defined HAVE_BIGDFT
2219    occopt_bigdft=dtset%occopt
2220    call wvl_occopt_abi2big(occopt_bigdft,occopt_bigdft,1)
2221 
2222 !Calculate occupations using BigDFT routine
2223    call evaltoocc(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, .false., &
2224 &   dtset%tsmear, wvl%wfs%ks%orbs,  occopt_bigdft)
2225 
2226 !Pass e_fermi from BigDFT object to ABINIT variable:
2227    energies%e_fermie = wvl%wfs%ks%orbs%efermi
2228 
2229 !Copy occupations from BigDFT to ABINIT variables
2230    call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs)
2231 #endif
2232 
2233    DBG_EXIT("COLL")
2234 
2235  end subroutine wvl_occ_bigdft