scf.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2014 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module scf_oct_m
  use batch_ops_oct_m
  use berry_oct_m
  use convergence_criterion_oct_m
  use criteria_factory_oct_m
  use debug_oct_m
  use density_oct_m
  use density_criterion_oct_m
  use eigensolver_oct_m
  use eigenval_criterion_oct_m
  use electron_space_oct_m
  use energy_calc_oct_m
  use energy_criterion_oct_m
  use forces_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use ks_potential_oct_m
  use lcao_oct_m
  use lda_u_oct_m
  use lda_u_io_oct_m
  use lda_u_mixer_oct_m
  use lalg_basic_oct_m
  use loct_oct_m
  use magnetic_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mix_oct_m
  use modelmb_exchange_syms_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use output_oct_m
  use output_low_oct_m
  use output_modelmb_oct_m
  use parser_oct_m
  use partial_charges_oct_m
  use profiling_oct_m
  use restart_oct_m
  use smear_oct_m
  use space_oct_m
  use species_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_io_oct_m
  use states_elec_restart_oct_m
  use stress_oct_m
  use symmetries_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use v_ks_oct_m
  use varinfo_oct_m
  use vdw_ts_oct_m
  use vtau_mixer_oct_m
  use walltimer_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
  use xc_f03_lib_m
  use xc_functional_oct_m
  use xc_interaction_oct_m
  use xc_oep_oct_m
  use xc_oep_photon_oct_m
 
  implicit none
 
  private
  public ::              &
    scf_t,               &
    scf_init,            &
    scf_mix_clear,       &
    scf_load,            &
    scf_start,           &
    scf_run,             &
    scf_iter,            &
    scf_iter_finish,     &
    scf_finish,          &
    scf_end,             &
    scf_state_info,      &
    scf_print_mem_use
 
  integer, public, parameter :: &
    VERB_NO      = 0,   &
    verb_compact = 1,   &
    verb_full    = 3
 
  type scf_t
    private
    integer, public :: max_iter
 
    real(real64), public :: lmm_r
 
    ! several convergence criteria
    logical :: conv_eigen_error
    logical :: check_conv
 
    integer :: mix_field
    logical :: lcao_restricted
    logical :: calc_force
    logical, public :: calc_stress
    logical :: calc_dipole
    logical :: calc_partial_charges
    type(mix_t) :: smix
    type(mixfield_t), pointer :: mixfield
    type(eigensolver_t) :: eigens
    integer :: mixdim1
    logical :: forced_finish
    type(lda_u_mixer_t) :: lda_u_mix
    type(vtau_mixer_t) :: vtau_mix
    type(berry_t) :: berry
    integer :: matvec
 
    type(restart_t), public :: restart_load, restart_dump
 
    type(criterion_list_t), public :: criterion_list
    real(real64) :: energy_in, energy_diff, abs_dens_diff, evsum_in, evsum_out, evsum_diff
 
    ! Variables needed to store information accross scf_start, scf_run, and scf_finish
    logical :: converged_current, converged_last
    integer :: verbosity_
    type(lcao_t) :: lcao
    real(real64), allocatable :: rhoout(:,:), rhoin(:,:)
    real(real64), allocatable :: vhxc_old(:,:)
    class(wfs_elec_t), allocatable :: psioutb(:, :)
    logical :: output_forces, calc_current, output_during_scf, finish
  end type scf_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine scf_init(scf, namespace, gr, ions, st, mc, hm, space)
    type(scf_t),              intent(inout) :: scf
    type(grid_t),             intent(in)    :: gr
    type(namespace_t),        intent(in)    :: namespace
    type(ions_t),             intent(in)    :: ions
    type(states_elec_t),      intent(in)    :: st
    type(multicomm_t),        intent(in)    :: mc
    type(hamiltonian_elec_t), intent(inout) :: hm
    class(space_t),           intent(in)    :: space
 
    real(real64) :: rmin
    integer :: mixdefault
    type(type_t) :: mix_type
    class(convergence_criterion_t), pointer    :: crit
    type(criterion_iterator_t) :: iter
    logical :: deactivate_oracle
 
    push_sub(scf_init)
 
    !%Variable MaximumIter
    !%Type integer
    !%Default 200
    !%Section SCF::Convergence
    !%Description
    !% Maximum number of SCF iterations. The code will stop even if convergence
    !% has not been achieved. -1 means unlimited.
    !% 0 means just do LCAO (or read from restart), compute the eigenvalues and energy,
    !% and stop, without updating the wavefunctions or density.
    !%
    !% If convergence criteria are set, the SCF loop will only stop once the criteria
    !% are fulfilled for two consecutive iterations.
    !%
    !% Note that this variable is also used in the section Calculation Modes::Unoccupied States,
    !% where it denotes the maximum number of calls of the eigensolver. In this context, the
    !% default value is 50.
    !%End
    call parse_variable(namespace, 'MaximumIter', 200, scf%max_iter)
 
    if (allocated(hm%vberry)) then
      call berry_init(scf%berry, namespace)
    end if
 
    !Create the list of convergence criteria
    call criteria_factory_init(scf%criterion_list, namespace, scf%check_conv)
    !Setting the pointers
    call iter%start(scf%criterion_list)
    do while (iter%has_next())
      crit => iter%get_next()
      select type (crit)
      type is (energy_criterion_t)
        call crit%set_pointers(scf%energy_diff, scf%energy_in)
      type is (density_criterion_t)
        call crit%set_pointers(scf%abs_dens_diff, st%qtot)
      type is (eigenval_criterion_t)
        call crit%set_pointers(scf%evsum_diff, scf%evsum_out)
      class default
        assert(.false.)
      end select
    end do
 
 
    if(.not. scf%check_conv .and. scf%max_iter < 0) then
      call messages_write("All convergence criteria are disabled. Octopus is cowardly refusing")
      call messages_new_line()
      call messages_write("to enter an infinite loop.")
      call messages_new_line()
      call messages_new_line()
      call messages_write("Please set one of the following variables to a positive value:")
      call messages_new_line()
      call messages_new_line()
      call messages_write(" | MaximumIter | ConvEnergy | ConvAbsDens | ConvRelDens |")
      call messages_new_line()
      call messages_write(" |  ConvAbsEv  | ConvRelEv  |")
      call messages_new_line()
      call messages_fatal(namespace=namespace)
    end if
 
    !%Variable ConvEigenError
    !%Type logical
    !%Default false
    !%Section SCF::Convergence
    !%Description
    !% If true, the calculation will not be considered converged unless all states have
    !% individual errors less than <tt>EigensolverTolerance</tt>.
    !% If <tt>ExtraStatesToConverge</tt> is set, the calculation will stop
    !% when all occupied states plus <tt>ExtraStatesToConverge</tt> extra states are converged.
    !%
    !% If this criterion is used, the SCF loop will only stop once it is
    !% fulfilled for two consecutive iterations.
    !%End
    call parse_variable(namespace, 'ConvEigenError', .false., scf%conv_eigen_error)
 
    if(scf%max_iter < 0) scf%max_iter = huge(scf%max_iter)
 
    call messages_obsolete_variable(namespace, 'What2Mix', 'MixField')
 
    ! now the eigensolver stuff
    deactivate_oracle = hm%theory_level == independent_particles
    call eigensolver_init(scf%eigens, namespace, gr, st, mc, space, deactivate_oracle)
 
    if(scf%eigens%es_type /= rs_evo) then
      !%Variable MixField
      !%Type integer
      !%Section SCF::Mixing
      !%Description
      !% Selects what should be mixed during the SCF cycle.  Note that
      !% currently the exact-exchange part of hybrid functionals is not
      !% mixed at all, which would require wavefunction-mixing, not yet
      !% implemented. This may lead to instabilities in the SCF cycle,
      !% so starting from a converged LDA/GGA calculation is recommended
      !% for hybrid functionals. The default depends on the <tt>TheoryLevel</tt>
      !% and the exchange-correlation potential used.
      !% This is not used in case of imaginary-time evolution.
      !%Option none 0
      !% No mixing is done. This is the default for independent
      !% particles.
      !%Option potential 1
      !% The Kohn-Sham potential is mixed. This is the default for other cases.
      !%Option density 2
      !% Mix the density.
      !%Option states 3
      !% (Experimental) Mix the states. In this case, the mixing is always linear.
      !%End
 
      mixdefault = option__mixfield__potential
      if(hm%theory_level == independent_particles) mixdefault = option__mixfield__none
 
      call parse_variable(namespace, 'MixField', mixdefault, scf%mix_field)
      if(.not.varinfo_valid_option('MixField', scf%mix_field)) call messages_input_error(namespace, 'MixField')
      call messages_print_var_option('MixField', scf%mix_field, "what to mix during SCF cycles", namespace=namespace)
 
      if (scf%mix_field == option__mixfield__potential .and. hm%theory_level == independent_particles) then
        call messages_write('Input: Cannot mix the potential for non-interacting particles.')
        call messages_fatal(namespace=namespace)
      end if
 
      if (scf%mix_field == option__mixfield__potential .and. hm%pcm%run_pcm) then
        call messages_write('Input: You have selected to mix the potential.', new_line = .true.)
        call messages_write('       This might produce convergence problems for solvated systems.', new_line = .true.)
        call messages_write('       Mix the Density instead.')
        call messages_warning(namespace=namespace)
      end if
 
      if(scf%mix_field == option__mixfield__density &
        .and. bitand(hm%xc%family, xc_family_oep + xc_family_mgga + xc_family_hyb_mgga + xc_family_nc_mgga) /= 0) then
 
        call messages_write('Input: You have selected to mix the density with OEP or MGGA XC functionals.', new_line = .true.)
        call messages_write('       This might produce convergence problems. Mix the potential instead.')
        call messages_warning(namespace=namespace)
      end if
 
      if(scf%mix_field == option__mixfield__states) then
        call messages_experimental('MixField = states', namespace=namespace)
      end if
 
      ! Handle mixing now...
      select case(scf%mix_field)
      case (option__mixfield__potential, option__mixfield__density)
        scf%mixdim1 = gr%np
      case(option__mixfield__states)
        ! we do not really need the mixer, except for the value of the mixing coefficient
        scf%mixdim1 = 1
      end select
 
      mix_type = type_float
 
      if (scf%mix_field /= option__mixfield__none) then
        call mix_init(scf%smix, namespace, space, gr%der, scf%mixdim1, st%d%nspin, func_type_ = mix_type)
      end if
 
      ! If we use DFT+U, we also have do mix it
      if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none ) then
        call lda_u_mixer_init(hm%lda_u, scf%lda_u_mix, st)
        call lda_u_mixer_init_auxmixer(hm%lda_u, namespace, scf%lda_u_mix, scf%smix, st)
      end if
 
      ! If we use tau-dependent MGGA, we need to mix vtau
      if(scf%mix_field == option__mixfield__potential) then
        call vtau_mixer_init_auxmixer(namespace, scf%vtau_mix, scf%smix, hm, gr%np, st%d%nspin)
      end if
 
      call mix_get_field(scf%smix, scf%mixfield)
    else
      scf%mix_field = option__mixfield__none
    end if
 
    !%Variable SCFinLCAO
    !%Type logical
    !%Default no
    !%Section SCF
    !%Description
    !% Performs the SCF cycle with the calculation restricted to the LCAO subspace.
    !% This may be useful for systems with convergence problems (first do a
    !% calculation within the LCAO subspace, then restart from that point for
    !% an unrestricted calculation).
    !%End
    call parse_variable(namespace, 'SCFinLCAO', .false., scf%lcao_restricted)
    if(scf%lcao_restricted) then
      call messages_experimental('SCFinLCAO', namespace=namespace)
      message(1) = 'Info: SCF restricted to LCAO subspace.'
      call messages_info(1, namespace=namespace)
 
      if(scf%conv_eigen_error) then
        message(1) = "ConvEigenError cannot be used with SCFinLCAO, since error is unknown."
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
 
    !%Variable SCFCalculateForces
    !%Type logical
    !%Section SCF
    !%Description
    !% This variable controls whether the forces on the ions are
    !% calculated at the end of a self-consistent iteration. The
    !% default is yes, unless the system only has user-defined
    !% species.
    !%End
    call parse_variable(namespace, 'SCFCalculateForces', .not. ions%only_user_def, scf%calc_force)
 
    if(scf%calc_force .and. gr%der%boundaries%spiralBC) then
      message(1) = 'Forces cannot be calculated when using spiral boundary conditions.'
      write(message(2),'(a)') 'Please use SCFCalculateForces = no.'
      call messages_fatal(2, namespace=namespace)
    end if
    if(scf%calc_force) then
      if (allocated(hm%ep%b_field) .or. allocated(hm%ep%a_static)) then
        write(message(1),'(a)') 'The forces are currently not properly calculated if static'
        write(message(2),'(a)') 'magnetic fields or static vector potentials are present.'
        write(message(3),'(a)') 'Please use SCFCalculateForces = no.'
        call messages_fatal(3, namespace=namespace)
      end if
    end if
 
    !%Variable SCFCalculateStress
    !%Type logical
    !%Section SCF
    !%Description
    !% This variable controls whether the stress on the lattice is
    !% calculated at the end of a self-consistent iteration. The
    !% default is no.
    !%End
    call parse_variable(namespace, 'SCFCalculateStress', .false. , scf%calc_stress)
 
    !%Variable SCFCalculateDipole
    !%Type logical
    !%Section SCF
    !%Description
    !% This variable controls whether the dipole is calculated at the
    !% end of a self-consistent iteration. For finite systems the
    !% default is yes. For periodic systems the default is no, unless
    !% an electric field is being applied in a periodic direction.
    !% The single-point Berry`s phase approximation is used for
    !% periodic directions. Ref:
    !% E Yaschenko, L Fu, L Resca, and R Resta, <i>Phys. Rev. B</i> <b>58</b>, 1222-1229 (1998).
    !%End
    call parse_variable(namespace, 'SCFCalculateDipole', .not. space%is_periodic(), scf%calc_dipole)
    if (allocated(hm%vberry)) scf%calc_dipole = .true.
 
    !%Variable SCFCalculatePartialCharges
    !%Type logical
    !%Default no
    !%Section SCF
    !%Description
    !% (Experimental) This variable controls whether partial charges
    !% are calculated at the end of a self-consistent iteration.
    !%End
    call parse_variable(namespace, 'SCFCalculatePartialCharges', .false., scf%calc_partial_charges)
    if (scf%calc_partial_charges) call messages_experimental('SCFCalculatePartialCharges', namespace=namespace)
 
    rmin = ions%min_distance()
 
    !%Variable LocalMagneticMomentsSphereRadius
    !%Type float
    !%Section Output
    !%Description
    !% The local magnetic moments are calculated by integrating the
    !% magnetization density in spheres centered around each atom.
    !% This variable controls the radius of the spheres.
    !% The default is half the minimum distance between two atoms
    !% in the input coordinates, or 100 a.u. if there is only one atom (for isolated systems).
    !%End
    call parse_variable(namespace, 'LocalMagneticMomentsSphereRadius', min(m_half*rmin, lmm_r_single_atom), scf%lmm_r, &
      unit=units_inp%length)
    ! this variable is also used in td/td_write.F90
 
    scf%forced_finish = .false.
 
    pop_sub(scf_init)
  end subroutine scf_init
 
 
  ! ---------------------------------------------------------
  subroutine scf_end(scf)
    type(scf_t),  intent(inout) :: scf
 
    class(convergence_criterion_t), pointer    :: crit
    type(criterion_iterator_t) :: iter
 
    push_sub(scf_end)
 
    call eigensolver_end(scf%eigens)
 
    if(scf%mix_field /= option__mixfield__none) call mix_end(scf%smix)
 
    nullify(scf%mixfield)
 
    if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none) then
      call lda_u_mixer_end(scf%lda_u_mix, scf%smix)
      call vtau_mixer_end(scf%vtau_mix, scf%smix)
    end if
 
    call iter%start(scf%criterion_list)
    do while (iter%has_next())
      crit => iter%get_next()
      safe_deallocate_p(crit)
    end do
 
    pop_sub(scf_end)
  end subroutine scf_end
 
 
  ! ---------------------------------------------------------
  subroutine scf_mix_clear(scf)
    type(scf_t), intent(inout) :: scf
 
    push_sub(scf_mix_clear)
 
    call mix_clear(scf%smix)
 
    if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none) then
      call lda_u_mixer_clear(scf%lda_u_mix, scf%smix)
      call vtau_mixer_clear(scf%vtau_mix, scf%smix)
    end if
 
    pop_sub(scf_mix_clear)
  end subroutine scf_mix_clear
 
  ! ---------------------------------------------------------
  subroutine scf_load(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, restart_load)
    type(scf_t),               intent(inout) :: scf
    type(namespace_t),         intent(in)    :: namespace
    type(electron_space_t),    intent(in)    :: space
    type(multicomm_t),         intent(in)    :: mc
    type(grid_t),              intent(inout) :: gr
    type(ions_t),              intent(in)    :: ions
    type(partner_list_t),      intent(in)    :: ext_partners
    type(states_elec_t),       intent(inout) :: st
    type(v_ks_t),              intent(inout) :: ks
    type(hamiltonian_elec_t),  intent(inout) :: hm
    type(restart_t),           intent(in)    :: restart_load
 
    integer :: ierr, is
 
    push_sub(scf_load)
 
    if (restart_has_flag(restart_load, restart_flag_rho)) then
      ! Load density and used it to recalculated the KS potential.
      call states_elec_load_rho(restart_load, space, st, gr, ierr)
      if (ierr /= 0) then
        message(1) = 'Unable to read density. Density will be calculated from states.'
        call messages_warning(1, namespace=namespace)
      else
        if (bitand(ks%xc_family, xc_family_oep) == 0) then
          call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
        else
          if (.not. restart_has_flag(restart_load, restart_flag_vhxc) .and. ks%oep%level /= oep_level_full) then
            call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
          end if
        end if
      end if
    end if
 
    if (restart_has_flag(restart_load, restart_flag_vhxc)) then
      call hm%ks_pot%load(restart_load, space, gr, ierr)
      if (ierr /= 0) then
        message(1) = 'Unable to read Vhxc. Vhxc will be calculated from states.'
        call messages_warning(1, namespace=namespace)
      else
        call hm%update(gr, namespace, space, ext_partners)
        if (bitand(ks%xc_family, xc_family_oep) /= 0) then
          if (ks%oep%level == oep_level_full) then
            do is = 1, st%d%nspin
              ks%oep%vxc(1:gr%np, is) = hm%ks_pot%vhxc(1:gr%np, is) - hm%ks_pot%vhartree(1:gr%np)
            end do
            call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
          end if
        end if
      end if
    end if
 
    if (restart_has_flag(restart_load, restart_flag_mix)) then
      if (scf%mix_field == option__mixfield__density .or. scf%mix_field == option__mixfield__potential) then
        call mix_load(namespace, restart_load, scf%smix, space, gr, ierr)
        if (ierr /= 0) then
          message(1) = "Unable to read mixing information. Mixing will start from scratch."
          call messages_warning(1, namespace=namespace)
        end if
      end if
    end if
 
    if(hm%lda_u_level /= dft_u_none) then
      call lda_u_load(restart_load, hm%lda_u, st, hm%energy%dft_u, ierr)
      if (ierr /= 0) then
        message(1) = "Unable to read DFT+U information. DFT+U data will be calculated from states."
        call messages_warning(1, namespace=namespace)
      end if
 
      ! As v_ks_calc has already been called, we need to update hm%energy%int_dft_u
      call v_ks_update_dftu_energy(ks, namespace, hm, st, hm%energy%int_dft_u)
    end if
 
    !TODO: Create a dedicated routine and call it from the initialize
 
!      if (present(outp) .and. mpi_grp_is_root(mpi_world)) then
!        call io_rm(STATIC_DIR //"info")
!      end if
!    end if
 
    pop_sub(scf_load)
  end subroutine scf_load
 
  ! ---------------------------------------------------------
  subroutine scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
    type(scf_t),               intent(inout) :: scf
    type(namespace_t),         intent(in)    :: namespace
    type(electron_space_t),    intent(in)    :: space
    type(grid_t),              intent(inout) :: gr
    type(ions_t),              intent(inout) :: ions
    type(states_elec_t),       intent(inout) :: st
    type(v_ks_t),              intent(inout) :: ks
    type(hamiltonian_elec_t),  intent(inout) :: hm
    type(output_t),  optional, intent(in)    :: outp
    integer,         optional, intent(in)    :: verbosity
 
    integer :: ib, iqn
 
    push_sub(scf_start)
 
    if(scf%forced_finish) then
      message(1) = "Previous clean stop, not doing SCF and quitting."
      call messages_fatal(1, only_root_writes = .true., namespace=namespace)
    end if
 
    if (.not. hm%is_hermitian()) then
      message(1) = "Trying to run a SCF calculation for a non-hermitian Hamiltonian. This is not supported."
      call messages_fatal(1, namespace=namespace)
    end if
 
    scf%verbosity_ = optional_default(verbosity, verb_full)
 
    scf%output_during_scf = .false.
    scf%output_forces = .false.
    scf%calc_current = .false.
 
    if (present(outp)) then
      ! if the user has activated output=stress but not SCFCalculateStress,
      ! we assume that is implied
      if (outp%what(option__output__stress)) then
        scf%calc_stress = .true.
      end if
 
      scf%output_during_scf = outp%duringscf
      scf%calc_current = output_needs_current(outp, states_are_real(st))
 
      if (outp%duringscf .and. outp%what(option__output__forces)) then
        scf%output_forces = .true.
      end if
    end if
 
    if(scf%lcao_restricted) then
      call lcao_init(scf%lcao, namespace, space, gr, ions, st, 1)
      if(.not. lcao_is_available(scf%lcao)) then
        message(1) = 'LCAO is not available. Cannot do SCF in LCAO.'
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    safe_allocate(scf%rhoout(1:gr%np, 1:st%d%nspin))
    safe_allocate(scf%rhoin (1:gr%np, 1:st%d%nspin))
 
    call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
    scf%rhoout = m_zero
 
    if (scf%calc_force .or. scf%output_forces) then
      !We store the Hxc potential for the contribution to the forces
      safe_allocate(scf%vhxc_old(1:gr%np, 1:st%d%nspin))
      call lalg_copy(gr%np, st%d%nspin, hm%ks_pot%vhxc, scf%vhxc_old)
    end if
 
 
    select case(scf%mix_field)
    case(option__mixfield__potential)
      call mixfield_set_vin(scf%mixfield, hm%ks_pot%vhxc)
      call vtau_mixer_set_vin(scf%vtau_mix, hm)
    case(option__mixfield__density)
      call mixfield_set_vin(scf%mixfield, scf%rhoin)
 
    case(option__mixfield__states)
 
      ! There is a ICE with foss2022a-serial. I am changing to allocate - NTD
      allocate(wfs_elec_t::scf%psioutb (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
 
      do iqn = st%d%kpt%start, st%d%kpt%end
        do ib = st%group%block_start, st%group%block_end
          call st%group%psib(ib, iqn)%copy_to(scf%psioutb(ib, iqn))
        end do
      end do
 
    end select
 
    call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
    ! If we use DFT+U, we also have do mix it
    if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none) then
      call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
    end if
 
    call create_convergence_file(static_dir, "convergence")
 
    if ( scf%verbosity_ /= verb_no ) then
      if(scf%max_iter > 0) then
        write(message(1),'(a)') 'Info: Starting SCF iteration.'
      else
        write(message(1),'(a)') 'Info: No SCF iterations will be done.'
        ! we cannot tell whether it is converged.
        scf%finish = .false.
      end if
      call messages_info(1, namespace=namespace)
    end if
 
    scf%converged_current = .false.
    scf%matvec = 0
 
    pop_sub(scf_start)
 
  contains
 
    ! -----------------------------------------------------
 
    subroutine create_convergence_file(dir, fname)
      character(len=*), intent(in) :: dir
      character(len=*), intent(in) :: fname
 
      integer :: iunit
      character(len=12) :: label
      if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes
        call io_mkdir(dir, namespace)
        iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write')
        write(iunit, '(a)', advance = 'no') '#iter energy           '
        label = 'energy_diff'
        write(iunit, '(1x,a)', advance = 'no') label
        label = 'abs_dens'
        write(iunit, '(1x,a)', advance = 'no') label
        label = 'rel_dens'
        write(iunit, '(1x,a)', advance = 'no') label
        label = 'abs_ev'
        write(iunit, '(1x,a)', advance = 'no') label
        label = 'rel_ev'
        write(iunit, '(1x,a)', advance = 'no') label
        if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
          .and. ks%theory_level /= generalized_kohn_sham_dft) then
          if (ks%oep%level == oep_level_full) then
            label = 'OEP norm2ss'
            write(iunit, '(1x,a)', advance = 'no') label
          end if
        end if
        write(iunit,'(a)') ''
        call io_close(iunit)
      end if
 
    end subroutine create_convergence_file
 
  end subroutine scf_start
 
  ! ---------------------------------------------------------
  subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
    verbosity, iters_done, restart_dump)
    type(scf_t),               intent(inout) :: scf
    type(namespace_t),         intent(in)    :: namespace
    type(electron_space_t),    intent(in)    :: space
    type(multicomm_t),         intent(in)    :: mc
    type(grid_t),              intent(inout) :: gr
    type(ions_t),              intent(inout) :: ions
    type(partner_list_t),      intent(in)    :: ext_partners
    type(states_elec_t),       intent(inout) :: st
    type(v_ks_t),              intent(inout) :: ks
    type(hamiltonian_elec_t),  intent(inout) :: hm
    type(output_t),  optional, intent(in)    :: outp
    integer,         optional, intent(in)    :: verbosity
    integer,         optional, intent(out)   :: iters_done
    type(restart_t), optional, intent(in)    :: restart_dump
 
    integer :: iter
    logical :: completed
 
    push_sub(scf_run)
 
    call scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
 
    ! SCF cycle
    do iter = 1, scf%max_iter
 
      call scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
        verbosity, iters_done, restart_dump)
 
      completed = scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
 
      if(scf%forced_finish .or. completed) then
        exit
      end if
    end do
 
    call scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
      verbosity, iters_done, restart_dump)
 
    pop_sub(scf_run)
  end subroutine scf_run
 
  ! ---------------------------------------------------------
  subroutine scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
    verbosity, iters_done, restart_dump)
    type(scf_t),               intent(inout) :: scf
    type(namespace_t),         intent(in)    :: namespace
    type(electron_space_t),    intent(in)    :: space
    type(multicomm_t),         intent(in)    :: mc
    type(grid_t),              intent(inout) :: gr
    type(ions_t),              intent(inout) :: ions
    type(partner_list_t),      intent(in)    :: ext_partners
    type(states_elec_t),       intent(inout) :: st
    type(v_ks_t),              intent(inout) :: ks
    type(hamiltonian_elec_t),  intent(inout) :: hm
    integer,                   intent(in)    :: iter
    type(output_t),  optional, intent(in)    :: outp
    integer,         optional, intent(in)    :: verbosity
    integer,         optional, intent(out)   :: iters_done
    type(restart_t), optional, intent(in)    :: restart_dump
 
    integer :: iqn, ib, ierr
    class(convergence_criterion_t), pointer    :: crit
    type(criterion_iterator_t) :: iterator
    logical :: is_crit_conv
    real(real64) :: etime, itime
 
    push_sub(scf_iter)
 
    call profiling_in("SCF_CYCLE")
 
    itime = loct_clock()
 
    ! this initialization seems redundant but avoids improper optimization at -O3 by PGI 7 on chum,
    ! which would cause a failure of testsuite/linear_response/04-vib_modes.03-vib_modes_fd.inp
    scf%eigens%converged = 0
 
    ! keep the information about the spectrum up to date, needed e.g. for Chebyshev expansion for imaginary time
    call hm%update_span(gr%spacing(1:space%dim), minval(st%eigenval(:, :)), namespace)
 
    !We update the quantities at the begining of the scf cycle
    if (iter == 1) then
      scf%evsum_in = states_elec_eigenvalues_sum(st)
    end if
    call iterator%start(scf%criterion_list)
    do while (iterator%has_next())
      crit => iterator%get_next()
      call scf_update_initial_quantity(scf, hm, crit)
    end do
 
    if (scf%calc_force .or. scf%output_forces) then
      !Used for computing the imperfect convegence contribution to the forces
      scf%vhxc_old(1:gr%np, 1:st%d%nspin) = hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin)
    end if
 
    if(scf%lcao_restricted) then
      call lcao_init_orbitals(scf%lcao, namespace, st, gr, ions)
      call lcao_wf(scf%lcao, st, gr, ions, hm, namespace)
    else
 
      !We check if the system is coupled with a partner that requires self-consistency
      !  if(hamiltonian_has_scf_partner(hm)) then
      if (allocated(hm%vberry)) then
        !In this case, v_Hxc is frozen and we do an internal SCF loop over the
        ! partners that require SCF
        ks%frozen_hxc = .true.
        ! call perform_scf_partners()
        call berry_perform_internal_scf(scf%berry, namespace, space, scf%eigens, gr, st, hm, iter, ks, ions, ext_partners)
        !and we unfreeze the potential once finished
        ks%frozen_hxc = .false.
      else
        scf%eigens%converged = 0
        call scf%eigens%run(namespace, gr, st, hm, iter)
      end if
    end if
 
    scf%matvec = scf%matvec + scf%eigens%matvec
 
    ! occupations
    call states_elec_fermi(st, namespace, gr)
    call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
 
    ! compute output density, potential (if needed) and eigenvalues sum
    call density_calc(st, gr, st%rho)
 
    call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoout)
 
    select case (scf%mix_field)
    case (option__mixfield__potential)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
      call mixfield_set_vout(scf%mixfield, hm%ks_pot%vhxc)
      call vtau_mixer_set_vout(scf%vtau_mix, hm)
    case (option__mixfield__density)
      call mixfield_set_vout(scf%mixfield, scf%rhoout)
    case(option__mixfield__states)
 
      do iqn = st%d%kpt%start, st%d%kpt%end
        do ib = st%group%block_start, st%group%block_end
          call st%group%psib(ib, iqn)%copy_data_to(gr%np, scf%psioutb(ib, iqn))
        end do
      end do
    end select
 
    if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none) then
      call lda_u_mixer_set_vout(hm%lda_u, scf%lda_u_mix)
    endif
 
    ! recalculate total energy
    call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit = -1)
 
    if (present(outp)) then
      ! compute forces only if requested
      if (outp%duringscf .and. outp%what_now(option__output__forces, iter)) then
        call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
      end if
    end if
 
    !We update the quantities at the end of the scf cycle
    call iterator%start(scf%criterion_list)
    do while (iterator%has_next())
      crit => iterator%get_next()
      call scf_update_diff_quantity(scf, hm, st, gr, scf%rhoout, scf%rhoin, crit)
    end do
 
    ! are we finished?
    scf%converged_last = scf%converged_current
 
    scf%converged_current = scf%check_conv .and. &
      (.not. scf%conv_eigen_error .or. all(scf%eigens%converged >= st%nst_conv))
    !Loop over the different criteria
    call iterator%start(scf%criterion_list)
    do while (iterator%has_next())
      crit => iterator%get_next()
      call crit%is_converged(is_crit_conv)
      scf%converged_current = scf%converged_current .and. is_crit_conv
    end do
 
    ! only finish if the convergence criteria are fulfilled in two
    ! consecutive iterations
    scf%finish = scf%converged_last .and. scf%converged_current
 
    etime = loct_clock() - itime
    call scf_write_iter(namespace)
 
    ! mixing
    select case (scf%mix_field)
    case (option__mixfield__density)
      ! mix input and output densities and compute new potential
      call mixing(namespace, scf%smix)
      call mixfield_get_vnew(scf%mixfield, st%rho)
      ! for spinors, having components 3 or 4 be negative is not unphysical
      if (minval(st%rho(1:gr%np, 1:st%d%spin_channels)) < -1e-6_real64) then
        write(message(1),*) 'Negative density after mixing. Minimum value = ', &
          minval(st%rho(1:gr%np, 1:st%d%spin_channels))
        call messages_warning(1, namespace=namespace)
      end if
      call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
    case (option__mixfield__potential)
      ! mix input and output potentials
      call mixing(namespace, scf%smix)
      call mixfield_get_vnew(scf%mixfield, hm%ks_pot%vhxc)
      call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
      call vtau_mixer_get_vnew(scf%vtau_mix, hm)
      call hamiltonian_elec_update_pot(hm, gr)
 
    case(option__mixfield__states)
 
      do iqn = st%d%kpt%start, st%d%kpt%end
        do ib = st%group%block_start, st%group%block_end
          call batch_scal(gr%np, m_one - mix_coefficient(scf%smix), st%group%psib(ib, iqn))
          call batch_axpy(gr%np, mix_coefficient(scf%smix), scf%psioutb(ib, iqn), st%group%psib(ib, iqn))
        end do
      end do
 
      call density_calc(st, gr, st%rho)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
 
    case (option__mixfield__none)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
    end select
 
 
    ! Are we asked to stop? (Whenever Fortran is ready for signals, this should go away)
    scf%forced_finish = clean_stop(mc%master_comm) .or. walltimer_alarm(mc%master_comm)
 
    if (scf%finish .and. st%modelmbparticles%nparticle > 0) then
      call modelmb_sym_all_states(space, gr, st)
    end if
 
    if (present(outp) .and. present(restart_dump)) then
      ! save restart information
 
      if ( (scf%finish .or. (modulo(iter, outp%restart_write_interval) == 0) &
        .or. iter == scf%max_iter .or. scf%forced_finish) ) then
 
        call states_elec_dump(restart_dump, space, st, gr, hm%kpoints, ierr, iter=iter)
        if (ierr /= 0) then
          message(1) = 'Unable to write states wavefunctions.'
          call messages_warning(1, namespace=namespace)
        end if
 
        call states_elec_dump_rho(restart_dump, space, st, gr, ierr, iter=iter)
        if (ierr /= 0) then
          message(1) = 'Unable to write density.'
          call messages_warning(1, namespace=namespace)
        end if
 
        if(hm%lda_u_level /= dft_u_none) then
          call lda_u_dump(restart_dump, namespace, hm%lda_u, st, gr, ierr)
          if (ierr /= 0) then
            message(1) = 'Unable to write DFT+U information.'
            call messages_warning(1, namespace=namespace)
          end if
        end if
 
        select case (scf%mix_field)
        case (option__mixfield__density)
          call mix_dump(namespace, restart_dump, scf%smix, space, gr, ierr)
          if (ierr /= 0) then
            message(1) = 'Unable to write mixing information.'
            call messages_warning(1, namespace=namespace)
          end if
        case (option__mixfield__potential)
          call hm%ks_pot%dump(restart_dump, space, gr, ierr)
          if (ierr /= 0) then
            message(1) = 'Unable to write Vhxc.'
            call messages_warning(1, namespace=namespace)
          end if
 
          call mix_dump(namespace, restart_dump, scf%smix, space, gr, ierr)
          if (ierr /= 0) then
            message(1) = 'Unable to write mixing information.'
            call messages_warning(1, namespace=namespace)
          end if
        end select
      end if
    end if
 
    call write_convergence_file(static_dir, "convergence")
 
    call profiling_out("SCF_CYCLE")
 
    pop_sub(scf_iter)
  contains
 
    ! ---------------------------------------------------------
    subroutine scf_write_iter(namespace)
      type(namespace_t), intent(in)  :: namespace
 
      character(len=50) :: str
      real(real64) :: dipole(1:space%dim)
 
      push_sub(scf_run.scf_write_iter)
 
      if ( scf%verbosity_ == verb_full ) then
 
        write(str, '(a,i5)') 'SCF CYCLE ITER #' ,iter
        call messages_print_with_emphasis(msg=trim(str), namespace=namespace)
        write(message(1),'(a,es15.8,2(a,es9.2))') ' etot  = ', units_from_atomic(units_out%energy, hm%energy%total), &
          ' abs_ev   = ', units_from_atomic(units_out%energy, scf%evsum_diff), &
          ' rel_ev   = ', scf%evsum_diff/(abs(scf%evsum_out)+1e-20)
        write(message(2),'(a,es15.2,2(a,es9.2))') &
          ' ediff = ', scf%energy_diff, ' abs_dens = ', scf%abs_dens_diff, &
          ' rel_dens = ', scf%abs_dens_diff/st%qtot
        call messages_info(2, namespace=namespace)
 
        if(.not.scf%lcao_restricted) then
          write(message(1),'(a,i0)') 'Matrix vector products: ', scf%eigens%matvec
          write(message(2),'(a,i0)') 'Converged eigenvectors: ', sum(scf%eigens%converged(1:st%nik))
          call messages_info(2, namespace=namespace)
          call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, scf%eigens%diff, compact = .true., namespace=namespace)
        else
          call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, compact = .true., namespace=namespace)
        end if
 
        if (allocated(hm%vberry)) then
          call calc_dipole(dipole, space, gr, st, ions)
          call write_dipole(st, hm, space, dipole, namespace=namespace)
        end if
 
        if(st%d%ispin > unpolarized) then
          call write_magnetic_moments(gr, st, ions, gr%der%boundaries, scf%lmm_r, namespace=namespace)
        end if
 
        if(hm%lda_u_level == dft_u_acbn0) then
          call lda_u_write_u(hm%lda_u, namespace=namespace)
          call lda_u_write_v(hm%lda_u, namespace=namespace)
        end if
 
        write(message(1),'(a)') ''
        write(message(2),'(a,i5,a,f14.2)') 'Elapsed time for SCF step ', iter,':', etime
        call messages_info(2, namespace=namespace)
 
        call scf_print_mem_use(namespace)
 
        call messages_print_with_emphasis(namespace=namespace)
 
      end if
 
      if ( scf%verbosity_ == verb_compact ) then
        write(message(1),'(a,i4,a,es15.8, a,es9.2, a, f7.1, a)') &
          'iter ', iter, &
          ' : etot ', units_from_atomic(units_out%energy, hm%energy%total), &
          ' : abs_dens', scf%abs_dens_diff, &
          ' : etime ', etime, 's'
        call messages_info(1, namespace=namespace)
      end if
 
      pop_sub(scf_run.scf_write_iter)
    end subroutine scf_write_iter
 
 
    ! -----------------------------------------------------
    subroutine write_convergence_file(dir, fname)
      character(len=*), intent(in) :: dir
      character(len=*), intent(in) :: fname
 
      integer :: iunit
 
      if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes
        call io_mkdir(dir, namespace)
        iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write', position='append')
        write(iunit, '(i5,es18.8)', advance = 'no') iter, units_from_atomic(units_out%energy, hm%energy%total)
        call iterator%start(scf%criterion_list)
        do while (iterator%has_next())
          crit => iterator%get_next()
          select type (crit)
          type is (energy_criterion_t)
            write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, crit%val_abs)
          type is (density_criterion_t)
            write(iunit, '(2es13.5)', advance = 'no') crit%val_abs, crit%val_rel
          type is (eigenval_criterion_t)
            write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, crit%val_abs)
            write(iunit, '(es13.5)', advance = 'no') crit%val_rel
          class default
            assert(.false.)
          end select
        end do
        if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
          .and. ks%theory_level /= generalized_kohn_sham_dft) then
          if (ks%oep%level == oep_level_full) then
            write(iunit, '(es13.5)', advance = 'no') ks%oep%norm2ss
          end if
        end if
        write(iunit,'(a)') ''
        call io_close(iunit)
      end if
    end subroutine write_convergence_file
 
  end subroutine scf_iter
 
  logical function scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, &
    verbosity, iters_done)  result(completed)
    type(scf_t),               intent(inout) :: scf
    type(namespace_t),         intent(in)    :: namespace
    type(electron_space_t),    intent(in)    :: space
    type(grid_t),              intent(inout) :: gr
    type(ions_t),              intent(inout) :: ions
    type(states_elec_t),       intent(inout) :: st
    type(v_ks_t),              intent(inout) :: ks
    type(hamiltonian_elec_t),  intent(inout) :: hm
    integer,                   intent(in)    :: iter
    type(output_t),  optional, intent(in)    :: outp
    integer,         optional, intent(in)    :: verbosity
    integer,         optional, intent(out)   :: iters_done
 
    character(len=MAX_PATH_LEN) :: dirname
    integer(int64) :: what_i
 
    push_sub(scf_iter_finish)
 
    completed = .false.
 
    if(scf%finish) then
      if(present(iters_done)) iters_done = iter
      if(scf%verbosity_ >= verb_compact) then
        write(message(1), '(a, i4, a)') 'Info: SCF converged in ', iter, ' iterations'
        write(message(2), '(a)')        ''
        call messages_info(2, namespace=namespace)
      end if
      completed = .true.
      pop_sub(scf_iter_finish)
      return
    end if
    if (present(outp)) then
      if (any(outp%what) .and. outp%duringscf) then
        do what_i = lbound(outp%what, 1), ubound(outp%what, 1)
          if (outp%what_now(what_i, iter)) then
            write(dirname,'(a,a,i4.4)') trim(outp%iter_dir),"scf.", iter
            call output_all(outp, namespace, space, dirname, gr, ions, iter, st, hm, ks)
            call output_modelmb(outp, namespace, space, dirname, gr, ions, iter, st)
            exit
          end if
        end do
      end if
    end if
 
    ! save information for the next iteration
    call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
 
    ! restart mixing
    if (scf%mix_field /= option__mixfield__none) then
      if (scf%smix%ns_restart > 0) then
        if (mod(iter, scf%smix%ns_restart) == 0) then
          message(1) = "Info: restarting mixing."
          call messages_info(1, namespace=namespace)
          call scf_mix_clear(scf)
        end if
      end if
    end if
 
    select case(scf%mix_field)
    case(option__mixfield__potential)
      call mixfield_set_vin(scf%mixfield, hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin))
      call vtau_mixer_set_vin(scf%vtau_mix, hm)
    case (option__mixfield__density)
      call mixfield_set_vin(scf%mixfield, scf%rhoin)
    end select
 
    !If we use LDA+U, we also have do mix it
    if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none) then
      call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
    end if
 
    ! check if debug mode should be enabled or disabled on the fly
    call io_debug_on_the_fly(namespace)
 
    pop_sub(scf_iter_finish)
  end function scf_iter_finish
 
  ! ---------------------------------------------------------
  subroutine scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
    verbosity, iters_done, restart_dump)
    type(scf_t),               intent(inout) :: scf
    type(namespace_t),         intent(in)    :: namespace
    type(electron_space_t),    intent(in)    :: space
    type(multicomm_t),         intent(in)    :: mc
    type(grid_t),              intent(inout) :: gr
    type(ions_t),              intent(inout) :: ions
    type(partner_list_t),      intent(in)    :: ext_partners
    type(states_elec_t),       intent(inout) :: st
    type(v_ks_t),              intent(inout) :: ks
    type(hamiltonian_elec_t),  intent(inout) :: hm
    integer,                   intent(in)    :: iter
    type(output_t),  optional, intent(in)    :: outp
    integer,         optional, intent(in)    :: verbosity
    integer,         optional, intent(out)   :: iters_done
    type(restart_t), optional, intent(in)    :: restart_dump
 
    integer :: iqn, ib
    class(convergence_criterion_t), pointer    :: crit
    type(criterion_iterator_t) :: iterator
 
 
    push_sub(scf_finish)
 
    if(scf%lcao_restricted) call lcao_end(scf%lcao)
 
    ! Compute the KS potential corresponding to the final density
    ! This is critical for getting consistent TD calculations
    if ((scf%max_iter > 0 .and. scf%mix_field == option__mixfield__potential) .or. scf%calc_current) then
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
        calc_current=scf%calc_current)
    end if
 
    select case(scf%mix_field)
    case(option__mixfield__states)
 
      do iqn = st%d%kpt%start, st%d%kpt%end
        do ib = st%group%block_start, st%group%block_end
          call scf%psioutb(ib, iqn)%end()
        end do
      end do
 
      ! There is a ICE with foss2022a-serial. I am changing to deallocate - NTD
      deallocate(scf%psioutb)
    end select
 
    safe_deallocate_a(scf%rhoout)
    safe_deallocate_a(scf%rhoin)
 
    if (scf%max_iter > 0 .and. any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted) then
      write(message(1),'(a)') 'Some of the states are not fully converged!'
      if (all(scf%eigens%converged >= st%nst_conv)) then
        write(message(2),'(a)') 'But all requested states to converge are converged.'
        call messages_info(2, namespace=namespace)
      else
        call messages_warning(1, namespace=namespace)
      end if
    end if
 
    if (.not.scf%finish) then
      write(message(1), '(a,i4,a)') 'SCF *not* converged after ', iter - 1, ' iterations.'
      call messages_warning(1, namespace=namespace)
    end if
 
    write(message(1), '(a,i10)') 'Info: Number of matrix-vector products: ', scf%matvec
    call messages_info(1)
 
    if (scf%calc_force) then
      call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
    end if
 
    if (scf%calc_stress) call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
 
    ! Update the eigenvalues, to match the KS potential that just got recomputed
    if (scf%mix_field == option__mixfield__potential) then
      call energy_calc_eigenvalues(namespace, hm, gr%der, st)
      call states_elec_fermi(st, namespace, gr)
    end if
 
    if(present(outp)) then
      ! output final information
      call scf_write_static(static_dir, "info")
      call output_all(outp, namespace, space, static_dir, gr, ions, -1, st, hm, ks)
      call output_modelmb(outp, namespace, space, static_dir, gr, ions, -1, st)
    end if
 
    if (space%is_periodic() .and. st%nik > st%d%nspin) then
      if (bitand(hm%kpoints%method, kpoints_path) /= 0) then
        call states_elec_write_bandstructure(static_dir, namespace, st%nst, st,  &
          ions, gr, hm%kpoints, hm%phase, vec_pot = hm%hm_base%uniform_vector_potential, &
          vec_pot_var = hm%hm_base%vector_potential)
      end if
    end if
 
    if (ks%vdw%vdw_correction == option__vdwcorrection__vdw_ts) then
      call vdw_ts_write_c6ab(ks%vdw%vdw_ts, ions, static_dir, 'c6ab_eff', namespace)
    end if
 
    safe_deallocate_a(scf%vhxc_old)
 
    pop_sub(scf_finish)
 
  contains
 
    ! ---------------------------------------------------------
    subroutine scf_write_static(dir, fname)
      character(len=*), intent(in) :: dir, fname
 
      integer :: iunit, iatom
      real(real64), allocatable :: hirshfeld_charges(:)
      real(real64) :: dipole(1:space%dim)
      real(real64) :: ex_virial
 
      push_sub(scf_run.scf_write_static)
 
      if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes
        call io_mkdir(dir, namespace)
        iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write')
 
        call grid_write_info(gr, iunit=iunit)
 
        call symmetries_write_info(gr%symm, space, iunit=iunit)
 
        if (space%is_periodic()) then
          call hm%kpoints%write_info(iunit=iunit)
          write(iunit,'(1x)')
        end if
 
        call v_ks_write_info(ks, iunit=iunit)
 
        ! scf information
        if(scf%finish) then
          write(iunit, '(a, i4, a)')'SCF converged in ', iter, ' iterations'
        else
          write(iunit, '(a)') 'SCF *not* converged!'
        end if
        write(iunit, '(1x)')
 
        if(any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted) then
          write(iunit,'(a)') 'Some of the states are not fully converged!'
          if (all(scf%eigens%converged >= st%nst_conv)) then
            write(iunit,'(a)') 'But all requested states to converge are converged.'
          end if
        end if
 
        call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, iunit=iunit)
        write(iunit, '(1x)')
 
        if (space%is_periodic()) then
          call states_elec_write_gaps(iunit, st, space)
          write(iunit, '(1x)')
        end if
 
        write(iunit, '(3a)') 'Energy [', trim(units_abbrev(units_out%energy)), ']:'
      else
        iunit = -1
      end if
 
      call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full = .true.)
 
      if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
      if(st%d%ispin > unpolarized) then
        call write_magnetic_moments(gr, st, ions, gr%der%boundaries, scf%lmm_r, iunit=iunit)
        if (mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
      end if
 
      if(st%d%ispin == spinors .and. space%dim == 3 .and. &
        (ks%theory_level == kohn_sham_dft .or. ks%theory_level == generalized_kohn_sham_dft) ) then
        call write_total_xc_torque(iunit, gr, hm%ks_pot%vxc, st)
        if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
      end if
 
      if(hm%lda_u_level == dft_u_acbn0) then
        call lda_u_write_u(hm%lda_u, iunit=iunit)
        call lda_u_write_v(hm%lda_u, iunit=iunit)
        if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
      end if
 
      if(scf%calc_dipole) then
        call calc_dipole(dipole, space, gr, st, ions)
        call write_dipole(st, hm, space, dipole, iunit=iunit)
      end if
 
      ! This only works when we do not have a correlation part
      if(ks%theory_level == kohn_sham_dft .and. &
        hm%xc%functional(func_c,1)%family == xc_family_none .and. st%d%ispin /= spinors) then
        call energy_calc_virial_ex(gr%der, hm%ks_pot%vxc, st, ex_virial)
 
        if (mpi_grp_is_root(mpi_world)) then
          write(iunit, '(3a)') 'Virial relation for exchange [', trim(units_abbrev(units_out%energy)), ']:'
          write(iunit,'(a,es14.6)') "Energy from the orbitals ", units_from_atomic(units_out%energy, hm%energy%exchange)
          write(iunit,'(a,es14.6)') "Energy from the potential (virial) ", units_from_atomic(units_out%energy, ex_virial)
          write(iunit, '(1x)')
        end if
      end if
 
      if(mpi_grp_is_root(mpi_world)) then
        if(scf%max_iter > 0) then
          write(iunit, '(a)') 'Convergence:'
          call iterator%start(scf%criterion_list)
          do while (iterator%has_next())
            crit => iterator%get_next()
            call crit%write_info(iunit)
          end do
          write(iunit,'(1x)')
        end if
        ! otherwise, these values are uninitialized, and unknown.
 
        if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
          .and. ks%theory_level /= generalized_kohn_sham_dft) then
          if ((ks%oep_photon%level == oep_level_full) .or. (ks%oep_photon%level == oep_level_kli)) then
            write(iunit, '(a)') 'Photon observables:'
            write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon number = ', ks%oep_photon%pt%number(1)
            write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon ex. = ', ks%oep_photon%pt%ex
            write(iunit,'(1x)')
          end if
        end if
 
        if (scf%calc_force) call forces_write_info(iunit, ions, dir, namespace)
 
        if (scf%calc_stress) then
          call output_stress(iunit, space%periodic_dim, st%stress_tensors, all_terms=.false.)
          call output_pressure(iunit, space%periodic_dim, st%stress_tensors%total)
        end if
 
      end if
 
      if(scf%calc_partial_charges) then
        safe_allocate(hirshfeld_charges(1:ions%natoms))
 
        call partial_charges_calculate(gr, st, ions, hirshfeld_charges)
 
        if(mpi_grp_is_root(mpi_world)) then
 
          write(iunit,'(a)') 'Partial ionic charges'
          write(iunit,'(a)') ' Ion                     Hirshfeld'
 
          do iatom = 1, ions%natoms
            write(iunit,'(i4,a10,f16.3)') iatom, trim(ions%atom(iatom)%species%get_label()), hirshfeld_charges(iatom)
 
          end do
 
        end if
 
        safe_deallocate_a(hirshfeld_charges)
 
      end if
 
      if(mpi_grp_is_root(mpi_world)) then
        call io_close(iunit)
      end if
 
      pop_sub(scf_run.scf_write_static)
    end subroutine scf_write_static
 
  end subroutine scf_finish
 
  ! ---------------------------------------------------------
  subroutine scf_state_info(namespace, st)
    type(namespace_t),    intent(in) :: namespace
    class(states_abst_t), intent(in) :: st
 
    push_sub(scf_state_info)
 
    if (states_are_real(st)) then
      call messages_write('Info: SCF using real wavefunctions.')
    else
      call messages_write('Info: SCF using complex wavefunctions.')
    end if
    call messages_info(namespace=namespace)
 
    pop_sub(scf_state_info)
 
  end subroutine scf_state_info
 
  ! ---------------------------------------------------------
  subroutine scf_print_mem_use(namespace)
    type(namespace_t),    intent(in) :: namespace
    real(real64) :: mem
    real(real64) :: mem_tmp
 
    push_sub(scf_print_mem_use)
 
    if(conf%report_memory) then
      mem = loct_get_memory_usage()/(1024.0_real64**2)
      call mpi_world%allreduce(mem, mem_tmp, 1, mpi_double_precision, mpi_sum)
      mem = mem_tmp
      write(message(1),'(a,f14.2)') 'Memory usage [Mbytes]     :', mem
      call messages_info(1, namespace=namespace)
    end if
 
    pop_sub(scf_print_mem_use)
  end subroutine scf_print_mem_use
 
  ! --------------------------------------------------------
  subroutine scf_update_initial_quantity(scf, hm, criterion)
    type(scf_t),                    intent(inout) :: scf
    type(hamiltonian_elec_t),       intent(in)    :: hm
    class(convergence_criterion_t), intent(in)    :: criterion
 
    push_sub(scf_update_initial_quantity)
 
    select type (criterion)
    type is (energy_criterion_t)
      scf%energy_in = hm%energy%total
    type is (density_criterion_t)
      !Do nothing here
    type is (eigenval_criterion_t)
      !Setting of the value is done in the scf_update_diff_quantity routine
    class default
      assert(.false.)
    end select
 
    pop_sub(scf_update_initial_quantity)
  end subroutine scf_update_initial_quantity
 
  ! --------------------------------------------------------
  subroutine scf_update_diff_quantity(scf, hm, st, gr, rhoout, rhoin, criterion)
    type(scf_t),                    intent(inout) :: scf
    type(hamiltonian_elec_t),       intent(in)    :: hm
    type(states_elec_t),            intent(in)    :: st
    type(grid_t),                   intent(in)    :: gr
    real(real64),                   intent(in)    :: rhoout(:,:), rhoin(:,:)
    class(convergence_criterion_t), intent(in)    :: criterion
 
    integer :: is
    real(real64), allocatable :: tmp(:)
 
    push_sub(scf_update_diff_quantity)
 
    select type (criterion)
    type is (energy_criterion_t)
      scf%energy_diff = abs(hm%energy%total - scf%energy_in)
 
    type is (density_criterion_t)
      scf%abs_dens_diff = m_zero
      safe_allocate(tmp(1:gr%np))
      do is = 1, st%d%nspin
        tmp(:) = abs(rhoin(1:gr%np, is) - rhoout(1:gr%np, is))
        scf%abs_dens_diff = scf%abs_dens_diff + dmf_integrate(gr, tmp)
      end do
      safe_deallocate_a(tmp)
 
    type is (eigenval_criterion_t)
      scf%evsum_out = states_elec_eigenvalues_sum(st)
      scf%evsum_diff = abs(scf%evsum_out - scf%evsum_in)
      scf%evsum_in = scf%evsum_out
 
    class default
      assert(.false.)
    end select
 
    pop_sub(scf_update_diff_quantity)
  end subroutine scf_update_diff_quantity
 
  ! ---------------------------------------------------------
  subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
    type(states_elec_t),         intent(in) :: st
    type(hamiltonian_elec_t),    intent(in) :: hm
    type(electron_space_t),      intent(in) :: space
    real(real64),                intent(in) :: dipole(:)
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    push_sub(write_dipole)
 
    if(mpi_grp_is_root(mpi_world)) then
      call output_dipole(dipole, space%dim, iunit=iunit, namespace=namespace)
 
      if (space%is_periodic()) then
        message(1) = "Defined only up to quantum of polarization (e * lattice vector)."
        message(2) = "Single-point Berry's phase method only accurate for large supercells."
        call messages_info(2, iunit=iunit, namespace=namespace)
 
        if (hm%kpoints%full%npoints > 1) then
          message(1) = &
            "WARNING: Single-point Berry's phase method for dipole should not be used when there is more than one k-point."
          message(2) = "Instead, finite differences on k-points (not yet implemented) are needed."
          call messages_info(2, iunit=iunit, namespace=namespace)
        end if
 
        if(.not. smear_is_semiconducting(st%smear)) then
          message(1) = "Single-point Berry's phase dipole calculation not correct without integer occupations."
          call messages_info(1, iunit=iunit, namespace=namespace)
        end if
      end if
 
      call messages_info(iunit=iunit, namespace=namespace)
    end if
 
    pop_sub(write_dipole)
  end subroutine write_dipole
 
 
end module scf_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: