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 comm_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 kpoints_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 photon_mode_oct_m
  use photon_mode_mf_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 varinfo_oct_m
  use x_slater_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_oep_oct_m
  use xc_sic_oct_m
  use xc_photons_oct_m
  use xc_vdw_oct_m
  use xc_oep_photon_oct_m
 
  ! from the dftd3 library
  use dftd3_api
 
  implicit none
 
  private
  public ::             &
    v_ks_t,             &
    v_ks_init,          &
    v_ks_end,           &
    v_ks_write_info,    &
    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
 
  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
    integer,                  public :: xc_photon = 0
    type(xc_t),               public :: xc
    type(xc_photons_t),       public :: xc_photons
    type(xc_oep_t),           public :: oep
    type(xc_oep_photon_t),    public :: oep_photon
    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(v_ks_calc_t)                :: calc
    logical                          :: calculate_current = .false.
    type(current_t)                  :: current_calculator
    logical                          :: include_td_field = .false.
    logical,                  public :: has_photons = .false.
    logical                          :: xc_photon_include_hartree = .true.
 
    real(real64),             public :: stress_xc_gga(3, 3)
    type(photon_mode_t), pointer, public :: pt => null()
    type(mf_t),               public :: pt_mx
  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>dft</tt>. When hybrid functionals are requested, through
    !% the <tt>XCFunctional</tt> variable, the default is
    !% <tt>hartree_fock</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.
    !%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, meta-GGA functionals, or DFT+U.
    !% 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.
    !%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%oep_photon%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()
    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_exchange_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
    !%Section Hamiltonian::XC
    !%Description
    !% Defines the exchange-correlation kernel. Only LDA kernels are available currently.
    !% The options are the same as <tt>XCFunctional</tt>.
    !% Note: the kernel is only needed for Casida, Sternheimer, or optimal-control calculations.
    !% Defaults:
    !% <br>1D: <tt>lda_x_1d + lda_c_1d_csc</tt>
    !% <br>2D: <tt>lda_x_2d + lda_c_2d_amgb</tt>
    !% <br>3D: <tt>lda_x + lda_c_pz_mod</tt>
    !%Option xc_functional -1
    !% The same functional defined by <tt>XCFunctional</tt>.
    !%End
    call parse_variable(namespace, 'XCKernel', default, 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')
 
    !%Variable XCPhotonFunctional
    !%Type integer
    !%Default 0
    !%Section Hamiltonian::XC
    !%Description
    !% Defines the exchange and correlation functionals to be used for the QEDFT
    !% description of the electron-photon system.
    !%Option none      0
    !% No functional is used
    !%Option photon_x_lda      10
    !% Exchange-only local density approcimation
    !%Option photon_xc_lda     11
    !% Exchange-correlation local density approcimation
    !%Option photon_x_wfn      20
    !% Exchange-only based on wave functions
    !%Option photon_xc_wfn     21
    !% Exchange-correlation based on wave functions
    !%End
 
    call parse_variable(namespace, 'XCPhotonFunctional', option__xcphotonfunctional__none, ks%xc_photon)
 
    !%Variable XCPhotonIncludeHartree
    !%Type logical
    !%Default yes
    !%Section Hamiltonian::XC
    !%Description
    !% Use the Hartree potential and energy in calculations
    !%End
 
    call parse_variable(namespace, 'XCPhotonIncludeHartree', .true., ks%xc_photon_include_hartree)
 
    if (.not. ks%xc_photon_include_hartree) then
      call messages_write('turn off hartree potential and energy')
      call messages_warning(namespace=namespace)
    end if
 
    ! 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 (accel_is_enabled() .and. (gr%parallel_in_domains .or. st%parallel_in_states .or. st%d%kpt%parallel)) then
        !At the moment this combination produces wrong results
        call messages_not_implemented("MGGA with energy functionals and CUDA+MPI")
      end if
 
      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) then
            call messages_not_implemented("Slater with k-points", 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%has_photons) .or. (ks%xc_photon /= 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 (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%xc_photon /= 0) then
      ! initilize the photon free variables
      call ks%xc_photons%init(namespace, ks%xc_photon , space, gr, st)
      ! remornalize the electron mass due to light-matter interaction; here we only deal with it in free space
      ks%oep_photon%level = oep_level_none
    end if
 
 
    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()
        class default
          xf = pseudo_exchange_unknown
          cf = pseudo_correlation_unknown
        end select
 
        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
 
      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()
 
    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)
 
    if (ks%xc_photon /= 0) then
      call ks%xc_photons%end()
    end if
 
    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_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)
    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.)
 
    call v_ks_calc_start(ks, namespace, space, hm, st, ions, hm%kpoints%latt, ext_partners, time, &
      calc_energy, calc_current_, 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)
 
    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, calc_current, 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),               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)    :: calc_current
    logical,                 optional, intent(in)    :: force_semilocal
 
    logical :: calc_current_
 
    push_sub(v_ks_calc_start)
 
    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))
 
    call profiling_in("KOHN_SHAM_CALC")
 
    assert(.not. ks%calc%calculating)
    ks%calc%calculating = .true.
 
    if (debug%info) then
      write(message(1), '(a)') 'Debug: Calculating Kohn-Sham potential.'
      call messages_info(1, namespace=namespace)
    end if
 
    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
      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 profiling_out("KOHN_SHAM_CALC")
      pop_sub(v_ks_calc_start)
      return
    end if
 
    safe_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
 
    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
 
    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
      safe_allocate(ks%calc%hf_st)
      call states_elec_copy(ks%calc%hf_st, st)
 
      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 OpenCL/CUDA. Please use domain parallelization')
          call messages_new_line()
          call messages_write("or disable acceleration using 'DisableAccel = yes'.")
          call messages_fatal(namespace=namespace)
        end if
        call states_elec_parallel_remote_access_start(ks%calc%hf_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%has_photons) .and. (ks%calc%time_present) .and. (ks%xc_photon == 0) ) then
      call mf_calc(ks%pt_mx, ks%gr, st, ions, ks%pt, 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
      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
 
      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
 
      integer :: ispin
 
      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
 
      safe_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, 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_gga, 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, 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_gga, 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
 
      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%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%gr, st, space, &
              ks%calc%energy%correlation, vxc = ks%calc%vxc)
 
          else
            if (ks%oep_photon%level /= oep_level_none) then
              if (states_are_real(st)) then
                call dxc_oep_photon_calc(ks%oep_photon, namespace, ks%xc, ks%gr, &
                  hm, st, space, ks%calc%energy%exchange, ks%calc%energy%correlation, vxc = ks%calc%vxc)
              else
                call zxc_oep_photon_calc(ks%oep_photon, namespace, ks%xc, ks%gr, &
                  hm, st, space, ks%calc%energy%exchange, ks%calc%energy%correlation, vxc = ks%calc%vxc)
              end if
              ks%calc%energy%photon_exchange = ks%oep_photon%pt%ex
            end if
          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%xc_photon /= 0) then
 
          call ks%xc_photons%v_ks(namespace, ks%calc%total_density, ks%calc%density, ks%gr, space, hm%psolver, hm%ep, st)
 
          ! add the photon-free px potential into the xc potential
          do ispin = 1, hm%d%spin_channels
            call lalg_axpy(ks%gr%np, m_one, ks%xc_photons%vpx(1:ks%gr%np), ks%calc%vxc(1:ks%gr%np, ispin) )
          end do
 
          ! photon-exchange energy
          ks%calc%energy%photon_exchange = ks%xc_photons%ex
        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 (ks%calc%calc_energy) then
        ! MGGA vtau contribution is done after copying vtau to hm%vtau
 
        if (hm%lda_u_level /= dft_u_none) then
          if (states_are_real(st)) then
            ks%calc%energy%int_dft_u = denergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
          else
            ks%calc%energy%int_dft_u = zenergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
          end if
        end if
      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
      hm%energy%intnvstatic = dmf_dotp(ks%gr, ks%calc%total_density, hm%v_static)
    else
      hm%energy%intnvstatic = m_zero
    end if
 
    if (ks%theory_level == independent_particles .or. abs(ks%sic%amaldi_factor) <= m_epsilon) then
 
      hm%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%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
          if (ks%gr%parallel_in_domains) 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%vxc)
        call move_alloc(ks%calc%vxc, hm%vxc)
 
        if (family_is_mgga_with_exc(hm%xc)) then
          call lalg_copy(ks%gr%np, hm%d%nspin, ks%calc%vtau, hm%vtau)
          call hm%hm_base%accel_copy_pot(ks%gr, hm%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%vxc = m_zero
      end if
 
      if (.not. ks%xc_photon_include_hartree) then
        hm%energy%hartree = m_zero
        hm%vhartree = m_zero
      end if
 
      ! Build Hartree + XC potential
 
      do ip = 1, ks%gr%np
        hm%vhxc(ip, 1) = hm%vxc(ip, 1) + hm%vhartree(ip)
      end do
      if (allocated(hm%vberry)) then
        do ip = 1, ks%gr%np
          hm%vhxc(ip, 1) = hm%vhxc(ip, 1) + hm%vberry(ip, 1)
        end do
      end if
 
      if (hm%d%ispin > unpolarized) then
        do ip = 1, ks%gr%np
          hm%vhxc(ip, 2) = hm%vxc(ip, 2) + hm%vhartree(ip)
        end do
        if (allocated(hm%vberry)) then
          do ip = 1, ks%gr%np
            hm%vhxc(ip, 2) = hm%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%vhxc(ip, ispin) = hm%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)) then
 
        ! swap the states object
        if (associated(hm%exxop%st)) then
          if (hm%exxop%st%parallel_in_states) 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
        if (associated(ks%calc%hf_st) .and. hm%exxop%useACE) then
          if (ks%calc%hf_st%parallel_in_states) call states_elec_parallel_remote_access_stop(ks%calc%hf_st)
        end if
 
        !At the moment this block is called before the reinit call. This way the LCAO does not call the
        !exchange operator.
        !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 (.not. optional_default(force_semilocal, .false.)) then
          if ((ks%theory_level == hartree_fock .or. ks%theory_level == rdmft &
            .or. (ks%theory_level == generalized_kohn_sham_dft &
            .and. family_is_hybrid(ks%xc))) .and. hm%exxop%useACE) then
            call xst%nullify()
            if (states_are_real(ks%calc%hf_st)) then
              call dexchange_operator_compute_potentials(hm%exxop, namespace, space, ks%gr, &
                ks%calc%hf_st, xst, hm%kpoints, exx_energy)
              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)
              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
        end if
 
        select case (ks%theory_level)
        case (generalized_kohn_sham_dft)
          if (family_is_hybrid(ks%xc)) then
            ! Add the energy only the ACE case. In the non-ACE case, the exchange energy is not correct
            hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
            call exchange_operator_reinit(hm%exxop, ks%xc%cam_omega, ks%xc%cam_alpha, ks%xc%cam_beta, ks%calc%hf_st)
          end if
        case (hartree_fock)
          ! Add the energy only the ACE case. In the non-ACE case, the exchange energy is not correct
          hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
          call exchange_operator_reinit(hm%exxop, ks%xc%cam_omega, ks%xc%cam_alpha, ks%xc%cam_beta, ks%calc%hf_st)
        case (hartree)
          call exchange_operator_reinit(hm%exxop, m_zero, m_one, m_zero, ks%calc%hf_st)
        case (rdmft)
          call exchange_operator_reinit(hm%exxop, m_zero, m_one, m_zero, ks%calc%hf_st)
        end select
      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_omega, ks%xc%cam_alpha, ks%xc%cam_beta)
      end if
    end if
 
    if (ks%has_photons .and. (ks%xc_photon == 0)) then
      if (associated(ks%pt_mx%vmf)) then
        forall(ip = 1:ks%gr%np) hm%vhxc(ip, 1) = hm%vhxc(ip, 1) + ks%pt_mx%vmf(ip)
        if (hm%d%ispin > unpolarized) then
          forall(ip = 1:ks%gr%np) hm%vhxc(ip, 2) = hm%vhxc(ip, 2) + ks%pt_mx%vmf(ip)
        end if
      end if
      hm%ep%photon_forces(1:space%dim) = ks%pt_mx%fmf(1:space%dim)
    end if
 
    if (ks%vdw%vdw_correction /= option__vdwcorrection__none) then
      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
 
  ! ---------------------------------------------------------
  !
  !
  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%vhartree, ks%calc%total_density)
    else
      ! The calculation was started by v_ks_calc_start.
      call dpoisson_solve_finish(hm%psolver, hm%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%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%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%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%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%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
 
end module v_ks_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: