hamiltonian_mxll.F90

Go to the documentation of this file.

!! Copyright (C) 2019 R. Jestaedt, H. Appel, F. Bonafe, M. Oliveira, N. Tancogne-Dejean
!!
!! 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_mxll_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use cube_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use energy_mxll_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_abst_oct_m
  use hamiltonian_elec_oct_m
  use, intrinsic :: iso_fortran_env
  use linear_medium_to_em_field_oct_m
  use math_oct_m
  use maxwell_boundary_op_oct_m
  use mesh_cube_parallel_map_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use nl_operator_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use states_elec_dim_oct_m
  use states_elec_oct_m
  use states_mxll_oct_m
 
  implicit none
 
  private
  public ::                                     &
    hamiltonian_mxll_t,                         &
    hamiltonian_mxll_init,                      &
    hamiltonian_mxll_end,                       &
    dhamiltonian_mxll_apply,                    &
    zhamiltonian_mxll_apply,                    &
    dhamiltonian_mxll_magnus_apply,             &
    zhamiltonian_mxll_magnus_apply,             &
    hamiltonian_mxll_apply_batch,               &
    hamiltonian_mxll_span,                      &
    hamiltonian_mxll_adjoint,                   &
    hamiltonian_mxll_not_adjoint,               &
    hamiltonian_mxll_hermitian,                 &
    hamiltonian_mxll_update,                    &
    hamiltonian_mxll_get_time,                  &
    hamiltonian_mxll_apply_packed,              &
    hamiltonian_mxll_apply_simple,              &
    mxll_update_pml_simple,                     &
    mxll_copy_pml_simple
 
  type, extends(hamiltonian_abst_t) :: hamiltonian_mxll_t
    integer                        :: dim
    logical :: adjoint = .false.
 
    real(real64) :: current_time
    logical :: apply_packed
 
    logical :: time_zero
 
    type(nl_operator_t), pointer   :: operators(:)
 
    type(bc_mxll_t)                :: bc
    type(derivatives_t), pointer, private :: der
    type(states_mxll_t), pointer   :: st
 
    integer                        :: rs_sign
 
    logical                        :: propagation_apply = .false.
 
    integer, pointer               :: rs_state_fft_map(:,:,:)
    integer, pointer               :: rs_state_fft_map_inv(:,:)
 
    logical                        :: mx_ma_coupling = .false.
    logical                        :: mx_ma_coupling_apply = .false.
    integer                        :: mx_ma_trans_field_calc_method
    logical                        :: mx_ma_trans_field_calc_corr = .false.
    integer                        :: mx_ma_coupling_points_number
    real(real64),   allocatable           :: mx_ma_coupling_points(:,:)
    integer, allocatable           :: mx_ma_coupling_points_map(:)
    integer                        :: mx_ma_coupling_order
    logical                        :: ma_mx_coupling       = .false.
    logical                        :: ma_mx_coupling_apply = .false.
 
    logical                        :: bc_add_ab_region  = .false.
    logical                        :: bc_zero           = .false.
    logical                        :: bc_constant       = .false.
    logical                        :: bc_mirror_pec     = .false.
    logical                        :: bc_mirror_pmc     = .false.
    logical                        :: bc_plane_waves    = .false.
    logical                        :: bc_medium         = .false.
 
    logical                        :: plane_waves                = .false.
    logical                        :: plane_waves_apply          = .false.
    logical                        :: spatial_constant           = .false.
    logical                        :: spatial_constant_apply     = .false.
    logical                        :: spatial_constant_propagate = .false.
 
    logical                        :: calc_medium_box = .false.
    type(single_medium_box_t), allocatable  :: medium_boxes(:)
    logical                         :: medium_boxes_initialized = .false.
 
    integer                        :: operator
    logical                        :: current_density_ext_flag = .false.
    logical                        :: current_density_from_medium = .false.
 
    type(energy_mxll_t)            :: energy
 
    logical                        :: cpml_hamiltonian = .false.
 
    logical                        :: diamag_current = .false.
    real(real64)                   :: c_factor
    real(real64)                   :: current_factor
 
    type(cube_t)                   :: cube
    type(mesh_cube_parallel_map_t) :: mesh_cube_map
 
  contains
    procedure :: update_span => hamiltonian_mxll_span
    procedure :: dapply => dhamiltonian_mxll_apply
    procedure :: zapply => zhamiltonian_mxll_apply
    procedure :: dmagnus_apply => dhamiltonian_mxll_magnus_apply
    procedure :: zmagnus_apply => zhamiltonian_mxll_magnus_apply
    procedure :: is_hermitian => hamiltonian_mxll_hermitian
  end type hamiltonian_mxll_t
 
  integer, public, parameter ::      &
    FARADAY_AMPERE              = 1, &
    faraday_ampere_medium       = 2, &
    mxll_simple                 = 3
 
contains
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_init(hm, namespace, gr, st)
    type(hamiltonian_mxll_t),                   intent(inout) :: hm
    type(namespace_t),                          intent(in)    :: namespace
    type(grid_t),                       target, intent(inout) :: gr
    type(states_mxll_t),                target, intent(inout) :: st
 
 
    push_sub(hamiltonian_mxll_init)
 
    call profiling_in('HAMILTONIAN_INIT')
 
    hm%dim = st%dim
    hm%st => st
 
    assert(associated(gr%der%lapl))
 
    hm%operators(1:3) => gr%der%grad(1:3) ! cross product for Maxwell calculation needs dimension >= 2
    hm%der => gr%der
    hm%rs_sign = st%rs_sign
 
    hm%mx_ma_coupling_apply = .false.
    hm%mx_ma_coupling  = .false.
    hm%ma_mx_coupling_apply = .false.
    hm%ma_mx_coupling  = .false.
 
    !%Variable MaxwellHamiltonianOperator
    !%Type integer
    !%Default faraday_ampere
    !%Section Maxwell
    !%Description
    !% With this variable the the Maxwell Hamiltonian operator can be selected
    !%Option faraday_ampere 1
    !% The propagation operation in vacuum with Spin 1 matrices without Gauss law condition.
    !%Option faraday_ampere_medium 2
    !% The propagation operation in medium with Spin 1 matrices without Gauss law condition
    !%Option simple 3
    !% A simpler implementation of the Hamiltonian, including PML. It does not use an extra
    !% transformation to (potentially) 6-element vectors, but uses directly the Riemann-Silberstein
    !% vector as a variable.
    !%End
    call parse_variable(namespace, 'MaxwellHamiltonianOperator', faraday_ampere, hm%operator)
 
    if (hm%operator == faraday_ampere_medium) then
      ! TODO: define this operator automatically if there is a linear medium system present
      hm%dim = 2*hm%dim
      hm%calc_medium_box = .true.
    end if
 
    !%Variable ExternalCurrent
    !%Type logical
    !%Default no
    !%Section Maxwell
    !%Description
    !% If an external current density will be used.
    !%End
    call parse_variable(namespace, 'ExternalCurrent', .false., hm%current_density_ext_flag)
 
    hm%plane_waves_apply = .false.
    hm%spatial_constant_apply = .false.
    hm%spatial_constant_propagate = .true. ! only used if spatially constant field is used
 
    hm%propagation_apply = .false.
 
    !%Variable SpeedOfLightFactor
    !%Type float
    !%Default 1.0
    !%Section Maxwell
    !%Description
    !% Fictitous factor to modify the speed of light in vacuum.
    !% Note: the proper way to handle linear media with a certain refractive index
    !% is by the user of the linear_medium system.
    !%End
    call parse_variable(namespace, 'SpeedOfLightFactor', m_one, hm%c_factor)
 
    !%Variable CurrentDensityFactor
    !%Type float
    !%Default 1.0
    !%Section Maxwell
    !%Description
    !% Fictitous factor to modify the current density coming from partner systems.
    !% Note: This factor does not affect the external current density prescribed by the
    !% <tt>UserDefinedMaxwellExternalCurrent</tt> block.
    !%End
    call parse_variable(namespace, 'CurrentDensityFactor', m_one, hm%current_factor)
 
    hm%rs_state_fft_map     => st%rs_state_fft_map
    hm%rs_state_fft_map_inv => st%rs_state_fft_map_inv
 
    call parse_variable(namespace, 'StatesPack', .true., hm%apply_packed)
 
    call parse_variable(namespace, 'TimeZero', .false., hm%time_zero)
    if (hm%time_zero) call messages_experimental('TimeZero')
 
    call hm%update_span(gr%spacing(1:gr%der%dim), m_zero, namespace)
 
    call profiling_out('HAMILTONIAN_INIT')
    pop_sub(hamiltonian_mxll_init)
  end subroutine hamiltonian_mxll_init
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_end(hm)
    type(hamiltonian_mxll_t), intent(inout) :: hm
 
    integer :: il
 
    push_sub(hamiltonian_mxll_end)
 
    call profiling_in("HAMILTONIAN_MXLL_END")
 
    nullify(hm%operators)
 
    call bc_mxll_end(hm%bc)
 
    if (allocated(hm%medium_boxes)) then
      do il = 1, size(hm%medium_boxes)
        call single_medium_box_end(hm%medium_boxes(il))
      end do
    end if
 
    call profiling_out("HAMILTONIAN_MXLL_END")
 
    pop_sub(hamiltonian_mxll_end)
  end subroutine hamiltonian_mxll_end
 
 
  ! ---------------------------------------------------------
  logical function hamiltonian_mxll_hermitian(hm)
    class(hamiltonian_mxll_t), intent(in) :: hm
 
    push_sub(hamiltonian_mxll_hermitian)
 
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      ! With PML, the Hamiltonian is not purely Hermitian
      hamiltonian_mxll_hermitian = .false.
    else
      hamiltonian_mxll_hermitian = .true.
    end if
 
    pop_sub(hamiltonian_mxll_hermitian)
  end function hamiltonian_mxll_hermitian
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_span(hm, delta, emin, namespace)
    class(hamiltonian_mxll_t), intent(inout) :: hm
    real(real64),              intent(in)    :: delta(:)
    real(real64),              intent(in)    :: emin
    type(namespace_t),         intent(in)    :: namespace
 
    integer :: i
    real(real64) :: emax, fd_factor
    real(real64), parameter :: fd_factors(4) = [1.0_real64, 1.3723_real64, 1.5861_real64, 1.7307_real64]
 
    push_sub(hamiltonian_mxll_span)
 
    ! the discretization of the gradient operator with finite differences of different order
    ! gives different prefactors, they can be obtained by Fourier transform of the stencil
    ! and taking the maximum of the resulting expressions
    if (hm%der%stencil_type == der_star .and. hm%der%order <= 4) then
      fd_factor = fd_factors(hm%der%order)
    else
      ! if we use a different stencil, use pi as an upper bound from the continuous solution
      fd_factor = m_pi
    end if
 
    emax = m_zero
    do i = 1, size(delta)
      emax = emax + m_one / delta(i)**2
    end do
    emax =  p_c * fd_factor * sqrt(emax/m_three)
 
    hm%spectral_middle_point = m_zero
    hm%spectral_half_span    = emax
 
    pop_sub(hamiltonian_mxll_span)
  end subroutine hamiltonian_mxll_span
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_adjoint(hm)
    type(hamiltonian_mxll_t), intent(inout) :: hm
 
    push_sub(hamiltonian_mxll_adjoint)
 
    if (.not. hm%adjoint) then
      hm%adjoint = .true.
    end if
 
    pop_sub(hamiltonian_mxll_adjoint)
  end subroutine hamiltonian_mxll_adjoint
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_not_adjoint(hm)
    type(hamiltonian_mxll_t), intent(inout) :: hm
 
    push_sub(hamiltonian_mxll_not_adjoint)
 
    if (hm%adjoint) then
      hm%adjoint = .false.
    end if
 
    pop_sub(hamiltonian_mxll_not_adjoint)
  end subroutine hamiltonian_mxll_not_adjoint
 
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_update(this, time)
    type(hamiltonian_mxll_t), intent(inout) :: this
    real(real64),   optional, intent(in)    :: time
 
    push_sub(hamiltonian_mxll_update)
 
    this%current_time = m_zero
    if (present(time)) this%current_time = time
 
    pop_sub(hamiltonian_mxll_update)
  end subroutine hamiltonian_mxll_update
 
  ! -----------------------------------------------------------------
 
  real(real64) function hamiltonian_mxll_get_time(this) result(time)
    type(hamiltonian_mxll_t),   intent(inout) :: this
 
    time = this%current_time
 
  end function hamiltonian_mxll_get_time
 
  ! -----------------------------------------------------------------
 
  logical pure function hamiltonian_mxll_apply_packed(this, mesh) result(apply)
    type(hamiltonian_mxll_t),   intent(in) :: this
    class(mesh_t),              intent(in) :: mesh
 
    apply = this%apply_packed
    if (mesh%use_curvilinear) apply = .false.
 
  end function hamiltonian_mxll_apply_packed
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_apply_batch(hm, namespace, der, psib, hpsib, time, terms, set_bc)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(namespace_t),        intent(in)    :: namespace
    type(derivatives_t),      intent(in)    :: der
    type(batch_t), target,    intent(inout) :: psib
    type(batch_t), target,    intent(inout) :: hpsib
    real(real64), optional,   intent(in)    :: time
    integer, optional,        intent(in)    :: terms
    logical, optional,        intent(in)    :: set_bc
 
    type(batch_t) :: gradb(der%dim)
    integer :: idir, ifield, field_dir, pml_dir, rs_sign
    integer :: ip, ip_in, il
    real(real64) :: pml_a, pml_b, pml_c
    complex(real64) :: pml_g, grad
    integer, parameter :: field_dirs(3, 2) = reshape([2, 3, 1, 3, 1, 2], [3, 2])
    real(real64) :: cc, aux_ep(3), aux_mu(3), sigma_e, sigma_m, ff_real(3), ff_imag(3), coeff_real, coeff_imag
    complex(real64) :: ff_plus(3), ff_minus(3), hpsi(6)
    integer, parameter :: sign_medium(6) = [1, 1, 1, -1, -1, -1]
    logical :: with_medium
    real(real64) :: p_c_
 
    push_sub(hamiltonian_mxll_apply_batch)
    call profiling_in("H_MXLL_APPLY_BATCH")
    assert(psib%status() == hpsib%status())
 
    assert(psib%nst == hpsib%nst)
    assert(hm%st%dim == 3)
    p_c_ = p_c * hm%c_factor
 
    !Not implemented at the moment
    assert(.not. present(terms))
    with_medium = hm%operator == faraday_ampere_medium
 
    if (present(time)) then
      if (abs(time - hm%current_time) > 1e-10_real64) then
        write(message(1),'(a)') 'hamiltonian_apply_batch time assertion failed.'
        write(message(2),'(a,f12.6,a,f12.6)') 'time = ', time, '; hm%current_time = ', hm%current_time
        call messages_fatal(2, namespace=namespace)
      end if
    end if
 
    if (hm%operator == mxll_simple) then
      call hamiltonian_mxll_apply_simple(hm, namespace, der%mesh, psib, hpsib, terms, set_bc)
      call profiling_out("H_MXLL_APPLY_BATCH")
      pop_sub(hamiltonian_mxll_apply_batch)
      return
    end if
 
    call zderivatives_batch_grad(der, psib, gradb, set_bc=set_bc)
 
    if (hm%cpml_hamiltonian) then
      call apply_pml_boundary(gradb)
    end if
 
    call zderivatives_batch_curl_from_gradient(der, hpsib, gradb)
 
    call scale_after_curl()
 
    if (hm%bc_constant .and. .not. with_medium) then
      call apply_constant_boundary()
    end if
 
    if (with_medium) then
      do idir = 1, 3
        if ((hm%bc%bc_type(idir) == mxll_bc_medium)) then
          call apply_medium_box(hm%bc%medium(idir))
        end if
      end do
 
      if (hm%calc_medium_box) then
        do il = 1, size(hm%medium_boxes)
          call apply_medium_box(hm%medium_boxes(il))
        end do
      end if
    end if
 
    do idir = 1, der%dim
      call gradb(idir)%end()
    end do
 
    call profiling_out("H_MXLL_APPLY_BATCH")
    pop_sub(hamiltonian_mxll_apply_batch)
 
  contains
    subroutine apply_pml_boundary(gradb)
      type(batch_t) :: gradb(:)
      type(accel_kernel_t), save :: ker_pml
      integer :: wgsize
      push_sub(hamiltonian_mxll_apply_batch.apply_pml_boundary)
      call profiling_in("APPLY_PML_BOUNDARY")
 
      if (with_medium) then
        rs_sign = 1
      else
        rs_sign = hm%rs_sign
      end if
      do idir = 1, der%dim
        call batch_scal(der%mesh%np, rs_sign * p_c_, gradb(idir))
      end do
 
      do pml_dir = 1, hm%st%dim
        if (hm%bc%bc_ab_type(pml_dir) == mxll_ab_cpml) then
          select case (gradb(pml_dir)%status())
          case (batch_not_packed)
            do ifield = 1, 2
              field_dir = field_dirs(pml_dir, ifield)
              !$omp parallel do private(ip, pml_c, pml_a, pml_b, pml_g, grad)
              do ip_in = 1, hm%bc%pml%points_number
                ip = hm%bc%pml%points_map(ip_in)
                pml_c = hm%bc%pml%c(ip_in, pml_dir)
                pml_a = hm%bc%pml%a(ip_in, pml_dir)
                pml_b = hm%bc%pml%b(ip_in, pml_dir)
                pml_g = hm%bc%pml%conv_plus(ip_in, pml_dir, field_dir)
                grad = gradb(pml_dir)%zff_linear(ip, field_dir)
                gradb(pml_dir)%zff_linear(ip, field_dir) = pml_c * ((m_one + pml_a)*grad/p_c_ &
                  + rs_sign * pml_b*pml_g)
                if (with_medium) then
                  pml_g = hm%bc%pml%conv_minus(ip_in, pml_dir, field_dir)
                  grad = gradb(pml_dir)%zff_linear(ip, field_dir+3)
                  gradb(pml_dir)%zff_linear(ip, field_dir+3) = pml_c * ((m_one + pml_a)*grad/p_c_ &
                    + rs_sign * pml_b*pml_g)
                end if
              end do
            end do
          case (batch_packed)
            !$omp parallel do private(ip, field_dir, pml_c, pml_a, pml_b, pml_g, grad, ifield)
            do ip_in = 1, hm%bc%pml%points_number
              ip = hm%bc%pml%points_map(ip_in)
              pml_c = hm%bc%pml%c(ip_in, pml_dir)
              pml_a = hm%bc%pml%a(ip_in, pml_dir)
              pml_b = hm%bc%pml%b(ip_in, pml_dir)
              do ifield = 1, 2
                field_dir = field_dirs(pml_dir, ifield)
                pml_g = hm%bc%pml%conv_plus(ip_in, pml_dir, field_dir)
                grad = gradb(pml_dir)%zff_pack(field_dir, ip)
                gradb(pml_dir)%zff_pack(field_dir, ip) = pml_c * ((m_one + pml_a)*grad/p_c_ &
                  + rs_sign * pml_b*pml_g)
                if (with_medium) then
                  pml_g = hm%bc%pml%conv_minus(ip_in, pml_dir, field_dir)
                  grad = gradb(pml_dir)%zff_pack(field_dir+3, ip)
                  gradb(pml_dir)%zff_pack(field_dir+3, ip) = pml_c * ((m_one + pml_a)*grad/p_c_ &
                    + rs_sign * pml_b*pml_g)
                end if
              end do
            end do
          case (batch_device_packed)
            call accel_kernel_start_call(ker_pml, 'pml.cl', 'pml_apply')
 
            if (with_medium) then
              call accel_set_kernel_arg(ker_pml, 0, 1_int32)
            else
              call accel_set_kernel_arg(ker_pml, 0, 0_int32)
            end if
            call accel_set_kernel_arg(ker_pml, 1, hm%bc%pml%points_number)
            call accel_set_kernel_arg(ker_pml, 2, pml_dir-1)
            call accel_set_kernel_arg(ker_pml, 3, p_c_)
            call accel_set_kernel_arg(ker_pml, 4, rs_sign)
            call accel_set_kernel_arg(ker_pml, 5, hm%bc%pml%buff_map)
            call accel_set_kernel_arg(ker_pml, 6, gradb(pml_dir)%ff_device)
            call accel_set_kernel_arg(ker_pml, 7, log2(int(gradb(pml_dir)%pack_size(1), int32)))
            call accel_set_kernel_arg(ker_pml, 8, hm%bc%pml%buff_a)
            call accel_set_kernel_arg(ker_pml, 9, hm%bc%pml%buff_b)
            call accel_set_kernel_arg(ker_pml, 10, hm%bc%pml%buff_c)
            call accel_set_kernel_arg(ker_pml, 11, hm%bc%pml%buff_conv_plus)
            call accel_set_kernel_arg(ker_pml, 12, hm%bc%pml%buff_conv_minus)
 
            wgsize = accel_max_workgroup_size()
            call accel_kernel_run(ker_pml, (/ accel_padded_size(hm%bc%pml%points_number) /), (/ wgsize /))
          end select
        end if
      end do
 
      if (accel_is_enabled()) then
        call accel_finish()
      end if
 
      call profiling_out("APPLY_PML_BOUNDARY")
      pop_sub(hamiltonian_mxll_apply_batch.apply_pml_boundary)
    end subroutine apply_pml_boundary
 
    subroutine scale_after_curl
      push_sub(hamiltonian_mxll_apply_batch.scale_after_curl)
      call profiling_in("SCALE_AFTER_CURL")
      if (.not. hm%cpml_hamiltonian) then
        ! if we do not need pml, scale after the curl because it is cheaper
        if (with_medium) then
          ! in case of a medium, multiply first 3 components with +, others with -
          call batch_scal(der%mesh%np, sign_medium * p_c_, hpsib)
        else
          call batch_scal(der%mesh%np, hm%rs_sign * p_c_, hpsib)
        end if
      else
        ! this is needed for PML computations with medium
        if (with_medium) then
          ! in case of a medium, multiply first 3 components with +, others with -
          call batch_scal(der%mesh%np, sign_medium * m_one, hpsib)
        end if
      end if
      call profiling_out("SCALE_AFTER_CURL")
      pop_sub(hamiltonian_mxll_apply_batch.scale_after_curl)
    end subroutine scale_after_curl
 
    subroutine apply_constant_boundary
      push_sub(hamiltonian_mxll_apply_batch.apply_constant_boundary)
      call profiling_in('APPLY_CONSTANT_BC')
      select case (hpsib%status())
      case (batch_not_packed)
        do field_dir = 1, hm%st%dim
          do ip_in = 1, hm%bc%constant_points_number
            ip = hm%bc%constant_points_map(ip_in)
            hpsib%zff_linear(ip, field_dir) = hm%st%rs_state_const(field_dir)
          end do
        end do
      case (batch_packed)
        do ip_in = 1, hm%bc%constant_points_number
          ip = hm%bc%constant_points_map(ip_in)
          do field_dir = 1, hm%st%dim
            hpsib%zff_pack(field_dir, ip) = hm%st%rs_state_const(field_dir)
          end do
        end do
      case (batch_device_packed)
        call messages_not_implemented("Maxwell constant boundary on GPU", namespace=namespace)
      end select
      call profiling_out('APPLY_CONSTANT_BC')
      pop_sub(hamiltonian_mxll_apply_batch.apply_constant_boundary)
    end subroutine apply_constant_boundary
 
    subroutine apply_medium_box(medium)
      type(single_medium_box_t),  intent(in) :: medium
 
      integer :: ifield
 
      push_sub(hamiltonian_mxll_apply_batch.apply_medium_box)
      call profiling_in("MEDIUM_BOX")
      assert(.not. medium%has_mapping)
      !$omp parallel do private(ip, cc, aux_ep, aux_mu, sigma_e, sigma_m, &
      !$omp ff_plus, ff_minus, hpsi, ff_real, ff_imag, ifield, coeff_real, coeff_imag)
      do ip = 1, medium%points_number
        cc          = medium%c(ip)/p_c
        if (abs(cc) <= m_epsilon) cycle
        aux_ep(1:3) = medium%aux_ep(ip, 1:3)
        aux_mu(1:3) = medium%aux_mu(ip, 1:3)
        sigma_e     = medium%sigma_e(ip)
        sigma_m     = medium%sigma_m(ip)
        select case (hpsib%status())
        case (batch_not_packed)
          ff_plus(1:3)  = psib%zff_linear(ip, 1:3)
          ff_minus(1:3) = psib%zff_linear(ip, 4:6)
          hpsi(1:6) = hpsib%zff_linear(ip, 1:6)
        case (batch_packed)
          ff_plus(1:3)  = psib%zff_pack(1:3, ip)
          ff_minus(1:3) = psib%zff_pack(4:6, ip)
          hpsi(1:6) = hpsib%zff_pack(1:6, ip)
        case (batch_device_packed)
          call messages_not_implemented("Maxwell Medium on GPU", namespace=namespace)
        end select
        ff_real = real(ff_plus+ff_minus, real64)
        ff_imag = aimag(ff_plus-ff_minus)
        aux_ep = dcross_product(aux_ep, ff_real)
        aux_mu = dcross_product(aux_mu, ff_imag)
        do ifield = 1, 3
          coeff_real = - cc * aux_ep(ifield) + sigma_m * ff_imag(ifield)
          coeff_imag = - cc * aux_mu(ifield) - sigma_e * ff_real(ifield)
          hpsi(ifield) = cc * hpsi(ifield) + cmplx(coeff_real, coeff_imag, real64)
          hpsi(ifield+3) = cc * hpsi(ifield+3) + cmplx(-coeff_real, coeff_imag, real64)
        end do
        select case (hpsib%status())
        case (batch_not_packed)
          hpsib%zff_linear(ip, 1:6) = hpsi(1:6)
        case (batch_packed)
          hpsib%zff_pack(1:6, ip) = hpsi(1:6)
        end select
      end do
      call profiling_out("MEDIUM_BOX")
      pop_sub(hamiltonian_mxll_apply_batch.apply_medium_box)
    end subroutine apply_medium_box
 
  end subroutine hamiltonian_mxll_apply_batch
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_mxll_apply_simple(hm, namespace, mesh, psib, hpsib, terms, set_bc)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(namespace_t),        intent(in)    :: namespace
    class(mesh_t),            intent(in)    :: mesh
    type(batch_t), target,    intent(inout) :: psib
    type(batch_t), target,    intent(inout) :: hpsib
    integer, optional,        intent(in)    :: terms
    logical, optional,        intent(in)    :: set_bc
 
    type(batch_t) :: gradb(hm%der%dim)
    integer :: idir
 
    push_sub(hamiltonian_mxll_apply_simple)
    call profiling_in("MXLL_HAMILTONIAN_SIMPLE")
 
    call zderivatives_batch_grad(hm%der, psib, gradb, set_bc=set_bc)
 
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      ! in the pml layer, the gradient is modified
      call mxll_apply_pml_simple(hm, gradb)
    end if
 
    ! compute the curl from the gradient
    call zderivatives_batch_curl_from_gradient(hm%der, hpsib, gradb)
 
    if (hm%calc_medium_box) then
      ! as the speed of light depends on space, hpsib is already scaled correctly
      call mxll_linear_medium_terms_simple(hm, hpsib)
    else
      ! scale rs_state_tmpb (prefactor of curl)
      call batch_scal(mesh%np, p_c*hm%c_factor, hpsib)
    end if
 
    do idir = 1, hm%der%dim
      call gradb(idir)%end()
    end do
 
    call profiling_out("MXLL_HAMILTONIAN_SIMPLE")
    pop_sub(hamiltonian_mxll_apply_simple)
  end subroutine hamiltonian_mxll_apply_simple
 
  ! ---------------------------------------------------------
  subroutine mxll_apply_pml_simple(hm, gradb)
    type(hamiltonian_mxll_t), target, intent(in)    :: hm
    type(batch_t),                    intent(inout) :: gradb(1:hm%st%dim)
 
    integer :: idir, jdir, ip_in, ip, wgsize
    type(accel_kernel_t), save :: ker_pml
 
 
    push_sub(mxll_apply_pml_simple)
    call profiling_in("APPLY_PML_SIMPLE")
    ! update pml values
    ! loop over gradient direction
    do jdir = 1, hm%st%dim
      select case (gradb(1)%status())
      case (batch_not_packed)
        ! loop over space direction
        do idir = 1, hm%st%dim
          if (idir == jdir) cycle
          do ip_in = 1, hm%bc%pml%points_number
            ip = hm%bc%pml%points_map(ip_in)
            ! update gradient to take pml into account when computing the curl
            gradb(jdir)%zff_linear(ip, idir) = gradb(jdir)%zff_linear(ip, idir) * (m_one + hm%bc%pml%c(ip_in, jdir)) + &
              hm%bc%pml%b(ip_in, jdir) * hm%bc%pml%conv_plus_old(ip_in, jdir, idir)
          end do
        end do
      case (batch_packed)
        do ip_in = 1, hm%bc%pml%points_number
          ip = hm%bc%pml%points_map(ip_in)
          ! loop over space direction
          do idir = 1, hm%st%dim
            if (idir == jdir) cycle
            ! update gradient to take pml into account when computing the curl
            gradb(jdir)%zff_pack(idir, ip) = gradb(jdir)%zff_pack(idir, ip) * (m_one + hm%bc%pml%c(ip_in, jdir)) +&
              hm%bc%pml%b(ip_in, jdir) * hm%bc%pml%conv_plus_old(ip_in, jdir, idir)
          end do
        end do
      case (batch_device_packed)
        call accel_kernel_start_call(ker_pml, 'pml.cl', 'pml_apply_new')
 
        call accel_set_kernel_arg(ker_pml, 0, hm%bc%pml%points_number)
        call accel_set_kernel_arg(ker_pml, 1, jdir-1)
        call accel_set_kernel_arg(ker_pml, 2, hm%bc%pml%buff_map)
        call accel_set_kernel_arg(ker_pml, 3, gradb(jdir)%ff_device)
        call accel_set_kernel_arg(ker_pml, 4, log2(int(gradb(jdir)%pack_size(1), int32)))
        call accel_set_kernel_arg(ker_pml, 5, hm%bc%pml%buff_b)
        call accel_set_kernel_arg(ker_pml, 6, hm%bc%pml%buff_c)
        call accel_set_kernel_arg(ker_pml, 7, hm%bc%pml%buff_conv_plus_old)
 
        wgsize = accel_max_workgroup_size()
        call accel_kernel_run(ker_pml, (/ accel_padded_size(hm%bc%pml%points_number) /), (/ wgsize /))
      end select
    end do
    call profiling_out("APPLY_PML_SIMPLE")
    pop_sub(mxll_apply_pml_simple)
  end subroutine mxll_apply_pml_simple
 
  ! ---------------------------------------------------------
  subroutine mxll_update_pml_simple(hm, rs_stateb)
    type(hamiltonian_mxll_t),intent(inout) :: hm
    type(batch_t),           intent(inout) :: rs_stateb
 
    integer :: wgsize, idir, jdir, ip, ip_in
    type(accel_kernel_t), save :: ker_pml_update
    type(batch_t) :: gradb(1:hm%st%dim)
 
    push_sub(mxll_update_pml_simple)
    call profiling_in("UPDATE_PML_SIMPLE")
 
    call zderivatives_batch_grad(hm%der, rs_stateb, gradb)
 
    do jdir = 1, hm%st%dim
      call rs_stateb%check_compatibility_with(gradb(jdir))
      select case (gradb(jdir)%status())
      case (batch_not_packed)
        ! loop over space direction
        do idir = 1, hm%st%dim
          if (idir == jdir) cycle
          do ip_in = 1, hm%bc%pml%points_number
            ip = hm%bc%pml%points_map(ip_in)
            hm%bc%pml%conv_plus(ip_in, jdir, idir) = hm%bc%pml%c(ip_in, jdir) * gradb(jdir)%zff_linear(ip, idir) +&
              hm%bc%pml%b(ip_in, jdir) * hm%bc%pml%conv_plus_old(ip_in, jdir, idir)
          end do
        end do
      case (batch_packed)
        do ip_in = 1, hm%bc%pml%points_number
          ip = hm%bc%pml%points_map(ip_in)
          ! loop over space direction
          do idir = 1, hm%st%dim
            if (idir == jdir) cycle
            hm%bc%pml%conv_plus(ip_in, jdir, idir) = hm%bc%pml%c(ip_in, jdir) * gradb(jdir)%zff_pack(idir, ip) +&
              hm%bc%pml%b(ip_in, jdir) * hm%bc%pml%conv_plus_old(ip_in, jdir, idir)
          end do
        end do
      case (batch_device_packed)
        call accel_kernel_start_call(ker_pml_update, 'pml.cl', 'pml_update_new')
 
        call accel_set_kernel_arg(ker_pml_update, 0, hm%bc%pml%points_number)
        call accel_set_kernel_arg(ker_pml_update, 1, jdir-1)
        call accel_set_kernel_arg(ker_pml_update, 2, hm%bc%pml%buff_map)
        call accel_set_kernel_arg(ker_pml_update, 3, gradb(jdir)%ff_device)
        call accel_set_kernel_arg(ker_pml_update, 4, log2(int(gradb(jdir)%pack_size(1), int32)))
        call accel_set_kernel_arg(ker_pml_update, 5, hm%bc%pml%buff_b)
        call accel_set_kernel_arg(ker_pml_update, 6, hm%bc%pml%buff_c)
        call accel_set_kernel_arg(ker_pml_update, 7, hm%bc%pml%buff_conv_plus)
        call accel_set_kernel_arg(ker_pml_update, 8, hm%bc%pml%buff_conv_plus_old)
 
        wgsize = accel_max_workgroup_size()
        call accel_kernel_run(ker_pml_update, (/ accel_padded_size(hm%bc%pml%points_number) /), (/ wgsize /))
      end select
      call gradb(jdir)%end()
    end do
    call profiling_out("UPDATE_PML_SIMPLE")
    pop_sub(mxll_update_pml_simple)
  end subroutine mxll_update_pml_simple
 
  ! ---------------------------------------------------------
  subroutine mxll_copy_pml_simple(hm, rs_stateb)
    type(hamiltonian_mxll_t),intent(inout) :: hm
    type(batch_t),           intent(inout) :: rs_stateb
 
    integer :: wgsize, idir, jdir, ip, ip_in
    type(accel_kernel_t), save :: ker_pml_copy
 
    push_sub(mxll_copy_pml_simple)
    call profiling_in("COPY_PML_SIMPLE")
 
    select case (rs_stateb%status())
    case (batch_packed, batch_not_packed)
      do jdir = 1, hm%st%dim
        ! loop over space direction
        do idir = 1, hm%st%dim
          if (idir == jdir) cycle
          do ip_in = 1, hm%bc%pml%points_number
            ip = hm%bc%pml%points_map(ip_in)
            hm%bc%pml%conv_plus_old(ip_in, jdir, idir) = hm%bc%pml%conv_plus(ip_in, jdir, idir)
          end do
        end do
      end do
    case (batch_device_packed)
      call accel_kernel_start_call(ker_pml_copy, 'pml.cl', 'pml_copy')
 
      call accel_set_kernel_arg(ker_pml_copy, 0, hm%bc%pml%points_number*9)
      call accel_set_kernel_arg(ker_pml_copy, 1, hm%bc%pml%buff_conv_plus)
      call accel_set_kernel_arg(ker_pml_copy, 2, hm%bc%pml%buff_conv_plus_old)
 
      wgsize = accel_max_workgroup_size()
      call accel_kernel_run(ker_pml_copy, (/ accel_padded_size(hm%bc%pml%points_number*9) /), (/ wgsize /))
    end select
    call profiling_out("COPY_PML_SIMPLE")
    pop_sub(mxll_copy_pml_simple)
  end subroutine mxll_copy_pml_simple
 
  ! ---------------------------------------------------------
  subroutine mxll_linear_medium_terms_simple(hm, rs_stateb)
    type(hamiltonian_mxll_t),intent(in)    :: hm
    type(batch_t),           intent(inout) :: rs_stateb
 
    integer :: il, ip
    logical :: need_to_pack
 
    push_sub(mxll_linear_medium_terms_simple)
    call profiling_in("LINEAR_MEDIUM_SIMPLE")
 
    do il = 1, size(hm%medium_boxes)
      assert(.not. hm%medium_boxes(il)%has_mapping)
      need_to_pack = .false.
      ! copy to the CPU for GPU runs
      if(rs_stateb%status() == batch_device_packed) then
        call rs_stateb%do_unpack(force=.true.)
        need_to_pack = .true.
      end if
      select case (rs_stateb%status())
      case (batch_not_packed)
        do ip = 1, hm%medium_boxes(il)%points_number
          if (abs(hm%medium_boxes(il)%c(ip)) <= m_epsilon) then
            ! Hamiltonian without medium
            rs_stateb%zff_linear(ip, 1:3) = p_c*hm%c_factor*rs_stateb%zff_linear(ip, 1:3)
          else
            ! Hamiltonian with medium terms
            rs_stateb%zff_linear(ip, 1:3) = hm%medium_boxes(il)%c(ip)*(rs_stateb%zff_linear(ip, 1:3) + &
              cmplx( &
              -dcross_product(hm%medium_boxes(il)%aux_ep(ip, 1:3)*m_two, &
              real(rs_stateb%zff_linear(ip, 1:3), real64)), &
              -dcross_product(hm%medium_boxes(il)%aux_mu(ip, 1:3)*m_two, &
              aimag(rs_stateb%zff_linear(ip, 1:3))), real64)) + &
              cmplx(&
              hm%medium_boxes(il)%sigma_m(ip)*aimag(rs_stateb%zff_linear(ip, 1:3)), &
              -hm%medium_boxes(il)%sigma_e(ip)*real(rs_stateb%zff_linear(ip, 1:3), real64), real64)
          end if
        end do
      case (batch_packed)
        do ip = 1, hm%medium_boxes(il)%points_number
          if (abs(hm%medium_boxes(il)%c(ip)) <= m_epsilon) then
            ! Hamiltonian without medium
            rs_stateb%zff_pack(1:3, ip) = p_c*hm%c_factor*rs_stateb%zff_pack(1:3, ip)
          else
            ! Hamiltonian with medium terms
            rs_stateb%zff_pack(1:3, ip) = hm%medium_boxes(il)%c(ip)*(rs_stateb%zff_pack(1:3, ip) + &
              cmplx( &
              -dcross_product(hm%medium_boxes(il)%aux_ep(ip, 1:3)*m_two, &
              real(rs_stateb%zff_pack(1:3, ip), real64)), &
              -dcross_product(hm%medium_boxes(il)%aux_mu(ip, 1:3)*m_two, &
              aimag(rs_stateb%zff_pack(1:3, ip))), real64)) + &
              cmplx(&
              hm%medium_boxes(il)%sigma_m(ip)*aimag(rs_stateb%zff_pack(1:3, ip)), &
              -hm%medium_boxes(il)%sigma_e(ip)*real(rs_stateb%zff_pack(1:3, ip), real64), real64)
          end if
        end do
      end select
      if(need_to_pack) then
        call rs_stateb%do_pack()
      end if
    end do
 
    call profiling_out("LINEAR_MEDIUM_SIMPLE")
    pop_sub(mxll_linear_medium_terms_simple)
  end subroutine mxll_linear_medium_terms_simple
 
  ! --------------------------------------------------------
  subroutine dhamiltonian_mxll_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
    class(hamiltonian_mxll_t), intent(in)    :: hm
    type(namespace_t),         intent(in)    :: namespace
    class(mesh_t),             intent(in)    :: mesh
    class(batch_t),    target, intent(inout) :: psib
    class(batch_t),    target, intent(inout) :: hpsib
    integer,         optional, intent(in)    :: terms
    logical,         optional, intent(in)    :: set_bc
 
    write(message(1),'(a)') 'dhamiltonian_mxll_apply not implemented (states are complex).'
    call messages_fatal(1, namespace=namespace)
 
  end subroutine dhamiltonian_mxll_apply
 
  ! ---------------------------------------------------------
  subroutine zhamiltonian_mxll_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
    class(hamiltonian_mxll_t),   intent(in)    :: hm
    type(namespace_t),           intent(in)    :: namespace
    class(mesh_t),               intent(in)    :: mesh
    class(batch_t),      target, intent(inout) :: psib
    class(batch_t),      target, intent(inout) :: hpsib
    integer,           optional, intent(in)    :: terms
    logical,           optional, intent(in)    :: set_bc
 
    complex(real64), allocatable :: rs_aux_in(:,:), rs_aux_out(:,:)
    integer :: ii
    logical :: on_gpu
 
    push_sub(zhamiltonian_mxll_apply)
 
    call profiling_in('ZHAMILTONIAN_MXLL_APPLY')
 
    on_gpu = psib%status() == batch_device_packed
    if (hm%operator == faraday_ampere_medium .and. on_gpu) then
      ! legacy code, keep for the moment
      safe_allocate(rs_aux_in(1:mesh%np_part, 1:hm%dim))
      safe_allocate(rs_aux_out(1:mesh%np, 1:hm%dim))
      call boundaries_set(hm%der%boundaries, mesh, psib)
      do ii = 1, hm%dim
        call batch_get_state(psib, ii,mesh%np_part, rs_aux_in(:, ii))
      end do
      ! This uses the old non-batch implementation
      call maxwell_hamiltonian_apply_fd(hm, hm%der, rs_aux_in, rs_aux_out)
      do ii = 1, hm%dim
        call batch_set_state(hpsib, ii, mesh%np, rs_aux_out(:, ii))
      end do
      safe_deallocate_a(rs_aux_in)
      safe_deallocate_a(rs_aux_out)
    else
      ! default branch, should be the only one at some point
      call hamiltonian_mxll_apply_batch(hm, namespace, hm%der, psib, hpsib, set_bc=set_bc)
    end if
 
    call profiling_out('ZHAMILTONIAN_MXLL_APPLY')
 
    pop_sub(zhamiltonian_mxll_apply)
  end subroutine zhamiltonian_mxll_apply
 
 
  ! ---------------------------------------------------------
  subroutine maxwell_hamiltonian_apply_fd(hm, der, psi, oppsi)
    type(hamiltonian_mxll_t),    intent(in)    :: hm
    type(derivatives_t),         intent(in)    :: der
    complex(real64), contiguous, intent(inout) :: psi(:,:)
    complex(real64),             intent(inout) :: oppsi(:,:)
 
    real(real64), pointer     :: mx_rho(:,:)
    complex(real64), allocatable :: tmp(:,:)
    complex(real64), pointer     :: kappa_psi(:,:)
    integer            :: np, np_part, ip, ip_in, rs_sign
    real(real64) :: P_c_
 
    push_sub(maxwell_hamiltonian_apply_fd)
 
    call profiling_in('MAXWELL_HAMILTONIAN_APPLY_FD')
 
    np = der%mesh%np
    np_part = der%mesh%np_part
    rs_sign = hm%rs_sign
    p_c_ = p_c * hm%c_factor
 
    select case (hm%operator)
 
      !=================================================================================================
      ! Maxwell Hamiltonian - Hamiltonian operation in vacuum via partial derivatives:
 
    case (faraday_ampere)
      call profiling_in('MXLL_HAM_APPLY_FD_FARADAY_AMP')
 
      safe_allocate(tmp(1:np,1:2))
      oppsi       = m_z0
 
      if (hm%diamag_current) then
        mx_rho    => hm%st%grid_rho
        kappa_psi => hm%st%kappa_psi
      end if
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector component 1 calculation:
      tmp = m_z0
      call zderivatives_partial(der, psi(:,3), tmp(:,1), 2, set_bc = .false.)
      call zderivatives_partial(der, psi(:,2), tmp(:,2), 3, set_bc = .false.)
      tmp = rs_sign * p_c_ * tmp
      call maxwell_pml_hamiltonian(hm, der, psi, 2, 3, tmp(:,1))
      call maxwell_pml_hamiltonian(hm, der, psi, 3, 2, tmp(:,2))
      oppsi(1:np,1) = ( tmp(1:np,1)-tmp(1:np,2))
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector component 2 calculation:
      tmp = m_z0
      call zderivatives_partial(der, psi(:,1), tmp(:,1), 3, set_bc = .false.)
      call zderivatives_partial(der, psi(:,3), tmp(:,2), 1, set_bc = .false.)
      tmp = rs_sign * p_c_ * tmp
      call maxwell_pml_hamiltonian(hm, der, psi, 3, 1, tmp(:,1))
      call maxwell_pml_hamiltonian(hm, der, psi, 1, 3, tmp(:,2))
      oppsi(1:np,2) = ( tmp(1:np,1)-tmp(1:np,2))
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector component 3 calculation:
      tmp = m_z0
      call zderivatives_partial(der, psi(:,2), tmp(:,1), 1, set_bc = .false.)
      call zderivatives_partial(der, psi(:,1), tmp(:,2), 2, set_bc = .false.)
      tmp = rs_sign * p_c_ * tmp
      call maxwell_pml_hamiltonian(hm, der, psi, 1, 2, tmp(:,1))
      call maxwell_pml_hamiltonian(hm, der, psi, 2, 1, tmp(:,2))
      oppsi(1:np,3) = ( tmp(1:np,1)-tmp(1:np,2))
 
      if (hm%bc_constant) then
        do ip_in = 1, hm%bc%constant_points_number
          ip = hm%bc%constant_points_map(ip_in)
          oppsi(ip,:) = hm%st%rs_state_const(:)
        end do
      end if
 
      safe_deallocate_a(tmp)
 
      call profiling_out('MXLL_HAM_APPLY_FD_FARADAY_AMP')
      !=================================================================================================
      ! Maxwell Hamiltonian - Hamiltonian operation in medium via partial derivatives:
 
    case (faraday_ampere_medium)
      call profiling_in('MXLL_HAM_APP_FAR_AMP_MED')
 
      safe_allocate(tmp(1:np,1:4))
      oppsi       = m_z0
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector component 1 calculation:
      tmp = m_z0
      call zderivatives_partial(der, psi(:,3), tmp(:,1), 2, set_bc = .false.)
      call zderivatives_partial(der, psi(:,2), tmp(:,3), 3, set_bc = .false.)
      call zderivatives_partial(der, psi(:,6), tmp(:,2), 2, set_bc = .false.)
      call zderivatives_partial(der, psi(:,5), tmp(:,4), 3, set_bc = .false.)
      tmp = p_c_ * tmp
      call maxwell_pml_hamiltonian_medium(hm, der, psi, 2, 3, tmp(:,1:2))
      call maxwell_pml_hamiltonian_medium(hm, der, psi, 3, 2, tmp(:,3:4))
      oppsi(1:np,1) =  (tmp(1:np,1)-tmp(1:np,3))
      oppsi(1:np,4) = -(tmp(1:np,2)-tmp(1:np,4))
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector component 2 calculation:
      tmp = m_z0
      call zderivatives_partial(der, psi(:,1), tmp(:,1), 3, set_bc = .false.)
      call zderivatives_partial(der, psi(:,3), tmp(:,3), 1, set_bc = .false.)
      call zderivatives_partial(der, psi(:,4), tmp(:,2), 3, set_bc = .false.)
      call zderivatives_partial(der, psi(:,6), tmp(:,4), 1, set_bc = .false.)
      tmp = p_c_ * tmp
      call maxwell_pml_hamiltonian_medium(hm, der, psi, 3, 1, tmp(:,1:2))
      call maxwell_pml_hamiltonian_medium(hm, der, psi, 1, 3, tmp(:,3:4))
      oppsi(1:np,2) =  (tmp(1:np,1)-tmp(1:np,3))
      oppsi(1:np,5) = -(tmp(1:np,2)-tmp(1:np,4))
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector component 3 calculation:
      tmp = m_z0
      call zderivatives_partial(der, psi(:,2), tmp(:,1), 1, set_bc = .false.)
      call zderivatives_partial(der, psi(:,1), tmp(:,3), 2, set_bc = .false.)
      call zderivatives_partial(der, psi(:,5), tmp(:,2), 1, set_bc = .false.)
      call zderivatives_partial(der, psi(:,4), tmp(:,4), 2, set_bc = .false.)
      tmp = p_c_ * tmp
      call maxwell_pml_hamiltonian_medium(hm, der, psi, 1, 2, tmp(:,1:2))
      call maxwell_pml_hamiltonian_medium(hm, der, psi, 2, 1, tmp(:,3:4))
      oppsi(1:np,3) =  (tmp(1:np,1)-tmp(1:np,3))
      oppsi(1:np,6) = -(tmp(1:np,2)-tmp(1:np,4))
 
 
      safe_deallocate_a(tmp)
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector calculation if medium boundaries is set:
      call maxwell_medium_boundaries_calculation(hm, psi, oppsi)
 
      !-----------------------------------------------------------------------------------------------
      ! Riemann-Silberstein vector calculation for medium boxes:
      call maxwell_medium_boxes_calculation(hm, der, psi, oppsi)
 
      call profiling_out('MXLL_HAM_APP_FAR_AMP_MED')
 
    end select
 
    call profiling_out('MAXWELL_HAMILTONIAN_APPLY_FD')
 
    pop_sub(maxwell_hamiltonian_apply_fd)
  end subroutine maxwell_hamiltonian_apply_fd
 
 
  ! ---------------------------------------------------------
  subroutine maxwell_pml_hamiltonian(hm, der, psi, dir1, dir2, tmp)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(derivatives_t),      intent(in)    :: der
    complex(real64),          intent(inout) :: psi(:,:)
    integer,                  intent(in)    :: dir1
    integer,                  intent(in)    :: dir2
    complex(real64),          intent(inout) :: tmp(:)
 
 
    push_sub(maxwell_pml_hamiltonian)
 
    call profiling_in('MAXWELL_PML_HAMILTONIAN')
 
    if ((hm%bc%bc_ab_type(dir1) == mxll_ab_cpml) .and. hm%cpml_hamiltonian) then
      call maxwell_pml_calculation_via_riemann_silberstein(hm, der, psi, dir1, dir2, tmp(:))
    end if
 
    call profiling_out('MAXWELL_PML_HAMILTONIAN')
 
    pop_sub(maxwell_pml_hamiltonian)
  end subroutine maxwell_pml_hamiltonian
 
  ! ---------------------------------------------------------
  subroutine maxwell_pml_hamiltonian_medium(hm, der, psi, dir1, dir2, tmp)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(derivatives_t),      intent(in)    :: der
    complex(real64),          intent(inout) :: psi(:,:)
    integer,                  intent(in)    :: dir1
    integer,                  intent(in)    :: dir2
    complex(real64),          intent(inout) :: tmp(:,:)
 
 
    push_sub(maxwell_pml_hamiltonian_medium)
 
    call profiling_in('MAXWELL_PML_HAMILTONIAN_MEDIUM')
 
    if ((hm%bc%bc_ab_type(dir1) == mxll_ab_cpml) .and. hm%cpml_hamiltonian) then
      call maxwell_pml_calculation_via_riemann_silberstein_medium(hm, der, psi, dir1, dir2, tmp(:,:))
    end if
 
    call profiling_out('MAXWELL_PML_HAMILTONIAN_MEDIUM')
 
    pop_sub(maxwell_pml_hamiltonian_medium)
  end subroutine maxwell_pml_hamiltonian_medium
 
  ! ---------------------------------------------------------
  subroutine maxwell_pml_calculation_via_riemann_silberstein(hm, der, psi, pml_dir, field_dir, pml)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(derivatives_t),      intent(in)    :: der
    integer,                  intent(in)    :: pml_dir
    complex(real64),          intent(inout) :: psi(:,:)
    integer,                  intent(in)    :: field_dir
    complex(real64),          intent(inout) :: pml(:)
 
    integer            :: ip, ip_in, rs_sign
    real(real64)       :: pml_c
    complex(real64), allocatable :: tmp_partial(:)
    complex(real64)    :: pml_a, pml_b, pml_g
 
    push_sub(maxwell_pml_calculation_via_riemann_silberstein)
 
    if (hm%cpml_hamiltonian) then
 
      rs_sign = hm%rs_sign
 
      safe_allocate(tmp_partial(1:der%mesh%np_part))
 
      call zderivatives_partial(der, psi(:,field_dir), tmp_partial(:), pml_dir, set_bc = .false.)
      do ip_in = 1, hm%bc%pml%points_number
        ip       = hm%bc%pml%points_map(ip_in)
        pml_c = hm%bc%pml%c(ip_in, pml_dir)
        pml_a = hm%bc%pml%a(ip_in, pml_dir)
        pml_b = hm%bc%pml%b(ip_in, pml_dir)
        pml_g = hm%bc%pml%conv_plus(ip_in, pml_dir, field_dir)
        pml(ip)  = rs_sign * pml_c * ( tmp_partial(ip) &
          +  real(pml_a, real64)  * real(tmp_partial(ip), real64)  &
          +  m_zi * aimag(pml_a) * aimag(tmp_partial(ip)) &
          +  real(pml_b, real64)  * real(pml_g, real64)  &
          +  m_zi * aimag(pml_b) * aimag(pml_g))
      end do
 
      safe_deallocate_a(tmp_partial)
    end if
 
    pop_sub(maxwell_pml_calculation_via_riemann_silberstein)
  end subroutine maxwell_pml_calculation_via_riemann_silberstein
 
 
  ! ---------------------------------------------------------
  subroutine maxwell_pml_calculation_via_riemann_silberstein_medium(hm, der, psi, pml_dir, field_dir, pml)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(derivatives_t),      intent(in)    :: der
    integer,                  intent(in)    :: pml_dir
    complex(real64),          intent(inout) :: psi(:,:)
    integer,                  intent(in)    :: field_dir
    complex(real64),          intent(inout) :: pml(:,:)
 
    integer            :: ip, ip_in, np
    real(real64)       :: pml_c(3)
    complex(real64), allocatable :: tmp_partial(:,:)
    complex(real64)    :: pml_a(3), pml_b(3), pml_g_p(3), pml_g_m(3)
 
    push_sub(maxwell_pml_calculation_via_riemann_silberstein_medium)
 
    if (hm%cpml_hamiltonian) then
 
      np = der%mesh%np
      safe_allocate(tmp_partial(1:np,1:2))
 
      call zderivatives_partial(der, psi(:,field_dir  ), tmp_partial(:,1), pml_dir, set_bc = .false.)
      call zderivatives_partial(der, psi(:,field_dir+3), tmp_partial(:,2), pml_dir, set_bc = .false.)
      do ip_in = 1, hm%bc%pml%points_number
        ip         = hm%bc%pml%points_map(ip_in)
        pml_c(1:3)   = hm%bc%pml%c(ip_in, 1:3)
        pml_a(1:3)   = hm%bc%pml%a(ip_in, 1:3)
        pml_b(1:3)   = hm%bc%pml%b(ip_in, 1:3)
        pml_g_p(1:3) = hm%bc%pml%conv_plus(ip_in, pml_dir, 1:3)
        pml_g_m(1:3) = hm%bc%pml%conv_minus(ip_in, pml_dir, 1:3)
        pml(ip, 1)   = pml_c(pml_dir) * tmp_partial(ip, 1) &
          + pml_c(pml_dir) * real(pml_a(pml_dir), real64) * real(tmp_partial(ip, 1), real64) &
          + m_zi * pml_c(pml_dir) * aimag(pml_a(pml_dir)) * aimag(tmp_partial(ip, 1)) &
          + pml_c(pml_dir) * real(pml_b(pml_dir), real64) * real(pml_g_p(field_dir), real64) &
          + m_zi * pml_c(pml_dir) * aimag(pml_b(pml_dir)) * aimag(pml_g_p(field_dir))
        pml(ip, 2)   = pml_c(pml_dir) * tmp_partial(ip, 2) &
          + pml_c(pml_dir) * real(pml_a(pml_dir), real64) * real(tmp_partial(ip, 2), real64) &
          + m_zi * pml_c(pml_dir) * aimag(pml_a(pml_dir)) * aimag(tmp_partial(ip, 2)) &
          + pml_c(pml_dir) * real(pml_b(pml_dir), real64) * real(pml_g_m(field_dir), real64) &
          + m_zi * pml_c(pml_dir) * aimag(pml_b(pml_dir)) * aimag(pml_g_m(field_dir))
      end do
 
    end if
 
    safe_deallocate_a(tmp_partial)
 
    pop_sub(maxwell_pml_calculation_via_riemann_silberstein_medium)
  end subroutine maxwell_pml_calculation_via_riemann_silberstein_medium
 
 
  ! ---------------------------------------------------------
  subroutine maxwell_medium_boundaries_calculation(hm, psi, oppsi)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    complex(real64),          intent(in)    :: psi(:,:)
    complex(real64),          intent(inout) :: oppsi(:,:)
 
    integer            :: ip, ip_in, idim
    real(real64)       :: cc, aux_ep(3), aux_mu(3), sigma_e, sigma_m
    complex(real64)    :: ff_plus(3), ff_minus(3)
 
    push_sub(maxwell_medium_boundaries_calculation)
 
    do idim = 1, 3
      if ((hm%bc%bc_type(idim) == mxll_bc_medium)) then
        do ip_in = 1, hm%bc%medium(idim)%points_number
          ip          = hm%bc%medium(idim)%points_map(ip_in)
          cc          = hm%bc%medium(idim)%c(ip_in)/p_c
          aux_ep(:)   = hm%bc%medium(idim)%aux_ep(ip_in, :)
          aux_mu(:)   = hm%bc%medium(idim)%aux_mu(ip_in, :)
          sigma_e     = hm%bc%medium(idim)%sigma_e(ip_in)
          sigma_m     = hm%bc%medium(idim)%sigma_m(ip_in)
          ff_plus(1)  = psi(ip, 1)
          ff_plus(2)  = psi(ip, 2)
          ff_plus(3)  = psi(ip, 3)
          ff_minus(1) = psi(ip, 4)
          ff_minus(2) = psi(ip, 5)
          ff_minus(3) = psi(ip, 6)
          aux_ep      = dcross_product(aux_ep,real(ff_plus+ff_minus, real64) )
          aux_mu      = dcross_product(aux_mu,aimag(ff_plus-ff_minus))
          oppsi(ip, 1) = oppsi(ip, 1)*cc                                         &
            - cc * aux_ep(1) - cc * m_zi * aux_mu(1)                  &
            - m_zi * sigma_e * real(ff_plus(1) + ff_minus(1), real64)          &
            - m_zi * sigma_m * m_zi * aimag(ff_plus(1) - ff_minus(1))
          oppsi(ip, 4) = oppsi(ip, 4)*cc                                         &
            + cc * aux_ep(1) - cc * m_zi * aux_mu(1)                  &
            - m_zi * sigma_e * real(ff_plus(1) + ff_minus(1), real64)          &
            + m_zi * sigma_m * m_zi * aimag(ff_plus(1) - ff_minus(1))
          oppsi(ip, 2) = oppsi(ip, 2)*cc                                         &
            - cc * aux_ep(2) - cc * m_zi * aux_mu(2)                  &
            - m_zi * sigma_e * real(ff_plus(2) + ff_minus(2), real64)          &
            - m_zi * sigma_m * m_zi * aimag(ff_plus(2) - ff_minus(2))
          oppsi(ip, 5) = oppsi(ip, 5)*cc                                         &
            + cc * aux_ep(2) - cc * m_zi * aux_mu(2)                  &
            - m_zi * sigma_e * real(ff_plus(2) + ff_minus(2), real64)          &
            + m_zi * sigma_m * m_zi * aimag(ff_plus(2) - ff_minus(2))
          oppsi(ip, 3) = oppsi(ip, 3)*cc                                         &
            - cc * aux_ep(3) - cc * m_zi * aux_mu(3)                  &
            - m_zi * sigma_e * real(ff_plus(3) + ff_minus(3), real64)          &
            - m_zi * sigma_m * m_zi * aimag(ff_plus(3) - ff_minus(3))
          oppsi(ip, 6) = oppsi(ip, 6)*cc                                         &
            + cc * aux_ep(3) - cc * m_zi * aux_mu(3)                  &
            - m_zi * sigma_e * real(ff_plus(3) + ff_minus(3), real64)          &
            + m_zi * sigma_m * m_zi * aimag(ff_plus(3) - ff_minus(3))
        end do
      end if
    end do
 
    pop_sub(maxwell_medium_boundaries_calculation)
  end subroutine maxwell_medium_boundaries_calculation
 
 
  ! ---------------------------------------------------------
  ! > Maxwell Hamiltonian including medium boxes
  subroutine maxwell_medium_boxes_calculation(hm, der, psi, oppsi)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(derivatives_t),      intent(in)    :: der
    complex(real64),          intent(in)    :: psi(:,:)
    complex(real64),          intent(inout) :: oppsi(:,:)
 
    integer            :: ip, il
    real(real64)       :: cc, aux_ep(3), aux_mu(3), sigma_e, sigma_m
    complex(real64)    :: ff_plus(3), ff_minus(3)
 
    push_sub(maxwell_medium_boxes_calculation)
 
    if (hm%calc_medium_box) then
      do il = 1, size(hm%medium_boxes)
        assert(.not. hm%medium_boxes(il)%has_mapping)
        do ip = 1, hm%medium_boxes(il)%points_number
          cc           = hm%medium_boxes(il)%c(ip)/p_c
          if (abs(cc) <= m_epsilon) cycle
          aux_ep(1:3)  = hm%medium_boxes(il)%aux_ep(ip, 1:3)
          aux_mu(1:3)  = hm%medium_boxes(il)%aux_mu(ip, 1:3)
          sigma_e      = hm%medium_boxes(il)%sigma_e(ip)
          sigma_m      = hm%medium_boxes(il)%sigma_m(ip)
          ff_plus(1)   = psi(ip, 1)
          ff_plus(2)   = psi(ip, 2)
          ff_plus(3)   = psi(ip, 3)
          ff_minus(1)  = psi(ip, 4)
          ff_minus(2)  = psi(ip, 5)
          ff_minus(3)  = psi(ip, 6)
          aux_ep       = dcross_product(aux_ep, real(ff_plus+ff_minus, real64) )
          aux_mu       = dcross_product(aux_mu, aimag(ff_plus-ff_minus))
          oppsi(ip, 1) = oppsi(ip,1)*cc                                          &
            - cc * aux_ep(1) - cc * m_zi * aux_mu(1)                  &
            - m_zi * sigma_e * real(ff_plus(1) + ff_minus(1), real64)          &
            - m_zi * sigma_m * m_zi * aimag(ff_plus(1) - ff_minus(1))
          oppsi(ip, 4) = oppsi(ip,4)*cc                                          &
            + cc * aux_ep(1) - cc * m_zi * aux_mu(1)                  &
            - m_zi * sigma_e * real(ff_plus(1) + ff_minus(1), real64)          &
            + m_zi * sigma_m * m_zi * aimag(ff_plus(1) - ff_minus(1))
          oppsi(ip, 2) = oppsi(ip,2)*cc                                          &
            - cc * aux_ep(2) - cc * m_zi * aux_mu(2)                  &
            - m_zi * sigma_e * real(ff_plus(2) + ff_minus(2), real64)          &
            - m_zi * sigma_m * m_zi * aimag(ff_plus(2) - ff_minus(2))
          oppsi(ip, 5) = oppsi(ip,5)*cc                                          &
            + cc * aux_ep(2) - cc * m_zi * aux_mu(2)                  &
            - m_zi * sigma_e * real(ff_plus(2) + ff_minus(2), real64)          &
            + m_zi * sigma_m * m_zi * aimag(ff_plus(2) - ff_minus(2))
          oppsi(ip, 3) = oppsi(ip,3)*cc                                          &
            - cc * aux_ep(3) - cc * m_zi * aux_mu(3)                  &
            - m_zi * sigma_e * real(ff_plus(3) + ff_minus(3), real64)          &
            - m_zi * sigma_m * m_zi * aimag(ff_plus(3) - ff_minus(3))
          oppsi(ip, 6) = oppsi(ip,6)*cc                                          &
            + cc * aux_ep(3) - cc * m_zi * aux_mu(3)                  &
            - m_zi * sigma_e * real(ff_plus(3) + ff_minus(3), real64)          &
            + m_zi * sigma_m * m_zi * aimag(ff_plus(3) - ff_minus(3))
        end do
      end do
    end if
 
    pop_sub(maxwell_medium_boxes_calculation)
  end subroutine maxwell_medium_boxes_calculation
 
  ! ---------------------------------------------------------
  subroutine dhamiltonian_mxll_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
    class(hamiltonian_mxll_t),   intent(in)    :: hm
    type(namespace_t),           intent(in)    :: namespace
    class(mesh_t),               intent(in)    :: mesh
    class(batch_t),              intent(inout) :: psib
    class(batch_t),              intent(inout) :: hpsib
    real(real64),                intent(in)    :: vmagnus(:, :, :)
 
    call messages_not_implemented ('dhamiltonian_mxll_magnus_apply', namespace=namespace)
 
  end subroutine dhamiltonian_mxll_magnus_apply
 
  ! ---------------------------------------------------------
  subroutine zhamiltonian_mxll_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
    class(hamiltonian_mxll_t),   intent(in)    :: hm
    type(namespace_t),           intent(in)    :: namespace
    class(mesh_t),               intent(in)    :: mesh
    class(batch_t),              intent(inout) :: psib
    class(batch_t),              intent(inout) :: hpsib
    real(real64),                intent(in)    :: vmagnus(:, :, :)
 
    call messages_not_implemented ('zhamiltonian_mxll_magnus_apply', namespace=namespace)
 
  end subroutine zhamiltonian_mxll_magnus_apply
 
end module hamiltonian_mxll_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: