v_ks.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 v_ks_oct_m
  use accel_oct_m
  use types_oct_m
  use current_oct_m
  use debug_oct_m
  use density_oct_m
  use derivatives_oct_m
  use energy_oct_m
  use energy_calc_oct_m
  use electron_space_oct_m
  use exchange_operator_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_elec_base_oct_m
  use interaction_partner_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use isdf_oct_m,        only: isdf_parallel_ace_compute_potentials => isdf_ace_compute_potentials
  use isdf_serial_oct_m, only: isdf_serial_ace_compute_potentials
  use kpoints_oct_m
  use ks_potential_oct_m
  use lalg_basic_oct_m
  use lda_u_oct_m
  use lattice_vectors_oct_m
  use magnetic_oct_m
  use magnetic_constrain_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use pseudo_oct_m
  use pseudopotential_oct_m
  use pcm_potential_oct_m
  use sort_oct_m
  use space_oct_m
  use species_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use states_elec_parallel_oct_m
  use sturm_liouville_oct_m
  use varinfo_oct_m
  use x_slater_oct_m
  use xc_cam_oct_m
  use xc_oct_m
  use xc_f03_lib_m
  use xc_fbe_oct_m
  use xc_functional_oct_m
  use xc_interaction_oct_m
  use xc_ks_inversion_oct_m
  use xc_noncollinear_oct_m
  use xc_oep_oct_m
  use xc_sic_oct_m
  use xc_vxc_oct_m
  use xc_vdw_oct_m
  use v_ks_photon_oct_m
 
  ! from the dftd3 library
  use dftd3_api
 
  implicit none
 
  private
  public ::             &
    v_ks_t,             &
    v_ks_photon_t,      &
    v_ks_init,          &
    v_ks_end,           &
    v_ks_write_info,    &
    v_ks_check_generalized_overlap, &
    v_ks_h_setup,       &
    v_ks_calc,          &
    v_ks_calc_t,        &
    v_ks_calc_start,    &
    v_ks_calc_finish,   &
    v_ks_freeze_hxc,    &
    v_ks_calculate_current, &
    v_ks_update_dftu_energy, &
    v_ks_is_orbital_dependent
 
  type v_ks_calc_t
    private
    logical                           :: calculating
    logical                           :: time_present
    real(real64)                      :: time
    real(real64),         allocatable :: density(:, :)
    logical                           :: total_density_alloc
    real(real64), pointer, contiguous :: total_density(:)
    type(energy_t),       allocatable :: energy
 
    type(states_elec_t),  pointer     :: hf_st
 
    real(real64),         allocatable :: vxc(:, :)
    real(real64),         allocatable :: vtau(:, :)
    real(real64),         allocatable :: axc(:, :, :)
    real(real64),         allocatable :: a_ind(:, :)
    real(real64),         allocatable :: b_ind(:, :)
    logical                           :: calc_energy
  end type v_ks_calc_t
 
  type v_ks_t
    private
    integer,                  public :: theory_level = -1
    logical,                  public :: frozen_hxc = .false. 
 
    integer,                  public :: xc_family = 0
    integer,                  public :: xc_flags = 0
    type(xc_t),               public :: xc
    type(xc_oep_t),           public :: oep
    type(xc_ks_inversion_t),  public :: ks_inversion
    type(xc_sic_t),           public :: sic
    type(xc_vdw_t),           public :: vdw
    type(grid_t), pointer,    public :: gr
    type(sturm_liouville_t),  public :: sl_solver
    type(v_ks_calc_t)                :: calc
    logical                          :: calculate_current = .false.
    type(current_t)                  :: current_calculator
    logical                          :: include_td_field = .false.
 
    real(real64),             public :: stress_xc(3, 3)
    type(v_ks_photon_t),      public :: v_ks_photons
  end type v_ks_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine v_ks_init(ks, namespace, gr, st, ions, mc, space, kpoints)
    type(v_ks_t),            intent(inout) :: ks
    type(namespace_t),       intent(in)    :: namespace
    type(grid_t),    target, intent(inout) :: gr
    type(states_elec_t),     intent(in)    :: st
    type(ions_t),            intent(inout) :: ions
    type(multicomm_t),       intent(in)    :: mc
    class(space_t),          intent(in)    :: space
    type(kpoints_t),         intent(in)    :: kpoints
 
    integer :: x_id, c_id, xk_id, ck_id, default, val
    logical :: parsed_theory_level, using_hartree_fock
    integer :: pseudo_x_functional, pseudo_c_functional
    integer :: oep_type
 
    push_sub(v_ks_init)
 
    ! We need to parse TheoryLevel and XCFunctional, this is
    ! complicated because they are interdependent.
 
    !%Variable TheoryLevel
    !%Type integer
    !%Section Hamiltonian
    !%Description
    !% The calculations can be run with different "theory levels" that
    !% control how electrons are simulated. The default is
    !% <tt>kohn_sham</tt>. When hybrid or meta-GGA functionals are
    !% requested, through the <tt>XCFunctional</tt> variable, the default
    !% is <tt>generalized_kohn_sham</tt>.
    !%Option independent_particles 2
    !% Particles will be considered as independent, <i>i.e.</i> as non-interacting.
    !% This mode is mainly used for testing purposes, as the code is usually
    !% much faster with <tt>independent_particles</tt>.
    !%Option hartree 1
    !% Calculation within the Hartree method (experimental). Note that, contrary to popular
    !% belief, the Hartree potential is self-interaction-free. Therefore, this run
    !% mode will not yield the same result as <tt>kohn-sham</tt> without exchange-correlation.
    !%Option hartree_fock 3
    !% This is the traditional Hartree-Fock scheme. Like the Hartree scheme, it is fully
    !% self-interaction-free.
    !%Option kohn_sham 4
    !% This is the default density-functional theory scheme. Note that you can also use
    !% hybrid functionals in this scheme, but they will be handled the "DFT" way, <i>i.e.</i>,
    !% solving the OEP equation. DFT+U (see <tt>DFTULevel</tt>) also runs within this
    !% scheme; it does not require <tt>generalized_kohn_sham</tt>.
    !%Option generalized_kohn_sham 5
    !% This is similar to the <tt>kohn-sham</tt> scheme, except that this allows for nonlocal operators.
    !% This is the default mode to run hybrid functionals or meta-GGA functionals.
    !% It can be more convenient to use <tt>kohn-sham</tt> DFT within the OEP scheme to get similar (but not the same) results.
    !% Note that within this scheme you can use a correlation functional, or a hybrid
    !% functional (see <tt>XCFunctional</tt>). In the latter case, you will be following the
    !% quantum-chemistry recipe to use hybrids.
    !% DFT+U (see <tt>DFTULevel</tt>) can be combined with any theory level and yields the
    !% same result under <tt>kohn_sham</tt> and <tt>generalized_kohn_sham</tt>; it only
    !% requires <tt>generalized_kohn_sham</tt> when combined with a hybrid functional, so
    !% that the exact-exchange part is treated as a nonlocal operator.
    !%Option rdmft 7
    !% (Experimental) Reduced Density Matrix functional theory.
    !%End
 
    ks%xc_family = xc_family_none
    ks%sic%level = sic_none
    ks%oep%level = oep_level_none
 
    ks%theory_level = kohn_sham_dft
    parsed_theory_level = .false.
 
    ! the user knows what he wants, give her that
    if (parse_is_defined(namespace, 'TheoryLevel')) then
      call parse_variable(namespace, 'TheoryLevel', kohn_sham_dft, ks%theory_level)
      if (.not. varinfo_valid_option('TheoryLevel', ks%theory_level)) call messages_input_error(namespace, 'TheoryLevel')
 
      parsed_theory_level = .true.
    end if
 
    ! parse the XC functional
 
    call get_functional_from_pseudos(pseudo_x_functional, pseudo_c_functional)
 
    default = 0
    if (ks%theory_level == kohn_sham_dft .or. ks%theory_level == generalized_kohn_sham_dft) then
      default = xc_get_default_functional(space%dim, pseudo_x_functional, pseudo_c_functional)
    end if
 
    if (.not. parse_is_defined(namespace, 'XCFunctional') &
      .and. (pseudo_x_functional /= pseudo_exchange_any .or. pseudo_c_functional /= pseudo_correlation_any)) then
      call messages_write('Info: the XCFunctional has been selected to match the pseudopotentials', new_line = .true.)
      call messages_write('      used in the calculation.')
      call messages_info(namespace=namespace)
    end if
 
    ! The description of this variable can be found in file src/xc/functionals_list.F90
    call parse_variable(namespace, 'XCFunctional', default, val)
 
    ! the first 3 digits of the number indicate the X functional and
    ! the next 3 the C functional.
    c_id = val / libxc_c_index
    x_id = val - c_id * libxc_c_index
 
    if ((x_id /= pseudo_x_functional .and. pseudo_x_functional /= pseudo_exchange_any) .or. &
      (c_id /= pseudo_c_functional .and. pseudo_c_functional /= pseudo_correlation_any)) then
      call messages_write('The XCFunctional that you selected does not match the one used', new_line = .true.)
      call messages_write('to generate the pseudopotentials.')
      call messages_warning(namespace=namespace)
    end if
 
    ! FIXME: we rarely need this. We should only parse when necessary.
 
    !%Variable XCKernel
    !%Type integer
    !%Default -1
    !%Section Hamiltonian::XC
    !%Description
    !% Defines the exchange-correlation kernel. LDA and GGA kernels are available;
    !% meta-GGAs and hybrids are not.
    !% The options are the same as <tt>XCFunctional</tt>.
    !% Note: the kernel is only needed for Casida, Sternheimer, or optimal-control calculations.
    !% Note: the second-order kernel needed for hyperpolarizabilities (<tt>EMHyperpol</tt>)
    !% is still restricted to LDA.
    !%Option xc_functional -1
    !% The same functional defined by <tt>XCFunctional</tt>. By default, this is the case.
    !%End
    call parse_variable(namespace, 'XCKernel', -1, val)
    if (-1 == val) then
      ck_id = c_id
      xk_id = x_id
    else
      ck_id = val / libxc_c_index
      xk_id = val - ck_id * libxc_c_index
    end if
 
    call messages_obsolete_variable(namespace, 'XFunctional', 'XCFunctional')
    call messages_obsolete_variable(namespace, 'CFunctional', 'XCFunctional')
 
    call ks%v_ks_photons%init(namespace)
 
    ! initialize XC modules
 
    ! This is a bit ugly, theory_level might not be generalized KS or HF now
    ! but it might become generalized KS or HF later. This is safe because it
    ! becomes generalized KS in the cases where the functional is hybrid
    ! and the ifs inside check for both conditions.
    using_hartree_fock = (ks%theory_level == hartree_fock) &
      .or. (ks%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(ks%xc))
    call xc_init(ks%xc, namespace, space%dim, space%periodic_dim, st%qtot, &
      x_id, c_id, xk_id, ck_id, hartree_fock = using_hartree_fock, ispin=st%d%ispin)
 
    ks%xc_family = ks%xc%family
    ks%xc_flags  = ks%xc%flags
 
    if (.not. parsed_theory_level) then
      default = kohn_sham_dft
 
      ! the functional is a hybrid, use Hartree-Fock as theory level by default
      if (family_is_hybrid(ks%xc) .or. family_is_mgga_with_exc(ks%xc)) then
        default = generalized_kohn_sham_dft
      end if
 
      ! In principle we do not need to parse. However we do it for consistency
      call parse_variable(namespace, 'TheoryLevel', default, ks%theory_level)
      if (.not. varinfo_valid_option('TheoryLevel', ks%theory_level)) call messages_input_error(namespace, 'TheoryLevel')
 
    end if
 
    ! In case we need OEP, we need to find which type of OEP it is
    oep_type = -1
    if (family_is_mgga_with_exc(ks%xc)) then
      call messages_experimental('MGGA energy functionals')
 
      if (ks%theory_level == kohn_sham_dft) then
        call messages_experimental("MGGA within the Kohn-Sham scheme")
        ks%xc_family = ior(ks%xc_family, xc_family_oep)
        oep_type = oep_type_mgga
      end if
    end if
 
    call messages_obsolete_variable(namespace, 'NonInteractingElectrons', 'TheoryLevel')
    call messages_obsolete_variable(namespace, 'HartreeFock', 'TheoryLevel')
 
    ! Due to how the code is made, we need to set this to have theory level other than DFT
    ! correct...
    ks%sic%amaldi_factor = m_one
 
    select case (ks%theory_level)
    case (independent_particles)
 
    case (hartree)
      call messages_experimental("Hartree theory level")
      if (space%periodic_dim == space%dim) then
        call messages_experimental("Hartree in fully periodic system")
      end if
      if (kpoints%full%npoints > 1) then
        call messages_not_implemented("Hartree with k-points", namespace=namespace)
      end if
 
    case (hartree_fock)
      if (kpoints%full%npoints > 1) then
        call messages_experimental("Hartree-Fock with k-points")
      end if
 
    case (generalized_kohn_sham_dft)
      if (kpoints%full%npoints > 1 .and. family_is_hybrid(ks%xc)) then
        call messages_experimental("Hybrid functionals with k-points")
      end if
 
    case (rdmft)
      call messages_experimental('RDMFT theory level')
 
    case (kohn_sham_dft)
 
      ! check for SIC
      if (bitand(ks%xc_family, xc_family_lda + xc_family_gga) /= 0) then
        call xc_sic_init(ks%sic, namespace, gr, st, mc, space)
      end if
 
      if (bitand(ks%xc_family, xc_family_oep) /= 0) then
        select case (ks%xc%functional(func_x,1)%id)
        case (xc_oep_x_slater)
          if (kpoints%reduced%npoints > 1 .and. st%d%ispin == spinors) then
            call messages_not_implemented("Slater with k-points and spinor wavefunctions", namespace=namespace)
          end if
          if (kpoints%use_symmetries) then
            call messages_not_implemented("Slater with k-points symmetries", namespace=namespace)
          end if
          ks%oep%level = oep_level_none
        case (xc_oep_x_fbe)
          if (kpoints%reduced%npoints > 1) then
            call messages_not_implemented("FBE functional with k-points", namespace=namespace)
          end if
          ks%oep%level = oep_level_none
        case default
          if((.not. ks%v_ks_photons%active()) .or. (ks%v_ks_photons%functional() /= 0)) then
            if(oep_type == -1) then ! Else we have a MGGA
              oep_type = oep_type_exx
            end if
            call xc_oep_init(ks%oep, namespace, gr, st, mc, space, oep_type)
          end if
        end select
      else
        ks%oep%level = oep_level_none
      end if
 
      if (bitand(ks%xc_family, xc_family_ks_inversion) /= 0) then
        call xc_ks_inversion_init(ks%ks_inversion, namespace, gr, ions, st, ks%xc, mc, space, kpoints)
      end if
 
    end select
 
    if (ks%theory_level /= kohn_sham_dft .and. parse_is_defined(namespace, "SICCorrection")) then
      message(1) = "SICCorrection can only be used with Kohn-Sham DFT"
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (st%d%ispin == spinors) then
      if (bitand(ks%xc_family, xc_family_mgga + xc_family_hyb_mgga) /= 0) then
        call messages_not_implemented("MGGA with spinors", namespace=namespace)
      end if
    end if
 
    ks%frozen_hxc = .false.
 
    call v_ks_write_info(ks, namespace=namespace)
 
    ks%gr => gr
    ks%calc%calculating = .false.
 
    !The value of ks%calculate_current is set to false or true by Output
    call current_init(ks%current_calculator, namespace)
 
    call ks%vdw%init(namespace, space, gr, ks%xc, ions, x_id, c_id)
    if (ks%vdw%vdw_correction /= option__vdwcorrection__none .and. ks%theory_level == rdmft) then
      message(1) = "VDWCorrection and RDMFT are not compatible"
      call messages_fatal(1, namespace=namespace)
    end if
    if (ks%vdw%vdw_correction /= option__vdwcorrection__none .and. ks%theory_level == independent_particles) then
      message(1) = "VDWCorrection and independent particles are not compatible"
      call messages_fatal(1, namespace=namespace)
    end if
 
    call ks%v_ks_photons%init_xc(namespace, space, gr, st)
 
    call sturm_liouville_init(ks%sl_solver, namespace, gr, space)
 
    ks%stress_xc = m_zero
 
    pop_sub(v_ks_init)
 
  contains
 
    subroutine get_functional_from_pseudos(x_functional, c_functional)
      integer, intent(out) :: x_functional
      integer, intent(out) :: c_functional
 
      integer :: xf, cf, ispecies
      logical :: warned_inconsistent
 
      x_functional = pseudo_exchange_any
      c_functional = pseudo_correlation_any
 
      warned_inconsistent = .false.
      do ispecies = 1, ions%nspecies
        select type(spec=>ions%species(ispecies)%s)
        class is(pseudopotential_t)
          xf = spec%x_functional()
          cf = spec%c_functional()
 
          if (xf == pseudo_exchange_unknown .or. cf == pseudo_correlation_unknown) then
            call messages_write("Unknown XC functional for species '"//trim(ions%species(ispecies)%s%get_label())//"'")
            call messages_warning(namespace=namespace)
            cycle
          end if
 
          if (x_functional == pseudo_exchange_any) then
            x_functional = xf
          else
            if (xf /= x_functional .and. .not. warned_inconsistent) then
              call messages_write('Inconsistent XC functional detected between species')
              call messages_warning(namespace=namespace)
              warned_inconsistent = .true.
            end if
          end if
 
          if (c_functional == pseudo_correlation_any) then
            c_functional = cf
          else
            if (cf /= c_functional .and. .not. warned_inconsistent) then
              call messages_write('Inconsistent XC functional detected between species')
              call messages_warning(namespace=namespace)
              warned_inconsistent = .true.
            end if
          end if
 
        class default
          xf = pseudo_exchange_unknown
          cf = pseudo_correlation_unknown
        end select
 
      end do
 
      assert(x_functional /= pseudo_exchange_unknown)
      assert(c_functional /= pseudo_correlation_unknown)
 
    end subroutine get_functional_from_pseudos
  end subroutine v_ks_init
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine v_ks_end(ks)
    type(v_ks_t),     intent(inout) :: ks
 
    push_sub(v_ks_end)
 
    call ks%vdw%end()
    call sturm_liouville_end(ks%sl_solver)
 
    select case (ks%theory_level)
    case (kohn_sham_dft)
      if (bitand(ks%xc_family, xc_family_ks_inversion) /= 0) then
        call xc_ks_inversion_end(ks%ks_inversion)
      end if
      if (bitand(ks%xc_family, xc_family_oep) /= 0) then
        call xc_oep_end(ks%oep)
      end if
      call xc_end(ks%xc)
    case (hartree_fock, generalized_kohn_sham_dft)
      call xc_end(ks%xc)
    end select
 
    call xc_sic_end(ks%sic)
 
    call ks%v_ks_photons%end()
 
    pop_sub(v_ks_end)
  end subroutine v_ks_end
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine v_ks_write_info(ks, iunit, namespace)
    type(v_ks_t),                intent(in) :: ks
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    push_sub(v_ks_write_info)
 
    call messages_print_with_emphasis(msg="Theory Level", iunit=iunit, namespace=namespace)
    call messages_print_var_option("TheoryLevel", ks%theory_level, iunit=iunit, namespace=namespace)
 
    select case (ks%theory_level)
    case (hartree_fock, generalized_kohn_sham_dft)
      call messages_info(iunit=iunit, namespace=namespace)
      call xc_write_info(ks%xc, iunit, namespace)
 
    case (kohn_sham_dft)
      call messages_info(iunit=iunit, namespace=namespace)
      call xc_write_info(ks%xc, iunit, namespace)
 
      call messages_info(iunit=iunit, namespace=namespace)
 
      call xc_sic_write_info(ks%sic, iunit, namespace)
      call xc_oep_write_info(ks%oep, iunit, namespace)
      call xc_ks_inversion_write_info(ks%ks_inversion, iunit, namespace)
 
    end select
 
    call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
 
    pop_sub(v_ks_write_info)
  end subroutine v_ks_write_info
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine v_ks_check_generalized_overlap(ks, st, hm, namespace)
    type(v_ks_t),             intent(in) :: ks
    type(states_elec_t),      intent(in) :: st
    type(hamiltonian_elec_t), intent(in) :: hm
    type(namespace_t),        intent(in) :: namespace
 
    push_sub(v_ks_check_generalized_overlap)
 
    if (st%requires_generalized_overlap) then
      if (st%d%ispin == spin_polarized) then
        call messages_not_implemented(&
          "Spin-polarized calculations with ultrasoft/PAW (the screened D_ij is not spin resolved)", &
          namespace=namespace)
      end if
      if (st%d%ispin == spinors) then
        call messages_not_implemented(&
          "Noncollinear (spinor) calculations with ultrasoft/PAW", namespace=namespace)
      end if
      if (family_is_mgga(ks%xc%family)) then
        call messages_not_implemented(&
          "meta-GGA with ultrasoft/PAW (the kinetic-energy-density augmentation is missing)", &
          namespace=namespace)
      end if
      if (ks%theory_level == hartree_fock .or. family_is_hybrid(ks%xc)) then
        call messages_not_implemented(&
          "Hartree-Fock / hybrid functionals with ultrasoft/PAW (exact-exchange augmentation is missing)", &
          namespace=namespace)
      end if
      if (ks%sic%level == sic_pz_oep .or. ks%oep%level /= oep_level_none) then
        call messages_not_implemented(&
          "PZ-SIC / OEP with ultrasoft/PAW (orbital-resolved augmentation is missing)", &
          namespace=namespace)
      end if
      if (hm%lda_u_level /= dft_u_none) then
        call messages_not_implemented(&
          "DFT+U with ultrasoft/PAW (the occupation matrix needs the S-augmented projection)", &
          namespace=namespace)
      end if
    end if
 
    pop_sub(v_ks_check_generalized_overlap)
  end subroutine v_ks_check_generalized_overlap
  ! ---------------------------------------------------------
 
 
  !----------------------------------------------------------
  subroutine v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, calc_eigenval, calc_current)
    type(namespace_t),        intent(in)    :: namespace
    type(electron_space_t),   intent(in)    :: space
    type(grid_t),             intent(in)    :: 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
    logical,        optional, intent(in)    :: calc_eigenval
    logical,        optional, intent(in)    :: calc_current
 
    integer, allocatable :: ind(:)
    integer :: ist, ik
    real(real64), allocatable :: copy_occ(:)
    logical :: calc_eigenval_
    logical :: calc_current_
 
    push_sub(v_ks_h_setup)
 
    calc_eigenval_ = optional_default(calc_eigenval, .true.)
    calc_current_ = optional_default(calc_current, .true.)
    call states_elec_fermi(st, namespace, gr)
    call density_calc(st, gr, st%rho, ep = hm%ep, vnl = hm%uspp_overlap, phase = hm%phase)
    call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
      calc_eigenval = calc_eigenval_, calc_current = calc_current_) ! get potentials
 
    if (st%restart_reorder_occs .and. .not. st%fromScratch) then
      message(1) = "Reordering occupations for restart."
      call messages_info(1, namespace=namespace)
 
      safe_allocate(ind(1:st%nst))
      safe_allocate(copy_occ(1:st%nst))
 
      do ik = 1, st%nik
        call sort(st%eigenval(:, ik), ind)
        copy_occ(1:st%nst) = st%occ(1:st%nst, ik)
        do ist = 1, st%nst
          st%occ(ist, ik) = copy_occ(ind(ist))
        end do
      end do
 
      safe_deallocate_a(ind)
      safe_deallocate_a(copy_occ)
    end if
 
    if (calc_eigenval_) call states_elec_fermi(st, namespace, gr) ! occupations
    call energy_calc_total(namespace, space, hm, gr, st, ext_partners)
 
    pop_sub(v_ks_h_setup)
  end subroutine v_ks_h_setup
 
  ! ---------------------------------------------------------
  subroutine v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
    calc_eigenval, time, calc_energy, calc_current, force_semilocal)
    type(v_ks_t),               intent(inout) :: ks
    type(namespace_t),          intent(in)    :: namespace
    type(electron_space_t),     intent(in)    :: space
    type(hamiltonian_elec_t),   intent(inout) :: hm
    type(states_elec_t),        intent(inout) :: st
    type(ions_t),               intent(in)    :: ions
    type(partner_list_t),       intent(in)    :: ext_partners
    logical,          optional, intent(in)    :: calc_eigenval
    real(real64),     optional, intent(in)    :: time
    logical,          optional, intent(in)    :: calc_energy
    logical,          optional, intent(in)    :: calc_current
    logical,          optional, intent(in)    :: force_semilocal
 
    logical :: calc_current_
 
    push_sub(v_ks_calc)
 
    calc_current_ = optional_default(calc_current, .true.)  &
      .and. (ks%calculate_current &
      .and. states_are_complex(st) &
      .or. hamiltonian_elec_needs_current(hm, states_are_real(st)))
 
    if (calc_current_) then
      call states_elec_allocate_current(st, space, ks%gr)
      call current_calculate(ks%current_calculator, namespace, ks%gr, hm, space, st)
    end if
 
    call v_ks_calc_start(ks, namespace, space, hm, st, ions, hm%kpoints%latt, ext_partners, time, &
      calc_energy, force_semilocal=force_semilocal)
    call v_ks_calc_finish(ks, hm, namespace, space, hm%kpoints%latt, st, &
      ext_partners, force_semilocal=force_semilocal)
 
    if (optional_default(calc_eigenval, .false.)) then
      call energy_calc_eigenvalues(namespace, hm, ks%gr%der, st)
    end if
 
    ! Update the magnetic constrain
    call magnetic_constrain_update(hm%magnetic_constrain, ks%gr, st%d, space, hm%kpoints%latt, ions%pos, st%rho)
    ! We add the potential to vxc, as this way the potential gets mixed together with vxc
    ! While this is not ideal, this is a simple practical solution
    if (hm%magnetic_constrain%level /= constrain_none) then
      call lalg_axpy(ks%gr%np, st%d%nspin, m_one, hm%magnetic_constrain%pot, hm%ks_pot%vhxc)
    end if
 
    pop_sub(v_ks_calc)
  end subroutine v_ks_calc
 
  ! ---------------------------------------------------------
 
  subroutine v_ks_calc_start(ks, namespace, space, hm, st, ions, latt, ext_partners, time, &
    calc_energy, force_semilocal)
    type(v_ks_t),              target, intent(inout) :: ks
    type(namespace_t),                 intent(in)    :: namespace
    class(space_t),                    intent(in)    :: space
    type(hamiltonian_elec_t),  target, intent(in)    :: hm
    type(states_elec_t),       target, intent(inout) :: st
    type(ions_t),                      intent(in)    :: ions
    type(lattice_vectors_t),           intent(in)    :: latt
    type(partner_list_t),              intent(in)    :: ext_partners
    real(real64),            optional, intent(in)    :: time
    logical,                 optional, intent(in)    :: calc_energy
    logical,                 optional, intent(in)    :: force_semilocal
 
    push_sub(v_ks_calc_start)
 
    call profiling_in("KOHN_SHAM_CALC")
 
    assert(.not. ks%calc%calculating)
    ks%calc%calculating = .true.
 
    write(message(1), '(a)') 'Debug: Calculating Kohn-Sham potential.'
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    ks%calc%time_present = present(time)
    ks%calc%time = optional_default(time, m_zero)
 
    ks%calc%calc_energy = optional_default(calc_energy, .true.)
 
    ! If the Hxc term is frozen, there is nothing more to do (WARNING: MISSING ks%calc%energy%intnvxc)
    if (ks%frozen_hxc) then
      call profiling_out("KOHN_SHAM_CALC")
      pop_sub(v_ks_calc_start)
      return
    end if
 
    allocate(ks%calc%energy)
 
    call energy_copy(hm%energy, ks%calc%energy)
 
    ks%calc%energy%intnvxc = m_zero
 
    nullify(ks%calc%total_density)
 
    if (ks%theory_level /= independent_particles .and. abs(ks%sic%amaldi_factor) > m_epsilon) then
 
      call calculate_density()
 
      if (poisson_is_async(hm%psolver)) then
        call dpoisson_solve_start(hm%psolver, ks%calc%total_density)
      end if
 
      if (ks%theory_level /= hartree .and. ks%theory_level /= rdmft) call v_a_xc(hm, force_semilocal)
    else
      ks%calc%total_density_alloc = .false.
    end if
 
    ! The exchange operator is computed from the states of the previous iteration
    ! This is done by copying the state object to ks%calc%hf_st
    ! For ACE, the states are the same in ks%calc%hf_st and st, as we compute the
    ! ACE potential in v_ks_finish, so the copy is not needed
    nullify(ks%calc%hf_st)
    if (ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
      .or. ks%theory_level == rdmft .or. (ks%theory_level == generalized_kohn_sham_dft &
      .and. family_is_hybrid(ks%xc))) then
 
      if (st%parallel_in_states) then
        if (accel_is_enabled()) then
          call messages_write('State parallelization of Hartree-Fock exchange is not supported')
          call messages_new_line()
          call messages_write('when running with GPUs. Please use domain parallelization')
          call messages_new_line()
          call messages_write("or disable acceleration using 'DisableAccel = yes'.")
          call messages_fatal(namespace=namespace)
        end if
      end if
 
      if (hm%exxop%useACE) then
        ks%calc%hf_st => st
      else
        safe_allocate(ks%calc%hf_st)
        call states_elec_copy(ks%calc%hf_st, st)
      end if
    end if
 
    ! Calculate the vector potential induced by the electronic current.
    ! WARNING: calculating the self-induced magnetic field here only makes
    ! sense if it is going to be used in the Hamiltonian, which does not happen
    ! now. Otherwise one could just calculate it at the end of the calculation.
    if (hm%self_induced_magnetic) then
      safe_allocate(ks%calc%a_ind(1:ks%gr%np_part, 1:space%dim))
      safe_allocate(ks%calc%b_ind(1:ks%gr%np_part, 1:space%dim))
      call magnetic_induced(namespace, ks%gr, st, hm%psolver, hm%kpoints, ks%calc%a_ind, ks%calc%b_ind)
    end if
 
    if ((ks%v_ks_photons%active()) .and. (ks%calc%time_present) .and. (ks%v_ks_photons%functional() == 0) ) then
      call ks%v_ks_photons%mf_calc(ks%gr, st, ions, time)
    end if
 
    ! if (ks%has_vibrations) then
    !   call vibrations_eph_coup(ks%vib, ks%gr, hm, ions, st)
    ! end if
 
    call profiling_out("KOHN_SHAM_CALC")
    pop_sub(v_ks_calc_start)
 
  contains
 
    subroutine calculate_density()
      integer :: ip
 
      push_sub(v_ks_calc_start.calculate_density)
 
      ! get density taking into account non-linear core corrections
      safe_allocate(ks%calc%density(1:ks%gr%np, 1:st%d%nspin))
      call states_elec_total_density(st, ks%gr, ks%calc%density)
 
      ! Amaldi correction on CPU
      if (ks%sic%level == sic_amaldi) then
        call lalg_scal(ks%gr%np, st%d%nspin, ks%sic%amaldi_factor, ks%calc%density)
      end if
 
      ! GPU counterpart of the CPU corrections above: the CPU path includes rho_core, frozen_rho,
      ! and Amaldi scaling directly into ks%calc%density before calling xc_get_vxc.
      ! On the GPU we cannot do the same to st%buff_density because:
      ! - it is in wavefunction storage layout (pnp, nspin) while libxc expects (spin_channels, np)
      ! - permanently modifying st%buff_density would corrupt the density for subsequent SCF steps.
      ! Instead we upload the correction arrays here so that xc_update_internal_quantities can apply
      ! them to dens_buff (already in libxc layout) via the xc_dens_apply_corrections kernel,
      ! immediately after the xc_dens_extract_block kernel runs.
      if (.not. ks%xc%xc_on_host .and. accel_buffer_is_allocated(st%buff_density)) then
        if (allocated(st%rho_core)) then
          call accel_create_buffer(ks%xc%quantities%buff_rho_core, &
            accel_mem_read_only, type_float, int(ks%gr%np, int64))
          call accel_write_buffer(ks%xc%quantities%buff_rho_core, &
            int(ks%gr%np, int64), st%rho_core)
        end if
        if (allocated(st%frozen_rho)) then
          call accel_create_buffer(ks%xc%quantities%buff_frozen_rho, &
            accel_mem_read_only, type_float, int(ks%gr%np, int64)*int(st%d%nspin, int64))
          call accel_write_buffer(ks%xc%quantities%buff_frozen_rho, &
            int(ks%gr%np, int64), int(st%d%nspin, int64), st%frozen_rho)
          ks%xc%quantities%frozen_rho_np = ks%gr%np
        end if
        if (ks%sic%level == sic_amaldi) then
          ks%xc%quantities%amaldi_factor = ks%sic%amaldi_factor
        end if
      end if
 
      nullify(ks%calc%total_density)
      if (allocated(st%rho_core) .or. hm%d%spin_channels > 1) then
        ks%calc%total_density_alloc = .true.
 
        safe_allocate(ks%calc%total_density(1:ks%gr%np))
 
        do ip = 1, ks%gr%np
          ks%calc%total_density(ip) = sum(ks%calc%density(ip, 1:hm%d%spin_channels))
        end do
 
        ! remove non-local core corrections
        if (allocated(st%rho_core)) then
          call lalg_axpy(ks%gr%np, -ks%sic%amaldi_factor, st%rho_core,  ks%calc%total_density)
        end if
      else
        ks%calc%total_density_alloc = .false.
        ks%calc%total_density => ks%calc%density(:, 1)
      end if
 
      pop_sub(v_ks_calc_start.calculate_density)
    end subroutine calculate_density
 
    ! ---------------------------------------------------------
    subroutine v_a_xc(hm, force_semilocal)
      type(hamiltonian_elec_t),  intent(in) :: hm
      logical, optional,         intent(in) :: force_semilocal
 
      push_sub(v_ks_calc_start.v_a_xc)
      call profiling_in("XC")
 
      ks%calc%energy%exchange = m_zero
      ks%calc%energy%correlation = m_zero
      ks%calc%energy%xc_j = m_zero
      ks%calc%energy%vdw = m_zero
 
      allocate(ks%calc%vxc(1:ks%gr%np, 1:st%d%nspin))
      ks%calc%vxc = m_zero
 
      if (family_is_mgga_with_exc(hm%xc)) then
        safe_allocate(ks%calc%vtau(1:ks%gr%np, 1:st%d%nspin))
        ks%calc%vtau = m_zero
      end if
 
      ! Get the *local* XC term
      if (ks%calc%calc_energy) then
        if (family_is_mgga_with_exc(hm%xc)) then
          call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, st%d%ispin, &
            latt%rcell_volume, ks%calc%vxc, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation, &
            deltaxc = ks%calc%energy%delta_xc, vtau = ks%calc%vtau, stress_xc=ks%stress_xc, &
            force_orbitalfree=force_semilocal)
        else
          call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, st%d%ispin, &
            latt%rcell_volume, ks%calc%vxc, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation, &
            deltaxc = ks%calc%energy%delta_xc, stress_xc=ks%stress_xc, force_orbitalfree=force_semilocal)
        end if
      else
        if (family_is_mgga_with_exc(hm%xc)) then
          call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, &
            st%d%ispin, latt%rcell_volume, ks%calc%vxc, vtau = ks%calc%vtau, stress_xc=ks%stress_xc, &
            force_orbitalfree=force_semilocal)
        else
          call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, &
            st%d%ispin, latt%rcell_volume, ks%calc%vxc, stress_xc=ks%stress_xc, force_orbitalfree=force_semilocal)
        end if
      end if
 
      !Noncollinear functionals
      if (bitand(hm%xc%family, xc_family_nc_lda + xc_family_nc_mgga) /= 0) then
        if (st%d%ispin /= spinors) then
          message(1) = "Noncollinear functionals can only be used with spinor wavefunctions."
          call messages_fatal(1)
        end if
 
        if (optional_default(force_semilocal, .false.)) then
          message(1) = "Cannot perform LCAO for noncollinear MGGAs."
          message(2) = "Please perform a LDA calculation first."
          call messages_fatal(2)
        end if
 
        if (ks%calc%calc_energy) then
          if (family_is_mgga_with_exc(hm%xc)) then
            call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, ks%calc%vxc, &
              vtau = ks%calc%vtau, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation)
          else
            call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, ks%calc%vxc, &
              ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation)
          end if
        else
          if (family_is_mgga_with_exc(hm%xc)) then
            call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, &
              ks%calc%vxc, vtau = ks%calc%vtau)
          else
            call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, ks%calc%vxc)
          end if
        end if
      end if
 
      call ks%vdw%calc(namespace, space, latt, ions%atom, ions%natoms, ions%pos, &
        ks%gr, st, ks%calc%energy%vdw, ks%calc%vxc)
 
      if (optional_default(force_semilocal, .false.)) then
        call profiling_out("XC")
        pop_sub(v_ks_calc_start.v_a_xc)
        return
      end if
 
      ! ADSIC correction
      if (ks%sic%level == sic_adsic) then
        if (family_is_mgga(hm%xc%family)) then
          call messages_not_implemented('ADSIC with MGGAs', namespace=namespace)
        end if
        if (ks%calc%calc_energy) then
          call xc_sic_calc_adsic(ks%sic, namespace, space, ks%gr, st, hm, ks%xc, ks%calc%density, &
            ks%calc%vxc, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation)
        else
          call xc_sic_calc_adsic(ks%sic, namespace, space, ks%gr, st, hm, ks%xc, ks%calc%density, &
            ks%calc%vxc)
        end if
      end if
      !PZ SIC is done in the finish routine as OEP full needs to update the Hamiltonian
 
      if (ks%theory_level == kohn_sham_dft) then
        ! The OEP family has to be handled specially
        ! Note that OEP is done in the finish state, as it requires updating the Hamiltonian and needs the new Hartre and vxc term
        if (bitand(ks%xc_family, xc_family_oep) /= 0 .or. family_is_mgga_with_exc(ks%xc)) then
 
          if (ks%xc%functional(func_x,1)%id == xc_oep_x_slater) then
            call x_slater_calc(namespace, ks%gr, space, hm%exxop, st, hm%kpoints, ks%calc%energy%exchange, &
              vxc = ks%calc%vxc)
          else if (ks%xc%functional(func_x,1)%id == xc_oep_x_fbe .or. ks%xc%functional(func_x,1)%id == xc_oep_x_fbe_sl) then
            call x_fbe_calc(ks%xc%functional(func_x,1)%id, namespace, hm%psolver, ks%sl_solver, ks%gr, st, space, &
              ks%calc%energy%exchange, vxc = ks%calc%vxc)
 
          else if (ks%xc%functional(func_c,1)%id == xc_lda_c_fbe_sl) then
 
            call fbe_c_lda_sl(namespace, hm%psolver, ks%sl_solver, ks%gr, st, space, ks%calc%energy%correlation, vxc = ks%calc%vxc)
 
          end if
 
        end if
 
        if (bitand(ks%xc_family, xc_family_ks_inversion) /= 0) then
          ! Also treat KS inversion separately (not part of libxc)
          call xc_ks_inversion_calc(ks%ks_inversion, namespace, space, ks%gr, hm, ext_partners, st, vxc = ks%calc%vxc, &
            time = ks%calc%time)
        end if
 
        ! compute the photon-free photon exchange potential and energy
        if (ks%v_ks_photons%functional() /= 0) then
          call ks%v_ks_photons%add_px(namespace, ks%calc%total_density, ks%gr, space, hm%psolver, st, &
            hm%d%spin_channels, ks%calc%vxc, ks%calc%energy%photon_exchange)
        end if
 
      end if
 
      if (ks%calc%calc_energy) then
        ! MGGA vtau contribution is done after copying vtau to hm%vtau
 
        call v_ks_update_dftu_energy(ks, namespace, hm, st, ks%calc%energy%int_dft_u)
      end if
 
      call profiling_out("XC")
      pop_sub(v_ks_calc_start.v_a_xc)
    end subroutine v_a_xc
 
  end subroutine v_ks_calc_start
  ! ---------------------------------------------------------
 
  subroutine v_ks_calc_finish(ks, hm, namespace, space, latt, st, ext_partners, force_semilocal)
    type(v_ks_t),     target, intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(lattice_vectors_t),  intent(in)    :: latt
    type(states_elec_t),      intent(inout) :: st
    type(partner_list_t),     intent(in)    :: ext_partners
    logical,       optional,  intent(in)    :: force_semilocal
 
    integer                           :: ip, ispin
    type(states_elec_t) :: xst
    real(real64) :: exx_energy
    real(real64) :: factor
 
    push_sub(v_ks_calc_finish)
 
    assert(ks%calc%calculating)
    ks%calc%calculating = .false.
 
    if (ks%frozen_hxc) then
      pop_sub(v_ks_calc_finish)
      return
    end if
 
    !change the pointer to the energy object
    safe_deallocate_a(hm%energy)
    call move_alloc(ks%calc%energy, hm%energy)
 
    if (hm%self_induced_magnetic) then
      hm%a_ind(1:ks%gr%np, 1:space%dim) = ks%calc%a_ind(1:ks%gr%np, 1:space%dim)
      hm%b_ind(1:ks%gr%np, 1:space%dim) = ks%calc%b_ind(1:ks%gr%np, 1:space%dim)
 
      safe_deallocate_a(ks%calc%a_ind)
      safe_deallocate_a(ks%calc%b_ind)
    end if
 
    if (allocated(hm%v_static)) then
      if (associated(ks%calc%total_density)) then
        hm%energy%intnvstatic = dmf_dotp(ks%gr, ks%calc%total_density, hm%v_static)
      else
        ! The total density is not computed for independent particles
        hm%energy%intnvstatic = m_zero
        do ispin = 1, hm%d%spin_channels
          hm%energy%intnvstatic = hm%energy%intnvstatic + dmf_dotp(ks%gr, st%rho(:, ispin), hm%v_static)
        end do
      end if
    else
      hm%energy%intnvstatic = m_zero
    end if
 
    if (ks%theory_level == independent_particles .or. abs(ks%sic%amaldi_factor) <= m_epsilon) then
 
      hm%ks_pot%vhxc = m_zero
      hm%energy%intnvxc     = m_zero
      hm%energy%hartree     = m_zero
      hm%energy%exchange    = m_zero
      hm%energy%exchange_hf = m_zero
      hm%energy%correlation = m_zero
    else
 
      hm%energy%hartree = m_zero
      call v_ks_hartree(namespace, ks, space, hm, ext_partners)
 
      if (.not. optional_default(force_semilocal, .false.)) then
        !PZ-SIC
        if(ks%sic%level == sic_pz_oep) then
          if (states_are_real(st)) then
            call dxc_oep_calc(ks%sic%oep, namespace, ks%xc, ks%gr, hm, st, space, &
              latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
          else
            call zxc_oep_calc(ks%sic%oep, namespace, ks%xc, ks%gr, hm, st, space, &
              latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
          end if
        end if
 
        ! OEP for exchange ad MGGAs (within Kohn-Sham DFT)
        if (ks%theory_level == kohn_sham_dft .and. ks%oep%level /= oep_level_none) then
          ! The OEP family has to be handled specially
          if (ks%xc%functional(func_x,1)%id == xc_oep_x .or. family_is_mgga_with_exc(ks%xc)) then
            if (states_are_real(st)) then
              call dxc_oep_calc(ks%oep, namespace, ks%xc, ks%gr, hm, st, space, &
                latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
            else
              call zxc_oep_calc(ks%oep, namespace, ks%xc, ks%gr, hm, st, space, &
                latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
            end if
          end if
        end if
      end if
 
      if (ks%theory_level == kohn_sham_dft) then
        call ks%v_ks_photons%oep_calc(namespace, ks%xc, ks%gr, hm, st, space, ks%calc%vxc)
      end if
 
 
      if (ks%calc%calc_energy) then
        ! Now we calculate Int[n vxc] = energy%intnvxc
        hm%energy%intnvxc = m_zero
 
        if (ks%theory_level /= independent_particles .and. ks%theory_level /= hartree .and. ks%theory_level /= rdmft) then
          do ispin = 1, hm%d%nspin
            if (ispin <= 2) then
              factor = m_one
            else
              factor = m_two
            end if
            hm%energy%intnvxc = hm%energy%intnvxc + &
              factor*dmf_dotp(ks%gr, st%rho(:, ispin), ks%calc%vxc(:, ispin), reduce = .false.)
          end do
          call ks%gr%allreduce(hm%energy%intnvxc)
        end if
      end if
 
 
      if (ks%theory_level /= hartree .and. ks%theory_level /= rdmft) then
        ! move allocation of vxc from ks%calc to hm
        safe_deallocate_a(hm%ks_pot%vxc)
        call move_alloc(ks%calc%vxc, hm%ks_pot%vxc)
 
        if (family_is_mgga_with_exc(hm%xc)) then
          call hm%ks_pot%set_vtau(ks%calc%vtau)
          safe_deallocate_a(ks%calc%vtau)
 
          ! We need to evaluate the energy after copying vtau to hm%vtau
          if (ks%theory_level == generalized_kohn_sham_dft .and. ks%calc%calc_energy) then
            ! MGGA vtau contribution
            if (states_are_real(st)) then
              hm%energy%intnvxc = hm%energy%intnvxc &
                + denergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_mgga)
            else
              hm%energy%intnvxc = hm%energy%intnvxc &
                + zenergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_mgga)
            end if
          end if
        end if
 
      else
        hm%ks_pot%vxc = m_zero
      end if
 
      if (.not. ks%v_ks_photons%includes_hartree()) then
        hm%energy%hartree = m_zero
        hm%ks_pot%vhartree = m_zero
      end if
 
      ! Build Hartree + XC potential
 
      do ip = 1, ks%gr%np
        hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vxc(ip, 1) + hm%ks_pot%vhartree(ip)
      end do
      if (allocated(hm%vberry)) then
        do ip = 1, ks%gr%np
          hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vhxc(ip, 1) + hm%vberry(ip, 1)
        end do
      end if
 
      if (hm%d%ispin > unpolarized) then
        do ip = 1, ks%gr%np
          hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vxc(ip, 2) + hm%ks_pot%vhartree(ip)
        end do
        if (allocated(hm%vberry)) then
          do ip = 1, ks%gr%np
            hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vhxc(ip, 2) + hm%vberry(ip, 2)
          end do
        end if
      end if
 
      if (hm%d%ispin == spinors) then
        do ispin=3, 4
          do ip = 1, ks%gr%np
            hm%ks_pot%vhxc(ip, ispin) = hm%ks_pot%vxc(ip, ispin)
          end do
        end do
      end if
 
      ! Note: this includes hybrids calculated with the Fock operator instead of OEP
      hm%energy%exchange_hf = m_zero
      if (ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
        .or. ks%theory_level == rdmft &
        .or. (ks%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(ks%xc))) then
 
        ! swap the states object
        if (.not. hm%exxop%useACE) then
          ! We also close the MPI remote memory access to the old object
          if (associated(hm%exxop%st)) then
            call states_elec_parallel_remote_access_stop(hm%exxop%st)
            call states_elec_end(hm%exxop%st)
            safe_deallocate_p(hm%exxop%st)
          end if
          ! We activate the MPI remote memory access for ks%calc%hf_st
          ! This allows to have all calls to exchange_operator_apply_standard to access
          ! the states over MPI
          call states_elec_parallel_remote_access_start(ks%calc%hf_st)
        end if
 
        ! The exchange operator will use ks%calc%hf_st
        ! For the ACE case, this is the same as st
        if (.not. optional_default(force_semilocal, .false.)) then
          select case (ks%theory_level)
          case (generalized_kohn_sham_dft)
            if (family_is_hybrid(ks%xc)) then
              call exchange_operator_reinit(hm%exxop, ks%xc%cam, ks%calc%hf_st)
            end if
          case (hartree_fock)
            call exchange_operator_reinit(hm%exxop, ks%xc%cam, ks%calc%hf_st)
          case (hartree, rdmft)
            call exchange_operator_reinit(hm%exxop, cam_exact_exchange, ks%calc%hf_st)
          end select
 
          !This should be changed and the CAM parameters should also be obtained from the restart information
          !Maybe the parameters should be mixed too.
          exx_energy = m_zero
          if (hm%exxop%useACE) then
            call xst%nullify()
            if (states_are_real(ks%calc%hf_st)) then
              ! TODO(Alex) Clean up nested if statements
              if (hm%exxop%with_isdf) then
                ! Find interpolation points from density, or read from file
                ! TODO(Alex) Issue #1195 Extend ISDF to spin-polarised systems
                call hm%exxop%isdf%get_interpolation_points(namespace, space, ks%gr, st%rho(1:ks%gr%np, 1))
                call isdf_ace_compute_potentials(hm%exxop, namespace, space, ks%gr, &
                  ks%calc%hf_st, xst, hm%kpoints)
              else
                call dexchange_operator_compute_potentials(hm%exxop, namespace, space, ks%gr, &
                  ks%calc%hf_st, xst, hm%kpoints)
              endif
              exx_energy = dexchange_operator_compute_ex(ks%gr, ks%calc%hf_st, xst)
              call dexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst)
            else
              call zexchange_operator_compute_potentials(hm%exxop, namespace, space, ks%gr, &
                ks%calc%hf_st, xst, hm%kpoints)
              exx_energy = zexchange_operator_compute_ex(ks%gr, ks%calc%hf_st, xst)
              if (hm%phase%is_allocated()) then
                call zexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst, hm%phase)
              else
                call zexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst)
              end if
            end if
            call states_elec_end(xst)
            exx_energy = exx_energy + hm%exxop%singul%energy
          end if
 
          ! Add the energy only the ACE case. In the non-ACE case, the singularity energy is added in energy_calc.F90
          select case (ks%theory_level)
          case (generalized_kohn_sham_dft)
            if (family_is_hybrid(ks%xc)) then
              hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
            end if
          case (hartree_fock)
            hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
          end select
        else
          ! If we ask for semilocal, we deactivate the exchange operator entirely
          call exchange_operator_reinit(hm%exxop, cam_null, ks%calc%hf_st)
        end if
      end if
 
    end if
 
    ! Because of the intent(in) in v_ks_calc_start, we need to update the parameters of hybrids for OEP
    ! here
    if (ks%theory_level == kohn_sham_dft .and. bitand(ks%xc_family, xc_family_oep) /= 0) then
      if (ks%xc%functional(func_x,1)%id /= xc_oep_x_slater .and. ks%xc%functional(func_x,1)%id /= xc_oep_x_fbe) then
        call exchange_operator_reinit(hm%exxop, ks%xc%cam)
      end if
    end if
 
    if (ks%v_ks_photons%active() .and. (ks%v_ks_photons%functional() == 0)) then
      call ks%v_ks_photons%add_mf_potential(ks%gr, hm%ks_pot%vhxc, hm%d%ispin, hm%ep%photon_forces(1:space%dim))
    end if
 
    if (ks%vdw%vdw_correction /= option__vdwcorrection__none) then
      assert(allocated(ks%vdw%forces))
      hm%ep%vdw_forces(:, :) = ks%vdw%forces(:, :)
      hm%ep%vdw_stress = ks%vdw%stress
      safe_deallocate_a(ks%vdw%forces)
    else
      hm%ep%vdw_forces = 0.0_real64
    end if
 
    if (ks%calc%time_present .or. hm%time_zero) then
      call hm%update(ks%gr, namespace, space, ext_partners, time = ks%calc%time)
    else
      call hamiltonian_elec_update_pot(hm, ks%gr)
    end if
 
 
    safe_deallocate_a(ks%calc%density)
    if (ks%calc%total_density_alloc) then
      safe_deallocate_p(ks%calc%total_density)
    end if
    nullify(ks%calc%total_density)
 
    pop_sub(v_ks_calc_finish)
  end subroutine v_ks_calc_finish
 
 
  !
  !! TODO(Alex) Once the implementation is finalised and benchmarked
  !! remove the serial version, and get rid of this routine.
  subroutine isdf_ace_compute_potentials(exxop, namespace, space, gr, hf_st, xst, kpoints)
    type(exchange_operator_t), intent(in ) :: exxop
    type(namespace_t),         intent(in ) :: namespace
    class(space_t),            intent(in ) :: space
    class(mesh_t),             intent(in ) :: gr
    type(states_elec_t),       intent(in ) :: hf_st
    type(kpoints_t),           intent(in ) :: kpoints
 
    type(states_elec_t),       intent(inout) :: xst
 
    if (exxop%isdf%use_serial) then
      call isdf_serial_ace_compute_potentials(exxop, namespace, space, gr, &
        hf_st, xst, kpoints)
    else
      call isdf_parallel_ace_compute_potentials(exxop, namespace, space, gr, &
        hf_st, xst, kpoints)
    endif
 
  end subroutine isdf_ace_compute_potentials
 
  ! ---------------------------------------------------------
  !
  !
  subroutine v_ks_hartree(namespace, ks, space, hm, ext_partners)
    type(namespace_t),                intent(in)    :: namespace
    type(v_ks_t),                     intent(inout) :: ks
    class(space_t),                   intent(in)    :: space
    type(hamiltonian_elec_t),         intent(inout) :: hm
    type(partner_list_t),             intent(in)    :: ext_partners
 
    push_sub(v_ks_hartree)
 
    if (.not. poisson_is_async(hm%psolver)) then
      ! solve the Poisson equation
      call dpoisson_solve(hm%psolver, namespace, hm%ks_pot%vhartree, ks%calc%total_density, reset=.false.)
    else
      ! The calculation was started by v_ks_calc_start.
      call dpoisson_solve_finish(hm%psolver, hm%ks_pot%vhartree)
    end if
 
    if (ks%calc%calc_energy) then
      ! Get the Hartree energy
      hm%energy%hartree = m_half*dmf_dotp(ks%gr, ks%calc%total_density, hm%ks_pot%vhartree)
    end if
 
    if(ks%calc%time_present) then
      if(hamiltonian_elec_has_kick(hm)) then
        call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
          ks%calc%total_density, hm%energy%pcm_corr, kick=hm%kick, time=ks%calc%time)
      else
        call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
          ks%calc%total_density, hm%energy%pcm_corr, time=ks%calc%time)
      end if
    else
      if(hamiltonian_elec_has_kick(hm)) then
        call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
          ks%calc%total_density, hm%energy%pcm_corr, kick=hm%kick)
      else
        call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
          ks%calc%total_density, hm%energy%pcm_corr)
      end if
    end if
 
    pop_sub(v_ks_hartree)
  end subroutine v_ks_hartree
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine v_ks_freeze_hxc(ks)
    type(v_ks_t), intent(inout) :: ks
 
    push_sub(v_ks_freeze_hxc)
 
    ks%frozen_hxc = .true.
 
    pop_sub(v_ks_freeze_hxc)
  end subroutine v_ks_freeze_hxc
  ! ---------------------------------------------------------
 
  subroutine v_ks_calculate_current(this, calc_cur)
    type(v_ks_t), intent(inout) :: this
    logical,      intent(in)    :: calc_cur
 
    push_sub(v_ks_calculate_current)
 
    this%calculate_current = calc_cur
 
    pop_sub(v_ks_calculate_current)
  end subroutine v_ks_calculate_current
 
  subroutine v_ks_update_dftu_energy(ks, namespace, hm, st, int_dft_u)
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(namespace_t),        intent(in)    :: namespace
    type(states_elec_t),      intent(inout) :: st
    real(real64),             intent(out)   :: int_dft_u
 
    int_dft_u = m_zero
    if (hm%lda_u_level == dft_u_none) return
 
    push_sub(v_ks_update_dftu_energy)
 
    if (states_are_real(st)) then
      int_dft_u = denergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
    else
      int_dft_u = zenergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
    end if
 
    pop_sub(v_ks_update_dftu_energy)
  end subroutine v_ks_update_dftu_energy
 
 
  logical function v_ks_is_orbital_dependent(ks)
    type(v_ks_t),  intent(in) :: ks
 
    push_sub(v_ks_is_orbital_dependent)
 
    v_ks_is_orbital_dependent = &
      ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
      .or. (ks%theory_level == kohn_sham_dft .and. xc_is_orbital_dependent(ks%xc)) &
      .or. (ks%theory_level == generalized_kohn_sham_dft .and. xc_is_orbital_dependent(ks%xc)) &
      .or. ks%sic%level == sic_pz_oep
 
    pop_sub(v_ks_is_orbital_dependent)
 
  end function v_ks_is_orbital_dependent
 
 
end module v_ks_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: