hamiltonian_elec.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2020 M. Marques, A. Castro, A. Rubio, G. Bertsch,
!!                         N. Tancogne-Dejean, M. Lueders
!!
!! 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 hamiltonian_elec_oct_m
  use absorbing_boundaries_oct_m
  use accel_oct_m
  use affine_coordinates_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use distributed_oct_m
  use energy_oct_m
  use electron_space_oct_m
  use exchange_operator_oct_m
  use external_potential_oct_m
  use hamiltonian_elec_base_oct_m
  use epot_oct_m
  use ext_partner_list_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_abst_oct_m
  use interaction_partner_oct_m
  use ion_electron_local_potential_oct_m
  use io_oct_m
  use ions_oct_m
  use kick_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use ks_potential_oct_m
  use lalg_basic_oct_m
  use lasers_oct_m
  use lattice_vectors_oct_m
  use lda_u_oct_m
  use linked_list_oct_m
  use magnetic_constrain_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use mxll_elec_coupling_oct_m
  use namespace_oct_m
  use nlcc_oct_m
  use nonlocal_pseudopotential_oct_m
  use oct_exchange_oct_m
  use output_low_oct_m
  use parser_oct_m
  use par_vec_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use projector_oct_m
  use pcm_oct_m
  use phase_oct_m
  use restart_oct_m
  use scissor_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 symmetries_oct_m
  use symm_op_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
  use xc_functional_oct_m
  use xc_oct_m
  use xc_f03_lib_m
  use xc_functional_oct_m
  use xc_interaction_oct_m
  use xc_photons_oct_m
  use zora_oct_m
 
  implicit none
 
  private
  public ::                               &
    hamiltonian_elec_t,                   &
    hamiltonian_elec_init,                &
    hamiltonian_elec_end,                 &
    dhamiltonian_elec_apply_single,       &
    zhamiltonian_elec_apply_single,       &
    zhamiltonian_elec_apply_all,          &
    dhamiltonian_elec_apply_batch,        &
    zhamiltonian_elec_apply_batch,        &
    hamiltonian_elec_diagonal,            &
    magnus,                               &
    dvmask,                               &
    zvmask,                               &
    hamiltonian_elec_inh_term,            &
    hamiltonian_elec_set_inh,             &
    hamiltonian_elec_remove_inh,          &
    hamiltonian_elec_adjoint,             &
    hamiltonian_elec_not_adjoint,         &
    hamiltonian_elec_epot_generate,       &
    hamiltonian_elec_needs_current,       &
    hamiltonian_elec_update_pot,          &
    hamiltonian_elec_update_with_ext_pot, &
    hamiltonian_elec_get_time,            &
    hamiltonian_elec_apply_packed,        &
    zhamiltonian_elec_apply_atom,         &
    hamiltonian_elec_has_kick,            &
    hamiltonian_elec_copy_and_set_phase
 
 
  type, extends(hamiltonian_abst_t) :: hamiltonian_elec_t
    ! Components are public by default
 
    type(space_t), private :: space
    type(states_elec_dim_t)  :: d
    type(hamiltonian_elec_base_t) :: hm_base
    type(phase_t) :: phase
    type(energy_t), allocatable  :: energy
    type(absorbing_boundaries_t) :: abs_boundaries
    type(ks_potential_t) :: ks_pot
    real(real64), allocatable :: vberry(:,:)
 
    type(derivatives_t), pointer, private :: der
 
    type(nonlocal_pseudopotential_t) :: vnl
 
    type(ions_t),     pointer :: ions
    real(real64) :: exx_coef
 
    type(poisson_t)          :: psolver
 
    logical :: self_induced_magnetic
    real(real64), allocatable :: a_ind(:, :)
    real(real64), allocatable :: b_ind(:, :)
 
    integer :: theory_level
    type(xc_t),             pointer :: xc
    type(xc_photons_t), pointer :: xc_photons
 
    type(epot_t) :: ep
    type(pcm_t)  :: pcm
 
    logical, private :: adjoint
 
    real(real64), private :: mass
 
    logical, private :: inh_term
    type(states_elec_t) :: inh_st
 
    type(oct_exchange_t) :: oct_exchange
 
    type(scissor_t) :: scissor
 
    real(real64) :: current_time
    logical, private :: is_applied_packed
 
    type(lda_u_t) :: lda_u
    integer       :: lda_u_level
 
    logical, public :: time_zero
 
    type(exchange_operator_t), public :: exxop
 
    type(kpoints_t), pointer, public :: kpoints => null()
 
    type(partner_list_t) :: external_potentials
    real(real64), allocatable, public  :: v_ext_pot(:)
    real(real64), allocatable, public  :: v_static(:)
 
    type(ion_electron_local_potential_t) :: v_ie_loc
    type(nlcc_t) :: nlcc
 
    type(magnetic_constrain_t) :: magnetic_constrain
 
    type(kick_t) :: kick
 
    type(mxll_coupling_t) :: mxll
    type(zora_t), pointer :: zora => null()
 
  contains
    procedure :: update => hamiltonian_elec_update
    procedure :: apply_packed => hamiltonian_elec_apply_packed
    procedure :: update_span => hamiltonian_elec_span
    procedure :: dapply => dhamiltonian_elec_apply
    procedure :: zapply => zhamiltonian_elec_apply
    procedure :: dmagnus_apply => dhamiltonian_elec_magnus_apply
    procedure :: zmagnus_apply => zhamiltonian_elec_magnus_apply
    procedure :: is_hermitian => hamiltonian_elec_hermitian
    procedure :: needs_mgga_term => hamiltonian_elec_needs_mgga_term
    procedure :: set_mass => hamiltonian_elec_set_mass
  end type hamiltonian_elec_t
 
  integer, public, parameter :: &
    LENGTH     = 1,             &
    velocity   = 2
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_init(hm, namespace, space, gr, ions, ext_partners, st, theory_level, xc, &
    mc, kpoints, need_exchange, xc_photons)
    type(hamiltonian_elec_t),           target, intent(inout) :: hm
    type(namespace_t),                          intent(in)    :: namespace
    class(space_t),                             intent(in)    :: space
    type(grid_t),                       target, intent(inout) :: gr
    type(ions_t),                       target, intent(inout) :: ions
    type(partner_list_t),                       intent(inout) :: ext_partners
    type(states_elec_t),                target, intent(inout) :: st
    integer,                                    intent(in)    :: theory_level
    type(xc_t),                         target, intent(in)    :: xc
    type(multicomm_t),                          intent(in)    :: mc
    type(kpoints_t),                    target, intent(in)    :: kpoints
    logical, optional,                          intent(in)    :: need_exchange
    type(xc_photons_t), optional,   target, intent(in)    :: xc_photons
 
 
    logical :: need_exchange_
    real(real64) :: rashba_coupling
 
    push_sub(hamiltonian_elec_init)
    call profiling_in('HAMILTONIAN_ELEC_INIT')
 
    ! make a couple of local copies
    hm%space = space
    hm%theory_level = theory_level
    call states_elec_dim_copy(hm%d, st%d)
 
    hm%kpoints => kpoints
 
    !%Variable ParticleMass
    !%Type float
    !%Default 1.0
    !%Section Hamiltonian
    !%Description
    !% It is possible to make calculations for a particle with a mass
    !% different from one (atomic unit of mass, or mass of the electron).
    !% This is useful to describe non-electronic systems, or for
    !% esoteric purposes.
    !%End
    call parse_variable(namespace, 'ParticleMass', m_one, hm%mass)
 
    !%Variable RashbaSpinOrbitCoupling
    !%Type float
    !%Default 0.0
    !%Section Hamiltonian
    !%Description
    !% (Experimental.) For systems described in 2D (electrons confined to 2D in semiconductor structures), one
    !% may add the Bychkov-Rashba spin-orbit coupling term [Bychkov and Rashba, <i>J. Phys. C: Solid
    !% State Phys.</i> <b>17</b>, 6031 (1984)]. This variable determines the strength
    !% of this perturbation, and has dimensions of energy times length.
    !%End
    call parse_variable(namespace, 'RashbaSpinOrbitCoupling', m_zero, rashba_coupling, units_inp%energy*units_inp%length)
    if (parse_is_defined(namespace, 'RashbaSpinOrbitCoupling')) then
      if (space%dim /= 2) then
        write(message(1),'(a)') 'Rashba spin-orbit coupling can only be used for two-dimensional systems.'
        call messages_fatal(1, namespace=namespace)
      end if
      call messages_experimental('RashbaSpinOrbitCoupling', namespace=namespace)
    end if
 
    call hm%hm_base%init(hm%d%nspin, hm%mass, rashba_coupling)
    call hm%vnl%init()
 
    assert(associated(gr%der%lapl))
    hm%hm_base%kinetic => gr%der%lapl
 
    safe_allocate(hm%energy)
 
    !Keep pointers to derivatives, geometry and xc
    hm%der => gr%der
    hm%ions => ions
    hm%xc => xc
 
    if(present(xc_photons)) then
      hm%xc_photons => xc_photons
    else
      hm%xc_photons => null()
    end if
 
    ! allocate potentials and density of the cores
    ! In the case of spinors, vxc_11 = hm%vxc(:, 1), vxc_22 = hm%vxc(:, 2), Re(vxc_12) = hm%vxc(:. 3);
    ! Im(vxc_12) = hm%vxc(:, 4)
    call hm%ks_pot%init(gr%der, gr%np, gr%np_part, hm%d%nspin, hm%theory_level, family_is_mgga_with_exc(hm%xc))
 
    !Initialize Poisson solvers
    call poisson_init(hm%psolver, namespace, space, gr%der, mc, gr%stencil, st%qtot)
 
    ! Initialize external potential
    call epot_init(hm%ep, namespace, gr, hm%ions, hm%psolver, hm%d%ispin, hm%xc%family, hm%kpoints)
    call kick_init(hm%kick, namespace, space, hm%kpoints, hm%d%ispin)
 
    hm%zora => zora_t(namespace, hm%der, hm%d, hm%ep, hm%mass)
 
    !Temporary construction of the ion-electron interactions
    call hm%v_ie_loc%init(gr, hm%psolver, hm%ions, namespace)
    if (hm%ep%nlcc) then
      call hm%nlcc%init(gr, hm%ions)
      safe_allocate(st%rho_core(1:gr%np))
      st%rho_core(:) = m_zero
    end if
 
    !Static magnetic field or rashba spin-orbit interaction requires complex wavefunctions
    if (parse_is_defined(namespace, 'StaticMagneticField') .or. list_has_gauge_field(ext_partners) .or. &
      parse_is_defined(namespace, 'RashbaSpinOrbitCoupling')) then
      call states_set_complex(st)
    end if
 
    !%Variable CalculateSelfInducedMagneticField
    !%Type logical
    !%Default no
    !%Section Hamiltonian
    !%Description
    !% The existence of an electronic current implies the creation of a self-induced magnetic
    !% field, which may in turn back-react on the system. Of course, a fully consistent treatment
    !% of this kind of effect should be done in QED theory, but we will attempt a first
    !% approximation to the problem by considering the lowest-order relativistic terms
    !% plugged into the normal Hamiltonian equations (spin-other-orbit coupling terms, etc.).
    !% For the moment being, none of this is done, but a first step is taken by calculating
    !% the induced magnetic field of a system that has a current, by considering the magnetostatic
    !% approximation and Biot-Savart law:
    !%
    !% <math> \nabla^2 \vec{A} + 4\pi\alpha \vec{J} = 0</math>
    !%
    !% <math> \vec{B} = \vec{\nabla} \times \vec{A}</math>
    !%
    !% If <tt>CalculateSelfInducedMagneticField</tt> is set to yes, this <i>B</i> field is
    !% calculated at the end of a <tt>gs</tt> calculation (nothing is done -- yet -- in the <tt>td</tt>case)
    !% and printed out, if the <tt>Output</tt> variable contains the <tt>potential</tt> keyword (the prefix
    !% of the output files is <tt>Bind</tt>).
    !%End
    call parse_variable(namespace, 'CalculateSelfInducedMagneticField', .false., hm%self_induced_magnetic)
    if (hm%self_induced_magnetic) then
      safe_allocate(hm%a_ind(1:gr%np_part, 1:space%dim))
      safe_allocate(hm%b_ind(1:gr%np_part, 1:space%dim))
 
      !(for dim = we could save some memory, but it is better to keep it simple)
    end if
 
    ! Absorbing boundaries
    call absorbing_boundaries_init(hm%abs_boundaries, namespace, space, gr)
 
    hm%inh_term = .false.
    call oct_exchange_remove(hm%oct_exchange)
 
    hm%adjoint = .false.
 
    call hm%phase%init(gr, hm%d%kpt, hm%kpoints, st%d, space)
 
    !%Variable DFTULevel
    !%Type integer
    !%Default no
    !%Section Hamiltonian::XC
    !%Description
    !% This variable selects which DFT+U expression is added to the Hamiltonian.
    !%Option dft_u_none 0
    !% No +U term is not applied.
    !%Option dft_u_empirical 1
    !% An empiricial Hubbard U is added on the orbitals specified in the block species
    !% with hubbard_l and hubbard_u
    !%Option dft_u_acbn0 2
    !% Octopus determines the effective U term using the
    !% ACBN0 functional as defined in PRX 5, 011006 (2015)
    !%End
    call parse_variable(namespace, 'DFTULevel', dft_u_none, hm%lda_u_level)
    call messages_print_var_option('DFTULevel', hm%lda_u_level, namespace=namespace)
    if (hm%lda_u_level /= dft_u_none) then
      call lda_u_init(hm%lda_u, namespace, space, hm%lda_u_level, gr, ions, st, mc, &
        hm%kpoints, hm%phase%is_allocated())
 
      !In the present implementation of DFT+U, in case of spinors, we have off-diagonal terms
      !in spin space which break the assumption of the generalized Bloch theorem
      if (kick_get_type(hm%kick) == kick_magnon_mode .and. gr%der%boundaries%spiral) then
        call messages_not_implemented("DFT+U with generalized Bloch theorem and magnon kick", namespace=namespace)
      end if
 
      ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
      if(hm%lda_u_level /= dft_u_none .and. hm%phase%is_allocated()) then
        call lda_u_build_phase_correction(hm%lda_u, space, hm%d, gr%der%boundaries, namespace, hm%kpoints)
      end if
    end if
 
    !%Variable HamiltonianApplyPacked
    !%Type logical
    !%Default yes
    !%Section Execution::Optimization
    !%Description
    !% If set to yes (the default), Octopus will 'pack' the
    !% wave-functions when operating with them. This might involve some
    !% additional copying but makes operations more efficient.
    !% See also the related <tt>StatesPack</tt> variable.
    !%End
    call parse_variable(namespace, 'HamiltonianApplyPacked', .true., hm%is_applied_packed)
 
    if (hm%theory_level == hartree_fock .and. st%parallel_in_states) then
      call messages_experimental('Hartree-Fock parallel in states', namespace=namespace)
    end if
 
    if (hm%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(hm%xc) &
      .and. st%parallel_in_states) then
      call messages_experimental('Hybrid functionals parallel in states', namespace=namespace)
    end if
 
 
    !%Variable TimeZero
    !%Type logical
    !%Default no
    !%Section Hamiltonian
    !%Description
    !% (Experimental) If set to yes, the ground state and other time
    !% dependent calculation will assume that they are done at time
    !% zero, so that all time depedent field at that time will be
    !% included.
    !%End
    call parse_variable(namespace, 'TimeZero', .false., hm%time_zero)
    if (hm%time_zero) call messages_experimental('TimeZero', namespace=namespace)
 
    !Cam parameters are irrelevant here and are updated later
    need_exchange_ = optional_default(need_exchange, .false.)
    if (hm%theory_level == hartree_fock .or. hm%theory_level == hartree &
      .or. hm%theory_level == rdmft .or. need_exchange_ .or. &
      (hm%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(hm%xc)) &
      .or. hm%xc%functional(func_x,1)%id == xc_oep_x_slater &
      .or. bitand(hm%xc%family, xc_family_oep) /= 0) then
      !We test Slater before OEP, as Slater is treated as OEP for the moment....
      if (hm%xc%functional(func_x,1)%id == xc_oep_x_slater) then
        call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
          hm%kpoints, m_zero, m_one, m_zero)
      else if (bitand(hm%xc%family, xc_family_oep) /= 0 .or. hm%theory_level == rdmft) then
        call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
          hm%kpoints, hm%xc%cam_omega, hm%xc%cam_alpha, hm%xc%cam_beta)
        if (hm%theory_level == rdmft) hm%exxop%useACE = .false.
      else
        call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
          hm%kpoints, m_zero, m_one, m_zero)
      end if
    end if
 
    if (hm%is_applied_packed .and. accel_is_enabled()) then
      ! Check if we can actually apply the hamiltonian packed
      if (gr%use_curvilinear) then
        if (accel_allow_cpu_only()) then
          hm%is_applied_packed = .false.
          call messages_write('Cannot use CUDA or OpenCL as curvilinear coordinates are used.')
          call messages_warning(namespace=namespace)
        else
          call messages_write('Cannot use CUDA or OpenCL as curvilinear coordinates are used.', new_line = .true.)
          call messages_write('Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
          call messages_fatal(namespace=namespace)
        end if
      end if
    end if
 
    !We are building the list of external potentials
    !This is done here at the moment, because we pass directly the mesh
    !TODO: Once the abstract Hamiltonian knows about an abstract basis, we might move this to the
    !      abstract Hamiltonian
    call load_external_potentials(hm%external_potentials, namespace)
 
    !Some checks which are electron specific, like k-points
    call external_potentials_checks()
 
    !At the moment we do only have static external potential, so we never update them
    call build_external_potentials()
 
    !Build the resulting interactions
    !TODO: This will be moved to the actual interactions
    call build_interactions()
 
    ! Constrained DFT for noncollinear magnetism
    if (hm%theory_level /= independent_particles) then
      call magnetic_constrain_init(hm%magnetic_constrain, namespace, gr, st%d, ions%natoms, ions%min_distance())
    end if
 
    ! init maxwell-electrons coupling
    call mxll_coupling_init(hm%mxll, st%d, gr, namespace, hm%mass)
 
    if (associated(hm%xc_photons)) then
      if (hm%xc_photons%wants_to_renormalize_mass()) then
        ! remornalize the electron mass due to light-matter interaction; here we only deal with it in free space
        call hm%set_mass(namespace, hm%xc_photons%get_renormalized_mass())
      end if
    end if
 
 
    call profiling_out('HAMILTONIAN_ELEC_INIT')
    pop_sub(hamiltonian_elec_init)
 
  contains
 
    ! ---------------------------------------------------------
    subroutine build_external_potentials()
      type(list_iterator_t) :: iter
      class(*), pointer :: potential
      integer :: iop
 
      push_sub(hamiltonian_elec_init.build_external_potentials)
 
      safe_allocate(hm%v_ext_pot(1:gr%np))
      hm%v_ext_pot(1:gr%np) = m_zero
 
      call iter%start(hm%external_potentials)
      do while (iter%has_next())
        potential => iter%get_next()
        select type (potential)
        class is (external_potential_t)
 
          call potential%allocate_memory(gr)
          call potential%calculate(namespace, gr, hm%psolver)
          !To preserve the old behavior, we are adding the various potentials
          !to the corresponding arrays
          select case (potential%type)
          case (external_pot_usdef, external_pot_from_file, external_pot_charge_density)
            call lalg_axpy(gr%np, m_one, potential%pot, hm%v_ext_pot)
 
          case (external_pot_static_bfield)
            if (.not. allocated(hm%ep%b_field)) then
              safe_allocate(hm%ep%b_field(1:3)) !Cannot be space%dim
              hm%ep%b_field(1:3) = m_zero
            end if
            hm%ep%b_field(1:3) = hm%ep%b_field(1:3) + potential%b_field(1:3)
 
            if (.not. allocated(hm%ep%a_static)) then
              safe_allocate(hm%ep%a_static(1:gr%np, 1:space%dim))
              hm%ep%a_static(1:gr%np, 1:space%dim) = m_zero
            end if
            call lalg_axpy(gr%np, space%dim, m_one, potential%a_static, hm%ep%a_static)
 
          case (external_pot_static_efield)
            if (.not. allocated(hm%ep%e_field)) then
              safe_allocate(hm%ep%e_field(1:space%dim))
              hm%ep%e_field(1:space%dim) = m_zero
            end if
            hm%ep%e_field(1:space%dim) = hm%ep%e_field(1:space%dim) + potential%e_field(1:space%dim)
 
            !In the fully periodic case, we use Berry phases
            if (space%periodic_dim < space%dim) then
              if (.not. allocated(hm%v_static)) then
                safe_allocate(hm%v_static(1:gr%np))
                hm%v_static(1:gr%np) = m_zero
              end if
              if (.not. allocated(hm%ep%v_ext)) then
                safe_allocate(hm%ep%v_ext(1:gr%np_part))
                hm%ep%v_ext(1:gr%np_part) = m_zero
              end if
              call lalg_axpy(gr%np, m_one, potential%pot, hm%v_static)
              call lalg_axpy(gr%np, m_one, potential%v_ext, hm%ep%v_ext)
            end if
 
            if (hm%kpoints%use_symmetries) then
              do iop = 1, symmetries_number(hm%kpoints%symm)
                if (iop == symmetries_identity_index(hm%kpoints%symm)) cycle
                if (.not. symm_op_invariant_cart(hm%kpoints%symm%ops(iop), hm%ep%e_field, 1e-5_real64)) then
                  message(1) = "The StaticElectricField breaks (at least) one of the symmetries used to reduce the k-points."
                  message(2) = "Set SymmetryBreakDir equal to StaticElectricField."
                  call messages_fatal(2, namespace=namespace)
                end if
              end do
            end if
 
          end select
          call potential%deallocate_memory()
 
        class default
          assert(.false.)
        end select
      end do
 
      pop_sub(hamiltonian_elec_init.build_external_potentials)
    end subroutine build_external_potentials
 
    ! ---------------------------------------------------------
    subroutine external_potentials_checks()
      type(list_iterator_t) :: iter
      class(*), pointer :: potential
 
      push_sub(hamiltonian_elec_init.external_potentials_checks)
 
      call iter%start(hm%external_potentials)
      do while (iter%has_next())
        potential => iter%get_next()
        select type (potential)
        class is (external_potential_t)
 
          if (potential%type == external_pot_static_efield .and. hm%kpoints%reduced%npoints > 1) then
            message(1) = "Applying StaticElectricField in a periodic direction is only accurate for large supercells."
            message(2) = "Single-point Berry phase is not appropriate when k-point sampling is needed."
            call messages_warning(2, namespace=namespace)
          end if
 
        class default
          assert(.false.)
        end select
      end do
 
      pop_sub(hamiltonian_elec_init.external_potentials_checks)
    end subroutine external_potentials_checks
 
 
    !The code in this routines needs to know about the external potentials.
    !This will be treated in the future by the interactions directly.
    subroutine build_interactions()
      logical :: external_potentials_present
      logical :: kick_present
 
      push_sub(hamiltonian_elec_init.build_interactions)
 
      if (allocated(hm%ep%e_field) .and. space%is_periodic() .and. .not. list_has_gauge_field(ext_partners)) then
        ! only need vberry if there is a field in a periodic direction
        ! and we are not setting a gauge field
        if (any(abs(hm%ep%e_field(1:space%periodic_dim)) > m_epsilon)) then
          safe_allocate(hm%vberry(1:gr%np, 1:hm%d%nspin))
          hm%vberry = m_zero
        end if
      end if
 
      external_potentials_present = epot_have_external_potentials(hm%ep) .or. &
        list_has_lasers(ext_partners) .or. allocated(hm%v_static)
 
 
      kick_present = hamiltonian_elec_has_kick(hm)
 
      call pcm_init(hm%pcm, namespace, space, ions, gr, st%qtot, st%val_charge, external_potentials_present, kick_present)
      if (hm%pcm%run_pcm) then
        if (hm%theory_level /= kohn_sham_dft) call messages_not_implemented("PCM for TheoryLevel /= kohn_sham", namespace=namespace)
      end if
 
      pop_sub(hamiltonian_elec_init.build_interactions)
 
    end subroutine build_interactions
 
 
  end subroutine hamiltonian_elec_init
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_end(hm)
    type(hamiltonian_elec_t), target, intent(inout) :: hm
 
    type(partner_iterator_t) :: iter
    class(interaction_partner_t), pointer :: potential
 
    push_sub(hamiltonian_elec_end)
 
    call hm%hm_base%end()
    call hm%vnl%end()
 
    call hm%phase%end()
 
    call hm%ks_pot%end()
    safe_deallocate_a(hm%vberry)
    safe_deallocate_a(hm%a_ind)
    safe_deallocate_a(hm%b_ind)
    safe_deallocate_a(hm%v_ext_pot)
 
    safe_deallocate_p(hm%zora)
 
    call poisson_end(hm%psolver)
 
    nullify(hm%xc)
 
    call kick_end(hm%kick)
    call epot_end(hm%ep)
    nullify(hm%ions)
 
    call absorbing_boundaries_end(hm%abs_boundaries)
 
    call states_elec_dim_end(hm%d)
 
    if (hm%scissor%apply) call scissor_end(hm%scissor)
 
    call exchange_operator_end(hm%exxop)
    call lda_u_end(hm%lda_u)
 
    safe_deallocate_a(hm%energy)
 
    if (hm%pcm%run_pcm) call pcm_end(hm%pcm)
 
    call hm%v_ie_loc%end()
    call hm%nlcc%end()
 
    call iter%start(hm%external_potentials)
    do while (iter%has_next())
      potential => iter%get_next()
      safe_deallocate_p(potential)
    end do
    call hm%external_potentials%empty()
    safe_deallocate_a(hm%v_static)
 
    call magnetic_constrain_end(hm%magnetic_constrain)
 
    call mxll_coupling_end(hm%mxll)
 
    pop_sub(hamiltonian_elec_end)
  end subroutine hamiltonian_elec_end
 
 
  ! ---------------------------------------------------------
  ! True if the Hamiltonian is Hermitian, false otherwise
  logical function hamiltonian_elec_hermitian(hm)
    class(hamiltonian_elec_t), intent(in) :: hm
 
    push_sub(hamiltonian_elec_hermitian)
    hamiltonian_elec_hermitian = .not.((hm%abs_boundaries%abtype == imaginary_absorbing) .or. &
      oct_exchange_enabled(hm%oct_exchange))
 
    pop_sub(hamiltonian_elec_hermitian)
  end function hamiltonian_elec_hermitian
 
  ! ---------------------------------------------------------
  ! True if the Hamiltonian needs to apply the mGGA term
  pure logical function hamiltonian_elec_needs_mgga_term(hm, terms)
    class(hamiltonian_elec_t), intent(in) :: hm
    integer,                   intent(in) :: terms
 
    hamiltonian_elec_needs_mgga_term = .false.
 
    if (bitand(term_mgga, terms) /= 0 .and. family_is_mgga_with_exc(hm%xc) &
      .and. hm%theory_level == generalized_kohn_sham_dft) then
      hamiltonian_elec_needs_mgga_term = .true.
    end if
 
    !For OEP
    if(terms == term_mgga .and. family_is_mgga_with_exc(hm%xc) .and. hm%theory_level == kohn_sham_dft) then
      hamiltonian_elec_needs_mgga_term = .true.
    end if
 
  end function hamiltonian_elec_needs_mgga_term
 
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_span(hm, delta, emin, namespace)
    class(hamiltonian_elec_t), intent(inout) :: hm
    real(real64),              intent(in)    :: delta(:)
    real(real64),              intent(in)    :: emin
    type(namespace_t),         intent(in)    :: namespace
 
    real(real64) :: emax
 
    push_sub(hamiltonian_elec_span)
 
    ! estimate maximum energy of discrete kinetic operator
    ! this neglects possible contributions from the non-local part of the pseudopotentials
    emax = derivatives_lapl_get_max_eigenvalue(hm%der)/m_two
 
    hm%spectral_middle_point = (emax + emin) / m_two
    hm%spectral_half_span    = (emax - emin) / m_two
 
    pop_sub(hamiltonian_elec_span)
  end subroutine hamiltonian_elec_span
 
 
  ! ---------------------------------------------------------
  pure logical function hamiltonian_elec_inh_term(hm) result(inh)
    type(hamiltonian_elec_t), intent(in) :: hm
 
    inh = hm%inh_term
  end function hamiltonian_elec_inh_term
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_set_inh(hm, st)
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(states_elec_t), intent(in)    :: st
 
    push_sub(hamiltonian_elec_set_inh)
 
    if (hm%inh_term) call states_elec_end(hm%inh_st)
    call states_elec_copy(hm%inh_st, st)
    hm%inh_term = .true.
 
    pop_sub(hamiltonian_elec_set_inh)
  end subroutine hamiltonian_elec_set_inh
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_remove_inh(hm)
    type(hamiltonian_elec_t), intent(inout) :: hm
 
    push_sub(hamiltonian_elec_remove_inh)
 
    if (hm%inh_term) then
      call states_elec_end(hm%inh_st)
      hm%inh_term = .false.
    end if
 
    pop_sub(hamiltonian_elec_remove_inh)
  end subroutine hamiltonian_elec_remove_inh
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_adjoint(hm)
    type(hamiltonian_elec_t), intent(inout) :: hm
 
    push_sub(hamiltonian_elec_adjoint)
 
    if (.not. hm%adjoint) then
      hm%adjoint = .true.
      if (hm%abs_boundaries%abtype == imaginary_absorbing) then
        hm%abs_boundaries%mf = -hm%abs_boundaries%mf
      end if
    end if
 
    pop_sub(hamiltonian_elec_adjoint)
  end subroutine hamiltonian_elec_adjoint
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_not_adjoint(hm)
    type(hamiltonian_elec_t), intent(inout) :: hm
 
    push_sub(hamiltonian_elec_not_adjoint)
 
    if (hm%adjoint) then
      hm%adjoint = .false.
      if (hm%abs_boundaries%abtype == imaginary_absorbing) then
        hm%abs_boundaries%mf = -hm%abs_boundaries%mf
      end if
    end if
 
    pop_sub(hamiltonian_elec_not_adjoint)
  end subroutine hamiltonian_elec_not_adjoint
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_update(this, mesh, namespace, space, ext_partners, time)
    class(hamiltonian_elec_t), intent(inout) :: this
    class(mesh_t),            intent(in)     :: mesh
    type(namespace_t),        intent(in)     :: namespace
    class(space_t),           intent(in)     :: space
    type(partner_list_t),     intent(in)     :: ext_partners
    real(real64), optional,          intent(in)     :: time
 
    integer :: ispin, ip, idir, iatom, ilaser
    real(real64)                 :: aa(space%dim), time_
    real(real64),    allocatable :: vp(:,:)
    type(lasers_t),      pointer :: lasers
    type(gauge_field_t), pointer :: gfield
    real(real64) :: am(space%dim)
 
    push_sub(hamiltonian_elec_update)
    call profiling_in("HAMILTONIAN_ELEC_UPDATE")
 
    this%current_time = m_zero
    if (present(time)) this%current_time = time
 
    time_ = optional_default(time, 0.0_real64)
 
    ! set everything to zero
    call this%hm_base%clear(mesh%np)
 
    ! alllocate the scalar potential for the xc, hartree and external potentials
    call this%hm_base%allocate_field(mesh, field_potential, &
      complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
 
    ! the lasers
    if (present(time) .or. this%time_zero) then
 
      lasers => list_get_lasers(ext_partners)
      if(associated(lasers)) then
        do ilaser = 1, lasers%no_lasers
          select case (laser_kind(lasers%lasers(ilaser)))
          case (e_field_scalar_potential, e_field_electric)
            do ispin = 1, this%d%spin_channels
              call laser_potential(lasers%lasers(ilaser), mesh,  &
                this%hm_base%potential(:, ispin), time_)
            end do
          case (e_field_magnetic)
            call this%hm_base%allocate_field(mesh, field_vector_potential + field_uniform_magnetic_field, &
              .false.)
            ! get the vector potential
            safe_allocate(vp(1:mesh%np, 1:space%dim))
            vp(1:mesh%np, 1:space%dim) = m_zero
            call laser_vector_potential(lasers%lasers(ilaser), mesh, vp, time_)
            !$omp parallel do private(idir) schedule(static)
            do ip = 1, mesh%np
              do idir = 1, space%dim
                this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + vp(ip, idir)/p_c
              end do
            end do
            ! and the magnetic field
            call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
            safe_deallocate_a(vp)
          case (e_field_vector_potential)
            call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
            ! get the uniform vector potential associated with a magnetic field
            aa = m_zero
            call laser_field(lasers%lasers(ilaser), aa, time_)
            this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
          end select
        end do
 
        if (lasers_with_nondipole_field(lasers)) then
          assert( allocated(this%hm_base%uniform_vector_potential))
          call lasers_nondipole_laser_field_step(lasers, am, time_)
          this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - am/p_c
        end if
      end if
 
      ! the gauge field
      gfield => list_get_gauge_field(ext_partners)
      if (associated(gfield)) then
        call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
        call gauge_field_get_vec_pot(gfield, aa)
        this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
      end if
 
      ! the electric field for a periodic system through the gauge field
      if (allocated(this%ep%e_field) .and. associated(gfield)) then
        this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
          this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
      end if
 
      ! add the photon-free mean-field vector potential
      if (associated(this%xc_photons)) then
        if(this%xc_photons%lpfmf .and. allocated(this%xc_photons%mf_vector_potential)) then
          call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
          ! here we put a minus sign in front of the mean field term to get the right answer (need to check the formula)
          this%hm_base%uniform_vector_potential(1:space%dim) = &
            this%hm_base%uniform_vector_potential(1:space%dim) - this%xc_photons%mf_vector_potential(1:space%dim)/p_c
        end if
      end if
 
    end if
 
    ! the vector potential of a static magnetic field
    if (allocated(this%ep%a_static)) then
      call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
      !ep%a_static contains 1/c A(r)
      !$omp parallel do private(idir) schedule(static)
      do ip = 1, mesh%np
        do idir = 1, space%dim
          this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
        end do
      end do
    end if
 
    ! add Maxwell coupling to Hamiltonian and sets the magnetic field for the Zeeman term added below
    call mxll_coupling_calc(this%mxll, this%hm_base, mesh, this%d, space)
 
    !The electric field was added to the KS potential
    call this%hm_base%accel_copy_pot(mesh)
 
    ! and the static magnetic field
    if (allocated(this%ep%b_field)) then
      call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
      do idir = 1, 3
        this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
      end do
    end if
 
    ! Combine the uniform and non-uniform fields and compute the Zeeman term
    call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
 
    ! This needs to be called at the end as the zeeman term enters the potential
    call hamiltonian_elec_update_pot(this, mesh, accumulate = .true.)
 
    if (this%mxll%test_equad) then
      call set_electric_quadrupole_pot(this%mxll, mesh)
    end if
 
    call build_phase()
 
    call profiling_out("HAMILTONIAN_ELEC_UPDATE")
    pop_sub(hamiltonian_elec_update)
 
  contains
 
    subroutine build_phase()
      integer :: ik, imat, nmat, max_npoints, offset
      integer :: ip
      integer :: iphase, nphase
 
      push_sub(hamiltonian_elec_update.build_phase)
 
      if ((.not. this%kpoints%gamma_only()) .or. allocated(this%hm_base%uniform_vector_potential)) then
 
        call profiling_in('UPDATE_PHASES')
        ! now regenerate the phases for the pseudopotentials
        do iatom = 1, this%ep%natoms
          call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints, &
            vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
        end do
 
        call profiling_out('UPDATE_PHASES')
      end if
 
      if (allocated(this%hm_base%uniform_vector_potential)) then
 
        call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
 
        ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
        if (this%lda_u_level /= dft_u_none) then
          call lda_u_build_phase_correction(this%lda_u, space, this%d, this%der%boundaries, namespace, this%kpoints, &
            vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
        end if
      end if
 
      max_npoints = this%vnl%max_npoints
      nmat = this%vnl%nprojector_matrices
 
 
      if (this%phase%is_allocated() .and. allocated(this%vnl%projector_matrices)) then
 
        nphase = 1
        if (this%der%boundaries%spiralBC) nphase = 3
 
        if (.not. allocated(this%vnl%projector_phases)) then
          safe_allocate(this%vnl%projector_phases(1:max_npoints, 1:nphase, nmat, this%d%kpt%start:this%d%kpt%end))
          if (accel_is_enabled()) then
            call accel_create_buffer(this%vnl%buff_projector_phases, accel_mem_read_only, &
              type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
            ! We need to save nphase, with which the array has been build,
            ! as the number might change throughout the run
            this%vnl%nphase = nphase
          end if
        end if
 
        offset = 0
        do ik = this%d%kpt%start, this%d%kpt%end
          do imat = 1, this%vnl%nprojector_matrices
            iatom = this%vnl%projector_to_atom(imat)
            do iphase = 1, nphase
              !$omp parallel do simd schedule(static)
              do ip = 1, this%vnl%projector_matrices(imat)%npoints
                this%vnl%projector_phases(ip, iphase, imat, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
              end do
 
              if (accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0) then
                call accel_write_buffer(this%vnl%buff_projector_phases, &
                  this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
                  offset = offset, async=.true.)
              end if
              offset = offset + this%vnl%projector_matrices(imat)%npoints
            end do
          end do
        end do
 
      end if
 
      call accel_finish()
 
      pop_sub(hamiltonian_elec_update.build_phase)
    end subroutine build_phase
 
  end subroutine hamiltonian_elec_update
 
 
  !----------------------------------------------------------------
  ! TODO: See Issue #1064
  subroutine hamiltonian_elec_update_pot(this, mesh, accumulate)
    type(hamiltonian_elec_t), intent(inout) :: this
    class(mesh_t),            intent(in)    :: mesh
    logical, optional,        intent(in)    :: accumulate
 
    integer :: ispin, ip
 
    push_sub(hamiltonian_elec_update_pot)
 
    ! By default we nullify first the result
    if (.not. optional_default(accumulate, .false.)) then
      !$omp parallel private(ip, ispin)
      do ispin = 1, this%d%nspin
        !$omp do simd schedule(static)
        do ip = 1, mesh%np
          this%hm_base%potential(ip, ispin) = m_zero
        end do
      end do
      !$omp end parallel
    end if
 
    !$omp parallel private(ip, ispin)
    do ispin = 1, this%d%nspin
      if (ispin <= 2) then
        !$omp do simd schedule(static)
        ! this%vhxc(ip, ispin) is added after the calculation of ZORA potential
        do ip = 1, mesh%np
          this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%ep%vpsl(ip) + this%v_ext_pot(ip)
        end do
 
        if (this%pcm%run_pcm) then
          if (this%pcm%solute) then
            !$omp do simd schedule(static)
            do ip = 1, mesh%np
              this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
                this%pcm%v_e_rs(ip) + this%pcm%v_n_rs(ip)
            end do
            !$omp end do simd nowait
          end if
          if (this%pcm%localf) then
            !$omp do simd schedule(static)
            do ip = 1, mesh%np
              this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
                this%pcm%v_ext_rs(ip)
            end do
            !$omp end do simd nowait
          end if
        end if
 
        if (this%abs_boundaries%abtype == imaginary_absorbing) then
          !$omp do simd schedule(static)
          do ip = 1, mesh%np
            this%hm_base%Impotential(ip, ispin) = this%hm_base%Impotential(ip, ispin) + this%abs_boundaries%mf(ip)
          end do
          !$omp end do simd nowait
        end if
      end if
    end do
    !$omp end parallel
 
    ! scalar relativistic ZORA contribution
    ! \boldsymbol{p} \frac{c^2}{2c^2 - V} \boldsymbol{p} \Phi^\mathrm{ZORA}
    if (this%ep%reltype == scalar_relativistic_zora .or. this%ep%reltype == fully_relativistic_zora) then
      call this%zora%update(this%der, this%hm_base%potential)
    end if
 
    !$omp parallel private(ip, ispin)
    do ispin = 1, this%d%nspin
      ! Adding Zeeman potential to hm_base%potential
      if (allocated(this%hm_base%zeeman_pot)) then
        !$omp do simd schedule(static)
        do ip = 1, mesh%np
          this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%hm_base%zeeman_pot(ip, ispin)
        end do
        !$omp end do simd nowait
      end if
 
      ! Adding Quadrupole potential from static E-field (test)
      if (this%mxll%test_equad) then
        !$omp do simd schedule(static)
        do ip = 1, mesh%np
          this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%mxll%e_quadrupole_pot(ip)
        end do
        !$omp end do simd
      end if
    end do
    !$omp end parallel
 
 
    ! Add the Hartree and KS potential
    call this%ks_pot%add_vhxc(this%hm_base%potential)
 
    call this%hm_base%accel_copy_pot(mesh)
 
    pop_sub(hamiltonian_elec_update_pot)
  end subroutine hamiltonian_elec_update_pot
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
    type(hamiltonian_elec_t), intent(inout) :: this
    type(namespace_t),        intent(in)    :: namespace
    class(electron_space_t),  intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(ions_t),     target, intent(inout) :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(states_elec_t),      intent(inout) :: st
    real(real64),   optional, intent(in)    :: time
 
    push_sub(hamiltonian_elec_epot_generate)
 
    this%ions => ions
    call epot_generate(this%ep, namespace, gr, this%ions, this%d)
 
    ! Interation terms are treated below
 
    ! First we add the static electric field
    if (allocated(this%ep%e_field) .and. space%periodic_dim < space%dim) then
      call lalg_axpy(gr%np, m_one, this%v_static, this%ep%vpsl)
    end if
 
    ! Here we need to pass this again, else test are failing.
    ! This is not a real problem, as the multisystem framework will indeed to this anyway
    this%v_ie_loc%atoms_dist => ions%atoms_dist
    this%v_ie_loc%atom => ions%atom
    call this%v_ie_loc%calculate()
 
    ! At the moment we need to add this to ep%vpsl, to keep the behavior of the code
    call lalg_axpy(gr%np, m_one, this%v_ie_loc%potential(:,1), this%ep%vpsl)
 
    ! Here we need to pass this again, else test are failing.
    ! This is not a real problem, as the multisystem framework will indeed to this anyway
    if (this%ep%nlcc) then
      this%nlcc%atoms_dist => ions%atoms_dist
      this%nlcc%atom => ions%atom
      call this%nlcc%calculate()
      call lalg_copy(gr%np, this%nlcc%density(:,1), st%rho_core)
    end if
 
    call this%vnl%build(space, gr, this%ep)
    call hamiltonian_elec_update(this, gr, namespace, space, ext_partners, time)
 
    ! Check if projectors are still compatible with apply_packed on GPU
    if (this%is_applied_packed .and. accel_is_enabled()) then
      if (this%ep%non_local .and. .not. this%vnl%apply_projector_matrices) then
        if (accel_allow_cpu_only()) then
          call messages_write('Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.')
          call messages_warning(namespace=namespace)
        else
          call messages_write('Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.',&
            new_line = .true.)
          call messages_write('Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
          call messages_fatal(namespace=namespace)
        end if
      end if
 
    end if
 
    if (this%pcm%run_pcm) then
      if (this%pcm%solute) then
        call pcm_calc_pot_rs(this%pcm, gr, this%psolver, ions = ions)
      end if
 
      !  Interpolation is needed, hence gr%np_part -> 1:gr%np
      if (this%pcm%localf .and. allocated(this%v_static)) then
        call pcm_calc_pot_rs(this%pcm, gr, this%psolver, v_ext = this%ep%v_ext(1:gr%np_part))
      end if
 
    end if
 
    call lda_u_update_basis(this%lda_u, space, gr, ions, st, this%psolver, namespace, this%kpoints, &
      this%phase%is_allocated())
 
    pop_sub(hamiltonian_elec_epot_generate)
  end subroutine hamiltonian_elec_epot_generate
 
  ! -----------------------------------------------------------------
 
  real(real64) function hamiltonian_elec_get_time(this) result(time)
    type(hamiltonian_elec_t),   intent(inout) :: this
 
    time = this%current_time
  end function hamiltonian_elec_get_time
 
  ! -----------------------------------------------------------------
 
  pure logical function hamiltonian_elec_apply_packed(this) result(apply)
    class(hamiltonian_elec_t),   intent(in) :: this
 
    apply = this%is_applied_packed
 
  end function hamiltonian_elec_apply_packed
 
 
  ! -----------------------------------------------------------------
  subroutine zhamiltonian_elec_apply_atom (hm, namespace, space, latt, species, pos, ia, mesh, psi, vpsi)
    type(hamiltonian_elec_t), intent(in)  :: hm
    type(namespace_t),        intent(in)  :: namespace
    class(space_t),           intent(in)  :: space
    type(lattice_vectors_t),  intent(in)  :: latt
    class(species_t),         intent(in)  :: species
    real(real64),             intent(in)  :: pos(1:space%dim)
    integer,                  intent(in)  :: ia
    class(mesh_t),            intent(in)  :: mesh
    complex(real64),          intent(in)  :: psi(:,:)
    complex(real64),          intent(out) :: vpsi(:,:)
 
    integer :: idim
    real(real64), allocatable :: vlocal(:)
    push_sub(zhamiltonian_elec_apply_atom)
 
    safe_allocate(vlocal(1:mesh%np_part))
    vlocal = m_zero
    call epot_local_potential(hm%ep, namespace, space, latt, mesh, species, pos, ia, vlocal)
 
    do idim = 1, hm%d%dim
      vpsi(1:mesh%np, idim) = vlocal(1:mesh%np) * psi(1:mesh%np, idim)
    end do
 
    safe_deallocate_a(vlocal)
    pop_sub(zhamiltonian_elec_apply_atom)
  end subroutine zhamiltonian_elec_apply_atom
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_update_with_ext_pot(this, mesh, space, ext_partners, time, mu)
    type(hamiltonian_elec_t), intent(inout) :: this
    class(space_t),           intent(in)    :: space
    class(mesh_t),            intent(in)    :: mesh
    type(partner_list_t),     intent(in)    :: ext_partners
    real(real64),             intent(in)    :: time(1:2)
    real(real64),             intent(in)    :: mu(1:2)
 
    integer :: ispin, ip, idir, iatom, ilaser, itime
    real(real64) :: aa(space%dim), time_
    real(real64), allocatable :: vp(:,:)
    real(real64), allocatable :: velectric(:)
    type(lasers_t), pointer :: lasers
    type(gauge_field_t), pointer :: gfield
 
    push_sub(hamiltonian_elec_update_with_ext_pot)
    call profiling_in("HAMILTONIAN_ELEC_UPDATE_EXT_POT")
 
    this%current_time = m_zero
    this%current_time = time(1)
 
    ! set everything to zero
    call this%hm_base%clear(mesh%np)
 
    ! the xc, hartree and external potentials
    call this%hm_base%allocate_field(mesh, field_potential, &
      complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
 
    do itime = 1, 2
      time_ = time(itime)
 
      lasers => list_get_lasers(ext_partners)
      if(associated(lasers)) then
        do ilaser = 1, lasers%no_lasers
          select case (laser_kind(lasers%lasers(ilaser)))
          case (e_field_scalar_potential, e_field_electric)
            safe_allocate(velectric(1:mesh%np))
            do ispin = 1, this%d%spin_channels
              velectric = m_zero
              call laser_potential(lasers%lasers(ilaser), mesh,  velectric, time_)
              !$omp parallel do simd schedule(static)
              do ip = 1, mesh%np
                this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + mu(itime) * velectric(ip)
              end do
            end do
            safe_deallocate_a(velectric)
          case (e_field_magnetic)
            call this%hm_base%allocate_field(mesh, field_vector_potential + field_uniform_magnetic_field, .false.)
            ! get the vector potential
            safe_allocate(vp(1:mesh%np, 1:space%dim))
            vp(1:mesh%np, 1:space%dim) = m_zero
            call laser_vector_potential(lasers%lasers(ilaser), mesh, vp, time_)
            do idir = 1, space%dim
              !$omp parallel do schedule(static)
              do ip = 1, mesh%np
                this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) &
                  - mu(itime) * vp(ip, idir)/p_c
              end do
            end do
            ! and the magnetic field
            call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
            safe_deallocate_a(vp)
          case (e_field_vector_potential)
            call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
            ! get the uniform vector potential associated with a magnetic field
            aa = m_zero
            call laser_field(lasers%lasers(ilaser), aa, time_)
            this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) &
              - mu(itime) * aa/p_c
          end select
        end do
      end if
 
      ! the gauge field
      gfield => list_get_gauge_field(ext_partners)
      if (associated(gfield)) then
        call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
        call gauge_field_get_vec_pot(gfield, aa)
        this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
      end if
 
      ! the electric field for a periodic system through the gauge field
      ! TODO: The condition is wrong here: the e_field should be in non-periodic dims as E field
      ! and as a gauge field in the periodic dim, unless we use a Bery phase, in which, we do not use it
      ! this way. But this is unrelated to the gauge field
      if (allocated(this%ep%e_field) .and. associated(gfield)) then
        this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
          this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
      end if
 
    end do
 
    ! the vector potential of a static magnetic field
    if (allocated(this%ep%a_static)) then
      call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
      !ep%a_static contains 1/c A(r)
      !$omp parallel do schedule(static) private(idir)
      do ip = 1, mesh%np
        do idir = 1, space%dim
          this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
        end do
      end do
    end if
 
    !The electric field is added to the KS potential
    call this%hm_base%accel_copy_pot(mesh)
 
    ! and the static magnetic field
    if (allocated(this%ep%b_field)) then
      call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
      do idir = 1, 3
        this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
      end do
    end if
 
    call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
 
    call hamiltonian_elec_update_pot(this, mesh)
 
    call build_phase()
 
    call profiling_out("HAMILTONIAN_ELEC_UPDATE_EXT_POT")
    pop_sub(hamiltonian_elec_update_with_ext_pot)
 
  contains
 
    subroutine build_phase()
      integer :: ik, imat, nmat, max_npoints, offset, iphase, nphase
 
      push_sub(hamiltonian_elec_update_with_ext_pot.build_phase)
 
      if ((.not. this%kpoints%gamma_only()) .or. allocated(this%hm_base%uniform_vector_potential)) then
 
        call profiling_in('UPDATE_PHASES')
        ! now regenerate the phases for the pseudopotentials
        do iatom = 1, this%ep%natoms
          call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints,  &
            vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
        end do
 
        call profiling_out('UPDATE_PHASES')
      end if
 
      if (allocated(this%hm_base%uniform_vector_potential)) then
        call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
      end if
 
      max_npoints = this%vnl%max_npoints
      nmat = this%vnl%nprojector_matrices
 
 
      if (this%phase%is_allocated() .and. allocated(this%vnl%projector_matrices)) then
 
        nphase = 1
        if (this%der%boundaries%spiralBC) nphase = 3
 
        if (.not. allocated(this%vnl%projector_phases)) then
          safe_allocate(this%vnl%projector_phases(1:max_npoints, nphase, nmat, this%d%kpt%start:this%d%kpt%end))
          if (accel_is_enabled()) then
            call accel_create_buffer(this%vnl%buff_projector_phases, accel_mem_read_only, &
              type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
          end if
        end if
 
        offset = 0
        do ik = this%d%kpt%start, this%d%kpt%end
          do imat = 1, this%vnl%nprojector_matrices
            iatom = this%vnl%projector_to_atom(imat)
            do iphase = 1, nphase
              !$omp parallel do schedule(static)
              do ip = 1, this%vnl%projector_matrices(imat)%npoints
                this%vnl%projector_phases(ip, imat, iphase, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
              end do
 
              if (accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0) then
                call accel_write_buffer(this%vnl%buff_projector_phases, &
                  this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
                  offset = offset)
              end if
              offset = offset + this%vnl%projector_matrices(imat)%npoints
            end do
          end do
        end do
 
      end if
 
      pop_sub(hamiltonian_elec_update_with_ext_pot.build_phase)
    end subroutine build_phase
 
  end subroutine hamiltonian_elec_update_with_ext_pot
 
  logical function hamiltonian_elec_needs_current(hm, states_are_real)
    type(hamiltonian_elec_t), intent(in) :: hm
    logical,                  intent(in) :: states_are_real
 
    hamiltonian_elec_needs_current = .false.
 
    if (hm%self_induced_magnetic) then
      if (.not. states_are_real) then
        hamiltonian_elec_needs_current = .true.
      else
        message(1) = 'No current density for real states since it is identically zero.'
        call messages_warning(1)
      end if
    end if
 
  end function hamiltonian_elec_needs_current
 
  ! ---------------------------------------------------------
  subroutine zhamiltonian_elec_apply_all(hm, namespace, mesh, st, hst)
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(namespace_t),        intent(in)    :: namespace
    class(mesh_t),            intent(in)    :: mesh
    type(states_elec_t),      intent(inout) :: st
    type(states_elec_t),      intent(inout) :: hst
 
    integer :: ik, ib, ist
    complex(real64), allocatable :: psi(:, :)
    complex(real64), allocatable :: psiall(:, :, :, :)
 
    push_sub(zhamiltonian_elec_apply_all)
 
    do ik = st%d%kpt%start, st%d%kpt%end
      do ib = st%group%block_start, st%group%block_end
        call zhamiltonian_elec_apply_batch(hm, namespace, mesh, st%group%psib(ib, ik), hst%group%psib(ib, ik))
      end do
    end do
 
    if (oct_exchange_enabled(hm%oct_exchange)) then
 
      safe_allocate(psiall(mesh%np_part, 1:hst%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
 
      call states_elec_get_state(st, mesh, psiall)
 
      call oct_exchange_prepare(hm%oct_exchange, mesh, psiall, hm%xc, hm%psolver, namespace)
 
      safe_deallocate_a(psiall)
 
      safe_allocate(psi(mesh%np_part, 1:hst%d%dim))
 
      do ik = 1, st%nik
        do ist = 1, st%nst
          call states_elec_get_state(hst, mesh, ist, ik, psi)
          call oct_exchange_operator(hm%oct_exchange, namespace, mesh, psi, ist, ik)
          call states_elec_set_state(hst, mesh, ist, ik, psi)
        end do
      end do
 
      safe_deallocate_a(psi)
 
    end if
 
    pop_sub(zhamiltonian_elec_apply_all)
  end subroutine zhamiltonian_elec_apply_all
 
 
  ! ---------------------------------------------------------
 
  subroutine magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
    type(hamiltonian_elec_t),    intent(in)    :: hm
    type(namespace_t),           intent(in)    :: namespace
    class(mesh_t),               intent(in)    :: mesh
    complex(real64), contiguous, intent(inout) :: psi(:,:)
    complex(real64), contiguous, intent(out)   :: hpsi(:,:)
    integer,                     intent(in)    :: ik
    real(real64),                intent(in)    :: vmagnus(:, :, :)
    logical, optional,           intent(in)    :: set_phase
 
    complex(real64), allocatable :: auxpsi(:, :), aux2psi(:, :)
    integer :: idim, ispin
 
    push_sub(magnus)
 
    ! We will assume, for the moment, no spinors.
    if (hm%d%dim /= 1) then
      call messages_not_implemented("Magnus with spinors", namespace=namespace)
    end if
 
    safe_allocate( auxpsi(1:mesh%np_part, 1:hm%d%dim))
    safe_allocate(aux2psi(1:mesh%np,      1:hm%d%dim))
 
    ispin = hm%d%get_spin_index(ik)
 
    ! Compute (T + Vnl)|psi> and store it
    call zhamiltonian_elec_apply_single(hm, namespace, mesh, psi, auxpsi, 1, ik, &
      terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
 
    ! H|psi>  =  (T + Vnl)|psi> + Vpsl|psi> + Vmagnus(t2)|psi> + Vborders
    do idim = 1, hm%d%dim
      call lalg_copy(mesh%np, auxpsi(:, idim), hpsi(:, idim))
      hpsi(1:mesh%np, idim) = hpsi(1:mesh%np, idim) + hm%ep%Vpsl(1:mesh%np)*psi(1:mesh%np,idim)
      call vborders(mesh, hm, psi(:, idim), hpsi(:, idim))
    end do
    hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + vmagnus(1:mesh%np, ispin, 2)*psi(1:mesh%np, 1)
 
    ! Add first term of the commutator:  - i Vmagnus(t1) (T + Vnl) |psi>
    hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) - m_zi*vmagnus(1:mesh%np, ispin, 1)*auxpsi(1:mesh%np, 1)
 
    ! Add second term of commutator:  i (T + Vnl) Vmagnus(t1) |psi>
    auxpsi(1:mesh%np, 1) = vmagnus(1:mesh%np, ispin, 1)*psi(1:mesh%np, 1)
    call zhamiltonian_elec_apply_single(hm, namespace, mesh, auxpsi, aux2psi, 1, ik, &
      terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
    hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + m_zi*aux2psi(1:mesh%np, 1)
 
    safe_deallocate_a(auxpsi)
    safe_deallocate_a(aux2psi)
    pop_sub(magnus)
  end subroutine magnus
 
  ! ---------------------------------------------------------
  subroutine vborders (mesh, hm, psi, hpsi)
    class(mesh_t),            intent(in)    :: mesh
    type(hamiltonian_elec_t), intent(in)    :: hm
    complex(real64),          intent(in)    :: psi(:)
    complex(real64),          intent(inout) :: hpsi(:)
 
    integer :: ip
 
    push_sub(vborders)
 
    if (hm%abs_boundaries%abtype == imaginary_absorbing) then
      do ip = 1, mesh%np
        hpsi(ip) = hpsi(ip) + m_zi*hm%abs_boundaries%mf(ip)*psi(ip)
      end do
    end if
 
    pop_sub(vborders)
  end subroutine vborders
 
  ! ---------------------------------------------------------
  logical function hamiltonian_elec_has_kick(hm)
    type(hamiltonian_elec_t), intent(in)  :: hm
 
    push_sub(hamiltonian_elec_has_kick)
 
    hamiltonian_elec_has_kick = (abs(hm%kick%delta_strength) >  m_epsilon)
 
    pop_sub(hamiltonian_elec_has_kick)
  end function hamiltonian_elec_has_kick
 
  !
  subroutine hamiltonian_elec_set_mass(this, namespace, mass)
    class(hamiltonian_elec_t) , intent(inout) :: this
    type(namespace_t),          intent(in)    :: namespace
    real(real64),               intent(in)    :: mass
 
    push_sub(hamiltonian_elec_set_mass)
 
    if (parse_is_defined(namespace, 'ParticleMass')) then
      message(1) = 'Attempting to redefine a non-unit electron mass'
      call messages_fatal(1)
    else
      this%mass = mass
    end if
 
    pop_sub(hamiltonian_elec_set_mass)
  end subroutine hamiltonian_elec_set_mass
 
  !----------------------------------------------------------
  subroutine hamiltonian_elec_copy_and_set_phase(hm, gr, kpt, psib, psib_with_phase)
    type(hamiltonian_elec_t), intent(in)  :: hm
    type(grid_t),             intent(in)  :: gr
    type(distributed_t),      intent(in)  :: kpt
    type(wfs_elec_t),         intent(in)  :: psib
    type(wfs_elec_t),         intent(out) :: psib_with_phase
 
    integer :: k_offset, n_boundary_points
 
    push_sub(hamiltonian_elec_copy_and_set_phase)
 
    call psib%copy_to(psib_with_phase)
    if (hm%phase%is_allocated()) then
      call hm%phase%apply_to(gr, gr%np, conjugate = .false., psib = psib_with_phase, src = psib, async=.true.)
      ! apply phase correction while setting boundary -> memory needs to be
      ! accessed only once
      k_offset = psib%ik - kpt%start
      n_boundary_points = int(gr%np_part - gr%np)
      call boundaries_set(gr%der%boundaries, gr, psib_with_phase, phase_correction = hm%phase%phase_corr(:, psib%ik), &
        buff_phase_corr = hm%phase%buff_phase_corr, offset=k_offset*n_boundary_points, async=.true.)
    else
      call psib%copy_data_to(gr%np, psib_with_phase)
      call boundaries_set(gr%der%boundaries, gr, psib_with_phase)
    end if
 
    call psib_with_phase%do_pack(copy = .true.)
 
    pop_sub(hamiltonian_elec_copy_and_set_phase)
  end subroutine hamiltonian_elec_copy_and_set_phase
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_diagonal (hm, mesh, diag, ik)
    type(hamiltonian_elec_t), intent(in)    :: hm
    class(mesh_t),            intent(in)    :: mesh
    real(real64), contiguous, intent(out)   :: diag(:,:)
    integer,                  intent(in)    :: ik
 
    integer :: idim, ip, ispin
 
    real(real64), allocatable  :: ldiag(:)
 
    push_sub(hamiltonian_elec_diagonal)
 
    safe_allocate(ldiag(1:mesh%np))
 
    diag = m_zero
 
    call derivatives_lapl_diag(hm%der, ldiag)
 
    do idim = 1, hm%d%dim
      diag(1:mesh%np, idim) = -m_half/hm%mass*ldiag(1:mesh%np)
    end do
 
    select case (hm%d%ispin)
 
    case (unpolarized, spin_polarized)
      ispin = hm%d%get_spin_index(ik)
      diag(1:mesh%np, 1) = diag(1:mesh%np, 1) + hm%ep%vpsl(1:mesh%np)
 
    case (spinors)
      do ip = 1, mesh%np
        diag(ip, 1) = diag(ip, 1) + hm%ep%vpsl(ip)
        diag(ip, 2) = diag(ip, 2) + hm%ep%vpsl(ip)
      end do
 
    end select
 
    call hm%ks_pot%add_vhxc(diag)
 
    pop_sub(hamiltonian_elec_diagonal)
  end subroutine hamiltonian_elec_diagonal
 
 
 
#include "undef.F90"
#include "real.F90"
#include "hamiltonian_elec_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "hamiltonian_elec_inc.F90"
 
end module hamiltonian_elec_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: