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)
  grcondft(3,natom)=grads of constrained DFT 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
  vxctau(nfftf,dtset%nspden,4*usekden)]=derivative of XC energy density wrt
      kinetic energy density (metaGGA cases)
  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 gred (hartree/bohr)
   gred(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_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE part of total energy.
   | e_nucdip(IN)=energy due to array of nuclear magnetic dipoles
   | 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.

SOURCE

2507 subroutine etotfor(atindx1,deltae,diffor,dtefield,dtset,&
2508 &  elast,electronpositron,energies,&
2509 &  etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,&
2510 &  grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfft,mpi_enreg,my_natom,nattyp,&
2511 &  nfft,ngfft,ngrvdw,nhat,nkxc,ntypat,nvresid,n1xccc,n3xccc,optene,optforces,optres,&
2512 &  pawang,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,red_ptot,psps,rhog,rhor,rmet,rprimd,&
2513 &  symrec,synlgr,ucvol,usepaw,vhartr,vpsp,vxc,vxctau,wvl,wvl_den,xccc3d,xred)
2514 
2515 !Arguments ------------------------------------
2516 !scalars
2517  integer,intent(in) :: my_natom,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nkxc,ntypat,optene,optforces
2518  integer,intent(in) :: optres,usepaw
2519  real(dp),intent(in) :: gsqcut
2520  real(dp),intent(inout) :: elast,ucvol
2521  real(dp),intent(out) :: deltae,diffor,etotal,maxfor
2522  type(MPI_type),intent(in) :: mpi_enreg
2523  type(efield_type),intent(in) :: dtefield
2524  type(dataset_type),intent(in) :: dtset
2525  type(electronpositron_type),pointer :: electronpositron
2526  type(energies_type),intent(inout) :: energies
2527  type(extfpmd_type),pointer,intent(inout) :: extfpmd
2528  type(pawang_type),intent(in) :: pawang
2529  type(pseudopotential_type),intent(in) :: psps
2530  type(wvl_internal_type), intent(in) :: wvl
2531  type(wvl_denspot_type), intent(inout) :: wvl_den
2532  type(fock_type),pointer, intent(inout) :: fock
2533 !arrays
2534  integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
2535  integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym)
2536  real(dp),intent(in) :: gmet(3,3),grchempottn(3,dtset%natom),grcondft(3,dtset%natom)
2537  real(dp),intent(in) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc)
2538  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),red_ptot(3)
2539  real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rmet(3,3)
2540  real(dp),intent(in) :: vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden)
2541  real(dp),intent(in) :: vxctau(nfft,dtset%nspden,4*dtset%usekden)
2542  real(dp),intent(in) :: xccc3d(n3xccc)
2543  real(dp),intent(inout) :: forold(3,dtset%natom),grnl(3*dtset%natom)
2544  real(dp),intent(inout) :: nhat(nfft,dtset%nspden*psps%usepaw)
2545  real(dp),intent(inout),target :: nvresid(nfft,dtset%nspden)
2546  real(dp),intent(inout) :: xred(3,dtset%natom)
2547  real(dp),intent(out) :: favg(3),gred(3,dtset%natom)
2548  real(dp),intent(inout) :: fcart(3,dtset%natom)
2549  real(dp),intent(inout) :: rprimd(3,3)
2550  real(dp),intent(out) :: gresid(3,dtset%natom),grhf(3,dtset%natom)
2551  real(dp),intent(inout) :: grxc(3,dtset%natom)
2552  real(dp),intent(out) :: synlgr(3,dtset%natom)
2553  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
2554  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
2555  type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw)
2556  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
2557 
2558 !Local variables-------------------------------
2559 !scalars
2560  integer :: comm_grid,dimnhat,ifft,ipositron,ispden,optgr,optgr2,option,optnc,optstr,optstr2,iir,jjr,kkr
2561  logical :: apply_residual
2562  real(dp) :: eenth,ucvol_
2563 !arrays
2564  real(dp),parameter :: k0(3)=(/zero,zero,zero/)
2565  real(dp) :: tsec(2),A(3,3),A1(3,3),A_new(3,3),efield_new(3)
2566  real(dp) :: dummy(0),nhat_dum(0,0)
2567  real(dp),allocatable :: vlocal(:,:)
2568  real(dp), ABI_CONTIGUOUS pointer :: resid(:,:)
2569 
2570 ! *********************************************************************
2571 
2572  call timab(80,1,tsec)
2573 
2574  ipositron=electronpositron_calctype(electronpositron)
2575 
2576  if (optene>-1) then
2577 
2578 !  Add the contribution of extfpmd to the entropy
2579    if(associated(extfpmd)) then
2580      energies%entropy=energies%entropy+extfpmd%entropy
2581    end if
2582 
2583 !  When the finite-temperature VG broadening scheme is used,
2584 !  the total entropy contribution "tsmear*entropy" has a meaning,
2585 !  and gather the two last terms of Eq.8 of VG paper
2586 !  Warning : might have to be changed for fixed moment calculations
2587    if(dtset%occopt>=3 .and. dtset%occopt<=8) then
2588      if (abs(dtset%tphysel) < tol10) then
2589        energies%e_entropy = - dtset%tsmear * energies%entropy
2590      else
2591        energies%e_entropy = - dtset%tphysel * energies%entropy
2592      end if
2593    else
2594      energies%e_entropy = zero
2595    end if
2596 
2597 !  Turn it into an electric enthalpy, refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]],
2598 !    the missing volume is added here
2599    energies%e_elecfield=zero
2600    if ((dtset%berryopt==4.or.dtset%berryopt==14).and.ipositron/=1) then
2601      energies%e_elecfield=-dot_product(dtset%red_efieldbar,red_ptot)  !!ebar_i p_i
2602      eenth=zero
2603      do iir=1,3
2604        do jjr=1,3
2605          eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)  !!g^{-1})_ij ebar_i ebar_j
2606        end do
2607      end do
2608      energies%e_elecfield=energies%e_elecfield-eenth*ucvol/(8._dp*pi)
2609    end if
2610 
2611 !  Turn it into an internal energy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]],
2612 !    but a little different: U=E_ks + (vol/8*pi) *  g^{-1})_ij ebar_i ebar_j
2613    if ((dtset%berryopt==6.or.dtset%berryopt==16).and.ipositron/=1) then
2614      energies%e_elecfield=zero
2615      eenth=zero
2616      do iir=1,3
2617        do jjr=1,3
2618          eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)  !! g^{-1})_ij ebar_i ebar_j
2619        end do
2620      end do
2621      energies%e_elecfield=energies%e_elecfield+eenth*ucvol/(8._dp*pi)
2622    end if
2623 
2624 !  Calculate internal energy and electric enthalpy for mixed BC case.
2625    if (dtset%berryopt==17.and.ipositron/=1) then
2626      energies%e_elecfield=zero
2627      A(:,:)=(four_pi/ucvol)*rmet(:,:)
2628      A1(:,:)=A(:,:) ; A_new(:,:)=A(:,:)
2629      efield_new(:)=dtset%red_efield(:)
2630      eenth=zero
2631      do kkr=1,3
2632        if (dtset%jfielddir(kkr)==1) then    ! fixed ebar direction
2633 !        step 1 add -ebar*p
2634          eenth=eenth-dtset%red_efieldbar(kkr)*red_ptot(kkr)
2635 !        step 2  chang to e_new (change e to ebar)
2636          efield_new(kkr)=dtset%red_efieldbar(kkr)
2637 !        step 3  chang matrix A to A1
2638          do iir=1,3
2639            do jjr=1,3
2640              if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr)
2641              if ((iir==kkr .and. jjr/=kkr) .or.  (iir/=kkr .and.  jjr==kkr)) &
2642 &             A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr)
2643              if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr)
2644            end do
2645          end do
2646          A(:,:)=A1(:,:) ; A_new(:,:)=A1(:,:)
2647        end if
2648      end do  ! end fo kkr
2649      do iir=1,3
2650        do jjr=1,3
2651          eenth= eenth+half*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr)
2652        end do
2653      end do
2654      energies%e_elecfield=energies%e_elecfield+eenth
2655    end if   ! berryopt==17
2656 
2657 !  Turn it into a magnetic enthalpy, by adding orbital electronic contribution
2658    energies%e_magfield = zero
2659 !  if (dtset%berryopt == 5 .and. ipositron/=1) then
2660 !  emag = dot_product(mag_cart,dtset%bfield)
2661 !  energies%e_magfield = emag
2662 !  end if
2663 
2664 !  Compute total (free)- energy by direct scheme
2665    if (optene==0) then
2666      etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc &
2667 &           + energies%e_localpsp + energies%e_corepsp &
2668 &           + energies%e_entropy + energies%e_elecfield &
2669 &           + energies%e_magfield + energies%e_nucdip &
2670 &           + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_hybcomp_v &
2671 &           + energies%e_constrained_dft + energies%e_ewald &
2672 &           + energies%e_chempot + energies%e_vdw_dftd
2673 
2674 !    +two*energies%e_fock-energies%e_fock0 ! The Fock energy is already included in the non-local one
2675 !    +energies%e_nlpsp_vfock - energies%e_fock0
2676 
2677 !    See similar section in m_energies.F90
2678 !    XG 20181025 This gives a variational energy in case of NCPP with all bands occupied - not yet for metals.
2679      if (usepaw==0) etotal = etotal + energies%e_nlpsp_vfock - energies%e_fock0
2680 !    XG 20181025 I was expecting the following to give also a variational energy in case of PAW, but this is not true.
2681 !    if (usepaw==1) etotal = etotal + energies%e_paw + energies%e_nlpsp_vfock - energies%e_fock0
2682 !    XG 20181025 So, the following is giving a non-variational expression ...
2683      if (usepaw==1) etotal = etotal + energies%e_paw + energies%e_fock
2684 
2685    end if
2686 
2687 !  Compute total (free) energy by double-counting scheme
2688    if (optene==1) then
2689      etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc &
2690 &     - energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc- energies%e_fock0 &
2691 &     + energies%e_entropy + energies%e_elecfield + energies%e_magfield &
2692 &     + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_constrained_dft
2693      etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
2694      if (usepaw/=0) etotal = etotal + energies%e_pawdc
2695    end if
2696 
2697 !  Additional stuff for electron-positron
2698    if (dtset%positron/=0) then
2699      if (ipositron==0) then
2700        energies%e_electronpositron  =zero
2701        energies%edc_electronpositron=zero
2702      else
2703        energies%e_electronpositron  =electronpositron%e_hartree+electronpositron%e_xc
2704        energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc
2705        if (usepaw==1) then
2706          energies%e_electronpositron  =energies%e_electronpositron  +electronpositron%e_paw
2707          energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc
2708        end if
2709      end if
2710      if (optene==0) electronpositron%e0=etotal
2711      if (optene==1) electronpositron%e0=etotal-energies%edc_electronpositron
2712      etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron
2713    end if
2714 
2715 !  Add the extfpmd energy contribution to the internal energy
2716    if(associated(extfpmd)) then
2717      energies%e_extfpmd=extfpmd%e_kinetic
2718      energies%edc_extfpmd=extfpmd%edc_kinetic
2719      if(optene==0) etotal=etotal+energies%e_extfpmd
2720      if(optene==1) etotal=etotal+energies%edc_extfpmd
2721    end if
2722 
2723 !  Compute energy residual
2724    deltae=etotal-elast
2725    elast=etotal
2726  end if !optene/=-1
2727 
2728  call timab(80,2,tsec)
2729 
2730 !------Compute forces-----------------------------------------------------
2731 
2732  if (optforces==1) then
2733 
2734 !  PAW: add gradients due to Dij derivatives to non-local term
2735    if (usepaw==1) then
2736      ABI_MALLOC(vlocal,(nfft,dtset%nspden))
2737      do ispden=1,min(dtset%nspden,2)
2738 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vlocal,vpsp,vxc)
2739        do ifft=1,nfft
2740          vlocal(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)
2741        end do
2742      end do
2743 
2744      if(dtset%nspden==4)then
2745        do ispden=3,4
2746 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vlocal,vxc)
2747          do ifft=1,nfft
2748            vlocal(ifft,ispden)=vxc(ifft,ispden)
2749          end do
2750        end do
2751      end if
2752      ucvol_=ucvol
2753 #if defined HAVE_BIGDFT
2754      if (dtset%usewvl==1) ucvol_=product(wvl_den%denspot%dpbox%hgrids)*real(product(wvl_den%denspot%dpbox%ndims),dp)
2755 #endif
2756      dimnhat=0;optgr=1;optgr2=0;optstr=0;optstr2=0
2757      comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl
2758      call pawgrnl(atindx1,dimnhat,dummy,1,dummy,grnl,gsqcut,mgfft,my_natom, &
2759 &     dtset%natom, nattyp,nfft,ngfft,nhat_dum,dummy,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,&
2760 &     pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred, &
2761 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=mpi_enreg%comm_fft,&
2762 &     comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,paral_kgb=mpi_enreg%paral_kgb)
2763      ABI_FREE(vlocal)
2764    end if
2765 
2766    apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. &
2767 &   abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5)
2768 
2769 !  If residual is a density residual (and forces from residual asked),
2770 !  has to convert it into a potential residual before calling forces routine
2771    if (apply_residual) then
2772      ABI_MALLOC(resid,(nfft,dtset%nspden))
2773      option=0; if (dtset%densfor_pred<0) option=1
2774      optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2
2775      call nres2vres(dtset,gsqcut,usepaw,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
2776 &     nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,&
2777 &     rhor,rprimd,usepaw,resid,xccc3d,xred,vxc)
2778    else
2779      resid => nvresid
2780    end if
2781    call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,gred,grchempottn,grcondft,gresid,grewtn,&
2782 &   grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfft,mpi_enreg,&
2783 &   n1xccc,n3xccc,nattyp,nfft,ngfft,ngrvdw,ntypat,pawrad,pawtab,&
2784 &   ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,dtset%usefock,resid,vxc,vxctau,wvl,wvl_den,xred,&
2785 &   electronpositron=electronpositron)
2786    if (apply_residual) then
2787      ABI_FREE(resid)
2788    end if
2789 
2790 !  Returned gred are full symmetrized gradients of Etotal
2791 !  wrt reduced coordinates xred, d(Etotal)/d(xred)
2792 !  Forces are contained in array fcart
2793 
2794  else   ! if optforces==0
2795    fcart=zero
2796    gred=zero
2797    favg=zero
2798    diffor=zero
2799    gresid=zero
2800    grhf=zero
2801    maxfor=zero
2802    synlgr=zero
2803  end if
2804 
2805  call timab(80,2,tsec)
2806 
2807 end subroutine etotfor

ABINIT/m_scfcv_core [ Modules ]

[ Top ] [ Modules ]

NAME

  m_scfcv_core

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 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 .

SOURCE

 16 #if defined HAVE_CONFIG_H
 17 #include "config.h"
 18 #endif
 19 
 20 #include "abi_common.h"
 21 
 22 ! nvtx related macro definition
 23 #include "nvtx_macros.h"
 24 
 25 module m_scfcv_core
 26 
 27  use defs_basis
 28  use defs_wvltypes
 29  use defs_rectypes
 30  use m_xmpi
 31  use m_abicore
 32  use m_wffile
 33  use m_rec
 34  use m_ab7_mixing
 35  use m_errors
 36  use m_efield
 37  use mod_prc_memory
 38  use m_nctk
 39  use m_hdr
 40  use m_xcdata
 41  use m_cgtools
 42  use m_dtfil
 43  use m_distribfft
 44  use m_extfpmd
 45 
 46  use m_nonlop,           only : nonlop_counter
 47  use defs_datatypes,     only : pseudopotential_type
 48  use defs_abitypes,      only : MPI_type
 49  use m_berryphase_new,   only : update_e_field_vars
 50  use m_dens,             only : constrained_dft_t, constrained_dft_ini, constrained_dft_free
 51  use m_time,             only : timab
 52  use m_fstrings,         only : int2char4, sjoin, itoa
 53  use m_symtk,            only : symmetrize_xred
 54  use m_geometry,         only : metric
 55  use m_fftcore,          only : getng, sphereboundary
 56  use m_time,             only : abi_wtime, sec2str
 57  use m_exit,             only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in
 58  use m_mpinfo,           only : destroy_mpi_enreg, iwrite_fftdatar, initmpi_seq, proc_distrb_cycle
 59  use m_ioarr,            only : fftdatar_write_from_hdr
 60  use m_results_gs ,      only : results_gs_type
 61  use m_scf_history,      only : scf_history_type, scf_history_init, scf_history_free
 62  use m_energies,         only : energies_type, energies_init, energies_copy
 63  use m_electronpositron, only : electronpositron_type, electronpositron_calctype
 64  use m_pawang,           only : pawang_type
 65  use m_pawrad,           only : pawrad_type
 66  use m_pawtab,           only : pawtab_type,pawtab_get_lsize
 67  use m_paw_an,           only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags
 68  use m_pawfgrtab,        only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
 69  use m_pawrhoij,         only : pawrhoij_type
 70  use m_pawcprj,          only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, &
 71 &                               pawcprj_free, pawcprj_axpby, pawcprj_put, pawcprj_getdim, pawcprj_reorder
 72  use m_pawdij,           only : pawdij, symdij
 73  use m_pawfgr,           only : pawfgr_type
 74  use m_paw_ij,           only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
 75  use m_paw_dmft,         only : paw_dmft_type
 76  use m_paw_nhat,         only : nhatgrid,wvl_nhatgrid,pawmknhat
 77  use m_paw_tools,        only : chkpawovlp
 78  use m_paw_denpot,       only : pawdenpot
 79  use m_paw_occupancies,  only : pawmkrhoij
 80  use m_paw_correlations, only : setnoccmmp,setrhoijpbe0
 81  use m_paw_mkrho,        only : pawmkrho
 82  use m_paw_uj,           only : pawuj_red, macro_uj_type
 83  use m_paw_dfpt,         only : pawgrnl
 84  use m_fock,             only : fock_type, fock_init, fock_destroy, fock_ACE_destroy, fock_common_destroy, &
 85                                 fock_BZ_destroy, fock_update_exc, fock_updatecwaveocc
 86  use m_gwls_hamiltonian, only : build_vxc
 87 #if defined HAVE_BIGDFT
 88  use BigDFT_API,         only : cprj_clean,cprj_paw_alloc
 89 #endif
 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 : make_vectornd,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 #if defined(HAVE_GPU_CUDA) && defined(HAVE_GPU_NVTX_V3)
121  use m_nvtx_data
122 #endif
123 
124  implicit none
125 
126  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
  itimes(2)=itime array, contain itime=itimes(1) and itimimage_gstate=itimes(2) from outer loops
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mcg=size of wave-functions array (cg) =mpw*my_nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =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.
  cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
  dtefield <type(efield_type)> = variables related to Berry phase
  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.

SOURCE

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

SOURCE

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