TABLE OF CONTENTS


ABINIT/etotfor [ Functions ]

[ Top ] [ Functions ]

NAME

 etotfor

FUNCTION

 This routine is called to compute the total energy and various parts of it.
 The routine computes -if requested- the forces.

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  dtefield <type(efield_type)> = variables related to Berry phase
  dtset <type(dataset_type)>=all input variables in this dataset
   | berryopt  = 4: electric field is on -> add the contribution of the
   |                - \Omega E.P term to the total energy
   |          /= 4: electric field is off
   | bfield = cartesian coordinates of magnetic field in atomic units
   | efield = cartesian coordinates of the electric field in atomic units
   | iatfix(3,natom)=1 for frozen atom along some direction, 0 for unfrozen
   | ionmov=governs the movement of atoms (see help file)
   | densfor_pred=governs the mixed electronic-atomic part of the preconditioner
   | natom=number of atoms in cell.
   | nconeq=number of atomic constraint equations
   | nspden=number of spin-density components
   | nsym=number of symmetry elements in space group
   | occopt=option for occupancies
   | prtvol=integer controlling volume of printed output
   | tsmear=smearing energy or temperature (if metal)
   | typat(natom)=type integer for each atom in cell
   | wtatcon(3,natom,nconeq)=weights for atomic constraints
  gmet(3,3)=metric tensor for G vecs (in bohr**-2)
  fock <type(fock_type)>= quantities to calculate Fock exact exchange
  grchempottn(3,natom)=grads of spatially-varying chemical potential energy (hartree)
  grewtn(3,natom)=grads of Ewald energy (hartree)
  grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree)
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  nattyp(ntypat)=number of atoms of each type
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
  nhat(nfft,nspden*usepaw)= -PAW only- compensation density
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  ntypat=number of types of atoms in unit cell.
  nvresid(nfft,nspden)=potential or density residual
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  n3xccc=dimension of the xccc3d array (0 or nfft).
  optene=option for the computation of total energy
         (-1=no computation; 0=direct scheme; 1=double-counting scheme)
  optforces=option for the computation of forces
  optres=0 if residual array (nvresid) contains the potential residual
        =1 if residual array (nvresid) contains the density residual
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information.
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhog(2,nfft)=array for Fourier transform of electron density
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symrec(3,3,nsym)=symmetry operations in reciprocal space
  ucvol=unit cell volume
  usepaw= 0 for non paw calculation; =1 for paw calculation
  vhartr(nfft)=array for holding Hartree potential
  vpsp(nfft)=array for holding local psp
  vxc(nfft,nspden)=array for holding XC potential
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  deltae=change in total energy between the previous and present SCF cycle
  etotal=total energy (hartree)
  ===== if optforces==1
   diffor=maximum absolute change in component of forces between present and previous SCF cycle.
   favg(3)=mean of fcart before correction for translational symmetry
   fcart(3,natom)=cartesian forces from fred (hartree/bohr)
   fred(3,natom)=symmetrized form of grtn (grads of Etot) (hartree)
   gresid(3,natom)=forces due to the residual of the density/potential
   grhf(3,natom)=Hellman-Feynman derivatives of the total energy
   grxc(3,natom)=d(Exc)/d(xred) derivatives (0 without core charges)
   maxfor=maximum absolute value of force
   synlgr(3,natom)=symmetrized form of grads of Enl (hartree)

SIDE EFFECTS

 Input/Output:
  elast=previous value of the energy,
        needed to compute deltae, then updated.
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  energies <type(energies_type)>=all part of total energy.
   | entropy(IN)=entropy due to the occupation number smearing (if metal)
   | e_localpsp(IN)=local psp energy (hartree)
   | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree)
   | e_chempot(IN)=energy from spatially varying chemical potential (hartree)
   | e_ewald(IN)=Ewald energy (hartree)
   | e_vdw_dftd(IN)=VdW DFT-D energy
   | e_hartree(IN)=Hartree part of total energy (hartree units)
   | e_corepsp(IN)=psp core-core energy
   | e_hybcomp_E0(IN)=energy compensation energy for the hybrid functionals at frozen density
   | e_hybcomp_v0(IN)=potential compensation energy for the hybrid functionals at frozen density
   | e_hybcomp_v (IN)=potential compensation energy for the hybrid functionals at self-consistent density
   | e_kinetic(IN)=kinetic energy part of total energy.
   | e_nonlocalpsp(IN)=nonlocal pseudopotential part of total energy.
   | e_xc(IN)=exchange-correlation energy (hartree)
   | e_xcdc(IN)=exchange-correlation double-counting energy (hartree)
   | e_paw(IN)=PAW spherical part energy
   | e_pawdc(IN)=PAW spherical part double-counting energy
   | e_elecfield(OUT)=the term of the energy functional that depends explicitely
   |                  on the electric field:  enefield = -ucvol*E*P
   | e_magfield(OUT)=the term of the energy functional that depends explicitely
   |                  on the magnetic field:  e_magfield = -ucvol*E*P
   | e_entropy(OUT)=entropy energy due to the occupation number smearing (if metal)
   |                this value is %entropy * dtset%tsmear (hartree).
  ===== if optforces==1
   forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
   grnl(3*natom)=gradients of Etot due to nonlocal contributions
                 Input for norm-conserving psps, output for PAW
  ===== if psps%usepaw==1
   pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
    (gradients of rhoij for each atom with respect to atomic positions are computed here)

NOTES

  In case of PAW calculations:
    All computations are done on the fine FFT grid.
    All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
    All arrays (densities/potentials...) are computed on this fine FFT grid.
  ! Developpers have to be careful when introducing others arrays:
      they have to be stored on the fine FFT grid.
  In case of norm-conserving calculations the FFT grid is the usual FFT grid.

PARENTS

      scfcv_core

CHILDREN

      forces,nres2vres,pawgrnl,timab

SOURCE

2426 subroutine etotfor(atindx1,deltae,diffor,dtefield,dtset,&
2427 &  elast,electronpositron,energies,&
2428 &  etotal,favg,fcart,fock,forold,fred,gmet,grchempottn,gresid,grewtn,grhf,grnl,grvdw,&
2429 &  grxc,gsqcut,indsym,kxc,maxfor,mgfft,mpi_enreg,my_natom,nattyp,&
2430 &  nfft,ngfft,ngrvdw,nhat,nkxc,ntypat,nvresid,n1xccc,n3xccc,optene,optforces,optres,&
2431 &  pawang,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,red_ptot,psps,rhog,rhor,rmet,rprimd,&
2432 &  symrec,synlgr,ucvol,usepaw,vhartr,vpsp,vxc,wvl,wvl_den,xccc3d,xred)
2433 
2434 
2435 !This section has been created automatically by the script Abilint (TD).
2436 !Do not modify the following lines by hand.
2437 #undef ABI_FUNC
2438 #define ABI_FUNC 'etotfor'
2439 !End of the abilint section
2440 
2441  implicit none
2442 
2443 !Arguments ------------------------------------
2444 !scalars
2445  integer,intent(in) :: my_natom,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nkxc,ntypat,optene,optforces
2446  integer,intent(in) :: optres,usepaw
2447  real(dp),intent(in) :: gsqcut
2448  real(dp),intent(inout) :: elast,ucvol
2449  real(dp),intent(out) :: deltae,diffor,etotal,maxfor
2450  type(MPI_type),intent(in) :: mpi_enreg
2451  type(efield_type),intent(in) :: dtefield
2452  type(dataset_type),intent(in) :: dtset
2453  type(electronpositron_type),pointer :: electronpositron
2454  type(energies_type),intent(inout) :: energies
2455  type(pawang_type),intent(in) :: pawang
2456  type(pseudopotential_type),intent(in) :: psps
2457  type(wvl_internal_type), intent(in) :: wvl
2458  type(wvl_denspot_type), intent(inout) :: wvl_den
2459  type(fock_type),pointer, intent(inout) :: fock
2460 !arrays
2461  integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
2462  integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym)
2463  real(dp),intent(in) :: gmet(3,3),grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc)
2464  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),red_ptot(3)
2465  real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rmet(3,3),rprimd(3,3)
2466  real(dp),intent(in) :: vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden)
2467  real(dp),intent(in) :: xccc3d(n3xccc)
2468  real(dp),intent(inout) :: forold(3,dtset%natom),grnl(3*dtset%natom)
2469  real(dp),intent(inout) :: nhat(nfft,dtset%nspden*psps%usepaw)
2470  real(dp),intent(inout),target :: nvresid(nfft,dtset%nspden)
2471  real(dp),intent(inout) :: xred(3,dtset%natom)
2472  real(dp),intent(out) :: favg(3),fred(3,dtset%natom)
2473  real(dp),intent(inout) :: fcart(3,dtset%natom)
2474  real(dp),intent(out) :: gresid(3,dtset%natom),grhf(3,dtset%natom)
2475  real(dp),intent(inout) :: grxc(3,dtset%natom)
2476  real(dp),intent(out) :: synlgr(3,dtset%natom)
2477  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
2478  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
2479  type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw)
2480  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
2481 
2482 !Local variables-------------------------------
2483 !scalars
2484  integer :: comm_grid,dimnhat,ifft,ipositron,ispden,optgr,optgr2,option,optnc,optstr,optstr2,iir,jjr,kkr
2485  logical :: apply_residual
2486  real(dp) :: eenth,ucvol_
2487 !arrays
2488  real(dp),parameter :: k0(3)=(/zero,zero,zero/)
2489  real(dp) :: tsec(2),A(3,3),A1(3,3),A_new(3,3),efield_new(3)
2490  real(dp) :: dummy(0),nhat_dum(0,0)
2491  real(dp),allocatable :: vlocal(:,:)
2492  real(dp), ABI_CONTIGUOUS pointer :: resid(:,:)
2493 
2494 ! *********************************************************************
2495 
2496  call timab(80,1,tsec)
2497 
2498  ipositron=electronpositron_calctype(electronpositron)
2499 
2500  if (optene>-1) then
2501 
2502 !  When the finite-temperature VG broadening scheme is used,
2503 !  the total entropy contribution "tsmear*entropy" has a meaning,
2504 !  and gather the two last terms of Eq.8 of VG paper
2505 !  Warning : might have to be changed for fixed moment calculations
2506    if(dtset%occopt>=3 .and. dtset%occopt<=8) then
2507      if (abs(dtset%tphysel) < tol10) then
2508        energies%e_entropy = - dtset%tsmear * energies%entropy
2509      else
2510        energies%e_entropy = - dtset%tphysel * energies%entropy
2511      end if
2512    else
2513      energies%e_entropy = zero
2514    end if
2515 
2516 !  Turn it into an electric enthalpy, refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]],
2517 !    the missing volume is added here
2518    energies%e_elecfield=zero
2519    if ((dtset%berryopt==4.or.dtset%berryopt==14).and.ipositron/=1) then
2520      energies%e_elecfield=-dot_product(dtset%red_efieldbar,red_ptot)  !!ebar_i p_i
2521      eenth=zero
2522      do iir=1,3
2523        do jjr=1,3
2524          eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)  !!g^{-1})_ij ebar_i ebar_j
2525        end do
2526      end do
2527      energies%e_elecfield=energies%e_elecfield-eenth*ucvol/(8._dp*pi)
2528    end if
2529 
2530 !  Turn it into an internal energy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]],
2531 !    but a little different: U=E_ks + (vol/8*pi) *  g^{-1})_ij ebar_i ebar_j
2532    if ((dtset%berryopt==6.or.dtset%berryopt==16).and.ipositron/=1) then
2533      energies%e_elecfield=zero
2534      eenth=zero
2535      do iir=1,3
2536        do jjr=1,3
2537          eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)  !! g^{-1})_ij ebar_i ebar_j
2538        end do
2539      end do
2540      energies%e_elecfield=energies%e_elecfield+eenth*ucvol/(8._dp*pi)
2541    end if
2542 
2543 !  Calculate internal energy and electric enthalpy for mixed BC case.
2544    if (dtset%berryopt==17.and.ipositron/=1) then
2545      energies%e_elecfield=zero
2546      A(:,:)=(four_pi/ucvol)*rmet(:,:)
2547      A1(:,:)=A(:,:) ; A_new(:,:)=A(:,:)
2548      efield_new(:)=dtset%red_efield(:)
2549      eenth=zero
2550      do kkr=1,3
2551        if (dtset%jfielddir(kkr)==1) then    ! fixed ebar direction
2552 !        step 1 add -ebar*p
2553          eenth=eenth-dtset%red_efieldbar(kkr)*red_ptot(kkr)
2554 !        step 2  chang to e_new (change e to ebar)
2555          efield_new(kkr)=dtset%red_efieldbar(kkr)
2556 !        step 3  chang matrix A to A1
2557          do iir=1,3
2558            do jjr=1,3
2559              if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr)
2560              if ((iir==kkr .and. jjr/=kkr) .or.  (iir/=kkr .and.  jjr==kkr)) &
2561 &             A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr)
2562              if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr)
2563            end do
2564          end do
2565          A(:,:)=A1(:,:) ; A_new(:,:)=A1(:,:)
2566        end if
2567      end do  ! end fo kkr
2568      do iir=1,3
2569        do jjr=1,3
2570          eenth= eenth+half*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr)
2571        end do
2572      end do
2573      energies%e_elecfield=energies%e_elecfield+eenth
2574    end if   ! berryopt==17
2575 
2576 !  Turn it into a magnetic enthalpy, by adding orbital electronic contribution
2577    energies%e_magfield = zero
2578 !  if (dtset%berryopt == 5 .and. ipositron/=1) then
2579 !  emag = dot_product(mag_cart,dtset%bfield)
2580 !  energies%e_magfield = emag
2581 !  end if
2582 
2583 !  Compute total (free)- energy by direct scheme
2584    if (optene==0) then
2585      etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc + &
2586 &     energies%e_localpsp + energies%e_corepsp + energies%e_fock+&
2587 &     energies%e_entropy + energies%e_elecfield + energies%e_magfield+&
2588 &     energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_hybcomp_v
2589      etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
2590      if (usepaw==0) etotal = etotal + energies%e_nonlocalpsp
2591      if (usepaw/=0) etotal = etotal + energies%e_paw
2592    end if
2593 
2594 !  Compute total (free) energy by double-counting scheme
2595    if (optene==1) then
2596      etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc &
2597 &     - energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc+ energies%e_fock- energies%e_fockdc &
2598 &     + energies%e_entropy + energies%e_elecfield + energies%e_magfield &
2599 &     + energies%e_hybcomp_E0 - energies%e_hybcomp_v0
2600      etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
2601      if (usepaw/=0) etotal = etotal + energies%e_pawdc
2602    end if
2603 
2604 !  Additional stuff for electron-positron
2605    if (dtset%positron/=0) then
2606      if (ipositron==0) then
2607        energies%e_electronpositron  =zero
2608        energies%edc_electronpositron=zero
2609      else
2610        energies%e_electronpositron  =electronpositron%e_hartree+electronpositron%e_xc
2611        energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc
2612        if (usepaw==1) then
2613          energies%e_electronpositron  =energies%e_electronpositron  +electronpositron%e_paw
2614          energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc
2615        end if
2616      end if
2617      if (optene==0) electronpositron%e0=etotal
2618      if (optene==1) electronpositron%e0=etotal-energies%edc_electronpositron
2619      etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron
2620    end if
2621 
2622 !  Compute energy residual
2623    deltae=etotal-elast
2624    elast=etotal
2625  end if !optene/=-1
2626 
2627  call timab(80,2,tsec)
2628 
2629 !------Compute forces-----------------------------------------------------
2630 
2631  if (optforces==1) then
2632 
2633 !  PAW: add gradients due to Dij derivatives to non-local term
2634    if (usepaw==1) then
2635      ABI_ALLOCATE(vlocal,(nfft,dtset%nspden))
2636      do ispden=1,min(dtset%nspden,2)
2637 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vlocal,vpsp,vxc)
2638        do ifft=1,nfft
2639          vlocal(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)
2640        end do
2641      end do
2642 
2643      if(dtset%nspden==4)then
2644        do ispden=3,4
2645 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vlocal,vxc)
2646          do ifft=1,nfft
2647            vlocal(ifft,ispden)=vxc(ifft,ispden)
2648          end do
2649        end do
2650      end if
2651      ucvol_=ucvol
2652 #if defined HAVE_BIGDFT
2653      if (dtset%usewvl==1) ucvol_=product(wvl_den%denspot%dpbox%hgrids)*real(product(wvl_den%denspot%dpbox%ndims),dp)
2654 #endif
2655      dimnhat=0;optgr=1;optgr2=0;optstr=0;optstr2=0
2656      comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl
2657      call pawgrnl(atindx1,dimnhat,dummy,1,dummy,grnl,gsqcut,mgfft,my_natom, &
2658 &     dtset%natom, nattyp,nfft,ngfft,nhat_dum,dummy,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,&
2659 &     pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred, &
2660 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=mpi_enreg%comm_fft,&
2661 &     comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,paral_kgb=mpi_enreg%paral_kgb)
2662      ABI_DEALLOCATE(vlocal)
2663    end if
2664 
2665    apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. &
2666 &   abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5)
2667 
2668 !  If residual is a density residual (and forces from residual asked),
2669 !  has to convert it into a potential residual before calling forces routine
2670    if (apply_residual) then
2671      ABI_ALLOCATE(resid,(nfft,dtset%nspden))
2672      option=0; if (dtset%densfor_pred<0) option=1
2673      optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2
2674      call nres2vres(dtset,gsqcut,usepaw,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
2675 &     nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,&
2676 &     rhor,rprimd,usepaw,resid,xccc3d,xred,vxc)
2677    else
2678      resid => nvresid
2679    end if
2680    call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,fred,grchempottn,gresid,grewtn,&
2681 &   grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfft,mpi_enreg,&
2682 &   n1xccc,n3xccc,nattyp,nfft,ngfft,ngrvdw,ntypat,&
2683 &   pawrad,pawtab,ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,dtset%usefock,resid,vxc,wvl,wvl_den,xred,&
2684 &   electronpositron=electronpositron)
2685    if (apply_residual) then
2686      ABI_DEALLOCATE(resid)
2687    end if
2688 
2689 !  Returned fred are full symmetrized gradients of Etotal
2690 !  wrt reduced coordinates xred, d(Etotal)/d(xred)
2691 !  Forces are contained in array fcart
2692 
2693  else   ! if optforces==0
2694    fcart=zero
2695    fred=zero
2696    favg=zero
2697    diffor=zero
2698    gresid=zero
2699    grhf=zero
2700    maxfor=zero
2701    synlgr=zero
2702  end if
2703 
2704  call timab(80,2,tsec)
2705 
2706 end subroutine etotfor

ABINIT/m_scfcv_core [ Modules ]

[ Top ] [ Modules ]

NAME

  m_scfcv_core

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

 21 #if defined HAVE_CONFIG_H
 22 #include "config.h"
 23 #endif
 24 
 25 #include "abi_common.h"
 26 
 27 module m_scfcv_core
 28 
 29  use defs_basis
 30  use defs_datatypes
 31  use defs_abitypes
 32  use defs_wvltypes
 33  use defs_rectypes
 34  use m_xmpi
 35  use m_abicore
 36  use m_wffile
 37  use m_rec
 38  use m_ab7_mixing
 39  use m_errors
 40  use m_efield
 41  use mod_prc_memory
 42  use m_nctk
 43  use m_hdr
 44  use m_xcdata
 45  use m_cgtools
 46 
 47  use m_berryphase_new,   only : update_e_field_vars
 48  use m_dtfil,            only : status
 49  use m_time,             only : timab
 50  use m_fstrings,         only : int2char4, sjoin
 51  use m_symtk,            only : symmetrize_xred
 52  use m_geometry,         only : metric
 53  use m_fftcore,          only : getng, sphereboundary
 54  use m_time,             only : abi_wtime, sec2str
 55  use m_exit,             only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in
 56  use m_mpinfo,           only : destroy_mpi_enreg, iwrite_fftdatar, initmpi_seq, proc_distrb_cycle
 57  use m_ioarr,            only : fftdatar_write_from_hdr
 58  use m_results_gs ,      only : results_gs_type
 59  use m_scf_history,      only : scf_history_type, scf_history_init, scf_history_free
 60  use m_energies,         only : energies_type, energies_init, energies_copy
 61  use m_electronpositron, only : electronpositron_type, electronpositron_calctype
 62  use m_pawang,           only : pawang_type
 63  use m_pawrad,           only : pawrad_type
 64  use m_pawtab,           only : pawtab_type,pawtab_get_lsize
 65  use m_paw_an,           only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags
 66  use m_pawfgrtab,        only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
 67  use m_pawrhoij,         only : pawrhoij_type
 68  use m_pawcprj,          only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, &
 69 &                               pawcprj_free, pawcprj_axpby, pawcprj_put, pawcprj_getdim, pawcprj_reorder
 70  use m_pawdij,           only : pawdij, symdij,pawdijhat
 71  use m_pawfgr,           only : pawfgr_type
 72  use m_paw_ij,           only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
 73  use m_paw_dmft,         only : paw_dmft_type
 74  use m_paw_nhat,         only : nhatgrid,wvl_nhatgrid,pawmknhat
 75  use m_paw_tools,        only : chkpawovlp
 76  use m_paw_denpot,       only : pawdenpot
 77  use m_paw_occupancies,  only : pawmkrhoij
 78  use m_paw_correlations, only : setnoccmmp,setrhoijpbe0
 79  use m_orbmag,           only : orbmag_type
 80  use m_paw_mkrho,        only : pawmkrho
 81  use m_paw_uj,           only : pawuj_red
 82  use m_paw_dfpt,         only : pawgrnl
 83  use m_fock,             only : fock_type, fock_init, fock_destroy, fock_ACE_destroy, fock_common_destroy, &
 84                                 fock_BZ_destroy, fock_update_exc, fock_updatecwaveocc
 85  use m_gwls_hamiltonian, only : build_vxc
 86 #if defined HAVE_BIGDFT
 87  use BigDFT_API,         only : cprj_clean,cprj_paw_alloc
 88 #endif
 89  use m_io_kss,           only : gshgg_mkncwrite
 90  use m_outxml,           only : out_resultsgs_XML, out_geometry_XML
 91  use m_kg,               only : getcut, getmpw, kpgio, getph
 92  use m_fft,              only : fourdp
 93  use m_vtorhorec,        only : first_rec, vtorhorec
 94  use m_vtorhotf,         only : vtorhotf
 95  use m_outscfcv,         only : outscfcv
 96  use m_afterscfloop,     only : afterscfloop
 97  use m_extraprho,        only : extraprho
 98  use m_spacepar,         only : setsym
 99  use m_newrho,           only : newrho
100  use m_newvtr,           only : newvtr
101  use m_vtorho,           only : vtorho
102  use m_setvtr,           only : setvtr
103  use m_mkrho,            only : mkrho
104  use m_rhotov,           only : rhotov
105  use m_forces,           only : fresid, forces
106  use m_dft_energy,       only : energy
107  use m_initylmg,         only : initylmg
108  use m_rhotoxc,          only : rhotoxc
109  use m_drivexc,          only : check_kxc
110  use m_odamix,           only : odamix
111  use m_common,           only : scprqt, prtene
112  use m_fourier_interpol, only : transgrid
113  use m_fock_getghc,      only : fock2ACE
114  use m_forstr,           only : nres2vres
115  use m_positron,         only : setup_positron
116  use m_cgprj,            only : ctocprj
117  use m_psolver,          only : psolver_rhohxc
118  use m_paw2wvl,          only : paw2wvl_ij, wvl_cprjreorder
119 
120  implicit none
121 
122  private

ABINIT/scfcv_core [ Functions ]

[ Top ] [ Functions ]

NAME

 scfcv_core

FUNCTION

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

INPUTS

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

OUTPUT

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

SIDE EFFECTS

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

NOTES

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

PARENTS

      m_scfcv

CHILDREN

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

SOURCE

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

ABINIT/wf_mixing [ Functions ]

[ Top ] [ Functions ]

NAME

 wf_mixing

FUNCTION

 Mixing of wavefunctions in the outer loop of a double loop SCF approach.
 Different algorithms are implemented, depending on the value of wfmixalg.

INPUTS

  atindx1(dtset%natom)=index table for atoms, inverse of atindx
  dtset <type(dataset_type)>=all input variables in this dataset
  istep=number of call the routine (usually the outer loop in the SCF double loop)
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of cprj array
  mpi_enreg=information about MPI parallelization
  nattyp(dtset%ntypat)=number of atoms of each type in cell.
  npwarr(nkpt)=number of planewaves in basis at this k point
  option= 1 for wf prediction for molecular dynamics; 0 otherwise.
  pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data

SIDE EFFECTS

  cg(2,mcg)= plane wave wavefunction coefficient
                          Value from previous SCF cycle is input and stored in some form
                          Extrapolated value is output
  cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
                          Value from previous SCF cycle is input and stored in some form
                          Extrapolated value is output
  scf_history_wf <type(scf_history_type)>=arrays obtained from previous SCF cycles

PARENTS

      scfcv_core

CHILDREN

      cgcprj_cholesky,dotprod_set_cgcprj,dotprodm_sumdiag_cgcprj
      lincom_cgcprj,pawcprj_alloc,pawcprj_axpby,pawcprj_free,pawcprj_get
      pawcprj_getdim,pawcprj_lincom,pawcprj_put,timab,xmpi_sum,zgesv

SOURCE

2748 subroutine wf_mixing(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
2749 & nattyp,npwarr,pawtab,scf_history_wf,option)
2750 
2751  !use m_scf_history
2752  use m_cgcprj,  only : dotprod_set_cgcprj, dotprodm_sumdiag_cgcprj, lincom_cgcprj, cgcprj_cholesky
2753 
2754 !This section has been created automatically by the script Abilint (TD).
2755 !Do not modify the following lines by hand.
2756 #undef ABI_FUNC
2757 #define ABI_FUNC 'wf_mixing'
2758 !End of the abilint section
2759 
2760  implicit none
2761 
2762 !Arguments ------------------------------------
2763 !scalars
2764  integer,intent(in) :: istep,mcg,mcprj,option
2765  type(MPI_type),intent(in) :: mpi_enreg
2766  type(dataset_type),intent(in) :: dtset
2767  type(scf_history_type),intent(inout) :: scf_history_wf
2768 !arrays
2769  integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat)
2770  integer,intent(in) :: npwarr(dtset%nkpt)
2771  real(dp), intent(inout) :: cg(2,mcg)
2772  type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj)
2773  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
2774 
2775 !Local variables-------------------------------
2776 !scalars
2777  integer :: hermitian
2778  integer :: ibdmix,ibdsp,ibg,ibg_hist,icg,icg_hist
2779  integer :: ierr,ikpt,indh,ind_biorthog,ind_biorthog_eff,ind_newwf,ind_residual,inplace
2780  integer :: iorder,iset2,isppol,istep_cycle,istep_new,istwf_k,kk,me_distrb,my_nspinor
2781  integer :: nband_k,nbdmix,npw_k,nset1,nset2,ntypat
2782  integer :: shift_set1,shift_set2,spaceComm_band,spare_mem,usepaw,wfmixalg
2783  real(dp) :: alpha,beta
2784  complex(dpc) :: sum_coeffs
2785 !arrays
2786  integer,allocatable :: ipiv(:),dimcprj(:)
2787  real(dp) :: tsec(2)
2788  real(dp),allocatable :: al(:,:),mmn(:,:,:)
2789  real(dp),allocatable :: dotprod_res(:,:,:),dotprod_res_k(:,:,:),res_mn(:,:,:),smn(:,:,:)
2790  complex(dpc),allocatable :: coeffs(:)
2791  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:)
2792 
2793 ! *************************************************************************
2794 
2795 !DEBUG
2796 !write(std_out,*)
2797 !write(std_out,*)' wf_mixing : enter, istep= ',istep
2798 !call flush(std_out)
2799 !write(std_out,*)' istep,scf_history_wf%alpha=',istep,scf_history_wf%alpha
2800 !write(std_out,*)' cg(1,1)=',cg(1,1)
2801 !write(std_out,*)' scf_history_wf%cg(1,1,1:5)=',scf_history_wf%cg(1,1,1:5)
2802 !ABI_ALLOCATE(cg_ref,(2,mcg))
2803 !cg_ref(:,:)=cg(:,:)
2804 !ABI_DATATYPE_ALLOCATE(cprj_ref,(dtset%natom,mcprj))
2805 !cprj_ref(:,:)=cprj(:,:)
2806 !      write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',&
2807 !&       scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)
2808 !       call flush(std_out)
2809 !ENDDEBUG
2810 
2811  if (istep==0) return
2812 
2813  ntypat=dtset%ntypat
2814  usepaw=dtset%usepaw
2815  wfmixalg=scf_history_wf%wfmixalg
2816  nbdmix=dtset%nbandhf
2817  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
2818  me_distrb=mpi_enreg%me_kpt
2819  spaceComm_band=xmpi_comm_self
2820 
2821  spare_mem=0
2822  if(scf_history_wf%history_size==wfmixalg-1)spare_mem=1
2823 
2824 !scf_history_wf%alpha contains dtset%wfmix
2825  alpha=scf_history_wf%alpha
2826  beta=one-scf_history_wf%alpha
2827  if(option==1) then
2828    wfmixalg=3
2829    beta=scf_history_wf%beta
2830    nbdmix=dtset%mband
2831  end if
2832  icg=0
2833  icg_hist=0
2834  ibg=0
2835  ibg_hist=0
2836 
2837 !Useful array
2838  ABI_ALLOCATE(dimcprj,(dtset%natom))
2839  if (usepaw==1) then
2840    iorder=0 ! There is no change of ordering in the mixing of wavefunctions
2841    call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
2842  end if
2843 
2844  if(istep==1)then
2845    do indh=1,scf_history_wf%history_size
2846      call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj)
2847    end do
2848  end if
2849 
2850  ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,my_nspinor*nbdmix))
2851  ABI_DATATYPE_ALLOCATE(cprj_kh,(dtset%natom,my_nspinor*nbdmix))
2852  if(usepaw==1) then
2853    call pawcprj_alloc(cprj_k,0,dimcprj)
2854    call pawcprj_alloc(cprj_kh,0,dimcprj)
2855  end if
2856  ABI_ALLOCATE(smn,(2,nbdmix,nbdmix))
2857  ABI_ALLOCATE(mmn,(2,nbdmix,nbdmix))
2858 
2859  if(wfmixalg>2)then
2860    nset1=1
2861    nset2=min(istep-1,wfmixalg-1)
2862    if(option==0) then
2863      ABI_ALLOCATE(dotprod_res_k,(2,1,nset2))
2864      ABI_ALLOCATE(dotprod_res,(2,1,nset2))
2865      ABI_ALLOCATE(res_mn,(2,wfmixalg-1,wfmixalg-1))
2866      dotprod_res=zero
2867      if(istep==1)then
2868        scf_history_wf%dotprod_sumdiag_cgcprj_ij=zero
2869      end if
2870    end if
2871  end if
2872 
2873 !Explanation for the index for the wavefunction stored in scf_history_wf
2874 !The reference is the cg+cprj output after the wf optimization at istep 1.
2875 !It comes as input to the present routine as cgcprj input at step 2, and is usually found at indh=1.
2876 
2877 !In the simple mixing case (wfmixalg==2), the reference is never stored, because it is used "on-the-fly" to biothogonalize the
2878 !previous input (that was stored in indh=1), then generate the next input, which is stored again in indh=1
2879 
2880 !When the storage is not spared:
2881 !- the values of indh from 2 to wfmixalg store the (computed here) biorthogonalized input cgcprj, then the residual
2882 !- the values of indh from wfmixalg+1 to 2*wfmixalg-1 store the biorthogonalized output cgcprj (coming as argument)
2883 
2884 !First step
2885  if (istep==1 .or. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) ) then
2886 
2887    indh=2   ! This input wavefunction is NOT the reference
2888    if(wfmixalg==2)indh=1 ! But this does not matter in the simple mixing case that has history_size=1
2889 
2890 !  Simply store the wavefunctions and cprj. However, nband_k might be different from nbandhf...
2891 !  LOOP OVER SPINS
2892    do isppol=1,dtset%nsppol
2893 
2894 !    BIG FAT k POINT LOOP
2895      do ikpt=1,dtset%nkpt
2896 
2897 !      Select k point to be treated by this proc
2898        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2899        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
2900 
2901        npw_k=npwarr(ikpt)
2902 
2903        scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
2904        if(usepaw==1) then
2905 !        scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix)
2906          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,dtset%mband,&
2907 &         dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,&
2908 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2909          call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
2910 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
2911 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
2912        end if
2913 
2914 !      Update the counters
2915        ibg=ibg+my_nspinor*nband_k
2916        ibg_hist=ibg_hist+my_nspinor*nbdmix
2917        icg=icg+my_nspinor*nband_k*npw_k
2918        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
2919 
2920      end do
2921    end do
2922 
2923  else
2924 !  From istep==2
2925 
2926 !  First part of the computation : biorthogonalization, and computation of the residual (possibly, prediction of the next input in the case of simple mixing)
2927 !  Index for the wavefunctions stored in scf_history_wf whose scalar products with the argument cgcprj will have to be computed.
2928    indh=1   ! This input wavefunction is the reference
2929    if(wfmixalg/=2 .and. istep==2)indh=2 ! except for istep=2 in the rmm-diis
2930 
2931    if(wfmixalg>2)then
2932 !    istep inside the cycle defined by wfmixalg, and next index. Then, indices of the wavefunction sets.
2933      istep_cycle=mod((istep-2),wfmixalg-1)
2934      istep_new=mod((istep-1),wfmixalg-1)
2935      ind_biorthog=1+wfmixalg+istep_cycle
2936      ind_residual=2+istep_cycle
2937      ind_newwf=2+istep_new
2938      shift_set1=ind_residual-1
2939      shift_set2=1
2940    end if
2941 
2942 !  LOOP OVER SPINS
2943    do isppol=1,dtset%nsppol
2944 
2945 !    BIG FAT k POINT LOOP
2946      do ikpt=1,dtset%nkpt
2947 
2948 !      Select k point to be treated by this proc
2949        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2950        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
2951 
2952        istwf_k=dtset%istwfk(ikpt)
2953        npw_k=npwarr(ikpt)
2954 
2955 !      Biorthogonalization
2956 
2957        if(usepaw==1) then
2958          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,dtset%mband,&
2959 &         dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,&
2960 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2961          call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
2962 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,my_nspinor,dtset%nsppol,0,&
2963 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2964        end if  !end usepaw=1
2965 
2966        hermitian=0
2967        if(wfmixalg==2 .or. istep==2)then
2968          call dotprod_set_cgcprj(atindx1,cg,scf_history_wf%cg(:,:,indh),cprj_k,cprj_kh,dimcprj,hermitian,&
2969 &         0,0,icg,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,&
2970 &         mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw)
2971        else
2972          call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,&
2973 &         0,0,icg_hist,icg,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,&
2974 &         mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw)
2975        end if
2976 
2977 !      Invert S matrix, that is NOT hermitian.
2978 !      Calculate M=S^-1
2979        mmn=zero
2980        do kk=1,nbdmix
2981          mmn(1,kk,kk)=one
2982        end do
2983 
2984        ABI_ALLOCATE(ipiv,(nbdmix))
2985 !      The smn is destroyed by the following inverse call
2986        call zgesv(nbdmix,nbdmix,smn,nbdmix,ipiv,mmn,nbdmix,ierr)
2987        ABI_DEALLOCATE(ipiv)
2988 !DEBUG
2989        if(ierr/=0)then
2990          MSG_ERROR(' The call to cgesv general inversion routine failed')
2991        end if
2992 !ENDDEBUG
2993 
2994 !      The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place
2995        if(wfmixalg==2 .or. istep==2)then
2996          inplace=1
2997          call lincom_cgcprj(mmn,scf_history_wf%cg(:,:,indh),cprj_kh,dimcprj,&
2998 &         icg_hist,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw)
2999        else
3000          inplace=0
3001          call lincom_cgcprj(mmn,cg,cprj_k,dimcprj,&
3002 &         icg,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,&
3003 &         cgout=scf_history_wf%cg(:,:,ind_biorthog),cprjout=scf_history_wf%cprj(:,:,ind_biorthog),icgout=icg_hist)
3004        end if
3005 
3006 !      The biorthogonalised set of wavefunctions is now stored at the proper place
3007 
3008 !      Finalize this first part of the computation, depending on the algorithm and the step.
3009 
3010        if(wfmixalg==2)then
3011 
3012 !        Wavefunction extrapolation, simple mixing case
3013 !        alpha contains dtset%wfmix, beta contains one-alpha
3014          cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=&
3015 &         alpha*cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)&
3016 &         +beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)
3017          if(usepaw==1) then
3018            do ibdmix=1,nbdmix
3019              call pawcprj_axpby(beta,alpha,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix))
3020            end do ! end loop on ibdmix
3021          end if
3022 
3023 !        Back to usual orthonormalization
3024          call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,&
3025 &         mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
3026 
3027 !        Store the newly extrapolated wavefunctions, orthonormalized, in scf_history_wf
3028          scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
3029          if(usepaw==1) then
3030            do ibdmix=1,nbdmix
3031              call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3032 &             nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3033 &             mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3034            end do ! end loop on ibdmix
3035          end if
3036 
3037        else  !  wfmixalg/=2
3038 !        RMM-DIIS
3039 
3040          if (istep==2)then
3041 !          Store the argument wf as the reference for all future steps, in scf_history_wf with index 1.
3042            scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
3043            if(usepaw==1) then
3044              do ibdmix=1,nbdmix
3045                call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,1),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3046 &               nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3047 &               mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3048              end do ! end loop on ibdmix
3049            end if
3050          end if
3051 
3052          ind_biorthog_eff=ind_biorthog
3053          if(istep==2)ind_biorthog_eff=1 ! The argument wf has not been stored in ind_biorthog
3054          if(option==0) then
3055 !          Compute the residual of the wavefunctions for this istep,
3056 !          that replaces the previously stored set of (biorthogonalized) input wavefunctions
3057            scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)=&
3058 &           scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)&
3059 &           -scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)
3060            if(usepaw==1) then
3061              do ibdmix=1,nbdmix
3062                call pawcprj_axpby(one,-one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff),&
3063 &               scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual))
3064              end do ! end loop on ibdmix
3065            end if
3066 
3067 !          Compute the new scalar products to fill the res_mn matrix
3068            call dotprodm_sumdiag_cgcprj(atindx1,scf_history_wf%cg,scf_history_wf%cprj,dimcprj,&
3069 &           ibg_hist,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcprj,dtset%mkmem,&
3070 &           mpi_enreg,scf_history_wf%history_size,dtset%natom,nattyp,nbdmix,npw_k,nset1,nset2,my_nspinor,dtset%nsppol,ntypat,&
3071 &           shift_set1,shift_set2,pawtab,dotprod_res_k,usepaw)
3072 
3073            dotprod_res=dotprod_res+dotprod_res_k
3074          end if
3075 !        scf_history_wf for index ind_biorthog will contain the extrapolated wavefunctions (and no more the output of the SCF loop).
3076          scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog)=&
3077 &         scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)+&
3078 &         (alpha-one)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)
3079          if(usepaw==1) then
3080            do ibdmix=1,nbdmix
3081              if(ind_biorthog/=ind_biorthog_eff)then
3082                scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog)=scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff)
3083              end if
3084              call pawcprj_axpby((alpha-one),one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual),&
3085 &             scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog))
3086            end do ! end loop on ibdmix
3087          end if
3088 
3089        end if
3090 
3091        ibg=ibg+my_nspinor*nband_k
3092        ibg_hist=ibg_hist+my_nspinor*nbdmix
3093        icg=icg+my_nspinor*nband_k*npw_k
3094        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
3095 
3096 !      End big k point loop
3097      end do
3098 !    End loop over spins
3099    end do
3100 
3101  end if ! istep>=2
3102 
3103  if(wfmixalg>2 .and. istep>1)then
3104 
3105 !DEBUG
3106 !      write(std_out,*)' '
3107 !      write(std_out,*)' Entering the residual minimisation part '
3108 !      write(std_out,*)' '
3109 !      call flush(std_out)
3110 !ENDDEBUG
3111 
3112    call timab(48,1,tsec)
3113    call xmpi_sum(dotprod_res,mpi_enreg%comm_kpt,ierr)
3114    call timab(48,2,tsec)
3115 
3116    scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set1,1+shift_set2:nset2+shift_set2)=dotprod_res(:,1,1:nset2)
3117    scf_history_wf%dotprod_sumdiag_cgcprj_ij(1,1+shift_set2:nset2+shift_set2,1+shift_set1)=dotprod_res(1,1,1:nset2)
3118    scf_history_wf%dotprod_sumdiag_cgcprj_ij(2,1+shift_set2:nset2+shift_set2,1+shift_set1)=-dotprod_res(2,1,1:nset2)
3119 
3120  end if ! wfmixalg>2 and istep>1
3121 
3122  if(wfmixalg>2 .and. istep>2)then
3123 
3124 !  Extract the relevant matrix R_mn
3125    res_mn(:,1:nset2,1:nset2)=&
3126 &   scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set2:nset2+shift_set2,1+shift_set2:nset2+shift_set2)
3127 
3128 !DEBUG
3129 !      write(std_out,*)' The matrix res_mn(:,1:nset2,1:nset2) is :'
3130 !      write(std_out,*)res_mn(:,1:nset2,1:nset2)
3131 !      call flush(std_out)
3132 !ENDDEBUG
3133 
3134 !  Solve R_mn \alpha_n = 1_m
3135    ABI_ALLOCATE(ipiv,(nset2))
3136    ABI_ALLOCATE(coeffs,(nset2))
3137    coeffs(:)=cone
3138 !  The res_mn is destroyed by the following inverse call
3139    call zgesv(nset2,1,res_mn,wfmixalg-1,ipiv,coeffs,nset2,ierr)
3140    ABI_DEALLOCATE(ipiv)
3141 !  The coefficients must sum to one
3142    sum_coeffs=sum(coeffs)
3143    coeffs=coeffs/sum_coeffs
3144 
3145 !DEBUG
3146 !      write(std_out,*)' The coefficients that minimize the residual have been found'
3147 !      write(std_out,*)' coeffs =',coeffs
3148 !      call flush(std_out)
3149 !ENDDEBUG
3150  end if ! wfmixalg>2 and istep>2
3151 
3152  if(wfmixalg>2 .and. istep>1)then
3153 
3154 !  Find the new "input" wavefunction, bi-orthogonalized, and store it replacing the adequate "old" input wavefunction.
3155 
3156    icg=0
3157    icg_hist=0
3158    ibg=0
3159    ibg_hist=0
3160    ABI_ALLOCATE(al,(2,nset2))
3161    if(istep>2)then
3162      do iset2=1,nset2
3163        al(1,iset2)=real(coeffs(iset2)) ; al(2,iset2)=aimag(coeffs(iset2))
3164      end do
3165    else
3166      al(1,1)=one ; al(2,1)=zero
3167    end if
3168 
3169 !DEBUG
3170 !      write(std_out,*)' Overload the coefficients, in order to simulate a simple mixing with wfmix '
3171 !      write(std_out,*)' Set al(1,ind_biorthog-3)=one, for ind_biorthog=',ind_biorthog
3172 !      write(std_out,*)' This will feed scf_history for set ind_biorthog-3+wfmixalg=',ind_biorthog-3+wfmixalg
3173 !      al(:,:)=zero
3174 !      al(1,ind_biorthog-3)=one
3175 !      call flush(std_out)
3176 !ENDDEBUG
3177 
3178 !  LOOP OVER SPINS
3179    do isppol=1,dtset%nsppol
3180 
3181 !    BIG FAT k POINT LOOP
3182      do ikpt=1,dtset%nkpt
3183 
3184 !      Select k point to be treated by this proc
3185        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
3186        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
3187 
3188        istwf_k=dtset%istwfk(ikpt)
3189        npw_k=npwarr(ikpt)
3190 
3191        if(istep>2)then
3192 !        Make the appropriate linear combination (from the extrapolated wfs)
3193          cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero
3194          do iset2=1,nset2
3195            cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)&
3196 &           +al(1,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)&
3197 &           -al(2,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)
3198            cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)&
3199 &           +al(1,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)&
3200 &           +al(2,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)
3201          end do
3202        else ! One needs a simple copy from the extrapolated wavefunctions
3203          cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1+wfmixalg)
3204        end if
3205 !      Note the storage in cprj_k. By the way, a simple copy might also be used in case istep=2.
3206        if(usepaw==1) then
3207          do ibdsp=1,my_nspinor*nbdmix
3208            call pawcprj_lincom(al,scf_history_wf%cprj(:,ibdsp,1+wfmixalg:nset2+wfmixalg),cprj_k(:,ibdsp:ibdsp),nset2)
3209          end do
3210        end if
3211 
3212 !      Store the newly extrapolated wavefunctions for this k point, still bi-orthonormalized, in scf_history_wf
3213        scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_newwf)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
3214        if(usepaw==1) then
3215          call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,ind_newwf),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3216 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3217 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3218        end if
3219 
3220 !      Back to usual orthonormalization for the cg and cprj_k
3221        call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,&
3222 &       mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
3223 
3224 !      Need to transfer cprj_k to cprj
3225        if(usepaw==1) then
3226          call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,&
3227 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3228 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3229        end if
3230 
3231        ibg=ibg+my_nspinor*nband_k
3232        ibg_hist=ibg_hist+my_nspinor*nbdmix
3233        icg=icg+my_nspinor*nband_k*npw_k
3234        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
3235 
3236      end do ! End big k point loop
3237    end do ! End loop over spins
3238 
3239    if(istep>2)then
3240      ABI_DEALLOCATE(coeffs)
3241    end if
3242    ABI_DEALLOCATE(al)
3243 
3244  end if ! wfmixalg>2 and istep>1
3245 
3246 !DEBUG
3247 ! write(std_out,*)' wf_mixing : exit '
3248 !      write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',&
3249 !&       scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)
3250 ! write(std_out,*)' cg(1:2,1:2)=',cg(1:2,1:2)
3251 ! write(std_out,*)' scf_history_wf%cg(1:2,1:2,1)=',scf_history_wf%cg(1:2,1:2,1)
3252 ! ABI_DEALLOCATE(cg_ref)
3253 ! ABI_DATATYPE_DEALLOCATE(cprj_ref)
3254 !ENDDEBUG
3255 
3256  if(usepaw==1) then
3257    call pawcprj_free(cprj_k)
3258    call pawcprj_free(cprj_kh)
3259  end if
3260  ABI_DATATYPE_DEALLOCATE(cprj_k)
3261  ABI_DATATYPE_DEALLOCATE(cprj_kh)
3262  ABI_DEALLOCATE(dimcprj)
3263  ABI_DEALLOCATE(mmn)
3264  ABI_DEALLOCATE(smn)
3265  if(wfmixalg>2.and.option==0)then
3266    ABI_DEALLOCATE(dotprod_res_k)
3267    ABI_DEALLOCATE(dotprod_res)
3268    ABI_DEALLOCATE(res_mn)
3269  end if
3270 
3271 end subroutine wf_mixing