diff --git a/src/Makefile.am b/src/Makefile.am
index 444d68d7da3d1c1c068a21fca7237163b6f446c8..9b57ede94f95e4a23135e630305d4e3f4cb17b69 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -55,6 +55,7 @@ basic_f_srcs =                                 \
 	basic/debug.F90                        \
 	basic/gdlib.F90                        \
 	basic/global.F90                       \
+	basic/global_h.F90                     \
 	basic/hardware.F90                     \
 	basic/heap.F90                         \
 	basic/iihash.F90                       \
@@ -466,9 +467,10 @@ hamiltonian_f_srcs =                                 \
 	hamiltonian/exchange_operator.F90            \
 	hamiltonian/xc_functional.F90                \
 	hamiltonian/gauge_field.F90                  \
-	hamiltonian/hamiltonian_abst.F90             \
+	hamiltonian/hamiltonian_abst_h.F90           \
 	hamiltonian/hamiltonian_elec_base.F90        \
 	hamiltonian/hamiltonian_elec.F90             \
+	hamiltonian/hamiltonian_elec_h.F90           \
 	hamiltonian/hgh_projector.F90                \
 	hamiltonian/hirshfeld.F90                    \
 	hamiltonian/ion_interaction.F90              \
@@ -551,6 +553,7 @@ multisystem_f_srcs =                                   \
 	multisystem/propagator_verlet.F90              \
 	multisystem/quantity.F90                       \
 	multisystem/system.F90                         \
+	multisystem/system_h.F90                       \
 	multisystem/system_factory_abst.F90
 
 multisystem_srcs = $(multisystem_f_srcs)
@@ -658,6 +661,7 @@ electrons_f_srcs =                                 \
 	electrons/stress.F90                       \
 	electrons/subspace.F90                     \
 	electrons/electrons.F90                    \
+	electrons/electrons_h.F90                  \
 	electrons/v_ks.F90                         \
 	electrons/x_fbe.F90                        \
 	electrons/x_slater.F90                     \
@@ -704,6 +708,7 @@ maxwell_f_srcs =                                 \
 	maxwell/external_densities.F90           \
 	maxwell/external_waves.F90               \
 	maxwell/hamiltonian_mxll.F90             \
+	maxwell/hamiltonian_mxll_h.F90           \
 	maxwell/propagator_mxll.F90              \
 	maxwell/dispersive_medium.F90            \
 	maxwell/linear_medium.F90                \
@@ -722,6 +727,7 @@ scf_f_srcs =                                 \
 	scf/criteria_factory.F90             \
 	scf/density_criterion.F90            \
 	scf/electrons_ground_state.F90       \
+	scf/electrons_ground_state_h.F90     \
 	scf/eigenval_criterion.F90           \
 	scf/energy_criterion.F90             \
 	scf/lcao.F90                         \
@@ -729,7 +735,9 @@ scf_f_srcs =                                 \
 	scf/mix.F90                          \
 	scf/mixing_preconditioner.F90        \
 	scf/rdmft.F90                        \
-	scf/scf.F90                          \
+	scf/scf_interface.F90                \
+	scf/scf_interface_h.F90              \
+	scf/scf_h.F90                        \
 	scf/unocc.F90
 
 scf_srcs = $(scf_f_srcs)
@@ -750,6 +758,7 @@ td_f_srcs =                                 \
 	td/propagator_base.F90              \
 	td/propagator_cn.F90                \
 	td/propagator_elec.F90              \
+	td/propagator_elec_h.F90            \
 	td/propagator_etrs.F90              \
 	td/propagator_expmid.F90            \
 	td/propagator_magnus.F90            \
@@ -757,7 +766,9 @@ td_f_srcs =                                 \
 	td/propagator_rk.F90                \
 	td/spectrum.F90                     \
 	td/td_calc.F90                      \
-	td/td.F90                           \
+	td/td_interface.F90                 \
+	td/td_interface_h.F90               \
+	td/td_h.F90                         \
 	td/td_write.F90                     \
 	td/td_write_low.F90                 \
 	td/propagation_ops_elec.F90
diff --git a/src/basic/CMakeLists.txt b/src/basic/CMakeLists.txt
index 96efa8109495d6eb432b7bd5321d845b121df6f3..fdf90ef09c7d5a330675c6be35322e4185cba258 100644
--- a/src/basic/CMakeLists.txt
+++ b/src/basic/CMakeLists.txt
@@ -17,6 +17,7 @@ target_sources(Octopus_lib PRIVATE
 		gdlib_f.c
 		getopt_f.c
 		global.F90
+		global_h.F90
 		hardware.F90
 		heap.F90
 		iihash.F90
diff --git a/src/basic/global.F90 b/src/basic/global.F90
index 1d610352feeae4e458fd0ea0dab83bbd2d2e6a34..9ffc7fb5fa39a6785998b7a820c650d72aeacfcf 100644
--- a/src/basic/global.F90
+++ b/src/basic/global.F90
@@ -18,197 +18,20 @@
 
 #include "global.h"
 
-module global_oct_m
-  use, intrinsic :: iso_fortran_env
+submodule (global_oct_m) impl
+  use global_oct_m
   use hardware_oct_m
   use loct_oct_m
-  use mpi_oct_m
   use varinfo_oct_m
 #ifdef HAVE_OPENMP
   use omp_lib
 #endif
 
   implicit none
-
-  private
-
-  !> Public types, variables and procedures.
-  public ::               &
-    conf_t,               &
-    global_init,          &
-    global_end,           &
-    init_octopus_globals, &
-    optional_default,     &
-    assert_die,           &
-    not_in_openmp,        &
-    operator(+),          &
-    bitand,               &
-    int32, int64,         &
-    real32, real64,       &
-    i4_to_i8,             &
-    i8_to_i4
-  ! Make these kind variables from kind_oct_m public here so that they are
-  ! available basically everywhere in the code. They still need to be in a
-  ! separate module because they are also needed in some low-level modules.
-
-  integer, public, parameter :: MAX_PATH_LEN=512
-  integer, public, parameter :: MAX_OUTPUT_TYPES=44
-
-  !> @brief Build configuration type
-  type conf_t
-    logical            :: devel_version              !< If true then allow unstable parts of the code
-    logical            :: report_memory
-    character(len=256) :: share = SHARE_DIR          !< Name of the share dir
-    character(len=256) :: git_commit  = GIT_COMMIT   !< hash of latest git commit
-    character(len=50)  :: config_time = BUILD_TIME   !< time octopus was configured
-    character(len=20)  :: version = PACKAGE_VERSION  !< version number
-    character(len=256) :: cc = CC                    !< C compiler
-    character(len=256) :: cxx = CXX                  !< C++ compiler
-    character(len=256) :: fc = FC                    !< Fortran compiler
-    ! Split flag definitions in case they don`t fit in one line, following preprocessing
-    character(len=256) :: cflags = &
-      CFLAGS //&
-      CFLAGS_EXTRA
-    character(len=256) :: cxxflags = &
-      CXXFLAGS //&
-      CXXFLAGS_EXTRA
-    character(len=256) :: fcflags = &
-      FCFLAGS //&
-      FCFLAGS_EXTRA
-    integer            :: target_states_block_size = -1
-  contains
-    procedure :: init => conf_init
-  end type conf_t
-
-  !> Global instance of Octopus configuration
-  type(conf_t),      public :: conf
-
-  real(real64), public, parameter :: R_SMALL = 1e-8_real64
-
-  !> Minimal distance between two distinguishable atoms
-  real(real64), public, parameter :: R_MIN_ATOM_DIST = 1e-3_real64
-
-  !> some mathematical constants
-  real(real64), public, parameter :: M_Pi        = 3.1415926535897932384626433832795029_real64
-  real(real64), public, parameter :: M_E         = 2.7182818284590452353602874713526625_real64
-  real(real64), public, parameter :: M_ZERO      = 0.0_real64
-  real(real64), public, parameter :: M_ONE       = 1.0_real64
-  real(real64), public, parameter :: M_TWO       = 2.0_real64
-  real(real64), public, parameter :: M_THREE     = 3.0_real64
-  real(real64), public, parameter :: M_FOUR      = 4.0_real64
-  real(real64), public, parameter :: M_FIVE      = 5.0_real64
-  real(real64), public, parameter :: M_HALF      = 0.5_real64
-  real(real64), public, parameter :: M_THIRD     = M_ONE/M_THREE
-  real(real64), public, parameter :: M_TWOTHIRD  = M_TWO/M_THREE
-  real(real64), public, parameter :: M_FOURTH    = M_ONE/M_FOUR
-  complex(real64), public, parameter :: M_z0        = (0.0_real64, 0.0_real64)
-  complex(real64), public, parameter :: M_z1        = (1.0_real64, 0.0_real64)
-  complex(real64), public, parameter :: M_z2        = (2.0_real64, 0.0_real64)
-  complex(real64), public, parameter :: M_z2I       = (0.0_real64, 2.0_real64)
-  complex(real64), public, parameter :: M_zI        = (0.0_real64, 1.0_real64)
-
-  real(real64), public, parameter :: M_EPSILON   =  epsilon(M_ONE)
-  real(real64), public, parameter :: M_TINY      =  tiny(M_ONE)
-  real(real64), public, parameter :: M_HUGE      =  huge(M_ONE)
-  real(real64), public, parameter :: M_MIN_EXP_ARG = -650_real64
-  real(real64), public, parameter :: M_MAX_EXP_ARG = 700_real64
-
-  !> Minimal occupation that is considered to be non-zero
-  real(real64), public, parameter :: M_MIN_OCC = 1.0e-10_real64
-  !> Minimal density that is considered to be non-zero
-  real(real64), public, parameter :: M_MIN_DENSITY = 1.0e-20_real64
-
-
-  !> some physical constants
-  real(real64), public, parameter :: P_a_B =  0.52917720859_real64
-  real(real64), public, parameter :: P_Ang =  M_ONE / P_a_B
-  real(real64), public, parameter :: P_Ry  =  13.60569193_real64
-  real(real64), public, parameter :: P_eV  =  M_ONE / P_Ry
-  real(real64), public, parameter :: P_Kb  =  8.617343e-5_real64/(M_TWO*P_Ry)  !< Boltzmann constant in Ha/K
-  real(real64), public, parameter :: P_c   =  137.035999679_real64
-  !< Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
-  real(real64), public, parameter :: P_g   =  2.00231930436118_real64
-  real(real64), public, parameter :: P_PROTON_CHARGE = -1.0_real64
-  real(real64), public, parameter :: P_ep  =  M_ONE/(M_FOUR*M_Pi)
-  real(real64), public, parameter :: P_mu  =  M_FOUR*M_PI/(P_c**2)
-
-  !> the standard input and output
-  integer, public :: stderr, stdin, stdout
-
-  !> global epoch time (time at startup)
-  integer, public :: s_epoch_sec, s_epoch_usec
-
-  !> The stack.
-  character(len=80), public          :: sub_stack(50)
-  real(real64), public               :: time_stack(50)
-  integer, public                    :: no_sub_stack = 0
-
-  !> Same for profiling mode.
-  logical, public :: in_profiling_mode = .false.
-
-  integer,    public :: global_alloc_err
-  integer(int64), public :: global_sizeof
-  character(len=100), public :: global_alloc_errmsg
-
-  ! The code directories should be defined here, and not hard coded in the Fortran files.
-  character(len=*), public, parameter :: GS_DIR = "gs/"
-  character(len=*), public, parameter :: TD_DIR = "td/"
-  character(len=*), public, parameter :: STATIC_DIR = "static/"
-  character(len=*), public, parameter :: EM_RESP_DIR = "em_resp/"
-  character(len=*), public, parameter :: EM_RESP_FD_DIR = "em_resp_fd/"
-  character(len=*), public, parameter :: KDOTP_DIR = "kdotp/"
-  character(len=*), public, parameter :: VIB_MODES_DIR = "vib_modes/"
-  character(len=*), public, parameter :: VDW_DIR = "vdw/"
-  character(len=*), public, parameter :: CASIDA_DIR = "casida/"
-  character(len=*), public, parameter :: OCT_DIR = "opt-control/"
-  character(len=*), public, parameter :: PCM_DIR = "pcm/"
-  character(len=*), public, parameter :: PARTITION_DIR = "partition/"
-
-  !> Alias MPI_COMM_UNDEFINED for the specific use case of initialising
-  !! Octopus utilities with no MPI support
-  integer, public, parameter :: SERIAL_DUMMY_COMM = MPI_COMM_UNDEFINED
-
-  ! End of declaration of public objects.
-  ! ---------------------------------------------------------
-
-  interface optional_default
-    module procedure doptional_default, zoptional_default, ioptional_default, loptional_default
-    module procedure looptional_default, soptional_default
-  end interface optional_default
-
-
-  !> This function is defined in messages.F90
-  interface
-    subroutine assert_die(s, f, l)
-      implicit none
-      character(len=*), intent(in) :: s, f
-      integer, intent(in) :: l
-    end subroutine assert_die
-  end interface
-
-  interface operator (+)
-    module procedure cat
-  end interface operator (+)
-
-  interface bitand
-    module procedure bitand48
-    module procedure bitand84
-    module procedure bitand88
-    module procedure bitand44
-  end interface bitand
-
-  interface i4_to_i8
-    module procedure i4_to_i8_0, i4_to_i8_1
-  end interface i4_to_i8
-
-  interface i8_to_i4
-    module procedure i8_to_i4_0, i8_to_i4_1
-  end interface i8_to_i4
-
 contains
 
   !> @brief Initialiser for conf_t
-  subroutine conf_init(this)
+  module subroutine conf_init(this)
     class(conf_t), intent(inout) :: this
 
     character(len=256) :: share
@@ -225,7 +48,7 @@ contains
   !! Main entry point for callers initialising Octopus.
   !! If a communicator is passed, no call is made to initialise MPI_COMM_WORLD.
   !! Else, Octopus initialises MPI_COMM_WORLD
-  subroutine global_init(communicator)
+  module subroutine global_init(communicator)
     integer, intent(in), optional :: communicator  !< Optional MPI communicator from caller
 
     integer :: comm
@@ -254,7 +77,7 @@ contains
   !! * Default CPU cache sizes
   !! * varinfo file, required for the parser
   !! * Configuration instance.
-  subroutine init_octopus_globals(comm)
+  module subroutine init_octopus_globals(comm)
     integer, intent(in) :: comm  !< MPI communicator. Can be a dummy value for serial apps.
 
     call mpi_grp_init(mpi_world, comm)
@@ -282,7 +105,7 @@ contains
 
 
   !> @brief Finalise parser varinfo file, and MPI
-  subroutine global_end()
+  module subroutine global_end()
 
     call varinfo_end()
     call mpi_mod_end()
@@ -290,9 +113,10 @@ contains
   end subroutine global_end
 
 
-  real(real64) pure function doptional_default(opt, def) result(val)
+  pure module function doptional_default(opt, def) result(val)
     real(real64), optional, intent(in) :: opt
     real(real64),    intent(in) :: def
+    real(real64) :: val
 
     val = def
     if (present(opt)) val = opt
@@ -300,9 +124,10 @@ contains
 
   !----------------------------------------------------------
 
-  complex(real64) pure function zoptional_default(opt, def) result(val)
+  pure module function zoptional_default(opt, def) result(val)
     complex(real64), optional, intent(in) :: opt
     complex(real64), intent(in) :: def
+    complex(real64) :: val
 
     val = def
     if (present(opt)) val = opt
@@ -310,9 +135,10 @@ contains
 
   !----------------------------------------------------------
 
-  integer pure function ioptional_default(opt, def) result(val)
+  pure module function ioptional_default(opt, def) result(val)
     integer, optional, intent(in) :: opt
     integer,           intent(in) :: def
+    integer :: val
 
     val = def
     if (present(opt)) val = opt
@@ -320,9 +146,10 @@ contains
 
   !----------------------------------------------------------
 
-  integer(int64) pure function loptional_default(opt, def) result(val)
+  pure module function loptional_default(opt, def) result(val)
     integer(int64), optional, intent(in) :: opt
     integer(int64),           intent(in) :: def
+    integer(int64) :: val
 
     val = def
     if (present(opt)) val = opt
@@ -330,9 +157,10 @@ contains
 
   !----------------------------------------------------------
 
-  logical pure function looptional_default(opt, def) result(val)
+  pure module function looptional_default(opt, def) result(val)
     logical, optional, intent(in) :: opt
     logical,           intent(in) :: def
+    logical :: val
 
     val = def
     if (present(opt)) val = opt
@@ -340,9 +168,10 @@ contains
 
   !----------------------------------------------------------
 
-  character(len=80) pure function soptional_default(opt, def) result(val)
+  pure module function soptional_default(opt, def) result(val)
     character(len=*), optional, intent(in) :: opt
     character(len=*),           intent(in) :: def
+    character(:), allocatable :: val
 
     val = def
     if (present(opt)) val = opt
@@ -350,105 +179,108 @@ contains
 
   !-----------------------------------------------------------
 
-  logical &
-#ifndef HAVE_OPENMP
-    pure &
-#endif
-    function not_in_openmp()
+  module function not_in_openmp() result(res)
+    logical :: res
 
 #ifdef HAVE_OPENMP
-    not_in_openmp = .not. omp_in_parallel()
+    res = .not. omp_in_parallel()
 #else
-    not_in_openmp = .true.
+    res = .true.
 #endif
 
   end function not_in_openmp
 
   !-----------------------------------------------------------
 
-  function cat(str1, str2)
+  module function cat(str1, str2) result(res)
     character(len=*), intent(in) :: str1
     character(len=*), intent(in) :: str2
+    character(len=len(str1) + len(str2)) :: res
 
-    character(len=len(str1) + len(str2)) :: cat
-    cat = str1//str2
+    res = str1//str2
 
   end function cat
 
   ! -----------------------------------------------------------
 
-  integer(int64) pure function bitand48(val1, val2)
+  pure module function bitand48(val1, val2) result(res)
     integer(int32), intent(in) :: val1
     integer(int64), intent(in) :: val2
+    integer(int64) :: res
 
-    bitand48 = iand(int(val1, int64), val2)
+    res = iand(int(val1, int64), val2)
 
   end function bitand48
 
   ! -----------------------------------------------------------
 
-  integer(int64) pure function bitand84(val1, val2)
+  pure module function bitand84(val1, val2) result(res)
     integer(int64), intent(in) :: val1
     integer(int32), intent(in) :: val2
+    integer(int64) :: res
 
-    bitand84 = iand(val1, int(val2, int64))
+    res = iand(val1, int(val2, int64))
 
   end function bitand84
 
   ! -----------------------------------------------------------
 
-  integer(int64) pure function bitand88(val1, val2)
+  pure module function bitand88(val1, val2) result(res)
     integer(int64), intent(in) :: val1
     integer(int64), intent(in) :: val2
+    integer(int64) :: res
 
-    bitand88 = iand(val1, val2)
+    res = iand(val1, val2)
 
   end function bitand88
 
   ! -----------------------------------------------------------
 
-  integer(int32) pure function bitand44(val1, val2)
+  pure module function bitand44(val1, val2) result(res)
     integer(int32), intent(in) :: val1
     integer(int32), intent(in) :: val2
+    integer(int32) :: res
 
-    bitand44 = iand(val1, val2)
+    res = iand(val1, val2)
 
   end function bitand44
 
   ! -----------------------------------------------------------
 
-  integer(int64) pure function i4_to_i8_0(ii)
+  pure module function i4_to_i8_0(ii) result(res)
     integer(int32), intent(in) :: ii
+    integer(int64) :: res
 
-    i4_to_i8_0 = int(ii, int64)
+    res = int(ii, int64)
   end function i4_to_i8_0
 
   ! -----------------------------------------------------------
 
-  integer(int32) pure function i8_to_i4_0(ii)
+  pure module function i8_to_i4_0(ii) result(res)
     integer(int64), intent(in) :: ii
+    integer(int32) :: res
 
-    i8_to_i4_0 = int(ii, int32)
+    res = int(ii, int32)
   end function i8_to_i4_0
 
   ! -----------------------------------------------------------
 
-  pure function i4_to_i8_1(ii)
+  pure module function i4_to_i8_1(ii) result(res)
     integer(int32), intent(in) :: ii(:)
-    integer(int64) :: i4_to_i8_1(lbound(ii, 1, kind=int64):ubound(ii, 1, kind=int64))
+    integer(int64) :: res(lbound(ii, 1, kind=int64):ubound(ii, 1, kind=int64))
 
-    i4_to_i8_1 = int(ii, int64)
+    res = int(ii, int64)
   end function i4_to_i8_1
 
   ! -----------------------------------------------------------
 
-  pure function i8_to_i4_1(ii)
+  pure module function i8_to_i4_1(ii) result(res)
     integer(int64), intent(in) :: ii(:)
-    integer(int32) :: i8_to_i4_1(lbound(ii, 1, kind=int64):ubound(ii, 1, kind=int64))
+    integer(int32) :: res(lbound(ii, 1, kind=int64):ubound(ii, 1, kind=int64))
 
-    i8_to_i4_1 = int(ii, int32)
+    res = int(ii, int32)
   end function i8_to_i4_1
-end module global_oct_m
+end submodule impl
 
 !! Local Variables:
 !! mode: f90
diff --git a/src/basic/global_h.F90 b/src/basic/global_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..fd28c70b26386fc06450df4436197493fc38a38a
--- /dev/null
+++ b/src/basic/global_h.F90
@@ -0,0 +1,290 @@
+#include "global.h"
+
+module global_oct_m
+  use mpi_oct_m
+  use, intrinsic :: iso_fortran_env
+
+  implicit none
+
+  private
+
+  !> Public types, variables and procedures.
+  public ::               &
+    conf_t,               &
+    global_init,          &
+    global_end,           &
+    init_octopus_globals, &
+    optional_default,     &
+    assert_die,           &
+    not_in_openmp,        &
+    operator(+),          &
+    bitand,               &
+    int32, int64,         &
+    real32, real64,       &
+    i4_to_i8,             &
+    i8_to_i4
+  ! Make these kind variables from kind_oct_m public here so that they are
+  ! available basically everywhere in the code. They still need to be in a
+  ! separate module because they are also needed in some low-level modules.
+
+  integer, public, parameter :: MAX_PATH_LEN=512
+  integer, public, parameter :: MAX_OUTPUT_TYPES=44
+
+  !> @brief Build configuration type
+  type conf_t
+    logical            :: devel_version              !< If true then allow unstable parts of the code
+    logical            :: report_memory
+    character(len=256) :: share = SHARE_DIR          !< Name of the share dir
+    character(len=256) :: git_commit  = GIT_COMMIT   !< hash of latest git commit
+    character(len=50)  :: config_time = BUILD_TIME   !< time octopus was configured
+    character(len=20)  :: version = PACKAGE_VERSION  !< version number
+    character(len=256) :: cc = CC                    !< C compiler
+    character(len=256) :: cxx = CXX                  !< C++ compiler
+    character(len=256) :: fc = FC                    !< Fortran compiler
+    ! Split flag definitions in case they don`t fit in one line, following preprocessing
+    character(len=256) :: cflags = &
+      CFLAGS //&
+      CFLAGS_EXTRA
+    character(len=256) :: cxxflags = &
+      CXXFLAGS //&
+      CXXFLAGS_EXTRA
+    character(len=256) :: fcflags = &
+      FCFLAGS //&
+      FCFLAGS_EXTRA
+    integer            :: target_states_block_size = -1
+  contains
+    procedure :: init => conf_init
+  end type conf_t
+
+  !> Global instance of Octopus configuration
+  type(conf_t),      public :: conf
+
+  real(real64), public, parameter :: R_SMALL = 1e-8_real64
+
+  !> Minimal distance between two distinguishable atoms
+  real(real64), public, parameter :: R_MIN_ATOM_DIST = 1e-3_real64
+
+  !> some mathematical constants
+  real(real64), public, parameter :: M_Pi        = 3.1415926535897932384626433832795029_real64
+  real(real64), public, parameter :: M_E         = 2.7182818284590452353602874713526625_real64
+  real(real64), public, parameter :: M_ZERO      = 0.0_real64
+  real(real64), public, parameter :: M_ONE       = 1.0_real64
+  real(real64), public, parameter :: M_TWO       = 2.0_real64
+  real(real64), public, parameter :: M_THREE     = 3.0_real64
+  real(real64), public, parameter :: M_FOUR      = 4.0_real64
+  real(real64), public, parameter :: M_FIVE      = 5.0_real64
+  real(real64), public, parameter :: M_HALF      = 0.5_real64
+  real(real64), public, parameter :: M_THIRD     = M_ONE/M_THREE
+  real(real64), public, parameter :: M_TWOTHIRD  = M_TWO/M_THREE
+  real(real64), public, parameter :: M_FOURTH    = M_ONE/M_FOUR
+  complex(real64), public, parameter :: M_z0        = (0.0_real64, 0.0_real64)
+  complex(real64), public, parameter :: M_z1        = (1.0_real64, 0.0_real64)
+  complex(real64), public, parameter :: M_z2        = (2.0_real64, 0.0_real64)
+  complex(real64), public, parameter :: M_z2I       = (0.0_real64, 2.0_real64)
+  complex(real64), public, parameter :: M_zI        = (0.0_real64, 1.0_real64)
+
+  real(real64), public, parameter :: M_EPSILON   =  epsilon(M_ONE)
+  real(real64), public, parameter :: M_TINY      =  tiny(M_ONE)
+  real(real64), public, parameter :: M_HUGE      =  huge(M_ONE)
+  real(real64), public, parameter :: M_MIN_EXP_ARG = -650_real64
+  real(real64), public, parameter :: M_MAX_EXP_ARG = 700_real64
+
+  !> Minimal occupation that is considered to be non-zero
+  real(real64), public, parameter :: M_MIN_OCC = 1.0e-10_real64
+  !> Minimal density that is considered to be non-zero
+  real(real64), public, parameter :: M_MIN_DENSITY = 1.0e-20_real64
+
+
+  !> some physical constants
+  real(real64), public, parameter :: P_a_B =  0.52917720859_real64
+  real(real64), public, parameter :: P_Ang =  M_ONE / P_a_B
+  real(real64), public, parameter :: P_Ry  =  13.60569193_real64
+  real(real64), public, parameter :: P_eV  =  M_ONE / P_Ry
+  real(real64), public, parameter :: P_Kb  =  8.617343e-5_real64/(M_TWO*P_Ry)  !< Boltzmann constant in Ha/K
+  real(real64), public, parameter :: P_c   =  137.035999679_real64
+  !< Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
+  real(real64), public, parameter :: P_g   =  2.00231930436118_real64
+  real(real64), public, parameter :: P_PROTON_CHARGE = -1.0_real64
+  real(real64), public, parameter :: P_ep  =  M_ONE/(M_FOUR*M_Pi)
+  real(real64), public, parameter :: P_mu  =  M_FOUR*M_PI/(P_c**2)
+
+  !> the standard input and output
+  integer, public :: stderr, stdin, stdout
+
+  !> global epoch time (time at startup)
+  integer, public :: s_epoch_sec, s_epoch_usec
+
+  !> The stack.
+  character(len=80), public          :: sub_stack(50)
+  real(real64), public               :: time_stack(50)
+  integer, public                    :: no_sub_stack = 0
+
+  !> Same for profiling mode.
+  logical, public :: in_profiling_mode = .false.
+
+  integer,    public :: global_alloc_err
+  integer(int64), public :: global_sizeof
+  character(len=100), public :: global_alloc_errmsg
+
+  ! The code directories should be defined here, and not hard coded in the Fortran files.
+  character(len=*), public, parameter :: GS_DIR = "gs/"
+  character(len=*), public, parameter :: TD_DIR = "td/"
+  character(len=*), public, parameter :: STATIC_DIR = "static/"
+  character(len=*), public, parameter :: EM_RESP_DIR = "em_resp/"
+  character(len=*), public, parameter :: EM_RESP_FD_DIR = "em_resp_fd/"
+  character(len=*), public, parameter :: KDOTP_DIR = "kdotp/"
+  character(len=*), public, parameter :: VIB_MODES_DIR = "vib_modes/"
+  character(len=*), public, parameter :: VDW_DIR = "vdw/"
+  character(len=*), public, parameter :: CASIDA_DIR = "casida/"
+  character(len=*), public, parameter :: OCT_DIR = "opt-control/"
+  character(len=*), public, parameter :: PCM_DIR = "pcm/"
+  character(len=*), public, parameter :: PARTITION_DIR = "partition/"
+
+  !> Alias MPI_COMM_UNDEFINED for the specific use case of initialising
+  !! Octopus utilities with no MPI support
+  integer, public, parameter :: SERIAL_DUMMY_COMM = MPI_COMM_UNDEFINED
+
+  ! End of declaration of public objects.
+  ! ---------------------------------------------------------
+
+  interface optional_default
+    module procedure doptional_default, zoptional_default, ioptional_default, loptional_default
+    module procedure looptional_default, soptional_default
+  end interface optional_default
+
+
+  !> This function is defined in messages.F90
+  interface
+    subroutine assert_die(s, f, l)
+      implicit none
+      character(len=*), intent(in) :: s, f
+      integer, intent(in) :: l
+    end subroutine assert_die
+  end interface
+
+  interface operator (+)
+    module procedure cat
+  end interface operator (+)
+
+  interface bitand
+    module procedure bitand48
+    module procedure bitand84
+    module procedure bitand88
+    module procedure bitand44
+  end interface bitand
+
+  interface i4_to_i8
+    module procedure i4_to_i8_0, i4_to_i8_1
+  end interface i4_to_i8
+
+  interface i8_to_i4
+    module procedure i8_to_i4_0, i8_to_i4_1
+  end interface i8_to_i4
+
+  interface
+    module subroutine conf_init(this)
+      class(conf_t), intent(inout) :: this
+    end subroutine conf_init
+
+    module subroutine global_init(communicator)
+      integer, intent(in), optional :: communicator
+    end subroutine global_init
+
+    module subroutine init_octopus_globals(comm)
+      integer, intent(in) :: comm
+    end subroutine init_octopus_globals
+
+    module subroutine global_end()
+    end subroutine global_end
+
+    pure module function doptional_default(opt, def) result(val)
+      real(real64), optional, intent(in) :: opt
+      real(real64),           intent(in) :: def
+      real(real64) :: val
+    end function doptional_default
+
+    pure module function zoptional_default(opt, def) result(val)
+      complex(real64), optional, intent(in) :: opt
+      complex(real64),           intent(in) :: def
+      complex(real64) :: val
+    end function zoptional_default
+
+    pure module function ioptional_default(opt, def) result(val)
+      integer, optional, intent(in) :: opt
+      integer,           intent(in) :: def
+      integer :: val
+    end function ioptional_default
+
+    pure module function loptional_default(opt, def) result(val)
+      integer(int64), optional, intent(in) :: opt
+      integer(int64),           intent(in) :: def
+      integer(int64) :: val
+    end function loptional_default
+
+    pure module function looptional_default(opt, def) result(val)
+      logical, optional, intent(in) :: opt
+      logical,           intent(in) :: def
+      logical :: val
+    end function looptional_default
+
+    pure module function soptional_default(opt, def) result(val)
+      character(len=*), optional, intent(in) :: opt
+      character(len=*),           intent(in) :: def
+      character(:), allocatable :: val
+    end function soptional_default
+
+    module function not_in_openmp() result(res)
+      logical :: res
+    end function not_in_openmp
+
+    module function cat(str1, str2) result(res)
+      character(len=*), intent(in) :: str1
+      character(len=*), intent(in) :: str2
+      character(len=len(str1) + len(str2)) :: res
+    end function cat
+
+    pure module function bitand48(val1, val2) result(res)
+      integer(int32), intent(in) :: val1
+      integer(int64), intent(in) :: val2
+      integer(int64) :: res
+    end function bitand48
+
+    pure module function bitand84(val1, val2) result(res)
+      integer(int64), intent(in) :: val1
+      integer(int32), intent(in) :: val2
+      integer(int64) :: res
+    end function bitand84
+
+    pure module function bitand88(val1, val2) result(res)
+      integer(int64), intent(in) :: val1
+      integer(int64), intent(in) :: val2
+      integer(int64) :: res
+    end function bitand88
+
+    pure module function bitand44(val1, val2) result(res)
+      integer(int32), intent(in) :: val1
+      integer(int32), intent(in) :: val2
+      integer(int32) :: res
+    end function bitand44
+
+    pure module function i4_to_i8_0(ii) result(res)
+      integer(int32), intent(in) :: ii
+      integer(int64) :: res
+    end function i4_to_i8_0
+
+    pure module function i8_to_i4_0(ii) result(res)
+      integer(int64), intent(in) :: ii
+      integer(int32) :: res
+    end function i8_to_i4_0
+
+    pure module function i4_to_i8_1(ii) result(res)
+      integer(int32), intent(in) :: ii(:)
+      integer(int64) :: res(lbound(ii, 1, kind=int64):ubound(ii, 1, kind=int64))
+    end function i4_to_i8_1
+
+    pure module function i8_to_i4_1(ii) result(res)
+      integer(int64), intent(in) :: ii(:)
+      integer(int32) :: res(lbound(ii, 1, kind=int64):ubound(ii, 1, kind=int64))
+    end function i8_to_i4_1
+  end interface
+end module global_oct_m
diff --git a/src/electrons/CMakeLists.txt b/src/electrons/CMakeLists.txt
index f87ba1d062cf6f881dbdba40f9c95555d29ce175..d602c56ffb27dad707310979de83214464f70749 100644
--- a/src/electrons/CMakeLists.txt
+++ b/src/electrons/CMakeLists.txt
@@ -9,6 +9,7 @@ target_sources(Octopus_lib PRIVATE
 		eigen_rmmdiis.F90
 		eigensolver.F90
 		electrons.F90
+		electrons_h.F90
 		electron_space.F90
 		elf.F90
 		energy_calc.F90
diff --git a/src/electrons/electrons.F90 b/src/electrons/electrons.F90
index 9d42228f18f7c0bd90b32899721173c0757f0d74..862d41260766c23d970f515923c0b65d0406e720 100644
--- a/src/electrons/electrons.F90
+++ b/src/electrons/electrons.F90
@@ -20,58 +20,39 @@
 
 #include "global.h"
 
-
-module electrons_oct_m
+submodule (electrons_oct_m) impl
+  use electrons_oct_m
   use accel_oct_m
   use absorbing_boundaries_oct_m
-  use algorithm_oct_m
-  use algorithm_factory_oct_m
   use calc_mode_par_oct_m
   use classical_particles_oct_m
-  use current_oct_m
   use current_to_mxll_field_oct_m
   use debug_oct_m
   use density_oct_m
-  use dipole_oct_m
-  use electron_space_oct_m
   use elf_oct_m
   use energy_calc_oct_m
   use ext_partner_list_oct_m
   use field_transfer_oct_m
   use forces_oct_m
-  use gauge_field_oct_m
   use global_oct_m
-  use grid_oct_m
-  use hamiltonian_elec_oct_m
   use hamiltonian_elec_base_oct_m
   use interaction_enum_oct_m
-  use interaction_oct_m
-  use interaction_partner_oct_m
-  use interaction_surrogate_oct_m
   use ion_dynamics_oct_m
-  use ions_oct_m
   use kick_oct_m
-  use kpoints_oct_m
   use lalg_basic_oct_m
   use lattice_vectors_oct_m
-  use lasers_oct_m
   use lda_u_oct_m
   use loct_oct_m
   use mesh_oct_m
   use messages_oct_m
   use modelmb_particles_oct_m
-  use mpi_oct_m
-  use multicomm_oct_m
   use mxll_e_field_to_matter_oct_m
   use mxll_b_field_to_matter_oct_m
   use mxll_vec_pot_to_matter_oct_m
   use mxll_elec_coupling_oct_m
-  use namespace_oct_m
   use output_oct_m
-  use output_low_oct_m
   use parser_oct_m
   use pes_oct_m
-  use photons_oct_m
   use photon_mode_oct_m
   use photon_mode_mf_oct_m
   use poisson_oct_m
@@ -87,100 +68,27 @@ module electrons_oct_m
   use profiling_oct_m
   use quantity_oct_m
   use regridding_oct_m
-  use scf_oct_m
+  use scf_interface_oct_m
   use space_oct_m
   use states_abst_oct_m
-  use states_elec_oct_m
   use states_elec_dim_oct_m
   use stress_oct_m
   use sort_oct_m
-  use system_oct_m
-  use td_oct_m
+  use td_interface_oct_m
   use td_write_oct_m
   use unit_system_oct_m
-  use v_ks_oct_m
   use xc_oct_m
   use xc_f03_lib_m
   use xc_oep_oct_m
-  use xc_interaction_oct_m
   use xc_oep_photon_oct_m
   use xc_functional_oct_m
 
   implicit none
 
-  private
-  public ::               &
-    electrons_t
-
-
-  !> @brief Class describing the electron system
-  !!
-  !! This class describes a system of electrons and ions.
-  !!
-  !! \todo move the ions into their own ions_t class.
-  type, extends(system_t) :: electrons_t
-    ! Components are public by default
-    type(electron_space_t)       :: space
-    class(ions_t),     pointer   :: ions => NULL() !< the ion component of the system
-    type(photons_t),   pointer   :: photons => null()
-    type(grid_t)                 :: gr    !< the mesh
-    type(states_elec_t)          :: st    !< the states
-    type(v_ks_t)                 :: ks    !< the Kohn-Sham potentials
-    type(output_t)               :: outp  !< the output
-    type(multicomm_t)            :: mc    !< index and domain communicators
-    type(hamiltonian_elec_t)     :: hm    !< the Hamiltonian
-    type(td_t)                   :: td    !< everything related to time propagation
-    type(current_t)              :: current_calculator
-    type(dipole_t)               :: dipole !< total dipole of electrons and ions
-    type(scf_t)                  :: scf   !< SCF for BOMD
-
-    type(kpoints_t) :: kpoints                   !< the k-points
-
-    logical :: generate_epot
-
-    type(states_elec_t)          :: st_copy  !< copy of the states
-
-    ! At the moment this is not treated as an external potential
-    class(lasers_t), pointer :: lasers => null()      !< lasers
-    class(gauge_field_t), pointer :: gfield => null()      !< gauge field
-
-    ! List with all the external partners
-    ! This will become a list of interactions in the future
-    type(partner_list_t) :: ext_partners
-
-    !TODO: have a list of self interactions
-    type(xc_interaction_t), pointer   :: xc_interaction => null()
-
-    logical :: ions_propagated = .false.
-  contains
-    procedure :: init_interaction => electrons_init_interaction
-    procedure :: init_parallelization => electrons_init_parallelization
-    procedure :: init_algorithm => electrons_init_algorithm
-    procedure :: initial_conditions => electrons_initial_conditions
-    procedure :: do_algorithmic_operation => electrons_do_algorithmic_operation
-    procedure :: is_tolerance_reached => electrons_is_tolerance_reached
-    procedure :: update_quantity => electrons_update_quantity
-    procedure :: init_interaction_as_partner => electrons_init_interaction_as_partner
-    procedure :: copy_quantities_to_interaction => electrons_copy_quantities_to_interaction
-    procedure :: output_start => electrons_output_start
-    procedure :: output_write => electrons_output_write
-    procedure :: output_finish => electrons_output_finish
-    procedure :: process_is_slave  => electrons_process_is_slave
-    procedure :: restart_write_data => electrons_restart_write_data
-    procedure :: restart_read_data => electrons_restart_read_data
-    procedure :: update_kinetic_energy => electrons_update_kinetic_energy
-    procedure :: propagation_start => electrons_propagation_start
-    final :: electrons_finalize
-  end type electrons_t
-
-  interface electrons_t
-    procedure electrons_constructor
-  end interface electrons_t
-
 contains
 
   !----------------------------------------------------------
-  function electrons_constructor(namespace, generate_epot) result(sys)
+  module function electrons_constructor(namespace, generate_epot) result(sys)
     class(electrons_t), pointer    :: sys
     type(namespace_t),  intent(in) :: namespace
     logical,  optional, intent(in) :: generate_epot
@@ -268,7 +176,7 @@ contains
   end function electrons_constructor
 
   ! ---------------------------------------------------------
-  subroutine electrons_init_interaction(this, interaction)
+  module subroutine electrons_init_interaction(this, interaction)
     class(electrons_t), target, intent(inout) :: this
     class(interaction_t),       intent(inout) :: interaction
 
@@ -331,7 +239,7 @@ contains
   end subroutine electrons_init_interaction
 
   ! ---------------------------------------------------------
-  subroutine electrons_init_parallelization(this, grp)
+  module subroutine electrons_init_parallelization(this, grp)
     class(electrons_t), intent(inout) :: this
     type(mpi_grp_t),    intent(in)    :: grp
 
@@ -524,7 +432,7 @@ contains
   end subroutine electrons_init_parallelization
 
   ! ---------------------------------------------------------
-  subroutine electrons_init_algorithm(this, factory)
+  module subroutine electrons_init_algorithm(this, factory)
     class(electrons_t),         intent(inout) :: this
     class(algorithm_factory_t), intent(in)    :: factory
 
@@ -550,7 +458,7 @@ contains
   end subroutine electrons_init_algorithm
 
   ! ---------------------------------------------------------
-  subroutine electrons_initial_conditions(this)
+  module subroutine electrons_initial_conditions(this)
     class(electrons_t), intent(inout) :: this
 
     PUSH_SUB(electrons_initial_conditions)
@@ -566,7 +474,7 @@ contains
   end subroutine electrons_initial_conditions
 
   ! ---------------------------------------------------------
-  subroutine electrons_propagation_start(this)
+  module subroutine electrons_propagation_start(this)
     class(electrons_t),      intent(inout) :: this
 
     PUSH_SUB(electrons_propagation_start)
@@ -575,16 +483,17 @@ contains
 
     ! additional initialization needed for electrons
     call td_init_with_wavefunctions(this%td, this%namespace, this%space, this%mc, this%gr, this%ions, &
-      this%ext_partners, this%st, this%ks, this%hm, this%outp, td_get_from_scratch(this%td))
+      this%ext_partners, this%st, this%ks, this%hm, this%outp, this%td%from_scratch)
 
     POP_SUB(electrons_propagation_start)
   end subroutine electrons_propagation_start
 
   ! ---------------------------------------------------------
-  logical function electrons_do_algorithmic_operation(this, operation, updated_quantities) result(done)
+  module function electrons_do_algorithmic_operation(this, operation, updated_quantities) result(done)
     class(electrons_t),             intent(inout) :: this
     class(algorithmic_operation_t), intent(in)    :: operation
     integer, allocatable,           intent(out)   :: updated_quantities(:)
+    logical :: done
 
     logical :: update_energy_
     type(gauge_field_t), pointer :: gfield
@@ -744,9 +653,10 @@ contains
   end function electrons_do_algorithmic_operation
 
   ! ---------------------------------------------------------
-  logical function electrons_is_tolerance_reached(this, tol) result(converged)
+  module function electrons_is_tolerance_reached(this, tol) result(converged)
     class(electrons_t), intent(in) :: this
     real(real64),       intent(in) :: tol
+    logical :: converged
 
     PUSH_SUB(electrons_is_tolerance_reached)
 
@@ -756,7 +666,7 @@ contains
   end function electrons_is_tolerance_reached
 
   ! ---------------------------------------------------------
-  subroutine electrons_update_quantity(this, iq)
+  module subroutine electrons_update_quantity(this, iq)
     class(electrons_t),   intent(inout) :: this
     integer,              intent(in)    :: iq
 
@@ -783,7 +693,7 @@ contains
   end subroutine electrons_update_quantity
 
   ! ---------------------------------------------------------
-  subroutine electrons_init_interaction_as_partner(partner, interaction)
+  module subroutine electrons_init_interaction_as_partner(partner, interaction)
     class(electrons_t),   intent(in)    :: partner
     class(interaction_surrogate_t), intent(inout) :: interaction
 
@@ -801,7 +711,7 @@ contains
   end subroutine electrons_init_interaction_as_partner
 
   ! ---------------------------------------------------------
-  subroutine electrons_copy_quantities_to_interaction(partner, interaction)
+  module subroutine electrons_copy_quantities_to_interaction(partner, interaction)
     class(electrons_t),   intent(inout) :: partner
     class(interaction_surrogate_t), intent(inout) :: interaction
 
@@ -823,7 +733,7 @@ contains
   end subroutine electrons_copy_quantities_to_interaction
 
   ! ---------------------------------------------------------
-  subroutine electrons_output_start(this)
+  module subroutine electrons_output_start(this)
     class(electrons_t), intent(inout) :: this
 
     PUSH_SUB(electrons_output_start)
@@ -832,7 +742,7 @@ contains
   end subroutine electrons_output_start
 
   ! ---------------------------------------------------------
-  subroutine electrons_output_write(this)
+  module subroutine electrons_output_write(this)
     class(electrons_t), intent(inout) :: this
 
     integer :: iter
@@ -858,7 +768,7 @@ contains
   end subroutine electrons_output_write
 
   ! ---------------------------------------------------------
-  subroutine electrons_output_finish(this)
+  module subroutine electrons_output_finish(this)
     class(electrons_t), intent(inout) :: this
 
     PUSH_SUB(electrons_output_finish)
@@ -867,8 +777,9 @@ contains
   end subroutine electrons_output_finish
 
   ! ---------------------------------------------------------
-  logical function electrons_process_is_slave(this) result(is_slave)
+  module function electrons_process_is_slave(this) result(is_slave)
     class(electrons_t), intent(in) :: this
+    logical :: is_slave
 
     PUSH_SUB(electrons_process_is_slave)
 
@@ -985,7 +896,7 @@ contains
   end subroutine electrons_exec_end_of_timestep_tasks
 
   ! ---------------------------------------------------------
-  subroutine electrons_restart_write_data(this)
+  module subroutine electrons_restart_write_data(this)
     class(electrons_t), intent(inout) :: this
 
     integer :: ierr
@@ -1014,8 +925,9 @@ contains
 
   ! ---------------------------------------------------------
   ! this function returns true if restart data could be read
-  logical function electrons_restart_read_data(this)
+  module function electrons_restart_read_data(this) result(res)
     class(electrons_t), intent(inout) :: this
+    logical :: res
 
     logical :: from_scratch
 
@@ -1030,10 +942,10 @@ contains
       call td_set_from_scratch(this%td, from_scratch)
       if (from_scratch) then
         ! restart data could not be loaded
-        electrons_restart_read_data = .false.
+        res = .false.
       else
         ! restart data could be loaded
-        electrons_restart_read_data = .true.
+        res = .true.
       end if
     end select
 
@@ -1042,7 +954,7 @@ contains
   end function electrons_restart_read_data
 
   !----------------------------------------------------------
-  subroutine electrons_update_kinetic_energy(this)
+  module subroutine electrons_update_kinetic_energy(this)
     class(electrons_t), intent(inout) :: this
 
     PUSH_SUB(electrons_update_kinetic_energy)
@@ -1099,12 +1011,9 @@ contains
   end subroutine get_fields_from_interaction
 
   !----------------------------------------------------------
-  subroutine electrons_finalize(sys)
+  module subroutine electrons_finalize(sys)
     type(electrons_t), intent(inout) :: sys
 
-    type(partner_iterator_t) :: iter
-    class(interaction_partner_t), pointer :: partner
-
     PUSH_SUB(electrons_finalize)
 
     if (associated(sys%algo)) then
@@ -1119,12 +1028,7 @@ contains
       call poisson_async_end(sys%hm%psolver, sys%mc)
     end if
 
-    call iter%start(sys%ext_partners)
-    do while (iter%has_next())
-      partner => iter%get_next()
-      SAFE_DEALLOCATE_P(partner)
-    end do
-    call sys%ext_partners%empty()
+    call deallocate_ext_partners()
 
     SAFE_DEALLOCATE_P(sys%xc_interaction)
 
@@ -1156,9 +1060,22 @@ contains
     call system_end(sys)
 
     POP_SUB(electrons_finalize)
+  contains
+    subroutine deallocate_ext_partners()
+
+      type(partner_iterator_t) :: iter
+      class(interaction_partner_t), pointer :: partner
+
+      call iter%start(sys%ext_partners)
+      do while (iter%has_next())
+        partner => iter%get_next()
+        SAFE_DEALLOCATE_P(partner)
+      end do
+      call sys%ext_partners%empty()
+    end subroutine deallocate_ext_partners
   end subroutine electrons_finalize
 
-end module electrons_oct_m
+end submodule impl
 
 !! Local Variables:
 !! mode: f90
diff --git a/src/electrons/electrons_h.F90 b/src/electrons/electrons_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..56d8d402d6e20c2b88f6faadfdd3702e976b4f9d
--- /dev/null
+++ b/src/electrons/electrons_h.F90
@@ -0,0 +1,198 @@
+module electrons_oct_m
+  use algorithm_factory_oct_m
+  use algorithm_oct_m
+  use current_oct_m
+  use dipole_oct_m
+  use electron_space_oct_m
+  use gauge_field_oct_m
+  use global_oct_m
+  use grid_oct_m
+  use hamiltonian_elec_oct_m
+  use interaction_oct_m
+  use interaction_partner_oct_m
+  use interaction_surrogate_oct_m
+  use ions_oct_m
+  use kpoints_oct_m
+  use lasers_oct_m
+  use mpi_oct_m
+  use multicomm_oct_m
+  use namespace_oct_m
+  use output_low_oct_m
+  use photons_oct_m
+  use scf_oct_m
+  use states_elec_oct_m
+  use system_oct_m
+  use td_oct_m
+  use v_ks_oct_m
+  use xc_interaction_oct_m
+
+  implicit none
+
+  private
+  public ::               &
+    electrons_t
+
+
+  !> @brief Class describing the electron system
+  !!
+  !! This class describes a system of electrons and ions.
+  !!
+  !! \todo move the ions into their own ions_t class.
+  type, extends(system_t) :: electrons_t
+    ! Components are public by default
+    type(electron_space_t)       :: space
+    class(ions_t),     pointer   :: ions => NULL() !< the ion component of the system
+    type(photons_t),   pointer   :: photons => null()
+    type(grid_t)                 :: gr    !< the mesh
+    type(states_elec_t)          :: st    !< the states
+    type(v_ks_t)                 :: ks    !< the Kohn-Sham potentials
+    type(output_t)               :: outp  !< the output
+    type(multicomm_t)            :: mc    !< index and domain communicators
+    type(hamiltonian_elec_t)     :: hm    !< the Hamiltonian
+    type(td_t)                   :: td    !< everything related to time propagation
+    type(current_t)              :: current_calculator
+    type(dipole_t)               :: dipole !< total dipole of electrons and ions
+    type(scf_t)                  :: scf   !< SCF for BOMD
+
+    type(kpoints_t) :: kpoints                   !< the k-points
+
+    logical :: generate_epot
+
+    type(states_elec_t)          :: st_copy  !< copy of the states
+
+    ! At the moment this is not treated as an external potential
+    class(lasers_t), pointer :: lasers => null()      !< lasers
+    class(gauge_field_t), pointer :: gfield => null()      !< gauge field
+
+    ! List with all the external partners
+    ! This will become a list of interactions in the future
+    type(partner_list_t) :: ext_partners
+
+    !TODO: have a list of self interactions
+    type(xc_interaction_t), pointer   :: xc_interaction => null()
+
+    logical :: ions_propagated = .false.
+  contains
+    procedure :: init_interaction => electrons_init_interaction
+    procedure :: init_parallelization => electrons_init_parallelization
+    procedure :: init_algorithm => electrons_init_algorithm
+    procedure :: initial_conditions => electrons_initial_conditions
+    procedure :: do_algorithmic_operation => electrons_do_algorithmic_operation
+    procedure :: is_tolerance_reached => electrons_is_tolerance_reached
+    procedure :: update_quantity => electrons_update_quantity
+    procedure :: init_interaction_as_partner => electrons_init_interaction_as_partner
+    procedure :: copy_quantities_to_interaction => electrons_copy_quantities_to_interaction
+    procedure :: output_start => electrons_output_start
+    procedure :: output_write => electrons_output_write
+    procedure :: output_finish => electrons_output_finish
+    procedure :: process_is_slave  => electrons_process_is_slave
+    procedure :: restart_write_data => electrons_restart_write_data
+    procedure :: restart_read_data => electrons_restart_read_data
+    procedure :: update_kinetic_energy => electrons_update_kinetic_energy
+    procedure :: propagation_start => electrons_propagation_start
+    final :: electrons_finalize
+  end type electrons_t
+
+  interface electrons_t
+    procedure electrons_constructor
+  end interface electrons_t
+
+  interface
+    module function electrons_constructor(namespace, generate_epot) result(sys)
+      type(namespace_t),  intent(in) :: namespace
+      logical,  optional, intent(in) :: generate_epot
+      class(electrons_t), pointer    :: sys
+    end function electrons_constructor
+
+    module subroutine electrons_init_interaction(this, interaction)
+      class(electrons_t), target, intent(inout) :: this
+      class(interaction_t),       intent(inout) :: interaction
+    end subroutine electrons_init_interaction
+
+    module subroutine electrons_init_parallelization(this, grp)
+      class(electrons_t), intent(inout) :: this
+      type(mpi_grp_t),    intent(in)    :: grp
+    end subroutine electrons_init_parallelization
+
+    module subroutine electrons_init_algorithm(this, factory)
+      class(electrons_t),         intent(inout) :: this
+      class(algorithm_factory_t), intent(in)    :: factory
+    end subroutine electrons_init_algorithm
+
+    module subroutine electrons_initial_conditions(this)
+      class(electrons_t), intent(inout) :: this
+    end subroutine electrons_initial_conditions
+
+    module subroutine electrons_propagation_start(this)
+      class(electrons_t),      intent(inout) :: this
+    end subroutine electrons_propagation_start
+
+    module function electrons_do_algorithmic_operation(this, operation, updated_quantities) result(done)
+      class(electrons_t),             intent(inout) :: this
+      class(algorithmic_operation_t), intent(in)    :: operation
+      integer, allocatable,           intent(out)   :: updated_quantities(:)
+      logical :: done
+    end function electrons_do_algorithmic_operation
+
+    module function electrons_is_tolerance_reached(this, tol) result(converged)
+      class(electrons_t), intent(in) :: this
+      real(real64),       intent(in) :: tol
+      logical :: converged
+    end function electrons_is_tolerance_reached
+
+    module subroutine electrons_update_quantity(this, iq)
+      class(electrons_t),   intent(inout) :: this
+      integer,              intent(in)    :: iq
+    end subroutine electrons_update_quantity
+
+    module subroutine electrons_update_exposed_quantity(partner, iq)
+      class(electrons_t), intent(inout) :: partner
+      integer,            intent(in)    :: iq
+    end subroutine electrons_update_exposed_quantity
+
+    module subroutine electrons_init_interaction_as_partner(partner, interaction)
+      class(electrons_t),             intent(in)    :: partner
+      class(interaction_surrogate_t), intent(inout) :: interaction
+    end subroutine electrons_init_interaction_as_partner
+
+    module subroutine electrons_copy_quantities_to_interaction(partner, interaction)
+      class(electrons_t),             intent(inout) :: partner
+      class(interaction_surrogate_t), intent(inout) :: interaction
+    end subroutine electrons_copy_quantities_to_interaction
+
+    module subroutine electrons_output_start(this)
+      class(electrons_t), intent(inout) :: this
+    end subroutine electrons_output_start
+
+    module subroutine electrons_output_write(this)
+      class(electrons_t), intent(inout) :: this
+    end subroutine electrons_output_write
+
+    module subroutine electrons_output_finish(this)
+      class(electrons_t), intent(inout) :: this
+    end subroutine electrons_output_finish
+
+    module function electrons_process_is_slave(this) result(is_slave)
+      class(electrons_t), intent(in) :: this
+      logical :: is_slave
+    end function electrons_process_is_slave
+
+    module subroutine electrons_restart_write_data(this)
+      class(electrons_t), intent(inout) :: this
+    end subroutine electrons_restart_write_data
+
+    module function electrons_restart_read_data(this) result(res)
+      class(electrons_t), intent(inout) :: this
+      logical :: res
+    end function electrons_restart_read_data
+
+    module subroutine electrons_update_kinetic_energy(this)
+      class(electrons_t), intent(inout) :: this
+    end subroutine electrons_update_kinetic_energy
+
+    module subroutine electrons_finalize(sys)
+      type(electrons_t), intent(inout) :: sys
+    end subroutine electrons_finalize
+  end interface
+
+end module electrons_oct_m
diff --git a/src/fdep/fortran_dependencies.pl b/src/fdep/fortran_dependencies.pl
index d9644ea6b71b40ae4bab6d619b7d396fa257f3e6..51bf106b3936241dafb06e18f14d0f87b29bc38f 100755
--- a/src/fdep/fortran_dependencies.pl
+++ b/src/fdep/fortran_dependencies.pl
@@ -17,8 +17,9 @@ my %files = ();
 # mode is the first argument: either mod or inc
 my $mode = shift;
 
-my $use_re = qr/^\s*use\s+(\S+)\s*$/;
-my $def_re = qr/^\s*(?:submodule|module)\s+(\S+)\s*$/;
+# Trick fdep to take `use *` and `submodule (*)` as dependency
+my $use_re = qr/^\s*(?:use|submodule)\s+\(?(\w+)\)?.*$/;
+my $def_re = qr/^\s*module\s+(\S+)\s*$/;
 my $inc_re = qr/^\s*(\S+)\s*$/;
 
 sub add_use {
diff --git a/src/hamiltonian/CMakeLists.txt b/src/hamiltonian/CMakeLists.txt
index 803f9218f56b7ad0b9301925139084d124d144da..b72bbeac074328b0c0cab2e069644931e6f5ed12 100644
--- a/src/hamiltonian/CMakeLists.txt
+++ b/src/hamiltonian/CMakeLists.txt
@@ -6,9 +6,10 @@ target_sources(Octopus_lib PRIVATE
 		exchange_operator.F90
 		ext_partner_list.F90
 		gauge_field.F90
-		hamiltonian_abst.F90
+		hamiltonian_abst_h.F90
 		hamiltonian_elec.F90
 		hamiltonian_elec_base.F90
+		hamiltonian_elec_h.F90
 		hgh_projector.F90
 		hirshfeld.F90
 		ion_interaction.F90
diff --git a/src/hamiltonian/hamiltonian_abst.F90 b/src/hamiltonian/hamiltonian_abst_h.F90
similarity index 97%
rename from src/hamiltonian/hamiltonian_abst.F90
rename to src/hamiltonian/hamiltonian_abst_h.F90
index 7ae72767d51559f80eb1829d22a0c4aa2a9849fc..8e3976d228bc0c6f0ba9c21ce7a18ed228306124 100644
--- a/src/hamiltonian/hamiltonian_abst.F90
+++ b/src/hamiltonian/hamiltonian_abst_h.F90
@@ -15,7 +15,6 @@
 !! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 !! 02110-1301, USA.
 !!
-#include "global.h"
 
 !> @brief This module defines an abstract class for Hamiltonians
 !!
@@ -31,11 +30,8 @@ module hamiltonian_abst_oct_m
 
   private
 
-  public ::                           &
-    hamiltonian_abst_t
-
   !> @brief The abstract Hamiltonian class defines a skeleton for specific implementations
-  type, abstract :: hamiltonian_abst_t
+  type, abstract, public :: hamiltonian_abst_t
     !> Spectral range
     real(real64) :: spectral_middle_point
     real(real64) :: spectral_half_span
diff --git a/src/hamiltonian/hamiltonian_elec.F90 b/src/hamiltonian/hamiltonian_elec.F90
index 1a9a97b196a572b55010034ee230b05fb7aa0af9..3db64544fca27bc5cf4bfe212a20a272bd539d66 100644
--- a/src/hamiltonian/hamiltonian_elec.F90
+++ b/src/hamiltonian/hamiltonian_elec.F90
@@ -19,228 +19,46 @@
 
 #include "global.h"
 
-module hamiltonian_elec_oct_m
-  use absorbing_boundaries_oct_m
+submodule (hamiltonian_elec_oct_m) impl
+  use hamiltonian_elec_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 external_potential_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 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 parser_oct_m
   use par_vec_oct_m
-  use poisson_oct_m
+  use parser_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 symmetries_oct_m
   use types_oct_m
   use unit_oct_m
   use unit_system_oct_m
-  use wfs_elec_oct_m
-  use xc_oct_m
   use xc_f03_lib_m
   use xc_functional_oct_m
   use xc_interaction_oct_m
-  use xc_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,        &
-    dhamiltonian_elec_diagonal,           &
-    zhamiltonian_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_dump_vhxc,           &
-    hamiltonian_elec_load_vhxc,           &
-    hamiltonian_elec_set_vhxc,            &
-    hamiltonian_elec_has_kick,            &
-    hamiltonian_elec_copy_and_set_phase
-
-
-  type, extends(hamiltonian_abst_t) :: hamiltonian_elec_t
-    ! Components are public by default
-
-    !> The Hamiltonian must know what are the "dimensions" of the spaces,
-    !! in order to be able to operate on the states.
-    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  !< absorbing boundaries
-    real(real64), allocatable :: vhartree(:) !< Hartree potential
-    real(real64), allocatable :: vxc(:,:)    !< XC potential
-    real(real64), allocatable :: vhxc(:,:)   !< XC potential + Hartree potential + Berry potential
-    real(real64), allocatable :: vtau(:,:)   !< Derivative of e_XC w.r.t. tau
-    real(real64), allocatable :: vberry(:,:) !< Berry phase potential from external e_field
-
-    type(derivatives_t), pointer, private :: der !< pointer to derivatives
-
-    type(nonlocal_pseudopotential_t) :: vnl  !< Nonlocal part of the pseudopotential
-
-    type(ions_t),     pointer :: ions
-    real(real64) :: exx_coef !< how much of EXX to mix
-
-    type(poisson_t)          :: psolver      !< Poisson solver
-
-    !> The self-induced vector potential and magnetic field
-    logical :: self_induced_magnetic
-    real(real64), allocatable :: a_ind(:, :)
-    real(real64), allocatable :: b_ind(:, :)
-
-    integer :: theory_level    !< copied from sys%ks
-    type(xc_t),             pointer :: xc    !< pointer to xc object
-    type(xc_photons_t), pointer :: xc_photons !< pointer to the xc_photons object
-
-    type(epot_t) :: ep         !< handles the external potential
-    type(pcm_t)  :: pcm        !< handles pcm variables
-
-    !> absorbing boundaries
-    logical, private :: adjoint
-
-    !> Mass of the particle (in most cases, mass = 1, electron mass)
-    real(real64), private :: mass
-
-    !> There may be an "inhomogeneous", "source", or "forcing" term (useful for the OCT formalism)
-    logical, private :: inh_term
-    type(states_elec_t) :: inh_st
-
-    !> There may also be a exchange-like term, similar to the one necessary for time-dependent
-    !! Hartree Fock, also useful only for the OCT equations
-    type(oct_exchange_t) :: oct_exchange
-
-    type(scissor_t) :: scissor
-
-    real(real64) :: current_time
-    logical, private :: is_applied_packed  !< This is initialized by the StatesPack variable.
-
-    !> For the DFT+U
-    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  !< List with all the external potentials
-    real(real64), allocatable, public  :: v_ext_pot(:)  !< the potential comming from external potentials
-    real(real64), allocatable, public  :: v_static(:)   !< static scalar potential
-
-    type(ion_electron_local_potential_t) :: v_ie_loc !< Ion-electron local potential interaction
-    type(nlcc_t) :: nlcc !< Ion-electron NLCC interaction
-
-    type(magnetic_constrain_t) :: magnetic_constrain
-
-    !> The possible kick
-    type(kick_t) :: kick
-
-    !> Maxwell-electrons coupling information
-    type(mxll_coupling_t) :: mxll
-    type(zora_t), pointer :: zora
-
-  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 :: set_mass => hamiltonian_elec_set_mass
-  end type hamiltonian_elec_t
-
-  integer, public, parameter :: &
-    LENGTH     = 1,             &
-    VELOCITY   = 2
-
-  integer, public, parameter ::        &
-    INDEPENDENT_PARTICLES = 2, &
-    HARTREE               = 1, &
-    HARTREE_FOCK          = 3, &
-    KOHN_SHAM_DFT         = 4, &
-    GENERALIZED_KOHN_SHAM_DFT = 5, &
-    RDMFT                 = 7
-
-
 contains
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_init(hm, namespace, space, gr, ions, ext_partners, st, theory_level, xc, &
+  module 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
@@ -687,7 +505,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_end(hm)
+  module subroutine hamiltonian_elec_end(hm)
     type(hamiltonian_elec_t), target, intent(inout) :: hm
 
     type(partner_iterator_t) :: iter
@@ -756,11 +574,12 @@ contains
 
   ! ---------------------------------------------------------
   ! True if the Hamiltonian is Hermitian, false otherwise
-  logical function hamiltonian_elec_hermitian(hm)
+  module function hamiltonian_elec_hermitian(hm) result(res)
     class(hamiltonian_elec_t), intent(in) :: hm
+    logical :: res
 
     PUSH_SUB(hamiltonian_elec_hermitian)
-    hamiltonian_elec_hermitian = .not.((hm%abs_boundaries%abtype == IMAGINARY_ABSORBING) .or. &
+    res = .not.((hm%abs_boundaries%abtype == IMAGINARY_ABSORBING) .or. &
       oct_exchange_enabled(hm%oct_exchange))
 
     POP_SUB(hamiltonian_elec_hermitian)
@@ -768,7 +587,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_span(hm, delta, emin, namespace)
+  module 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
@@ -805,15 +624,16 @@ contains
 
 
   ! ---------------------------------------------------------
-  pure logical function hamiltonian_elec_inh_term(hm) result(inh)
+  pure module function hamiltonian_elec_inh_term(hm) result(inh)
     type(hamiltonian_elec_t), intent(in) :: hm
+    logical :: inh
 
     inh = hm%inh_term
   end function hamiltonian_elec_inh_term
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_set_inh(hm, st)
+  module subroutine hamiltonian_elec_set_inh(hm, st)
     type(hamiltonian_elec_t), intent(inout) :: hm
     type(states_elec_t), intent(in)    :: st
 
@@ -828,7 +648,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_remove_inh(hm)
+  module subroutine hamiltonian_elec_remove_inh(hm)
     type(hamiltonian_elec_t), intent(inout) :: hm
 
     PUSH_SUB(hamiltonian_elec_remove_inh)
@@ -842,7 +662,7 @@ contains
   end subroutine hamiltonian_elec_remove_inh
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_adjoint(hm)
+  module subroutine hamiltonian_elec_adjoint(hm)
     type(hamiltonian_elec_t), intent(inout) :: hm
 
     PUSH_SUB(hamiltonian_elec_adjoint)
@@ -859,7 +679,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_not_adjoint(hm)
+  module subroutine hamiltonian_elec_not_adjoint(hm)
     type(hamiltonian_elec_t), intent(inout) :: hm
 
     PUSH_SUB(hamiltonian_elec_not_adjoint)
@@ -877,7 +697,7 @@ contains
 
   ! ---------------------------------------------------------
   !> (re-)build the Hamiltonian for the next application:
-  subroutine hamiltonian_elec_update(this, mesh, namespace, space, ext_partners, time)
+  module 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
@@ -1108,7 +928,7 @@ contains
   !>@brief Update the KS potential of the electronic Hamiltonian
   !!
   ! TODO: See Issue #1064
-  subroutine hamiltonian_elec_update_pot(this, mesh, accumulate)
+  module 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
@@ -1228,7 +1048,7 @@ contains
   end subroutine hamiltonian_elec_update_pot
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
+  module 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
@@ -1311,16 +1131,18 @@ contains
 
   ! -----------------------------------------------------------------
 
-  real(real64) function hamiltonian_elec_get_time(this) result(time)
+  module function hamiltonian_elec_get_time(this) result(time)
     type(hamiltonian_elec_t),   intent(inout) :: this
+    real(real64) :: time
 
     time = this%current_time
   end function hamiltonian_elec_get_time
 
   ! -----------------------------------------------------------------
 
-  pure logical function hamiltonian_elec_apply_packed(this) result(apply)
+  pure module function hamiltonian_elec_apply_packed(this) result(apply)
     class(hamiltonian_elec_t),   intent(in) :: this
+    logical :: apply
 
     apply = this%is_applied_packed
 
@@ -1328,7 +1150,7 @@ contains
 
 
   ! -----------------------------------------------------------------
-  subroutine zhamiltonian_elec_apply_atom (hm, namespace, space, latt, species, pos, ia, mesh, psi, vpsi)
+  module 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
@@ -1358,7 +1180,7 @@ contains
 
 
   ! -----------------------------------------------------------------
-  subroutine hamiltonian_elec_dump_vhxc(restart, hm, space, mesh, ierr)
+  module subroutine hamiltonian_elec_dump_vhxc(restart, hm, space, mesh, ierr)
     type(restart_t),          intent(in)  :: restart
     type(hamiltonian_elec_t), intent(in)  :: hm
     class(space_t),           intent(in)  :: space
@@ -1454,7 +1276,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_load_vhxc(restart, hm, space, mesh, ierr)
+  module subroutine hamiltonian_elec_load_vhxc(restart, hm, space, mesh, ierr)
     type(restart_t),          intent(in)    :: restart
     type(hamiltonian_elec_t), intent(inout) :: hm
     class(space_t),           intent(in)    :: space
@@ -1524,7 +1346,7 @@ contains
   !! CFM4 propagator. It updates the Hamiltonian by considering a
   !! weighted sum of the external potentials at times time(1) and time(2),
   !! weighted by alpha(1) and alpha(2).
-  subroutine hamiltonian_elec_update_with_ext_pot(this, mesh, space, ext_partners, time, mu)
+  module 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
@@ -1720,7 +1542,7 @@ contains
   end subroutine hamiltonian_elec_update_with_ext_pot
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_elec_set_vhxc(hm, mesh, vold, vold_tau)
+  module subroutine hamiltonian_elec_set_vhxc(hm, mesh, vold, vold_tau)
     type(hamiltonian_elec_t), intent(inout)  :: hm
     class(mesh_t),            intent(in)     :: mesh
     real(real64),             intent(in)     :: vold(:, :)
@@ -1736,15 +1558,16 @@ contains
     POP_SUB(hamiltonian_elec_set_vhxc)
   end subroutine hamiltonian_elec_set_vhxc
 
-  logical function hamiltonian_elec_needs_current(hm, states_are_real)
+  module function hamiltonian_elec_needs_current(hm, states_are_real) result(res)
     type(hamiltonian_elec_t), intent(in) :: hm
     logical,                  intent(in) :: states_are_real
+    logical :: res
 
-    hamiltonian_elec_needs_current = .false.
+    res = .false.
 
     if (hm%self_induced_magnetic) then
       if (.not. states_are_real) then
-        hamiltonian_elec_needs_current = .true.
+        res = .true.
       else
         message(1) = 'No current density for real states since it is identically zero.'
         call messages_warning(1)
@@ -1754,7 +1577,7 @@ contains
   end function hamiltonian_elec_needs_current
 
   ! ---------------------------------------------------------
-  subroutine zhamiltonian_elec_apply_all(hm, namespace, mesh, st, hst)
+  module 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
@@ -1803,7 +1626,7 @@ contains
 
   ! ---------------------------------------------------------
 
-  subroutine magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
+  module 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
@@ -1855,7 +1678,7 @@ contains
   end subroutine magnus
 
   ! ---------------------------------------------------------
-  subroutine vborders (mesh, hm, psi, hpsi)
+  module subroutine vborders (mesh, hm, psi, hpsi)
     class(mesh_t),            intent(in)    :: mesh
     type(hamiltonian_elec_t), intent(in)    :: hm
     complex(real64),          intent(in)    :: psi(:)
@@ -1875,19 +1698,20 @@ contains
   end subroutine vborders
 
   ! ---------------------------------------------------------
-  logical function hamiltonian_elec_has_kick(hm)
+  module function hamiltonian_elec_has_kick(hm) result(res)
     type(hamiltonian_elec_t), intent(in)  :: hm
+    logical :: res
 
     PUSH_SUB(hamiltonian_elec_has_kick)
 
-    hamiltonian_elec_has_kick = (abs(hm%kick%delta_strength) >  M_EPSILON)
+    res = (abs(hm%kick%delta_strength) >  M_EPSILON)
 
     POP_SUB(hamiltonian_elec_has_kick)
   end function hamiltonian_elec_has_kick
 
   !> set the effective electron mass, checking whether it was previously redefined.
   !
-  subroutine hamiltonian_elec_set_mass(this, namespace, mass)
+  module 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
@@ -1911,7 +1735,7 @@ contains
   !! If no phase is defined, a packed copy of psib is returned
   !!
   !! TODO: This should should probably belong to wfs_elec_t, but cannot due to circular dependencies
-  subroutine hamiltonian_elec_copy_and_set_phase(hm, gr, kpt, psib, psib_with_phase)
+  module 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              !< k-point distribution
@@ -1949,7 +1773,7 @@ contains
 #include "complex.F90"
 #include "hamiltonian_elec_inc.F90"
 
-end module hamiltonian_elec_oct_m
+end submodule impl
 
 !! Local Variables:
 !! mode: f90
diff --git a/src/hamiltonian/hamiltonian_elec_h.F90 b/src/hamiltonian/hamiltonian_elec_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..bf15740e4fca4516a3d54c67eecd86ddd5292c8f
--- /dev/null
+++ b/src/hamiltonian/hamiltonian_elec_h.F90
@@ -0,0 +1,532 @@
+module hamiltonian_elec_oct_m
+  use absorbing_boundaries_oct_m
+  use batch_oct_m
+  use derivatives_oct_m
+  use distributed_oct_m
+  use electron_space_oct_m
+  use energy_oct_m
+  use epot_oct_m
+  use exchange_operator_oct_m
+  use global_oct_m
+  use grid_oct_m
+  use hamiltonian_abst_oct_m
+  use hamiltonian_elec_base_oct_m
+  use interaction_partner_oct_m
+  use ion_electron_local_potential_oct_m
+  use ions_oct_m
+  use kick_oct_m
+  use kpoints_oct_m
+  use lattice_vectors_oct_m
+  use lda_u_oct_m
+  use magnetic_constrain_oct_m
+  use mesh_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 pcm_oct_m
+  use phase_oct_m
+  use poisson_oct_m
+  use restart_oct_m
+  use scissor_oct_m
+  use space_oct_m
+  use species_oct_m
+  use states_elec_dim_oct_m
+  use states_elec_oct_m
+  use wfs_elec_oct_m
+  use xc_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,        &
+    dhamiltonian_elec_diagonal,           &
+    zhamiltonian_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_dump_vhxc,           &
+    hamiltonian_elec_load_vhxc,           &
+    hamiltonian_elec_set_vhxc,            &
+    hamiltonian_elec_has_kick,            &
+    hamiltonian_elec_copy_and_set_phase
+
+
+  type, extends(hamiltonian_abst_t) :: hamiltonian_elec_t
+    ! Components are public by default
+
+    !> The Hamiltonian must know what are the "dimensions" of the spaces,
+    !! in order to be able to operate on the states.
+    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  !< absorbing boundaries
+    real(real64), allocatable :: vhartree(:) !< Hartree potential
+    real(real64), allocatable :: vxc(:,:)    !< XC potential
+    real(real64), allocatable :: vhxc(:,:)   !< XC potential + Hartree potential + Berry potential
+    real(real64), allocatable :: vtau(:,:)   !< Derivative of e_XC w.r.t. tau
+    real(real64), allocatable :: vberry(:,:) !< Berry phase potential from external e_field
+
+    type(derivatives_t), pointer, private :: der !< pointer to derivatives
+
+    type(nonlocal_pseudopotential_t) :: vnl  !< Nonlocal part of the pseudopotential
+
+    type(ions_t),     pointer :: ions
+    real(real64) :: exx_coef !< how much of EXX to mix
+
+    type(poisson_t)          :: psolver      !< Poisson solver
+
+    !> The self-induced vector potential and magnetic field
+    logical :: self_induced_magnetic
+    real(real64), allocatable :: a_ind(:, :)
+    real(real64), allocatable :: b_ind(:, :)
+
+    integer :: theory_level    !< copied from sys%ks
+    type(xc_t),             pointer :: xc    !< pointer to xc object
+    type(xc_photons_t), pointer :: xc_photons !< pointer to the xc_photons object
+
+    type(epot_t) :: ep         !< handles the external potential
+    type(pcm_t)  :: pcm        !< handles pcm variables
+
+    !> absorbing boundaries
+    logical, private :: adjoint
+
+    !> Mass of the particle (in most cases, mass = 1, electron mass)
+    real(real64), private :: mass
+
+    !> There may be an "inhomogeneous", "source", or "forcing" term (useful for the OCT formalism)
+    logical, private :: inh_term
+    type(states_elec_t) :: inh_st
+
+    !> There may also be a exchange-like term, similar to the one necessary for time-dependent
+    !! Hartree Fock, also useful only for the OCT equations
+    type(oct_exchange_t) :: oct_exchange
+
+    type(scissor_t) :: scissor
+
+    real(real64) :: current_time
+    logical, private :: is_applied_packed  !< This is initialized by the StatesPack variable.
+
+    !> For the DFT+U
+    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  !< List with all the external potentials
+    real(real64), allocatable, public  :: v_ext_pot(:)  !< the potential comming from external potentials
+    real(real64), allocatable, public  :: v_static(:)   !< static scalar potential
+
+    type(ion_electron_local_potential_t) :: v_ie_loc !< Ion-electron local potential interaction
+    type(nlcc_t) :: nlcc !< Ion-electron NLCC interaction
+
+    type(magnetic_constrain_t) :: magnetic_constrain
+
+    !> The possible kick
+    type(kick_t) :: kick
+
+    !> Maxwell-electrons coupling information
+    type(mxll_coupling_t) :: mxll
+    type(zora_t), pointer :: zora
+
+  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 :: set_mass => hamiltonian_elec_set_mass
+  end type hamiltonian_elec_t
+
+  integer, public, parameter :: &
+    LENGTH     = 1,             &
+    VELOCITY   = 2
+
+  integer, public, parameter ::        &
+    INDEPENDENT_PARTICLES = 2, &
+    HARTREE               = 1, &
+    HARTREE_FOCK          = 3, &
+    KOHN_SHAM_DFT         = 4, &
+    GENERALIZED_KOHN_SHAM_DFT = 5, &
+    RDMFT                 = 7
+
+  interface
+    module 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
+    end subroutine hamiltonian_elec_init
+
+    module subroutine hamiltonian_elec_end(hm)
+      type(hamiltonian_elec_t), target, intent(inout) :: hm
+    end subroutine hamiltonian_elec_end
+
+    module function hamiltonian_elec_hermitian(hm) result(res)
+      class(hamiltonian_elec_t), intent(in) :: hm
+      logical :: res
+    end function hamiltonian_elec_hermitian
+
+    module 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
+    end subroutine hamiltonian_elec_span
+
+    pure module function hamiltonian_elec_inh_term(hm) result(inh)
+      type(hamiltonian_elec_t), intent(in) :: hm
+      logical :: inh
+    end function hamiltonian_elec_inh_term
+
+    module subroutine hamiltonian_elec_set_inh(hm, st)
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      type(states_elec_t), intent(in)    :: st
+    end subroutine hamiltonian_elec_set_inh
+
+    module subroutine hamiltonian_elec_remove_inh(hm)
+      type(hamiltonian_elec_t), intent(inout) :: hm
+    end subroutine hamiltonian_elec_remove_inh
+
+    module subroutine hamiltonian_elec_adjoint(hm)
+      type(hamiltonian_elec_t), intent(inout) :: hm
+    end subroutine hamiltonian_elec_adjoint
+
+    module subroutine hamiltonian_elec_not_adjoint(hm)
+      type(hamiltonian_elec_t), intent(inout) :: hm
+    end subroutine hamiltonian_elec_not_adjoint
+
+    module 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
+    end subroutine hamiltonian_elec_update
+
+    module 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
+    end subroutine hamiltonian_elec_update_pot
+
+    module 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
+    end subroutine hamiltonian_elec_epot_generate
+
+    module function hamiltonian_elec_get_time(this) result(time)
+      type(hamiltonian_elec_t),   intent(inout) :: this
+      real(real64) :: time
+    end function hamiltonian_elec_get_time
+
+    pure module function hamiltonian_elec_apply_packed(this) result(apply)
+      class(hamiltonian_elec_t),   intent(in) :: this
+      logical :: apply
+    end function hamiltonian_elec_apply_packed
+
+    module 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(:,:)  !< (gr%np_part, hm%d%dim)
+      complex(real64),          intent(out) :: vpsi(:,:) !< (gr%np, hm%d%dim)
+    end subroutine zhamiltonian_elec_apply_atom
+
+    module subroutine hamiltonian_elec_dump_vhxc(restart, hm, space, mesh, ierr)
+      type(restart_t),          intent(in)  :: restart
+      type(hamiltonian_elec_t), intent(in)  :: hm
+      class(space_t),           intent(in)  :: space
+      class(mesh_t),            intent(in)  :: mesh
+      integer,                  intent(out) :: ierr
+    end subroutine hamiltonian_elec_dump_vhxc
+
+    module subroutine hamiltonian_elec_load_vhxc(restart, hm, space, mesh, ierr)
+      type(restart_t),          intent(in)    :: restart
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      class(space_t),           intent(in)    :: space
+      class(mesh_t),            intent(in)    :: mesh
+      integer,                  intent(out)   :: ierr
+    end subroutine hamiltonian_elec_load_vhxc
+
+    module 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)
+    end subroutine hamiltonian_elec_update_with_ext_pot
+
+    module subroutine hamiltonian_elec_set_vhxc(hm, mesh, vold, vold_tau)
+      type(hamiltonian_elec_t), intent(inout)  :: hm
+      class(mesh_t),            intent(in)     :: mesh
+      real(real64),             intent(in)     :: vold(:, :)
+      real(real64), optional,   intent(in)     :: vold_tau(:, :)
+    end subroutine hamiltonian_elec_set_vhxc
+
+    module function hamiltonian_elec_needs_current(hm, states_are_real) result(res)
+      type(hamiltonian_elec_t), intent(in) :: hm
+      logical,                  intent(in) :: states_are_real
+      logical :: res
+    end function hamiltonian_elec_needs_current
+
+    module 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
+    end subroutine zhamiltonian_elec_apply_all
+
+    module 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),          intent(inout) :: psi(:,:)
+      complex(real64),          intent(out)   :: hpsi(:,:)
+      integer,                  intent(in)    :: ik
+      real(real64),             intent(in)    :: vmagnus(:, :, :)
+      logical, optional,        intent(in)    :: set_phase
+    end subroutine magnus
+
+    module 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(:)
+    end subroutine vborders
+
+    module function hamiltonian_elec_has_kick(hm) result(res)
+      type(hamiltonian_elec_t), intent(in)  :: hm
+      logical :: res
+    end function hamiltonian_elec_has_kick
+
+    module subroutine dhamiltonian_elec_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
+      class(hamiltonian_elec_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
+    end subroutine dhamiltonian_elec_apply
+
+    module subroutine dhamiltonian_elec_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
+      class(hamiltonian_elec_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(:, :, :)
+    end subroutine dhamiltonian_elec_magnus_apply
+
+    module subroutine dhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
+      type(hamiltonian_elec_t),    intent(in)    :: hm
+      type(namespace_t),           intent(in)    :: namespace
+      class(mesh_t),               intent(in)    :: mesh
+      type(wfs_elec_t),    target, intent(inout) :: psib
+      type(wfs_elec_t),    target, intent(inout) :: hpsib
+      integer,           optional, intent(in)    :: terms
+      logical,           optional, intent(in)    :: set_bc
+    end subroutine dhamiltonian_elec_apply_batch
+
+    module subroutine dhamiltonian_elec_external(this, mesh, psib, vpsib)
+      type(hamiltonian_elec_t),    intent(in)    :: this
+      class(mesh_t),               intent(in)    :: mesh
+      type(wfs_elec_t),            intent(in)    :: psib
+      type(wfs_elec_t),            intent(inout) :: vpsib
+    end subroutine dhamiltonian_elec_external
+
+    module subroutine dhamiltonian_elec_apply_single(hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
+      type(hamiltonian_elec_t),         intent(in)    :: hm
+      type(namespace_t),                intent(in)    :: namespace
+      class(mesh_t),                    intent(in)    :: mesh
+      integer,                          intent(in)    :: ist
+      integer,                          intent(in)    :: ik
+      real(real64), contiguous, target, intent(inout) :: psi(:,:)
+      real(real64), contiguous, target, intent(inout) :: hpsi(:,:)
+      integer, optional,                intent(in)    :: terms
+      logical, optional,                intent(in)    :: set_bc
+      logical, optional,                intent(in)    :: set_phase
+    end subroutine dhamiltonian_elec_apply_single
+
+    module subroutine dhamiltonian_elec_magnus_apply_batch(hm, namespace, mesh, psib, hpsib, vmagnus)
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      type(namespace_t),        intent(in)    :: namespace
+      class(mesh_t),            intent(in)    :: mesh
+      type(wfs_elec_t),         intent(inout) :: psib
+      type(wfs_elec_t),         intent(inout) :: hpsib
+      real(real64),             intent(in)    :: vmagnus(:, :, :)
+    end subroutine dhamiltonian_elec_magnus_apply_batch
+
+    module subroutine dh_mgga_terms(hm, mesh, psib, hpsib, ghost_update)
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      class(mesh_t),            intent(in)    :: mesh
+      type(wfs_elec_t),         intent(inout) :: psib
+      type(wfs_elec_t),         intent(inout) :: hpsib
+      logical,                  intent(in)    :: ghost_update
+    end subroutine dh_mgga_terms
+
+    module subroutine dvmask(mesh, hm, st)
+      class(mesh_t),            intent(in)    :: mesh
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      type(states_elec_t),      intent(inout) :: st
+    end subroutine dvmask
+
+    module subroutine dhamiltonian_elec_diagonal(hm, mesh, diag, ik)
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      class(mesh_t),            intent(in)    :: mesh
+      real(real64),             intent(out)   :: diag(:,:)
+      integer,                  intent(in)    :: ik
+    end subroutine dhamiltonian_elec_diagonal
+
+    module subroutine zhamiltonian_elec_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
+      class(hamiltonian_elec_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
+    end subroutine zhamiltonian_elec_apply
+
+    module subroutine zhamiltonian_elec_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
+      class(hamiltonian_elec_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(:, :, :)
+    end subroutine zhamiltonian_elec_magnus_apply
+
+    module subroutine zhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
+      type(hamiltonian_elec_t),    intent(in)    :: hm
+      type(namespace_t),           intent(in)    :: namespace
+      class(mesh_t),               intent(in)    :: mesh
+      type(wfs_elec_t),    target, intent(inout) :: psib
+      type(wfs_elec_t),    target, intent(inout) :: hpsib
+      integer,           optional, intent(in)    :: terms
+      logical,           optional, intent(in)    :: set_bc
+    end subroutine zhamiltonian_elec_apply_batch
+
+    module subroutine zhamiltonian_elec_external(this, mesh, psib, vpsib)
+      type(hamiltonian_elec_t),    intent(in)    :: this
+      class(mesh_t),               intent(in)    :: mesh
+      type(wfs_elec_t),            intent(in)    :: psib
+      type(wfs_elec_t),            intent(inout) :: vpsib
+    end subroutine zhamiltonian_elec_external
+
+    module subroutine zhamiltonian_elec_apply_single(hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
+      type(hamiltonian_elec_t),            intent(in)    :: hm
+      type(namespace_t),                   intent(in)    :: namespace
+      class(mesh_t),                       intent(in)    :: mesh
+      integer,                             intent(in)    :: ist
+      integer,                             intent(in)    :: ik
+      complex(real64), contiguous, target, intent(inout) :: psi(:,:)
+      complex(real64), contiguous, target, intent(inout) :: hpsi(:,:)
+      integer, optional,                   intent(in)    :: terms
+      logical, optional,                   intent(in)    :: set_bc
+      logical, optional,                   intent(in)    :: set_phase
+    end subroutine zhamiltonian_elec_apply_single
+
+    module subroutine zhamiltonian_elec_magnus_apply_batch(hm, namespace, mesh, psib, hpsib, vmagnus)
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      type(namespace_t),        intent(in)    :: namespace
+      class(mesh_t),            intent(in)    :: mesh
+      type(wfs_elec_t),         intent(inout) :: psib
+      type(wfs_elec_t),         intent(inout) :: hpsib
+      real(real64),             intent(in)    :: vmagnus(:, :, :)
+    end subroutine zhamiltonian_elec_magnus_apply_batch
+
+    module subroutine zh_mgga_terms(hm, mesh, psib, hpsib, ghost_update)
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      class(mesh_t),            intent(in)    :: mesh
+      type(wfs_elec_t),         intent(inout) :: psib
+      type(wfs_elec_t),         intent(inout) :: hpsib
+      logical,                  intent(in)    :: ghost_update
+    end subroutine zh_mgga_terms
+
+    module subroutine zvmask(mesh, hm, st)
+      class(mesh_t),            intent(in)    :: mesh
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      type(states_elec_t),      intent(inout) :: st
+    end subroutine zvmask
+
+    module subroutine zhamiltonian_elec_diagonal(hm, mesh, diag, ik)
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      class(mesh_t),            intent(in)    :: mesh
+      complex(real64),          intent(out)   :: diag(:,:)
+      integer,                  intent(in)    :: ik
+    end subroutine zhamiltonian_elec_diagonal
+
+    module 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
+    end subroutine hamiltonian_elec_set_mass
+
+    module 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              !< k-point distribution
+      type(wfs_elec_t),         intent(in)  :: psib             !< Batched wave functions
+      type(wfs_elec_t),         intent(out) :: psib_with_phase  !< Batched wave functions with phase applied
+    end subroutine hamiltonian_elec_copy_and_set_phase
+  end interface
+end module hamiltonian_elec_oct_m
diff --git a/src/hamiltonian/hamiltonian_elec_inc.F90 b/src/hamiltonian/hamiltonian_elec_inc.F90
index 898297ed27644dd98681a9468c415ac010cc7de1..a71db78c8cef6363028a002deaea29cd10c42aa6 100644
--- a/src/hamiltonian/hamiltonian_elec_inc.F90
+++ b/src/hamiltonian/hamiltonian_elec_inc.F90
@@ -17,7 +17,7 @@
 !!
 
 ! ---------------------------------------------------------
-subroutine X(hamiltonian_elec_apply) (hm, namespace, mesh, psib, hpsib, terms, set_bc)
+module subroutine X(hamiltonian_elec_apply) (hm, namespace, mesh, psib, hpsib, terms, set_bc)
   class(hamiltonian_elec_t),   intent(in)    :: hm
   type(namespace_t),           intent(in)    :: namespace
   class(mesh_t),               intent(in)    :: mesh
@@ -46,7 +46,7 @@ subroutine X(hamiltonian_elec_apply) (hm, namespace, mesh, psib, hpsib, terms, s
 end subroutine X(hamiltonian_elec_apply)
 
 ! ---------------------------------------------------------
-subroutine X(hamiltonian_elec_magnus_apply) (hm, namespace, mesh, psib, hpsib, vmagnus)
+module subroutine X(hamiltonian_elec_magnus_apply) (hm, namespace, mesh, psib, hpsib, vmagnus)
   class(hamiltonian_elec_t),   intent(in)    :: hm
   type(namespace_t),           intent(in)    :: namespace
   class(mesh_t),               intent(in)    :: mesh
@@ -74,7 +74,7 @@ subroutine X(hamiltonian_elec_magnus_apply) (hm, namespace, mesh, psib, hpsib, v
 end subroutine X(hamiltonian_elec_magnus_apply)
 
 ! ---------------------------------------------------------
-subroutine X(hamiltonian_elec_apply_batch) (hm, namespace, mesh, psib, hpsib, terms, set_bc)
+module subroutine X(hamiltonian_elec_apply_batch) (hm, namespace, mesh, psib, hpsib, terms, set_bc)
   type(hamiltonian_elec_t),    intent(in)    :: hm
   type(namespace_t),           intent(in)    :: namespace
   class(mesh_t),               intent(in)    :: mesh
@@ -283,7 +283,7 @@ end subroutine X(hamiltonian_elec_apply_batch)
 
 ! ---------------------------------------------------------
 
-subroutine X(hamiltonian_elec_external)(this, mesh, psib, vpsib)
+module subroutine X(hamiltonian_elec_external)(this, mesh, psib, vpsib)
   type(hamiltonian_elec_t),    intent(in)    :: this
   class(mesh_t),               intent(in)    :: mesh
   type(wfs_elec_t),            intent(in)    :: psib
@@ -332,7 +332,7 @@ end subroutine X(hamiltonian_elec_external)
 
 ! ---------------------------------------------------------
 
-subroutine X(hamiltonian_elec_apply_single) (hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
+module subroutine X(hamiltonian_elec_apply_single) (hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
   type(hamiltonian_elec_t), intent(in)    :: hm
   type(namespace_t),        intent(in)    :: namespace
   class(mesh_t),            intent(in)    :: mesh
@@ -366,7 +366,7 @@ subroutine X(hamiltonian_elec_apply_single) (hm, namespace, mesh, psi, hpsi, ist
 end subroutine X(hamiltonian_elec_apply_single)
 
 
-subroutine X(hamiltonian_elec_magnus_apply_batch) (hm, namespace, mesh, psib, hpsib, vmagnus)
+module subroutine X(hamiltonian_elec_magnus_apply_batch) (hm, namespace, mesh, psib, hpsib, vmagnus)
   type(hamiltonian_elec_t), intent(in)    :: hm
   type(namespace_t),        intent(in)    :: namespace
   class(mesh_t),            intent(in)    :: mesh
@@ -418,7 +418,7 @@ subroutine X(hamiltonian_elec_magnus_apply_batch) (hm, namespace, mesh, psib, hp
 end subroutine X(hamiltonian_elec_magnus_apply_batch)
 
 ! ---------------------------------------------------------
-subroutine X(h_mgga_terms) (hm, mesh, psib, hpsib, ghost_update)
+module subroutine X(h_mgga_terms) (hm, mesh, psib, hpsib, ghost_update)
   type(hamiltonian_elec_t), intent(in)    :: hm
   class(mesh_t),            intent(in)    :: mesh
   type(wfs_elec_t),         intent(inout) :: psib
@@ -482,7 +482,7 @@ end subroutine X(h_mgga_terms)
 
 
 ! ---------------------------------------------------------
-subroutine X(vmask) (mesh, hm, st)
+module subroutine X(vmask) (mesh, hm, st)
   class(mesh_t),       intent(in)    :: mesh
   type(hamiltonian_elec_t), intent(in)    :: hm
   type(states_elec_t), intent(inout) :: st
@@ -513,7 +513,7 @@ end subroutine X(vmask)
 
 
 ! ---------------------------------------------------------
-subroutine X(hamiltonian_elec_diagonal) (hm, mesh, diag, ik)
+module subroutine X(hamiltonian_elec_diagonal) (hm, mesh, diag, ik)
   type(hamiltonian_elec_t), intent(in)    :: hm
   class(mesh_t),            intent(in)    :: mesh
   R_TYPE,              intent(out)   :: diag(:,:) !< hpsi(gr%mesh%np, hm%d%dim)
diff --git a/src/main/geom_opt.F90 b/src/main/geom_opt.F90
index 51a13c1c99111e381938a438b82180784f58a224..1249dddaa6613e0aced4412d46820da745538228 100644
--- a/src/main/geom_opt.F90
+++ b/src/main/geom_opt.F90
@@ -46,6 +46,7 @@ module geom_opt_oct_m
   use profiling_oct_m
   use read_coords_oct_m
   use restart_oct_m
+  use scf_interface_oct_m
   use scf_oct_m
   use species_oct_m
   use states_elec_oct_m
diff --git a/src/main/phonons_fd.F90 b/src/main/phonons_fd.F90
index 6eff8e1c605b7f044413a3398290107ca79a7c0e..92ad4e2715759a9cf757e319ac1bc615294bd95e 100644
--- a/src/main/phonons_fd.F90
+++ b/src/main/phonons_fd.F90
@@ -39,6 +39,7 @@ module phonons_fd_oct_m
   use parser_oct_m
   use profiling_oct_m
   use restart_oct_m
+  use scf_interface_oct_m
   use scf_oct_m
   use space_oct_m
   use states_elec_oct_m
diff --git a/src/main/run.F90 b/src/main/run.F90
index 708273bdf962a85e4652618a7cdc00460fa0939a..dc15b4ee4e41e76c8cfdea751f99bce031021469 100644
--- a/src/main/run.F90
+++ b/src/main/run.F90
@@ -55,7 +55,7 @@ module run_oct_m
   use static_pol_oct_m
   use system_factory_oct_m
   use system_oct_m
-  use td_oct_m
+  use td_interface_oct_m
   use test_oct_m
   use time_dependent_oct_m
   use unit_system_oct_m
diff --git a/src/main/static_pol.F90 b/src/main/static_pol.F90
index ac404ef5b56f96b30d3a4d4dbfa4b5017c3511eb..4af6a1b65b66f5dec36a0dbe6aa7e9b21d584861 100644
--- a/src/main/static_pol.F90
+++ b/src/main/static_pol.F90
@@ -39,6 +39,7 @@ module static_pol_oct_m
   use parser_oct_m
   use profiling_oct_m
   use restart_oct_m
+  use scf_interface_oct_m
   use scf_oct_m
   use space_oct_m
   use states_abst_oct_m
diff --git a/src/main/time_dependent.F90 b/src/main/time_dependent.F90
index 556c733da5d67338f83c60ee3b6fc0e7310985a0..468f14985f732d3dafcb3e3a5f20397dca19755a 100644
--- a/src/main/time_dependent.F90
+++ b/src/main/time_dependent.F90
@@ -30,7 +30,7 @@ module time_dependent_oct_m
   use profiling_oct_m
   use restart_oct_m
   use system_oct_m
-  use td_oct_m
+  use td_interface_oct_m
   use walltimer_oct_m
 
   implicit none
diff --git a/src/maxwell/CMakeLists.txt b/src/maxwell/CMakeLists.txt
index 74a639c24ed05817facf9ae366cbd2854541c665..6a6ffcb669da746dae6914a708a7189de862bf56 100644
--- a/src/maxwell/CMakeLists.txt
+++ b/src/maxwell/CMakeLists.txt
@@ -4,6 +4,7 @@ target_sources(Octopus_lib PRIVATE
 		external_densities.F90
                 external_waves.F90
 		hamiltonian_mxll.F90
+		hamiltonian_mxll_h.F90
 		linear_medium.F90
 		maxwell.F90
 		maxwell_boundary_op.F90
diff --git a/src/maxwell/hamiltonian_mxll.F90 b/src/maxwell/hamiltonian_mxll.F90
index 599ff2c43b1ef0b6fa626246c46145b3ee2eef11..62774d567de65985d84af6fca49c75377f763ce2 100644
--- a/src/maxwell/hamiltonian_mxll.F90
+++ b/src/maxwell/hamiltonian_mxll.F90
@@ -18,146 +18,28 @@
 
 #include "global.h"
 
-module hamiltonian_mxll_oct_m
+submodule (hamiltonian_mxll_oct_m) impl
+  use 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
-    !> absorbing boundaries
-    logical :: adjoint = .false.
-
-    real(real64) :: current_time
-    logical :: apply_packed  !< This is initialized by the StatesPack variable.
-
-    logical :: time_zero
-
-    type(nl_operator_t), pointer   :: operators(:)
-
-    type(bc_mxll_t)                :: bc
-    type(derivatives_t), pointer, private :: der !< pointer to derivatives
-    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_periodic       = .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.
-
-    !> maxwell hamiltonian_mxll
-    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
 
   ! ---------------------------------------------------------
   !> Initializing the Maxwell Hamiltonian
-  subroutine hamiltonian_mxll_init(hm, namespace, gr, st)
+  module 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
@@ -258,7 +140,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_mxll_end(hm)
+  module subroutine hamiltonian_mxll_end(hm)
     type(hamiltonian_mxll_t), intent(inout) :: hm
 
     integer :: il
@@ -284,16 +166,17 @@ contains
 
 
   ! ---------------------------------------------------------
-  logical function hamiltonian_mxll_hermitian(hm)
+  module function hamiltonian_mxll_hermitian(hm) result(res)
     class(hamiltonian_mxll_t), intent(in) :: hm
+    logical :: res
 
     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.
+      res = .false.
     else
-      hamiltonian_mxll_hermitian = .true.
+      res = .true.
     end if
 
     POP_SUB(hamiltonian_mxll_hermitian)
@@ -301,7 +184,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_mxll_span(hm, delta, emin, namespace)
+  module 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
@@ -337,7 +220,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_mxll_adjoint(hm)
+  module subroutine hamiltonian_mxll_adjoint(hm)
     type(hamiltonian_mxll_t), intent(inout) :: hm
 
     PUSH_SUB(hamiltonian_mxll_adjoint)
@@ -351,7 +234,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_mxll_not_adjoint(hm)
+  module subroutine hamiltonian_mxll_not_adjoint(hm)
     type(hamiltonian_mxll_t), intent(inout) :: hm
 
     PUSH_SUB(hamiltonian_mxll_not_adjoint)
@@ -366,7 +249,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Maxwell Hamiltonian update (here only the time is updated, can maybe be added to another routine)
-  subroutine hamiltonian_mxll_update(this, time)
+  module subroutine hamiltonian_mxll_update(this, time)
     type(hamiltonian_mxll_t), intent(inout) :: this
     real(real64),   optional, intent(in)    :: time
 
@@ -380,8 +263,9 @@ contains
 
   ! -----------------------------------------------------------------
 
-  real(real64) function hamiltonian_mxll_get_time(this) result(time)
+  module function hamiltonian_mxll_get_time(this) result(time)
     type(hamiltonian_mxll_t),   intent(inout) :: this
+    real(real64) :: time
 
     time = this%current_time
 
@@ -389,9 +273,10 @@ contains
 
   ! -----------------------------------------------------------------
 
-  logical pure function hamiltonian_mxll_apply_packed(this, mesh) result(apply)
+  pure module function hamiltonian_mxll_apply_packed(this, mesh) result(apply)
     type(hamiltonian_mxll_t),   intent(in) :: this
     class(mesh_t),              intent(in) :: mesh
+    logical :: apply
 
     apply = this%apply_packed
     if (mesh%use_curvilinear) apply = .false.
@@ -399,7 +284,7 @@ contains
   end function hamiltonian_mxll_apply_packed
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_mxll_apply_batch(hm, namespace, der, psib, hpsib, time, terms, set_bc)
+  module 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
@@ -680,7 +565,7 @@ contains
   end subroutine hamiltonian_mxll_apply_batch
 
   ! ---------------------------------------------------------
-  subroutine hamiltonian_mxll_apply_simple(hm, namespace, mesh, psib, hpsib, terms, set_bc)
+  module 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
@@ -722,7 +607,7 @@ contains
   end subroutine hamiltonian_mxll_apply_simple
 
   ! ---------------------------------------------------------
-  subroutine mxll_apply_pml_simple(hm, gradb)
+  module 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)
 
@@ -779,7 +664,7 @@ contains
   end subroutine mxll_apply_pml_simple
 
   ! ---------------------------------------------------------
-  subroutine mxll_update_pml_simple(hm, rs_stateb)
+  module subroutine mxll_update_pml_simple(hm, rs_stateb)
     type(hamiltonian_mxll_t),intent(inout) :: hm
     type(batch_t),           intent(inout) :: rs_stateb
 
@@ -838,7 +723,7 @@ contains
   end subroutine mxll_update_pml_simple
 
   ! ---------------------------------------------------------
-  subroutine mxll_copy_pml_simple(hm, rs_stateb)
+  module subroutine mxll_copy_pml_simple(hm, rs_stateb)
     type(hamiltonian_mxll_t),intent(inout) :: hm
     type(batch_t),           intent(inout) :: rs_stateb
 
@@ -875,7 +760,7 @@ contains
   end subroutine mxll_copy_pml_simple
 
   ! ---------------------------------------------------------
-  subroutine mxll_linear_medium_terms_simple(hm, rs_stateb)
+  module subroutine mxll_linear_medium_terms_simple(hm, rs_stateb)
     type(hamiltonian_mxll_t),intent(in)    :: hm
     type(batch_t),           intent(inout) :: rs_stateb
 
@@ -942,7 +827,7 @@ contains
 
   ! --------------------------------------------------------
   !> Apply hamiltonian to real states (not possible)
-  subroutine dhamiltonian_mxll_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
+  module 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
@@ -958,7 +843,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Applying the Maxwell Hamiltonian on Maxwell states
-  subroutine zhamiltonian_mxll_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
+  module 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
@@ -1004,7 +889,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Applying the Maxwell Hamiltonian on Maxwell states with finite difference
-  subroutine maxwell_hamiltonian_apply_fd(hm, der, psi, oppsi)
+  module subroutine maxwell_hamiltonian_apply_fd(hm, der, psi, oppsi)
     type(hamiltonian_mxll_t), intent(in)    :: hm
     type(derivatives_t),      intent(in)    :: der
     complex(real64),          intent(inout) :: psi(:,:)
@@ -1152,7 +1037,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Maxwell Hamiltonian is updated for the PML calculation
-  subroutine maxwell_pml_hamiltonian(hm, der, psi, dir1, dir2, tmp)
+  module 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(:,:)
@@ -1176,7 +1061,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Maxwell Hamiltonian is updated for the PML calculation
-  subroutine maxwell_pml_hamiltonian_medium(hm, der, psi, dir1, dir2, tmp)
+  module 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(:,:)
@@ -1200,7 +1085,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Maxwell Hamiltonian is updated for the PML calculation via Riemann-Silberstein vector
-  subroutine maxwell_pml_calculation_via_riemann_silberstein(hm, der, psi, pml_dir, field_dir, pml)
+  module 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
@@ -1245,7 +1130,7 @@ contains
   ! ---------------------------------------------------------
   !> Maxwell Hamiltonian is updated for the PML calculation via Riemann-Silberstein
   !> vector with medium inside the box
-  subroutine maxwell_pml_calculation_via_riemann_silberstein_medium(hm, der, psi, pml_dir, field_dir, pml)
+  module 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
@@ -1296,7 +1181,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Maxwell Hamiltonian for medium boundaries
-  subroutine maxwell_medium_boundaries_calculation(hm, psi, oppsi)
+  module 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(:,:)
@@ -1358,7 +1243,7 @@ contains
 
   ! ---------------------------------------------------------
   ! > Maxwell Hamiltonian including medium boxes
-  subroutine maxwell_medium_boxes_calculation(hm, der, psi, oppsi)
+  module 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(:,:)
@@ -1421,7 +1306,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Maxwell hamiltonian Magnus (not implemented)
-  subroutine dhamiltonian_mxll_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
+  module 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
@@ -1435,7 +1320,7 @@ contains
 
   ! ---------------------------------------------------------
   !> Maxwell hamiltonian Magnus (not implemented)
-  subroutine zhamiltonian_mxll_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
+  module 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
@@ -1447,7 +1332,7 @@ contains
 
   end subroutine zhamiltonian_mxll_magnus_apply
 
-end module hamiltonian_mxll_oct_m
+end submodule impl
 
 !! Local Variables:
 !! mode: f90
diff --git a/src/maxwell/hamiltonian_mxll_h.F90 b/src/maxwell/hamiltonian_mxll_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..e4b90ec990729ed9675c4b08bf4f2cad25f86eaa
--- /dev/null
+++ b/src/maxwell/hamiltonian_mxll_h.F90
@@ -0,0 +1,306 @@
+module hamiltonian_mxll_oct_m
+  use batch_oct_m
+  use cube_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 linear_medium_to_em_field_oct_m
+  use maxwell_boundary_op_oct_m
+  use mesh_cube_parallel_map_oct_m
+  use mesh_oct_m
+  use namespace_oct_m
+  use nl_operator_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
+    !> absorbing boundaries
+    logical :: adjoint = .false.
+
+    real(real64) :: current_time
+    logical :: apply_packed  !< This is initialized by the StatesPack variable.
+
+    logical :: time_zero
+
+    type(nl_operator_t), pointer   :: operators(:)
+
+    type(bc_mxll_t)                :: bc
+    type(derivatives_t), pointer, private :: der !< pointer to derivatives
+    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_periodic       = .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.
+
+    !> maxwell hamiltonian_mxll
+    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
+
+  interface
+    module 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
+    end subroutine hamiltonian_mxll_init
+
+    module subroutine hamiltonian_mxll_end(hm)
+      type(hamiltonian_mxll_t), intent(inout) :: hm
+    end subroutine hamiltonian_mxll_end
+
+    module function hamiltonian_mxll_hermitian(hm) result(res)
+      class(hamiltonian_mxll_t), intent(in) :: hm
+      logical :: res
+    end function hamiltonian_mxll_hermitian
+
+    module 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
+    end subroutine hamiltonian_mxll_span
+
+    module subroutine hamiltonian_mxll_adjoint(hm)
+      type(hamiltonian_mxll_t), intent(inout) :: hm
+    end subroutine hamiltonian_mxll_adjoint
+
+    module subroutine hamiltonian_mxll_not_adjoint(hm)
+      type(hamiltonian_mxll_t), intent(inout) :: hm
+    end subroutine hamiltonian_mxll_not_adjoint
+
+    module subroutine hamiltonian_mxll_update(this, time)
+      type(hamiltonian_mxll_t), intent(inout) :: this
+      real(real64),   optional, intent(in)    :: time
+    end subroutine hamiltonian_mxll_update
+
+    module function hamiltonian_mxll_get_time(this) result(time)
+      type(hamiltonian_mxll_t),   intent(inout) :: this
+      real(real64) :: time
+    end function hamiltonian_mxll_get_time
+
+    pure module function hamiltonian_mxll_apply_packed(this, mesh) result(apply)
+      type(hamiltonian_mxll_t),   intent(in) :: this
+      class(mesh_t),              intent(in) :: mesh
+      logical :: apply
+    end function hamiltonian_mxll_apply_packed
+
+    module 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
+    end subroutine hamiltonian_mxll_apply_batch
+
+    module 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
+    end subroutine hamiltonian_mxll_apply_simple
+
+    module 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)
+    end subroutine mxll_apply_pml_simple
+
+    module subroutine mxll_copy_pml_simple(hm, rs_stateb)
+      type(hamiltonian_mxll_t),intent(inout) :: hm
+      type(batch_t),           intent(inout) :: rs_stateb
+    end subroutine mxll_copy_pml_simple
+
+    module subroutine mxll_update_pml_simple(hm, rs_stateb)
+      type(hamiltonian_mxll_t),intent(inout) :: hm
+      type(batch_t),           intent(inout) :: rs_stateb
+    end subroutine mxll_update_pml_simple
+
+    module subroutine mxll_linear_medium_terms_simple(hm, rs_stateb)
+      type(hamiltonian_mxll_t),intent(in)    :: hm
+      type(batch_t),           intent(inout) :: rs_stateb
+    end subroutine mxll_linear_medium_terms_simple
+
+    module 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
+    end subroutine dhamiltonian_mxll_apply
+
+    module 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
+    end subroutine zhamiltonian_mxll_apply
+
+    module subroutine maxwell_hamiltonian_apply_fd(hm, der, psi, oppsi)
+      type(hamiltonian_mxll_t), intent(in)    :: hm
+      type(derivatives_t),      intent(in)    :: der
+      complex(real64),          intent(inout) :: psi(:,:)
+      complex(real64),          intent(inout) :: oppsi(:,:)
+    end subroutine maxwell_hamiltonian_apply_fd
+
+    module 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(:)
+    end subroutine maxwell_pml_hamiltonian
+
+    module 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(:,:)
+    end subroutine maxwell_pml_hamiltonian_medium
+
+    module 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(:)
+    end subroutine maxwell_pml_calculation_via_riemann_silberstein
+
+    module 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(:,:)
+    end subroutine maxwell_pml_calculation_via_riemann_silberstein_medium
+
+    module 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(:,:)
+    end subroutine maxwell_medium_boundaries_calculation
+
+    module 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(:,:)
+    end subroutine maxwell_medium_boxes_calculation
+
+    module 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(:, :, :)
+    end subroutine dhamiltonian_mxll_magnus_apply
+
+    module 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(:, :, :)
+    end subroutine zhamiltonian_mxll_magnus_apply
+  end interface
+end module hamiltonian_mxll_oct_m
diff --git a/src/multisystem/CMakeLists.txt b/src/multisystem/CMakeLists.txt
index 51bc99876c3c4ac2002168b06c21330b4df1d421..1c30e995de14b2474fe208f1d5fe2a8019f4d5ee 100644
--- a/src/multisystem/CMakeLists.txt
+++ b/src/multisystem/CMakeLists.txt
@@ -28,4 +28,5 @@ target_sources(Octopus_lib PRIVATE
 		quantity.F90
 		system.F90
 		system_factory_abst.F90
+		system_h.F90
 		)
diff --git a/src/multisystem/system.F90 b/src/multisystem/system.F90
index 7cd8c44b7da7f511c729981770b398acc3d20e10..2c8faaa7de2425a65a4c96c5d71454c587bb2da7 100644
--- a/src/multisystem/system.F90
+++ b/src/multisystem/system.F90
@@ -20,23 +20,14 @@
 
 #include "global.h"
 
-!> This module implements the abstract system type.
-!!
-module system_oct_m
-  use algorithm_oct_m
-  use algorithm_factory_oct_m
+submodule (system_oct_m) impl
+  use system_oct_m
   use debug_oct_m
   use ghost_interaction_oct_m
   use global_oct_m
-  use interactions_factory_abst_oct_m
-  use interaction_partner_oct_m
-  use interaction_oct_m
-  use iteration_counter_oct_m
   use messages_oct_m
-  use mpi_oct_m
   use namespace_oct_m
   use multisystem_debug_oct_m
-  use linked_list_oct_m
   use parser_oct_m
   use profiling_oct_m
   use quantity_oct_m
@@ -44,181 +35,6 @@ module system_oct_m
   use unit_system_oct_m
   use varinfo_oct_m
   implicit none
-
-  private
-  public ::                                    &
-    system_t,                                  &
-    system_execute_algorithm,                  &
-    system_init_parallelization,               &
-    system_init_algorithm,                     &
-    system_init_iteration_counters,            &
-    system_reset_iteration_counters,           &
-    system_create_interactions,                &
-    system_propagation_start,                  &
-    system_propagation_finish,                 &
-    system_restart_read,                       &
-    system_restart_write,                      &
-    system_update_potential_energy,            &
-    system_update_total_energy,                &
-    system_end,                                &
-    system_list_t,                             &
-    system_iterator_t
-
-  type :: barrier_t
-    logical :: active
-    real(real64) :: target_time
-  end type barrier_t
-
-  integer, parameter, public :: &
-    NUMBER_BARRIERS = 1,        &
-    BARRIER_RESTART = 1
-
-  !> @brief Abstract class for systems
-  !!
-  !! All explicit systems are derived from this class.
-  type, extends(interaction_partner_t), abstract :: system_t
-    private
-    type(iteration_counter_t),    public :: iteration
-    class(algorithm_t),  pointer, public :: algo => null()
-
-    integer, allocatable, public :: supported_interactions(:)
-    type(interaction_list_t), public :: interactions !< List with all the interactions of this system
-
-    type(mpi_grp_t), public :: grp  !< mpi group for this system
-
-    type(barrier_t) :: barrier(NUMBER_BARRIERS)
-    real(real64), public :: kinetic_energy     !< Energy not from interactions, like the kinetic energy
-    real(real64), public :: potential_energy   !< Energy from the interactions with external systems
-    real(real64), public :: internal_energy    !< Energy from the interactions with itself and for containers the kinetic energy of its constituents
-    real(real64), public :: total_energy       !< Sum of internal, external, and self energy
-
-  contains
-    procedure :: execute_algorithm =>  system_execute_algorithm !< @copydoc system_oct_m::system_execute_algorithm
-    procedure :: reset_iteration_counters => system_reset_iteration_counters !< @copydoc system_oct_m::system_reset_iteration_counters
-    procedure :: init_algorithm => system_init_algorithm !< @copydoc system_oct_m::system_init_algorithm
-    procedure :: algorithm_finished => system_algorithm_finished !< @copydoc system_oct_m::system_algorithm_finished
-    procedure :: init_iteration_counters => system_init_iteration_counters !< @copydoc system_oct_m::system_init_iteration_counters
-    procedure :: create_interactions => system_create_interactions !< @copydoc system_oct_m::system_create_interactions
-    procedure :: init_parallelization => system_init_parallelization !< @copydoc system_oct_m::system_init_parallelization
-    procedure :: update_couplings => system_update_couplings  !< @copydoc system_oct_m::system_update_couplings
-    procedure :: update_interactions => system_update_interactions !< @copydoc system_oct_m::system_update_interactions
-    procedure :: update_interactions_start => system_update_interactions_start !< @copydoc system_oct_m::system_update_interactions_start
-    procedure :: update_interactions_finish => system_update_interactions_finish !< @copydoc system_oct_m::system_update_interactions_finish
-    procedure :: propagation_start => system_propagation_start !< @copydoc system_oct_m::system_propagation_start
-    procedure :: propagation_finish => system_propagation_finish !< @copydoc system_oct_m::system_propagation_finish
-    procedure :: iteration_info => system_iteration_info !< @copydoc system_oct_m::system_iteration_info
-    procedure :: restart_write => system_restart_write !< @copydoc system_oct_m::system_restart_write
-    procedure :: restart_read => system_restart_read !< @copydoc system_oct_m::system_restart_read
-    procedure :: output_start => system_output_start !< @copydoc system_oct_m::system_output_start
-    procedure :: output_write => system_output_write !< @copydoc system_oct_m::system_output_write
-    procedure :: output_finish => system_output_finish !< @copydoc system_oct_m::system_output_finish
-    procedure :: process_is_slave => system_process_is_slave !< @copydoc system_oct_m::system_process_is_slave
-    procedure :: start_barrier => system_start_barrier !< @copydoc system_oct_m::system_start_barrier
-    procedure :: end_barrier => system_end_barrier !< @copydoc system_oct_m::system_end_barrier
-    procedure :: arrived_at_barrier => system_arrived_at_barrier !< @copydoc system_oct_m::system_arrived_at_barrier
-    procedure :: arrived_at_any_barrier => system_arrived_at_any_barrier !< @copydoc system_oct_m::system_arrived_at_any_barrier
-    procedure :: update_potential_energy => system_update_potential_energy !< @copydoc system_oct_m::system_update_potential_energy
-    procedure :: update_internal_energy => system_update_internal_energy !< @copydoc system_oct_m::system_update_internal_energy
-    procedure :: update_total_energy => system_update_total_energy !< @copydoc system_oct_m::system_update_total_energy
-    procedure(system_init_interaction),          deferred :: init_interaction !< @copydoc system_oct_m::system_init_interaction
-    procedure(system_initial_conditions),        deferred :: initial_conditions !< @copydoc system_oct_m::system_initial_conditions
-    procedure(system_do_algorithmic_operation),  deferred :: do_algorithmic_operation !< @copydoc system_oct_m::system_do_algorithmic_operation
-    procedure(system_is_tolerance_reached),      deferred :: is_tolerance_reached !< @copydoc system_oct_m::system_is_tolerance_reached
-    procedure(system_restart_write_data),        deferred :: restart_write_data !< @copydoc system_oct_m::system_restart_write_data
-    procedure(system_restart_read_data),         deferred :: restart_read_data !< @copydoc system_oct_m::system_restart_read_data
-    procedure(system_update_kinetic_energy),     deferred :: update_kinetic_energy !< @copydoc system_oct_m::system_update_kinetic_energy
-  end type system_t
-
-  abstract interface
-
-    ! ---------------------------------------------------------
-    !> @brief initialize a given interaction of the system
-    subroutine system_init_interaction(this, interaction)
-      import system_t
-      import interaction_t
-      class(system_t), target, intent(inout) :: this
-      class(interaction_t), intent(inout) :: interaction
-    end subroutine system_init_interaction
-
-    ! ---------------------------------------------------------
-    !> set initial conditions for a system
-    subroutine system_initial_conditions(this)
-      import system_t
-      class(system_t), intent(inout) :: this
-    end subroutine system_initial_conditions
-
-    ! ---------------------------------------------------------
-    !> @brief Execute one operation that is part of a larger algorithm.  Returns true
-    !! if the operation was successfully executed, false otherwise.
-    !!
-    !! Unsuccessful operations can occur, e.g. of quantities from an interaction
-    !! are required, but the interaction is still behind in terms of the iteration counters.
-    !!
-    !! On output, the routine should also provide a list quantities that were
-    !! updated. If no quantitiy was updated, then the corresponding array should
-    !! be left unallocated.
-    logical function system_do_algorithmic_operation(this, operation, updated_quantities) result(done)
-      import system_t
-      import algorithmic_operation_t
-      class(system_t),                intent(inout) :: this
-      class(algorithmic_operation_t), intent(in)    :: operation
-      integer, allocatable,           intent(out)   :: updated_quantities(:)
-    end function system_do_algorithmic_operation
-
-    ! ---------------------------------------------------------
-    !> @brief check whether a system has reached a given tolerance
-    logical function system_is_tolerance_reached(this, tol)
-      use, intrinsic :: iso_fortran_env
-      import system_t
-      class(system_t), intent(in) :: this
-      real(real64),    intent(in) :: tol
-    end function system_is_tolerance_reached
-
-    ! ---------------------------------------------------------
-    !> @brief For some algorithms it might be necessary to store the status of a system at a given algorithmic step
-    !!
-    !! This should be implemented by each system in this routine.
-    subroutine system_store_current_status(this)
-      import system_t
-      class(system_t), intent(inout) :: this
-    end subroutine system_store_current_status
-
-    ! ---------------------------------------------------------
-    subroutine system_restart_write_data(this)
-      import system_t
-      class(system_t), intent(inout) :: this
-    end subroutine system_restart_write_data
-
-    ! ---------------------------------------------------------
-    ! this function returns true if restart data could be read
-    logical function system_restart_read_data(this)
-      import system_t
-      class(system_t), intent(inout) :: this
-    end function system_restart_read_data
-    subroutine system_update_kinetic_energy(this)
-      import system_t
-      class(system_t),      intent(inout) :: this
-    end subroutine system_update_kinetic_energy
-
-  end interface
-
-  !> @brief These classes extends the list and list iterator to create a system list.
-  !!
-  !! Since a list of systems is also a list of interaction partners, the system
-  !! list is an extension of the partner list.
-  type, extends(partner_list_t) :: system_list_t
-    private
-  contains
-    procedure :: add => system_list_add_node !< @copydoc system_oct_m::system_list_add_node
-    procedure :: contains => system_list_contains !< @copydoc system_oct_m::system_list_contains
-  end type system_list_t
-
-  type, extends(linked_list_iterator_t) :: system_iterator_t
-    private
-  contains
-    procedure :: get_next => system_iterator_get_next !< @copydoc system_oct_m::system_iterator_get_next
-  end type system_iterator_t
-
 contains
 
   ! ---------------------------------------------------------
@@ -234,7 +50,7 @@ contains
   !! The couplings update is always considered a barrier, even if the update was
   !! successful. This is to allow other system to also update their couplings
   !! with this system before it moves on to the next operations.
-  subroutine system_execute_algorithm(this)
+  module subroutine system_execute_algorithm(this)
     class(system_t),     intent(inout) :: this
 
     type(algorithmic_operation_t) :: operation
@@ -352,7 +168,7 @@ contains
   end subroutine system_execute_algorithm
 
   ! ---------------------------------------------------------
-  subroutine system_reset_iteration_counters(this, accumulated_iterations)
+  module subroutine system_reset_iteration_counters(this, accumulated_iterations)
     class(system_t),      intent(inout) :: this
     integer,              intent(in)    :: accumulated_iterations
 
@@ -405,7 +221,7 @@ contains
   !! and all available partners. Any class overriding this method must make sure
   !! ghost interactions are properly created or the framework might not work
   !! correctly.
-  recursive subroutine system_create_interactions(this, interaction_factory, available_partners)
+  recursive module subroutine system_create_interactions(this, interaction_factory, available_partners)
     class(system_t),                       intent(inout) :: this                !< system for which interactions are created.
     class(interactions_factory_abst_t),    intent(in)    :: interaction_factory !< factory that creates the actual interactions
     class(partner_list_t),         target, intent(in)    :: available_partners  !< a list of available partners for the given system.
@@ -556,8 +372,9 @@ contains
   !! This function loops over all interactions and the corresponding interaction partners
   !! and attempts to update their couplings to the requested iteration. It returns true if all
   !! couplings have been successfully updated.
-  logical function system_update_couplings(this) result(all_updated)
+  module function system_update_couplings(this) result(all_updated)
     class(system_t),      intent(inout) :: this
+    logical :: all_updated
 
     class(interaction_t),              pointer :: interaction
     type(interaction_iterator_t) :: iter
@@ -593,7 +410,7 @@ contains
   !!
   !! First we try to update the systems own quantities required for the interaction,
   !! and then try to update the interaction itself.
-  subroutine system_update_interactions(this)
+  module subroutine system_update_interactions(this)
     class(system_t),      intent(inout) :: this
 
     integer :: iq, q_id, n_quantities
@@ -666,7 +483,7 @@ contains
   end subroutine system_update_interactions
 
   ! ---------------------------------------------------------
-  subroutine system_update_interactions_start(this)
+  module subroutine system_update_interactions_start(this)
     class(system_t), intent(inout) :: this
 
     PUSH_SUB(system_update_interactions_start)
@@ -678,7 +495,7 @@ contains
   end subroutine system_update_interactions_start
 
   ! ---------------------------------------------------------
-  subroutine system_update_interactions_finish(this)
+  module subroutine system_update_interactions_finish(this)
     class(system_t), intent(inout) :: this
 
     PUSH_SUB(system_update_interactions_finish)
@@ -690,7 +507,7 @@ contains
   end subroutine system_update_interactions_finish
 
   ! ---------------------------------------------------------
-  subroutine system_restart_write(this)
+  module subroutine system_restart_write(this)
     class(system_t), intent(inout) :: this
 
     logical :: restart_write
@@ -729,8 +546,9 @@ contains
 
   ! ---------------------------------------------------------
   ! this function returns true if restart data could be read
-  logical function system_restart_read(this)
+  module function system_restart_read(this) result(res)
     class(system_t), intent(inout) :: this
+    logical :: res
 
     type(interaction_iterator_t) :: iter
     class(interaction_t),              pointer :: interaction
@@ -740,19 +558,19 @@ contains
 
     ! do some generic restart steps here
     ! read iteration data
-    system_restart_read = this%iteration%restart_read('restart_iteration_system', this%namespace)
-    system_restart_read = system_restart_read .and. &
+    res = this%iteration%restart_read('restart_iteration_system', this%namespace)
+    res = res .and. &
       this%algo%iteration%restart_read('restart_iteration_propagator', this%namespace)
     call iter%start(this%interactions)
     do while (iter%has_next())
       interaction => iter%get_next()
-      system_restart_read = system_restart_read .and. interaction%restart_read(this%namespace)
+      res = res .and. interaction%restart_read(this%namespace)
       ! reduce by one because of the first UPDATE_INTERACTIONS
       interaction%iteration = interaction%iteration - 1
     end do
     do ii = 1, MAX_QUANTITIES
       if (this%quantities(ii)%required) then
-        system_restart_read = system_restart_read .and. &
+        res = res .and. &
           this%quantities(ii)%iteration%restart_read('restart_iteration_quantity_'//trim(QUANTITY_LABEL(ii)), &
           this%namespace)
       end if
@@ -762,9 +580,9 @@ contains
       end if
     end do
     ! the following call is delegated to the corresponding system
-    system_restart_read = system_restart_read .and. this%restart_read_data()
+    res = res .and. this%restart_read_data()
 
-    if (system_restart_read) then
+    if (res) then
       message(1) = "Successfully read restart data for system "//trim(this%namespace%get())
       call messages_info(1, namespace=this%namespace)
     end if
@@ -773,7 +591,7 @@ contains
   end function system_restart_read
 
   ! ---------------------------------------------------------
-  subroutine system_output_start(this)
+  module subroutine system_output_start(this)
     class(system_t), intent(inout) :: this
 
     PUSH_SUB(system_output_start)
@@ -785,7 +603,7 @@ contains
   end subroutine system_output_start
 
   ! ---------------------------------------------------------
-  subroutine system_output_write(this)
+  module subroutine system_output_write(this)
     class(system_t), intent(inout) :: this
 
     PUSH_SUB(system_output_write)
@@ -797,7 +615,7 @@ contains
   end subroutine system_output_write
 
   ! ---------------------------------------------------------
-  subroutine system_output_finish(this)
+  module subroutine system_output_finish(this)
     class(system_t), intent(inout) :: this
 
     PUSH_SUB(system_output_finish)
@@ -809,7 +627,7 @@ contains
   end subroutine system_output_finish
 
   ! ---------------------------------------------------------
-  subroutine system_init_algorithm(this, factory)
+  module subroutine system_init_algorithm(this, factory)
     class(system_t),            intent(inout) :: this
     class(algorithm_factory_t), intent(in)    :: factory
 
@@ -832,7 +650,7 @@ contains
   end subroutine system_init_algorithm
 
   ! ---------------------------------------------------------------------------------------
-  recursive function system_algorithm_finished(this) result(finished)
+  recursive module function system_algorithm_finished(this) result(finished)
     class(system_t),       intent(in) :: this
     logical :: finished
 
@@ -847,7 +665,7 @@ contains
   !! before the algorithm iteration counter. This is necessary, as the interactions and on-demand quantities
   !! first need to be updated.
   !
-  subroutine system_init_iteration_counters(this)
+  module subroutine system_init_iteration_counters(this)
     class(system_t),            intent(inout) :: this
 
     type(interaction_iterator_t) :: iter
@@ -880,7 +698,7 @@ contains
   end subroutine system_init_iteration_counters
 
   ! ---------------------------------------------------------
-  subroutine system_propagation_start(this)
+  module subroutine system_propagation_start(this)
     class(system_t),      intent(inout) :: this
 
     logical :: all_updated
@@ -934,7 +752,7 @@ contains
   end subroutine system_propagation_start
 
   ! ---------------------------------------------------------
-  subroutine system_propagation_finish(this)
+  module subroutine system_propagation_finish(this)
     class(system_t),      intent(inout) :: this
     type(event_handle_t) :: debug_handle
 
@@ -967,7 +785,7 @@ contains
   end subroutine system_propagation_finish
 
   ! ---------------------------------------------------------
-  subroutine system_iteration_info(this)
+  module subroutine system_iteration_info(this)
     class(system_t), intent(in) :: this
 
     real(real64) :: energy
@@ -996,19 +814,20 @@ contains
   end subroutine system_iteration_info
 
   ! ---------------------------------------------------------
-  logical function system_process_is_slave(this)
+  module function system_process_is_slave(this) result(res)
     class(system_t), intent(in) :: this
+    logical :: res
 
     PUSH_SUB(system_process_is_slave)
 
     ! By default an MPI process is not a slave
-    system_process_is_slave = .false.
+    res = .false.
 
     POP_SUB(system_process_is_slave)
   end function system_process_is_slave
 
   ! ---------------------------------------------------------
-  subroutine system_end(this)
+  module subroutine system_end(this)
     class(system_t), intent(inout) :: this
 
     type(interaction_iterator_t) :: iter
@@ -1033,7 +852,7 @@ contains
   end subroutine system_end
 
   ! ---------------------------------------------------------
-  subroutine system_list_add_node(this, partner)
+  module subroutine system_list_add_node(this, partner)
     class(system_list_t)         :: this
     class(interaction_partner_t), target :: partner
 
@@ -1050,9 +869,10 @@ contains
   end subroutine system_list_add_node
 
   ! ---------------------------------------------------------
-  recursive logical function system_list_contains(this, partner) result(contains)
+  recursive module function system_list_contains(this, partner) result(contains)
     class(system_list_t)         :: this
     class(interaction_partner_t), target :: partner
+    logical :: contains
 
     type(partner_iterator_t)  :: iterator
     class(interaction_partner_t),  pointer :: system
@@ -1078,7 +898,7 @@ contains
   end function system_list_contains
 
   ! ---------------------------------------------------------
-  function system_iterator_get_next(this) result(system)
+  module function system_iterator_get_next(this) result(system)
     class(system_iterator_t), intent(inout) :: this
     class(system_t),          pointer       :: system
 
@@ -1098,7 +918,7 @@ contains
   !> Basic functionality: copy the MPI group.
   !! This function needs to be implemented by extended types
   !! that need more initialization for their parallelization.
-  subroutine system_init_parallelization(this, grp)
+  module subroutine system_init_parallelization(this, grp)
     class(system_t), intent(inout) :: this
     type(mpi_grp_t), intent(in)    :: grp
 
@@ -1113,7 +933,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine system_start_barrier(this, target_time, barrier_index)
+  module subroutine system_start_barrier(this, target_time, barrier_index)
     class(system_t), intent(inout) :: this
     real(real64),    intent(in)    :: target_time
     integer,         intent(in)    :: barrier_index
@@ -1127,7 +947,7 @@ contains
   end subroutine system_start_barrier
 
   ! ---------------------------------------------------------
-  subroutine system_end_barrier(this, barrier_index)
+  module subroutine system_end_barrier(this, barrier_index)
     class(system_t), intent(inout) :: this
     integer,         intent(in)    :: barrier_index
 
@@ -1140,19 +960,20 @@ contains
   end subroutine system_end_barrier
 
   ! ---------------------------------------------------------
-  logical function system_arrived_at_barrier(this, barrier_index)
+  module function system_arrived_at_barrier(this, barrier_index) result(res)
     class(system_t), intent(inout) :: this
     integer,         intent(in)    :: barrier_index
+    logical :: res
 
     type(iteration_counter_t) :: iteration
 
     PUSH_SUB(system_arrived_at_barrier)
 
-    system_arrived_at_barrier = .false.
+    res = .false.
     if (this%barrier(barrier_index)%active) then
       iteration = this%iteration + 1
       if (iteration%value() > this%barrier(barrier_index)%target_time) then
-        system_arrived_at_barrier = .true.
+        res = .true.
       end if
     end if
 
@@ -1160,16 +981,17 @@ contains
   end function system_arrived_at_barrier
 
   ! ---------------------------------------------------------
-  logical function system_arrived_at_any_barrier(this)
+  module function system_arrived_at_any_barrier(this) result(res)
     class(system_t), intent(inout) :: this
+    logical :: res
 
     integer :: ii
 
     PUSH_SUB(system_arrived_at_any_barrier)
 
-    system_arrived_at_any_barrier = .false.
+    res = .false.
     do ii = 1, NUMBER_BARRIERS
-      system_arrived_at_any_barrier = system_arrived_at_any_barrier &
+      res = res &
         .or. this%arrived_at_barrier(ii)
     end do
 
@@ -1182,7 +1004,7 @@ contains
   !! The potential energy is defined as the sum of all energies
   !! arising from interactions with external systems.
   !! (Note that multisystems override this function)
-  subroutine system_update_potential_energy(this)
+  module subroutine system_update_potential_energy(this)
     class(system_t), intent(inout) :: this
 
     type(interaction_iterator_t) :: iter
@@ -1209,7 +1031,7 @@ contains
   !! The internal energy is defined as the sum of all energies
   !! arising from intra-interactions and the entropy terms (if available).
   !! (Note that multisystems override this function)
-  subroutine system_update_internal_energy(this)
+  module subroutine system_update_internal_energy(this)
     class(system_t), intent(inout) :: this
 
     type(interaction_iterator_t) :: iter
@@ -1234,7 +1056,7 @@ contains
   !> Calculate the total energy of the system.
   !! The total energy is defined as the sum of
   !! the kinetic, the internal and the potential energies.
-  subroutine system_update_total_energy(this)
+  module subroutine system_update_total_energy(this)
     class(system_t), intent(inout) :: this
 
     PUSH_SUB(system_update_total_energy)
@@ -1253,7 +1075,7 @@ contains
     POP_SUB(system_update_total_energy)
   end subroutine system_update_total_energy
 
-end module system_oct_m
+end submodule impl
 
 !! Local Variables:
 !! mode: f90
diff --git a/src/multisystem/system_h.F90 b/src/multisystem/system_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..709db1a6bcd4b34110d34e8fe51dd4f572e28e63
--- /dev/null
+++ b/src/multisystem/system_h.F90
@@ -0,0 +1,334 @@
+!> This module implements the abstract system type.
+!!
+module system_oct_m
+  use algorithm_factory_oct_m
+  use algorithm_oct_m
+  use global_oct_m
+  use interaction_oct_m
+  use interaction_partner_oct_m
+  use interactions_factory_abst_oct_m
+  use iteration_counter_oct_m
+  use linked_list_oct_m
+  use mpi_oct_m
+  implicit none
+
+  private
+  public ::                                    &
+    system_t,                                  &
+    system_execute_algorithm,                  &
+    system_init_parallelization,               &
+    system_init_algorithm,                     &
+    system_init_iteration_counters,            &
+    system_reset_iteration_counters,           &
+    system_create_interactions,                &
+    system_propagation_start,                  &
+    system_propagation_finish,                 &
+    system_restart_read,                       &
+    system_restart_write,                      &
+    system_update_potential_energy,            &
+    system_update_total_energy,                &
+    system_end,                                &
+    system_list_t,                             &
+    system_iterator_t
+
+  type :: barrier_t
+    logical :: active
+    real(real64) :: target_time
+  end type barrier_t
+
+  integer, parameter, public :: &
+    NUMBER_BARRIERS = 1,        &
+    BARRIER_RESTART = 1
+
+  !> @brief Abstract class for systems
+  !!
+  !! All explicit systems are derived from this class.
+  type, extends(interaction_partner_t), abstract :: system_t
+    private
+    type(iteration_counter_t),    public :: iteration
+    class(algorithm_t),  pointer, public :: algo => null()
+
+    integer, allocatable, public :: supported_interactions(:)
+    type(interaction_list_t), public :: interactions !< List with all the interactions of this system
+
+    type(mpi_grp_t), public :: grp  !< mpi group for this system
+
+    type(barrier_t) :: barrier(NUMBER_BARRIERS)
+    real(real64), public :: kinetic_energy     !< Energy not from interactions, like the kinetic energy
+    real(real64), public :: potential_energy   !< Energy from the interactions with external systems
+    real(real64), public :: internal_energy    !< Energy from the interactions with itself and for containers the kinetic energy of its constituents
+    real(real64), public :: total_energy       !< Sum of internal, external, and self energy
+
+  contains
+    procedure :: execute_algorithm =>  system_execute_algorithm !< @copydoc system_oct_m::system_execute_algorithm
+    procedure :: reset_iteration_counters => system_reset_iteration_counters !< @copydoc system_oct_m::system_reset_iteration_counters
+    procedure :: init_algorithm => system_init_algorithm !< @copydoc system_oct_m::system_init_algorithm
+    procedure :: algorithm_finished => system_algorithm_finished !< @copydoc system_oct_m::system_algorithm_finished
+    procedure :: init_iteration_counters => system_init_iteration_counters !< @copydoc system_oct_m::system_init_iteration_counters
+    procedure :: create_interactions => system_create_interactions !< @copydoc system_oct_m::system_create_interactions
+    procedure :: init_parallelization => system_init_parallelization !< @copydoc system_oct_m::system_init_parallelization
+    procedure :: update_couplings => system_update_couplings  !< @copydoc system_oct_m::system_update_couplings
+    procedure :: update_interactions => system_update_interactions !< @copydoc system_oct_m::system_update_interactions
+    procedure :: update_interactions_start => system_update_interactions_start !< @copydoc system_oct_m::system_update_interactions_start
+    procedure :: update_interactions_finish => system_update_interactions_finish !< @copydoc system_oct_m::system_update_interactions_finish
+    procedure :: propagation_start => system_propagation_start !< @copydoc system_oct_m::system_propagation_start
+    procedure :: propagation_finish => system_propagation_finish !< @copydoc system_oct_m::system_propagation_finish
+    procedure :: iteration_info => system_iteration_info !< @copydoc system_oct_m::system_iteration_info
+    procedure :: restart_write => system_restart_write !< @copydoc system_oct_m::system_restart_write
+    procedure :: restart_read => system_restart_read !< @copydoc system_oct_m::system_restart_read
+    procedure :: output_start => system_output_start !< @copydoc system_oct_m::system_output_start
+    procedure :: output_write => system_output_write !< @copydoc system_oct_m::system_output_write
+    procedure :: output_finish => system_output_finish !< @copydoc system_oct_m::system_output_finish
+    procedure :: process_is_slave => system_process_is_slave !< @copydoc system_oct_m::system_process_is_slave
+    procedure :: start_barrier => system_start_barrier !< @copydoc system_oct_m::system_start_barrier
+    procedure :: end_barrier => system_end_barrier !< @copydoc system_oct_m::system_end_barrier
+    procedure :: arrived_at_barrier => system_arrived_at_barrier !< @copydoc system_oct_m::system_arrived_at_barrier
+    procedure :: arrived_at_any_barrier => system_arrived_at_any_barrier !< @copydoc system_oct_m::system_arrived_at_any_barrier
+    procedure :: update_potential_energy => system_update_potential_energy !< @copydoc system_oct_m::system_update_potential_energy
+    procedure :: update_internal_energy => system_update_internal_energy !< @copydoc system_oct_m::system_update_internal_energy
+    procedure :: update_total_energy => system_update_total_energy !< @copydoc system_oct_m::system_update_total_energy
+    procedure(system_init_interaction),          deferred :: init_interaction !< @copydoc system_oct_m::system_init_interaction
+    procedure(system_initial_conditions),        deferred :: initial_conditions !< @copydoc system_oct_m::system_initial_conditions
+    procedure(system_do_algorithmic_operation),  deferred :: do_algorithmic_operation !< @copydoc system_oct_m::system_do_algorithmic_operation
+    procedure(system_is_tolerance_reached),      deferred :: is_tolerance_reached !< @copydoc system_oct_m::system_is_tolerance_reached
+    procedure(system_restart_write_data),        deferred :: restart_write_data !< @copydoc system_oct_m::system_restart_write_data
+    procedure(system_restart_read_data),         deferred :: restart_read_data !< @copydoc system_oct_m::system_restart_read_data
+    procedure(system_update_kinetic_energy),     deferred :: update_kinetic_energy !< @copydoc system_oct_m::system_update_kinetic_energy
+  end type system_t
+
+  abstract interface
+
+    ! ---------------------------------------------------------
+    !> @brief initialize a given interaction of the system
+    subroutine system_init_interaction(this, interaction)
+      import system_t
+      import interaction_t
+      class(system_t), target, intent(inout) :: this
+      class(interaction_t), intent(inout) :: interaction
+    end subroutine system_init_interaction
+
+    ! ---------------------------------------------------------
+    !> set initial conditions for a system
+    subroutine system_initial_conditions(this)
+      import system_t
+      class(system_t), intent(inout) :: this
+    end subroutine system_initial_conditions
+
+    ! ---------------------------------------------------------
+    !> @brief Execute one operation that is part of a larger algorithm.  Returns true
+    !! if the operation was successfully executed, false otherwise.
+    !!
+    !! Unsuccessful operations can occur, e.g. of quantities from an interaction
+    !! are required, but the interaction is still behind in terms of the iteration counters.
+    !!
+    !! On output, the routine should also provide a list quantities that were
+    !! updated. If no quantitiy was updated, then the corresponding array should
+    !! be left unallocated.
+    logical function system_do_algorithmic_operation(this, operation, updated_quantities) result(done)
+      import system_t
+      import algorithmic_operation_t
+      class(system_t),                intent(inout) :: this
+      class(algorithmic_operation_t), intent(in)    :: operation
+      integer, allocatable,           intent(out)   :: updated_quantities(:)
+    end function system_do_algorithmic_operation
+
+    ! ---------------------------------------------------------
+    !> @brief check whether a system has reached a given tolerance
+    logical function system_is_tolerance_reached(this, tol)
+      use, intrinsic :: iso_fortran_env
+      import system_t
+      class(system_t), intent(in) :: this
+      real(real64),    intent(in) :: tol
+    end function system_is_tolerance_reached
+
+    ! ---------------------------------------------------------
+    !> @brief For some algorithms it might be necessary to store the status of a system at a given algorithmic step
+    !!
+    !! This should be implemented by each system in this routine.
+    subroutine system_store_current_status(this)
+      import system_t
+      class(system_t), intent(inout) :: this
+    end subroutine system_store_current_status
+
+    ! ---------------------------------------------------------
+    subroutine system_restart_write_data(this)
+      import system_t
+      class(system_t), intent(inout) :: this
+    end subroutine system_restart_write_data
+
+    ! ---------------------------------------------------------
+    ! this function returns true if restart data could be read
+    logical function system_restart_read_data(this)
+      import system_t
+      class(system_t), intent(inout) :: this
+    end function system_restart_read_data
+    subroutine system_update_kinetic_energy(this)
+      import system_t
+      class(system_t),      intent(inout) :: this
+    end subroutine system_update_kinetic_energy
+
+  end interface
+
+  !> @brief These classes extends the list and list iterator to create a system list.
+  !!
+  !! Since a list of systems is also a list of interaction partners, the system
+  !! list is an extension of the partner list.
+  type, extends(partner_list_t) :: system_list_t
+    private
+  contains
+    procedure :: add => system_list_add_node !< @copydoc system_oct_m::system_list_add_node
+    procedure :: contains => system_list_contains !< @copydoc system_oct_m::system_list_contains
+  end type system_list_t
+
+  type, extends(linked_list_iterator_t) :: system_iterator_t
+    private
+  contains
+    procedure :: get_next => system_iterator_get_next !< @copydoc system_oct_m::system_iterator_get_next
+  end type system_iterator_t
+
+  ! Subroutine/Functions
+  interface
+    module subroutine system_execute_algorithm(this)
+      class(system_t),     intent(inout) :: this
+    end subroutine system_execute_algorithm
+
+    module subroutine system_reset_iteration_counters(this, accumulated_iterations)
+      class(system_t),      intent(inout) :: this
+      integer,              intent(in)    :: accumulated_iterations
+    end subroutine system_reset_iteration_counters
+
+    recursive module subroutine system_create_interactions(this, interaction_factory, available_partners)
+      class(system_t),                       intent(inout) :: this
+      class(interactions_factory_abst_t),    intent(in)    :: interaction_factory
+      class(partner_list_t),         target, intent(in)    :: available_partners
+    end subroutine system_create_interactions
+
+    module function system_update_couplings(this) result(all_updated)
+      class(system_t),      intent(inout) :: this
+      logical :: all_updated
+    end function system_update_couplings
+
+    module subroutine system_update_interactions(this)
+      class(system_t),      intent(inout) :: this
+    end subroutine system_update_interactions
+
+    module subroutine system_update_interactions_start(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_update_interactions_start
+
+    module subroutine system_update_interactions_finish(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_update_interactions_finish
+
+    module subroutine system_restart_write(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_restart_write
+
+    module function system_restart_read(this) result(res)
+      class(system_t), intent(inout) :: this
+      logical :: res
+    end function system_restart_read
+
+    module subroutine system_output_start(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_output_start
+
+    module subroutine system_output_write(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_output_write
+
+    module subroutine system_output_finish(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_output_finish
+
+    module subroutine system_init_algorithm(this, factory)
+      class(system_t),            intent(inout) :: this
+      class(algorithm_factory_t), intent(in)    :: factory
+    end subroutine system_init_algorithm
+
+    recursive module function system_algorithm_finished(this) result(finished)
+      class(system_t),       intent(in) :: this
+      logical :: finished
+    end function system_algorithm_finished
+
+    module subroutine system_init_iteration_counters(this)
+      class(system_t),            intent(inout) :: this
+    end subroutine system_init_iteration_counters
+
+    module subroutine system_propagation_start(this)
+      class(system_t),      intent(inout) :: this
+    end subroutine system_propagation_start
+
+    module subroutine system_propagation_finish(this)
+      class(system_t),      intent(inout) :: this
+    end subroutine system_propagation_finish
+
+    module subroutine system_iteration_info(this)
+      class(system_t), intent(in) :: this
+    end subroutine system_iteration_info
+
+    module function system_process_is_slave(this) result(res)
+      class(system_t), intent(in) :: this
+      logical :: res
+    end function system_process_is_slave
+
+    module subroutine system_end(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_end
+
+    module subroutine system_list_add_node(this, partner)
+      class(system_list_t)         :: this
+      class(interaction_partner_t), target :: partner
+    end subroutine system_list_add_node
+
+    recursive module function system_list_contains(this, partner) result(contains)
+      class(system_list_t)         :: this
+      class(interaction_partner_t), target :: partner
+      logical contains
+    end function system_list_contains
+
+    module function system_iterator_get_next(this) result(system)
+      class(system_iterator_t), intent(inout) :: this
+      class(system_t),          pointer       :: system
+    end function system_iterator_get_next
+
+    module subroutine system_init_parallelization(this, grp)
+      class(system_t), intent(inout) :: this
+      type(mpi_grp_t), intent(in)    :: grp
+    end subroutine system_init_parallelization
+
+    module subroutine system_start_barrier(this, target_time, barrier_index)
+      class(system_t), intent(inout) :: this
+      real(real64),           intent(in)    :: target_time
+      integer,         intent(in)    :: barrier_index
+    end subroutine system_start_barrier
+
+    module subroutine system_end_barrier(this, barrier_index)
+      class(system_t), intent(inout) :: this
+      integer,         intent(in)    :: barrier_index
+    end subroutine system_end_barrier
+
+    module function system_arrived_at_barrier(this, barrier_index) result(res)
+      class(system_t), intent(inout) :: this
+      integer,         intent(in)    :: barrier_index
+      logical :: res
+    end function system_arrived_at_barrier
+
+    module function system_arrived_at_any_barrier(this) result(res)
+      class(system_t), intent(inout) :: this
+      logical :: res
+    end function system_arrived_at_any_barrier
+
+    module subroutine system_update_potential_energy(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_update_potential_energy
+
+    module subroutine system_update_internal_energy(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_update_internal_energy
+
+    module subroutine system_update_total_energy(this)
+      class(system_t), intent(inout) :: this
+    end subroutine system_update_total_energy
+  end interface
+end module system_oct_m
diff --git a/src/opt_control/opt_control.F90 b/src/opt_control/opt_control.F90
index 4530b3d3c7459f66b4cef85fadf7fe4f5ffb901c..2953a3d09b623888bfa60ffc7781b7d6617962cf 100644
--- a/src/opt_control/opt_control.F90
+++ b/src/opt_control/opt_control.F90
@@ -58,6 +58,7 @@ module opt_control_oct_m
   use electrons_oct_m
   use target_oct_m
   use td_oct_m
+  use td_interface_oct_m
 
   implicit none
 
diff --git a/src/scf/CMakeLists.txt b/src/scf/CMakeLists.txt
index 75f20d3a9894a4e9f6d690f335bc8fd849d8c657..1bb8cb1b2a712a49d1909bac7f1e8fa99133cbb7 100644
--- a/src/scf/CMakeLists.txt
+++ b/src/scf/CMakeLists.txt
@@ -3,13 +3,16 @@ target_sources(Octopus_lib PRIVATE
 		density_criterion.F90
 		eigenval_criterion.F90
 		electrons_ground_state.F90
+		electrons_ground_state_h.F90
 		energy_criterion.F90
 		lcao.F90
 		lda_u_mixer.F90
 		mix.F90
 		mixing_preconditioner.F90
 		rdmft.F90
-		scf.F90
+		scf_h.F90
+		scf_interface.F90
+		scf_interface_h.F90
 		unocc.F90
 		)
 ## Unused sources
diff --git a/src/scf/electrons_ground_state.F90 b/src/scf/electrons_ground_state.F90
index 0e608c8e9c595ad1f257bbac67f108ea7ccf7012..c4a7c70fb1cd07bcd9b2377b6792780ac1c5afc0 100644
--- a/src/scf/electrons_ground_state.F90
+++ b/src/scf/electrons_ground_state.F90
@@ -18,42 +18,31 @@
 
 #include "global.h"
 
-module electrons_ground_state_oct_m
+submodule (electrons_ground_state_oct_m) impl
+  use electrons_ground_state_oct_m
   use debug_oct_m
-  use electron_space_oct_m
   use global_oct_m
-  use grid_oct_m
   use hamiltonian_elec_oct_m
-  use interaction_partner_oct_m
   use io_function_oct_m
-  use ions_oct_m
   use lcao_oct_m
   use math_oct_m
   use mesh_oct_m
   use messages_oct_m
-  use multicomm_oct_m
-  use namespace_oct_m
-  use output_low_oct_m
   use pcm_oct_m
   use rdmft_oct_m
   use restart_oct_m
+  use scf_interface_oct_m
   use scf_oct_m
   use space_oct_m
   use states_abst_oct_m
-  use states_elec_oct_m
   use states_elec_restart_oct_m
-  use v_ks_oct_m
 
   implicit none
 
-  private
-  public ::                       &
-    electrons_ground_state_run
-
 contains
 
   ! ---------------------------------------------------------
-  subroutine electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, fromScratch)
+  module subroutine electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, fromScratch)
     type(namespace_t),        intent(in)    :: namespace
     type(multicomm_t),        intent(in)    :: mc
     type(grid_t),             intent(inout) :: gr
@@ -188,7 +177,7 @@ contains
     POP_SUB(ground_state_run_legacy)
   end subroutine electrons_ground_state_run
 
-end module electrons_ground_state_oct_m
+end submodule impl
 
 !! Local Variables:
 !! mode: f90
diff --git a/src/scf/electrons_ground_state_h.F90 b/src/scf/electrons_ground_state_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..da0b3ac4e1fd2b0a9b0aaa4c486f43b30331f388
--- /dev/null
+++ b/src/scf/electrons_ground_state_h.F90
@@ -0,0 +1,34 @@
+module electrons_ground_state_oct_m
+  use electron_space_oct_m
+  use grid_oct_m
+  use hamiltonian_elec_oct_m
+  use interaction_partner_oct_m
+  use ions_oct_m
+  use multicomm_oct_m
+  use namespace_oct_m
+  use output_low_oct_m
+  use states_elec_oct_m
+  use v_ks_oct_m
+
+  implicit none
+
+  private
+  public ::                       &
+    electrons_ground_state_run
+
+  interface
+    module subroutine electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, fromScratch)
+      type(namespace_t),        intent(in)    :: namespace
+      type(multicomm_t),        intent(in)    :: mc
+      type(grid_t),             intent(inout) :: gr
+      type(ions_t),             intent(inout) :: ions
+      type(partner_list_t),     intent(in)    :: ext_partners
+      type(states_elec_t),      intent(inout) :: st
+      type(v_ks_t),             intent(inout) :: ks
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      type(output_t),           intent(in)    :: outp
+      type(electron_space_t),   intent(in)    :: space
+      logical,                  intent(inout) :: fromScratch
+    end subroutine electrons_ground_state_run
+  end interface
+end module electrons_ground_state_oct_m
diff --git a/src/scf/scf_h.F90 b/src/scf/scf_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..920eda90d05d1a4d61edf9b7d73c80dd3256c7a5
--- /dev/null
+++ b/src/scf/scf_h.F90
@@ -0,0 +1,47 @@
+module scf_oct_m
+  use berry_oct_m
+  use convergence_criterion_oct_m
+  use eigensolver_oct_m
+  use global_oct_m
+  use lda_u_mixer_oct_m
+  use mix_oct_m
+
+  implicit none
+
+  private
+
+  integer, public, parameter :: &
+    VERB_NO      = 0,   &
+    VERB_COMPACT = 1,   &
+    VERB_FULL    = 3
+
+  !> some variables used for the SCF cycle
+  type, public :: scf_t
+    private
+    integer, public :: max_iter   !< maximum number of SCF iterations
+
+    real(real64), public :: lmm_r
+
+    ! several convergence criteria
+    logical, public :: conv_eigen_error
+    logical, public :: check_conv
+
+    integer, public :: mix_field
+    logical, public :: lcao_restricted
+    logical, public :: calc_force
+    logical, public :: calc_stress
+    logical, public :: calc_dipole
+    logical, public :: calc_partial_charges
+    type(mix_t), public :: smix
+    type(mixfield_t), public, pointer :: mixfield
+    type(eigensolver_t), public :: eigens
+    integer, public :: mixdim1
+    logical, public :: forced_finish !< remember if 'touch stop' was triggered earlier.
+    type(lda_u_mixer_t), public :: lda_u_mix
+    type(berry_t), public :: berry
+    integer, public :: matvec !< number matrix-vector products
+
+    type(criterion_list_t), public :: criterion_list
+    real(real64), public :: energy_in, energy_diff, abs_dens_diff, evsum_in, evsum_out, evsum_diff
+  end type scf_t
+end module scf_oct_m
diff --git a/src/scf/scf.F90 b/src/scf/scf_interface.F90
similarity index 96%
rename from src/scf/scf.F90
rename to src/scf/scf_interface.F90
index 2b8eb87c97bdd708c05545cb2037c9e9bc6a662b..0c6fa7cdd4a6bb44c7005495caaf6e2c4e66c4c7 100644
--- a/src/scf/scf.F90
+++ b/src/scf/scf_interface.F90
@@ -18,56 +18,43 @@
 
 #include "global.h"
 
-module scf_oct_m
+submodule (scf_interface_oct_m) impl
+  use scf_interface_oct_m
   use batch_ops_oct_m
   use berry_oct_m
   use convergence_criterion_oct_m
   use criteria_factory_oct_m
   use debug_oct_m
-  use density_oct_m
   use density_criterion_oct_m
+  use density_oct_m
   use eigensolver_oct_m
   use eigenval_criterion_oct_m
-  use electron_space_oct_m
   use energy_calc_oct_m
   use energy_criterion_oct_m
   use forces_oct_m
-  use global_oct_m
-  use grid_oct_m
-  use hamiltonian_elec_oct_m
-  use interaction_partner_oct_m
   use io_oct_m
-  use ions_oct_m
-  use, intrinsic :: iso_fortran_env
   use kpoints_oct_m
+  use lalg_basic_oct_m
   use lcao_oct_m
-  use lda_u_oct_m
   use lda_u_io_oct_m
   use lda_u_mixer_oct_m
-  use lalg_basic_oct_m
+  use lda_u_oct_m
   use loct_oct_m
   use magnetic_oct_m
   use math_oct_m
-  use mesh_oct_m
   use mesh_function_oct_m
+  use mesh_oct_m
   use messages_oct_m
   use mix_oct_m
   use modelmb_exchange_syms_oct_m
   use mpi_oct_m
-  use multicomm_oct_m
-  use namespace_oct_m
-  use output_oct_m
-  use output_low_oct_m
   use output_modelmb_oct_m
+  use output_oct_m
   use parser_oct_m
   use partial_charges_oct_m
   use profiling_oct_m
-  use restart_oct_m
   use smear_oct_m
-  use space_oct_m
   use species_oct_m
-  use states_abst_oct_m
-  use states_elec_oct_m
   use states_elec_io_oct_m
   use states_elec_restart_oct_m
   use stress_oct_m
@@ -76,68 +63,21 @@ module scf_oct_m
   use unit_oct_m
   use unit_system_oct_m
   use utils_oct_m
-  use v_ks_oct_m
   use varinfo_oct_m
   use vdw_ts_oct_m
   use walltimer_oct_m
   use wfs_elec_oct_m
-  use xc_oct_m
   use xc_f03_lib_m
   use xc_interaction_oct_m
+  use xc_oct_m
   use xc_oep_oct_m
   use xc_oep_photon_oct_m
 
   implicit none
-
-  private
-  public ::              &
-    scf_t,               &
-    scf_init,            &
-    scf_mix_clear,       &
-    scf_run,             &
-    scf_end,             &
-    scf_state_info,      &
-    scf_print_mem_use
-
-  integer, public, parameter :: &
-    VERB_NO      = 0,   &
-    VERB_COMPACT = 1,   &
-    VERB_FULL    = 3
-
-  !> some variables used for the SCF cycle
-  type scf_t
-    private
-    integer, public :: max_iter   !< maximum number of SCF iterations
-
-    real(real64), public :: lmm_r
-
-    ! several convergence criteria
-    logical :: conv_eigen_error
-    logical :: check_conv
-
-    integer :: mix_field
-    logical :: lcao_restricted
-    logical :: calc_force
-    logical, public :: calc_stress
-    logical :: calc_dipole
-    logical :: calc_partial_charges
-    type(mix_t) :: smix
-    type(mixfield_t), pointer :: mixfield
-    type(eigensolver_t) :: eigens
-    integer :: mixdim1
-    logical :: forced_finish !< remember if 'touch stop' was triggered earlier.
-    type(lda_u_mixer_t) :: lda_u_mix
-    type(berry_t) :: berry
-    integer :: matvec !< number matrix-vector products
-
-    type(criterion_list_t), public :: criterion_list
-    real(real64) :: energy_in, energy_diff, abs_dens_diff, evsum_in, evsum_out, evsum_diff
-  end type scf_t
-
 contains
 
   ! ---------------------------------------------------------
-  subroutine scf_init(scf, namespace, gr, ions, st, mc, hm, space)
+  module subroutine scf_init(scf, namespace, gr, ions, st, mc, hm, space)
     type(scf_t),              intent(inout) :: scf
     type(grid_t),             intent(in)    :: gr
     type(namespace_t),        intent(in)    :: namespace
@@ -432,7 +372,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine scf_end(scf)
+  module subroutine scf_end(scf)
     type(scf_t),  intent(inout) :: scf
 
     class(convergence_criterion_t), pointer    :: crit
@@ -459,7 +399,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine scf_mix_clear(scf)
+  module subroutine scf_mix_clear(scf)
     type(scf_t), intent(inout) :: scf
 
     PUSH_SUB(scf_mix_clear)
@@ -473,7 +413,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
+  module subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
     verbosity, iters_done, restart_load, restart_dump)
     type(scf_t),               intent(inout) :: scf !< self consistent cycle
     type(namespace_t),         intent(in)    :: namespace
@@ -1319,7 +1259,7 @@ contains
   end subroutine scf_run
 
   ! ---------------------------------------------------------
-  subroutine scf_state_info(namespace, st)
+  module subroutine scf_state_info(namespace, st)
     type(namespace_t),    intent(in) :: namespace
     class(states_abst_t), intent(in) :: st
 
@@ -1337,7 +1277,7 @@ contains
   end subroutine scf_state_info
 
   ! ---------------------------------------------------------
-  subroutine scf_print_mem_use(namespace)
+  module subroutine scf_print_mem_use(namespace)
     type(namespace_t),    intent(in) :: namespace
     real(real64) :: mem
     real(real64) :: mem_tmp
@@ -1419,7 +1359,7 @@ contains
   end subroutine scf_update_diff_quantity
 
 
-end module scf_oct_m
+end submodule impl
 
 
 !! Local Variables:
diff --git a/src/scf/scf_interface_h.F90 b/src/scf/scf_interface_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..6da78e0c57808103dce4ab4983672f4ec98bd91b
--- /dev/null
+++ b/src/scf/scf_interface_h.F90
@@ -0,0 +1,77 @@
+module scf_interface_oct_m
+  use electron_space_oct_m
+  use global_oct_m
+  use grid_oct_m
+  use hamiltonian_elec_oct_m
+  use interaction_partner_oct_m
+  use ions_oct_m
+  use multicomm_oct_m
+  use namespace_oct_m
+  use output_low_oct_m
+  use restart_oct_m
+  use scf_oct_m
+  use space_oct_m
+  use states_abst_oct_m
+  use states_elec_oct_m
+  use v_ks_oct_m
+
+  implicit none
+
+  private
+  public ::              &
+    scf_init,            &
+    scf_mix_clear,       &
+    scf_run,             &
+    scf_end,             &
+    scf_state_info,      &
+    scf_print_mem_use
+
+  interface
+    module subroutine scf_init(scf, namespace, gr, ions, st, mc, hm, space)
+      type(scf_t),              intent(inout) :: scf
+      type(grid_t),             intent(in)    :: gr
+      type(namespace_t),        intent(in)    :: namespace
+      type(ions_t),             intent(in)    :: ions
+      type(states_elec_t),      intent(in)    :: st
+      type(multicomm_t),        intent(in)    :: mc
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      class(space_t),           intent(in)    :: space
+    end subroutine scf_init
+
+    module subroutine scf_mix_clear(scf)
+      type(scf_t), intent(inout) :: scf
+    end subroutine scf_mix_clear
+
+    module subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
+      verbosity, iters_done, restart_load, restart_dump)
+      type(scf_t),               intent(inout) :: scf !< self consistent cycle
+      type(namespace_t),         intent(in)    :: namespace
+      type(electron_space_t),    intent(in)    :: space
+      type(multicomm_t),         intent(in)    :: mc
+      type(grid_t),              intent(inout) :: gr !< grid
+      type(ions_t),              intent(inout) :: ions !< geometry
+      type(partner_list_t),      intent(in)    :: ext_partners
+      type(states_elec_t),       intent(inout) :: st !< States
+      type(v_ks_t),              intent(inout) :: ks !< Kohn-Sham
+      type(hamiltonian_elec_t),  intent(inout) :: hm !< Hamiltonian
+      type(output_t),  optional, intent(in)    :: outp
+      integer,         optional, intent(in)    :: verbosity
+      integer,         optional, intent(out)   :: iters_done
+      type(restart_t), optional, intent(in)    :: restart_load
+      type(restart_t), optional, intent(in)    :: restart_dump
+    end subroutine scf_run
+
+    module subroutine scf_end(scf)
+      type(scf_t),  intent(inout) :: scf
+    end subroutine scf_end
+
+    module subroutine scf_state_info(namespace, st)
+      type(namespace_t),    intent(in) :: namespace
+      class(states_abst_t), intent(in) :: st
+    end subroutine scf_state_info
+
+    module subroutine scf_print_mem_use(namespace)
+      type(namespace_t),    intent(in) :: namespace
+    end subroutine scf_print_mem_use
+  end interface
+end module scf_interface_oct_m
diff --git a/src/scf/unocc.F90 b/src/scf/unocc.F90
index 6373a84cd622003b3778379334594cf38fe87650..9318ddf10363635688525ccd6ac8a88c7be9d110 100644
--- a/src/scf/unocc.F90
+++ b/src/scf/unocc.F90
@@ -40,6 +40,7 @@ module unocc_oct_m
   use parser_oct_m
   use profiling_oct_m
   use restart_oct_m
+  use scf_interface_oct_m
   use scf_oct_m
   use space_oct_m
   use states_abst_oct_m
diff --git a/src/td/CMakeLists.txt b/src/td/CMakeLists.txt
index 7d266368fd1e19334a5d24526599af309d6f61cf..819f0ca6484e51fb94c066423a30c4e765e13be8 100644
--- a/src/td/CMakeLists.txt
+++ b/src/td/CMakeLists.txt
@@ -9,14 +9,17 @@ target_sources(Octopus_lib PRIVATE
 		propagator_base.F90
 		propagator_cn.F90
 		propagator_elec.F90
+		propagator_elec_h.F90
 		propagator_etrs.F90
 		propagator_expmid.F90
 		propagator_magnus.F90
 		propagator_qoct.F90
 		propagator_rk.F90
 		spectrum.F90
-		td.F90
 		td_calc.F90
+		td_h.F90
+		td_interface.F90
+		td_interface_h.F90
 		td_write.F90
 		td_write_low.F90
 		)
diff --git a/src/td/propagator_elec.F90 b/src/td/propagator_elec.F90
index 68f9f50b48976855ee3ff72042633c7d96817ca2..20cd7a8fc2a471c5f7460d23fbb11a95650fcc86 100644
--- a/src/td/propagator_elec.F90
+++ b/src/td/propagator_elec.F90
@@ -18,33 +18,21 @@
 
 #include "global.h"
 
-module propagator_elec_oct_m
+submodule (propagator_elec_oct_m) impl
+  use propagator_elec_oct_m
   use debug_oct_m
-  use electron_space_oct_m
   use energy_calc_oct_m
   use exponential_oct_m
   use ext_partner_list_oct_m
   use forces_oct_m
   use gauge_field_oct_m
-  use grid_oct_m
-  use global_oct_m
-  use hamiltonian_elec_oct_m
-  use interaction_partner_oct_m
-  use ion_dynamics_oct_m
-  use ions_oct_m
-  use, intrinsic :: iso_fortran_env
   use lasers_oct_m
   use lda_u_oct_m
   use parser_oct_m
   use mesh_function_oct_m
   use messages_oct_m
-  use multicomm_oct_m
-  use namespace_oct_m
-  use opt_control_state_oct_m
-  use output_low_oct_m
   use potential_interpolation_oct_m
   use profiling_oct_m
-  use propagator_base_oct_m
   use propagator_cn_oct_m
   use propagator_etrs_oct_m
   use propagator_expmid_oct_m
@@ -52,34 +40,17 @@ module propagator_elec_oct_m
   use propagator_qoct_oct_m
   use propagator_rk_oct_m
   use propagator_verlet_oct_m
-  use scf_oct_m
+  use scf_interface_oct_m
   use sparskit_oct_m
-  use space_oct_m
-  use states_elec_oct_m
   use stress_oct_m
-  use td_write_oct_m
-  use v_ks_oct_m
   use varinfo_oct_m
   use xc_oct_m
-
   implicit none
 
-  private
-  public ::                         &
-    propagator_elec_init,                &
-    propagator_elec_end,                 &
-    propagator_elec_copy,                &
-    propagator_elec_run_zero_iter,       &
-    propagator_elec_dt,                  &
-    propagator_elec_set_scf_prop,        &
-    propagator_elec_remove_scf_prop,     &
-    propagator_elec_ions_are_propagated, &
-    propagator_elec_dt_bo
-
 contains
 
   ! ---------------------------------------------------------
-  subroutine propagator_elec_copy(tro, tri)
+  module subroutine propagator_elec_copy(tro, tri)
     type(propagator_base_t), intent(inout) :: tro
     type(propagator_base_t), intent(in)    :: tri
 
@@ -117,7 +88,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine propagator_elec_init(gr, namespace, st, tr, have_fields, family_is_mgga_with_exc)
+  module subroutine propagator_elec_init(gr, namespace, st, tr, have_fields, family_is_mgga_with_exc)
     type(grid_t),            intent(in)    :: gr
     type(namespace_t),       intent(in)    :: namespace
     type(states_elec_t),     intent(in)    :: st
@@ -389,7 +360,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine propagator_elec_set_scf_prop(tr, threshold)
+  module subroutine propagator_elec_set_scf_prop(tr, threshold)
     type(propagator_base_t), intent(inout) :: tr
     real(real64), optional,         intent(in)    :: threshold
 
@@ -406,7 +377,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine propagator_elec_remove_scf_prop(tr)
+  module subroutine propagator_elec_remove_scf_prop(tr)
     type(propagator_base_t), intent(inout) :: tr
 
     PUSH_SUB(propagator_elec_remove_scf_prop)
@@ -419,7 +390,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine propagator_elec_end(tr)
+  module subroutine propagator_elec_end(tr)
     type(propagator_base_t), intent(inout) :: tr
 
     PUSH_SUB(propagator_elec_end)
@@ -442,7 +413,7 @@ contains
 
 
   ! ---------------------------------------------------------
-  subroutine propagator_elec_run_zero_iter(hm, gr, tr)
+  module subroutine propagator_elec_run_zero_iter(hm, gr, tr)
     type(hamiltonian_elec_t),  intent(in)    :: hm
     type(grid_t),              intent(in)    :: gr
     type(propagator_base_t),   intent(inout) :: tr
@@ -459,7 +430,7 @@ contains
   !> Propagates st from time - dt to t.
   !! If dt<0, it propagates *backwards* from t+|dt| to t
   ! ---------------------------------------------------------
-  subroutine propagator_elec_dt(ks, namespace, space, hm, gr, st, tr, time, dt, nt, &
+  module subroutine propagator_elec_dt(ks, namespace, space, hm, gr, st, tr, time, dt, nt, &
     ions_dyn, ions, ext_partners, outp, write_handler, scsteps, update_energy, qcchi)
     type(v_ks_t),                        target, intent(inout) :: ks
     type(namespace_t),                           intent(in)    :: namespace
@@ -614,8 +585,9 @@ contains
 
 
   ! ---------------------------------------------------------
-  logical pure function propagator_elec_ions_are_propagated(tr) result(propagated)
+  pure module function propagator_elec_ions_are_propagated(tr) result(propagated)
     type(propagator_base_t), intent(in) :: tr
+    logical :: propagated
 
     select case (tr%method)
     case (PROP_ETRS, PROP_AETRS, PROP_CAETRS, PROP_EXPLICIT_RUNGE_KUTTA4)
@@ -630,7 +602,8 @@ contains
 
   ! ---------------------------------------------------------
 
-  subroutine propagator_elec_dt_bo(scf, namespace, space, gr, ks, st, hm, ions, ext_partners, mc, outp, iter, dt, ions_dyn, scsteps)
+  module subroutine propagator_elec_dt_bo(scf, namespace, space, gr, ks, st, hm, ions, ext_partners, &
+    mc, outp, iter, dt, ions_dyn, scsteps)
     type(scf_t),              intent(inout) :: scf
     type(namespace_t),        intent(in)    :: namespace
     type(electron_space_t),   intent(in)    :: space
@@ -692,7 +665,7 @@ contains
     POP_SUB(propagator_elec_dt_bo)
   end subroutine propagator_elec_dt_bo
 
-end module propagator_elec_oct_m
+end submodule impl
 
 
 !! Local Variables:
diff --git a/src/td/propagator_elec_h.F90 b/src/td/propagator_elec_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..bc6db9574c9e604ce51c5cb46e70364698bd7113
--- /dev/null
+++ b/src/td/propagator_elec_h.F90
@@ -0,0 +1,114 @@
+module propagator_elec_oct_m
+  use electron_space_oct_m
+  use global_oct_m
+  use grid_oct_m
+  use hamiltonian_elec_oct_m
+  use interaction_partner_oct_m
+  use ion_dynamics_oct_m
+  use ions_oct_m
+  use multicomm_oct_m
+  use namespace_oct_m
+  use opt_control_state_oct_m
+  use output_low_oct_m
+  use propagator_base_oct_m
+  use scf_oct_m
+  use space_oct_m
+  use states_elec_oct_m
+  use td_write_oct_m
+  use v_ks_oct_m
+
+  implicit none
+
+  private
+  public ::                         &
+    propagator_elec_init,                &
+    propagator_elec_end,                 &
+    propagator_elec_copy,                &
+    propagator_elec_run_zero_iter,       &
+    propagator_elec_dt,                  &
+    propagator_elec_set_scf_prop,        &
+    propagator_elec_remove_scf_prop,     &
+    propagator_elec_ions_are_propagated, &
+    propagator_elec_dt_bo
+
+  interface
+    module subroutine propagator_elec_copy(tro, tri)
+      type(propagator_base_t), intent(inout) :: tro
+      type(propagator_base_t), intent(in)    :: tri
+    end subroutine propagator_elec_copy
+
+    module subroutine propagator_elec_init(gr, namespace, st, tr, have_fields, family_is_mgga_with_exc)
+      type(grid_t),            intent(in)    :: gr
+      type(namespace_t),       intent(in)    :: namespace
+      type(states_elec_t),     intent(in)    :: st
+      type(propagator_base_t), intent(inout) :: tr
+      logical,                 intent(in)    :: have_fields
+      logical,                 intent(in)    :: family_is_mgga_with_exc
+    end subroutine propagator_elec_init
+
+    module subroutine propagator_elec_set_scf_prop(tr, threshold)
+      type(propagator_base_t), intent(inout) :: tr
+      real(real64), optional,  intent(in)    :: threshold
+    end subroutine propagator_elec_set_scf_prop
+
+    module subroutine propagator_elec_remove_scf_prop(tr)
+      type(propagator_base_t), intent(inout) :: tr
+    end subroutine propagator_elec_remove_scf_prop
+
+    module subroutine propagator_elec_end(tr)
+      type(propagator_base_t), intent(inout) :: tr
+    end subroutine propagator_elec_end
+
+    module subroutine propagator_elec_run_zero_iter(hm, gr, tr)
+      type(hamiltonian_elec_t),  intent(in)    :: hm
+      type(grid_t),              intent(in)    :: gr
+      type(propagator_base_t),   intent(inout) :: tr
+    end subroutine propagator_elec_run_zero_iter
+
+    module subroutine propagator_elec_dt(ks, namespace, space, hm, gr, st, tr, time, dt, nt, &
+      ions_dyn, ions, ext_partners, outp, write_handler, scsteps, update_energy, qcchi)
+      type(v_ks_t),                        target, intent(inout) :: ks
+      type(namespace_t),                           intent(in)    :: namespace
+      type(electron_space_t),                      intent(in)    :: space
+      type(hamiltonian_elec_t),            target, intent(inout) :: hm
+      type(grid_t),                        target, intent(inout) :: gr
+      type(states_elec_t),                 target, intent(inout) :: st
+      type(propagator_base_t),             target, intent(inout) :: tr
+      real(real64),                                intent(in)    :: time
+      real(real64),                                intent(in)    :: dt
+      integer,                                     intent(in)    :: nt
+      type(ion_dynamics_t),                        intent(inout) :: ions_dyn
+      type(ions_t),                                intent(inout) :: ions
+      type(partner_list_t),                        intent(in)    :: ext_partners
+      type(output_t),                              intent(in)    :: outp
+      type(td_write_t),                            intent(in)    :: write_handler
+      integer,                   optional,         intent(out)   :: scsteps
+      logical,                   optional,         intent(in)    :: update_energy
+      type(opt_control_state_t), optional, target, intent(inout) :: qcchi
+    end subroutine propagator_elec_dt
+
+    pure module function propagator_elec_ions_are_propagated(tr) result(propagated)
+      type(propagator_base_t), intent(in) :: tr
+      logical :: propagated
+    end function propagator_elec_ions_are_propagated
+
+    module subroutine propagator_elec_dt_bo(scf, namespace, space, gr, ks, st, hm, ions, ext_partners, &
+      mc, outp, iter, dt, ions_dyn, scsteps)
+      type(scf_t),              intent(inout) :: scf
+      type(namespace_t),        intent(in)    :: namespace
+      type(electron_space_t),   intent(in)    :: space
+      type(grid_t),             intent(inout) :: gr
+      type(v_ks_t),             intent(inout) :: ks
+      type(states_elec_t),      intent(inout) :: st
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      type(ions_t),             intent(inout) :: ions
+      type(partner_list_t),     intent(in)    :: ext_partners
+      type(multicomm_t),        intent(inout) :: mc    !< index and domain communicators
+      type(output_t),           intent(inout) :: outp
+      integer,                  intent(in)    :: iter
+      real(real64),             intent(in)    :: dt
+      type(ion_dynamics_t),     intent(inout) :: ions_dyn
+      integer,                  intent(inout) :: scsteps
+    end subroutine propagator_elec_dt_bo
+  end interface
+end module propagator_elec_oct_m
diff --git a/src/td/td_h.F90 b/src/td/td_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..1867f3a53bd8afe573c60ef66070ad2e75b64477
--- /dev/null
+++ b/src/td/td_h.F90
@@ -0,0 +1,56 @@
+module td_oct_m
+  use electron_space_oct_m
+  use global_oct_m
+  use grid_oct_m
+  use hamiltonian_elec_oct_m
+  use interaction_partner_oct_m
+  use ion_dynamics_oct_m
+  use ions_oct_m
+  use multicomm_oct_m
+  use namespace_oct_m
+  use output_low_oct_m
+  use pes_oct_m
+  use propagator_base_oct_m
+  use restart_oct_m
+  use scf_oct_m
+  use space_oct_m
+  use states_elec_oct_m
+  use td_write_oct_m
+  use v_ks_oct_m
+
+  implicit none
+
+  private
+
+  !> Parameters.
+  integer, parameter, public :: &
+    EHRENFEST = 1,   &
+    BO        = 2
+
+  type, public :: td_t
+    private
+    type(propagator_base_t), public :: tr             !< contains the details of the time-evolution
+    type(scf_t),             public :: scf
+    type(ion_dynamics_t),    public :: ions_dyn
+    real(real64),            public :: dt             !< time step
+    integer,                 public :: max_iter       !< maximum number of iterations to perform
+    integer,                 public :: iter           !< the actual iteration
+    logical,                 public :: recalculate_gs !< Recalculate ground-state along the evolution.
+
+    type(pes_t),             public :: pesv
+
+    integer,                 public :: dynamics
+    integer,                 public :: energy_update_iter
+    real(real64),            public :: scissor
+
+    logical                         :: freeze_occ
+    logical                         :: freeze_u
+    integer,                 public :: freeze_orbitals
+
+    logical,                 public :: from_scratch = .false.
+
+    type(td_write_t), public        :: write_handler
+    type(restart_t),         public :: restart_load
+    type(restart_t),         public :: restart_dump
+  end type td_t
+end module td_oct_m
diff --git a/src/td/td.F90 b/src/td/td_interface.F90
similarity index 93%
rename from src/td/td.F90
rename to src/td/td_interface.F90
index 7f267f0a56d110f35865a393e6ae1c845d1bb52d..2a3e95e6ce969dfff29b12939a5e4f98fe55db77 100644
--- a/src/td/td.F90
+++ b/src/td/td_interface.F90
@@ -18,43 +18,35 @@
 
 #include "global.h"
 
-module td_oct_m
+submodule (td_interface_oct_m) impl
+  use td_interface_oct_m
   use absorbing_boundaries_oct_m
   use boundaries_oct_m
   use calc_mode_par_oct_m
-  use current_oct_m
   use classical_particle_oct_m
+  use current_oct_m
   use debug_oct_m
   use density_oct_m
-  use energy_calc_oct_m
   use electrons_ground_state_oct_m
-  use electron_space_oct_m
+  use energy_calc_oct_m
   use epot_oct_m
   use ext_partner_list_oct_m
   use forces_oct_m
   use gauge_field_oct_m
   use global_oct_m
-  use grid_oct_m
-  use hamiltonian_elec_oct_m
-  use interaction_partner_oct_m
   use io_oct_m
   use ion_dynamics_oct_m
-  use ions_oct_m
   use kick_oct_m
-  use, intrinsic :: iso_fortran_env
   use lasers_oct_m
-  use lda_u_oct_m
   use lda_u_io_oct_m
+  use lda_u_oct_m
   use linked_list_oct_m
   use loct_oct_m
   use maxwell_boundary_op_oct_m
   use mesh_oct_m
   use messages_oct_m
   use mpi_oct_m
-  use multicomm_oct_m
-  use namespace_oct_m
   use output_oct_m
-  use output_low_oct_m
   use parser_oct_m
   use pes_oct_m
   use photon_mode_mf_oct_m
@@ -62,84 +54,28 @@ module td_oct_m
   use poisson_oct_m
   use potential_interpolation_oct_m
   use profiling_oct_m
-  use propagator_oct_m
   use propagator_elec_oct_m
-  use propagator_base_oct_m
+  use propagator_oct_m
   use restart_oct_m
+  use scf_interface_oct_m
   use scf_oct_m
   use scissor_oct_m
-  use space_oct_m
   use states_abst_oct_m
-  use states_elec_oct_m
   use states_elec_restart_oct_m
   use stress_oct_m
   use td_write_oct_m
   use types_oct_m
   use unit_oct_m
   use unit_system_oct_m
-  use v_ks_oct_m
   use varinfo_oct_m
   use walltimer_oct_m
   use xc_oct_m
 
   implicit none
 
-  private
-  public ::                    &
-    td_t,                      &
-    td_run,                    &
-    td_run_init,               &
-    td_init,                   &
-    td_init_run,               &
-    td_end,                    &
-    td_end_run,                &
-    td_write_iter,             &
-    td_check_point,            &
-    td_dump,                   &
-    td_allocate_wavefunctions, &
-    td_init_gaugefield,        &
-    td_load_restart_from_gs,   &
-    td_load_restart_from_td,   &
-    td_init_with_wavefunctions,&
-    td_get_from_scratch,       &
-    td_set_from_scratch
-
-  !> Parameters.
-  integer, parameter, public :: &
-    EHRENFEST = 1,   &
-    BO        = 2
-
-  type td_t
-    private
-    type(propagator_base_t), public :: tr             !< contains the details of the time-evolution
-    type(scf_t),             public :: scf
-    type(ion_dynamics_t),    public :: ions_dyn
-    real(real64),            public :: dt             !< time step
-    integer,                 public :: max_iter       !< maximum number of iterations to perform
-    integer,                 public :: iter           !< the actual iteration
-    logical,                 public :: recalculate_gs !< Recalculate ground-state along the evolution.
-
-    type(pes_t),             public :: pesv
-
-    integer,                 public :: dynamics
-    integer,                 public :: energy_update_iter
-    real(real64)                    :: scissor
-
-    logical                         :: freeze_occ
-    logical                         :: freeze_u
-    integer                         :: freeze_orbitals
-
-    logical                         :: from_scratch = .false.
-
-    type(td_write_t), public        :: write_handler
-    type(restart_t)                 :: restart_load
-    type(restart_t)                 :: restart_dump
-  end type td_t
-
-
 contains
 
-  subroutine td_run_init()
+  module subroutine td_run_init()
 
     PUSH_SUB(td_run_init)
 
@@ -150,7 +86,7 @@ contains
 
   ! ---------------------------------------------------------
 
-  subroutine td_init(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp)
+  module subroutine td_init(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp)
     type(td_t),               intent(inout) :: td
     type(namespace_t),        intent(in)    :: namespace
     class(space_t),           intent(in)    :: space
@@ -396,7 +332,7 @@ contains
   end subroutine td_init
 
   ! ---------------------------------------------------------
-  subroutine td_init_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
+  module subroutine td_init_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
     type(td_t),               intent(inout) :: td
     type(namespace_t),        intent(in)    :: namespace
     type(multicomm_t),        intent(inout) :: mc
@@ -448,7 +384,7 @@ contains
   end subroutine td_init_run
 
   ! ---------------------------------------------------------
-  subroutine td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
+  module subroutine td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
     type(td_t),               intent(inout) :: td
     type(namespace_t),        intent(in)    :: namespace
     type(multicomm_t),        intent(inout) :: mc
@@ -482,7 +418,7 @@ contains
   end subroutine td_allocate_wavefunctions
 
   ! ---------------------------------------------------------
-  subroutine td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
+  module subroutine td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
     type(td_t),               intent(inout) :: td
     type(namespace_t),        intent(in)    :: namespace
     type(grid_t),             intent(inout) :: gr
@@ -513,7 +449,7 @@ contains
   end subroutine td_init_gaugefield
 
   ! ---------------------------------------------------------
-  subroutine td_end(td)
+  module subroutine td_end(td)
     type(td_t),               intent(inout) :: td
 
     PUSH_SUB(td_end)
@@ -528,7 +464,7 @@ contains
   end subroutine td_end
 
   ! ---------------------------------------------------------
-  subroutine td_end_run(td, st, hm)
+  module subroutine td_end_run(td, st, hm)
     type(td_t),               intent(inout) :: td
     type(states_elec_t),      intent(inout) :: st
     type(hamiltonian_elec_t), intent(inout) :: hm
@@ -548,7 +484,7 @@ contains
   end subroutine td_end_run
 
   ! ---------------------------------------------------------
-  subroutine td_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
+  module subroutine td_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
     type(td_t),               intent(inout) :: td
     type(namespace_t),        intent(in)    :: namespace
     type(multicomm_t),        intent(inout) :: mc
@@ -649,7 +585,7 @@ contains
   end subroutine td_print_header
 
   ! ---------------------------------------------------------
-  subroutine td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
+  module subroutine td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
     iter, scsteps, etime, stopping, from_scratch)
     type(td_t),               intent(inout) :: td
     type(namespace_t),        intent(in)    :: namespace
@@ -745,7 +681,7 @@ contains
   end subroutine td_update_elapsed_time
 
   ! ---------------------------------------------------------
-  subroutine td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
+  module subroutine td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
     type(td_t),                  intent(inout) :: td
     type(namespace_t),           intent(in)    :: namespace
     type(electron_space_t),      intent(in)    :: space
@@ -1021,7 +957,7 @@ contains
   end subroutine td_init_ions_and_forces
 
   ! ---------------------------------------------------------
-  subroutine td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
+  module subroutine td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
     type(td_t),                  intent(inout) :: td
     type(namespace_t),           intent(in)    :: namespace
     class(space_t),              intent(in)    :: space
@@ -1057,7 +993,7 @@ contains
   end subroutine td_load_restart_from_td
 
   ! ---------------------------------------------------------
-  subroutine td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
+  module subroutine td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
     type(td_t),                  intent(inout) :: td
     type(namespace_t),           intent(in)    :: namespace
     class(space_t),              intent(in)    :: space
@@ -1189,7 +1125,7 @@ contains
   end subroutine td_read_coordinates
 
   ! ---------------------------------------------------------
-  subroutine td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
+  module subroutine td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
     type(td_t),               intent(in)  :: td
     type(namespace_t),        intent(in)  :: namespace
     class(space_t),           intent(in)  :: space
@@ -1413,18 +1349,19 @@ contains
   end subroutine td_load_frozen
 
   ! ---------------------------------------------------------
-  logical function td_get_from_scratch(td)
+  module function td_get_from_scratch(td) result(res)
     type(td_t), intent(in) :: td
+    logical :: res
 
     PUSH_SUB(td_get_from_scratch)
 
-    td_get_from_scratch = td%from_scratch
+    res = td%from_scratch
 
     POP_SUB(td_get_from_scratch)
   end function td_get_from_scratch
 
   ! ---------------------------------------------------------
-  subroutine td_set_from_scratch(td, from_scratch)
+  module subroutine td_set_from_scratch(td, from_scratch)
     type(td_t), intent(inout) :: td
     logical,    intent(in)    :: from_scratch
 
@@ -1434,7 +1371,7 @@ contains
 
     POP_SUB(td_set_from_scratch)
   end subroutine td_set_from_scratch
-end module td_oct_m
+end submodule impl
 
 !! Local Variables:
 !! mode: f90
diff --git a/src/td/td_interface_h.F90 b/src/td/td_interface_h.F90
new file mode 100644
index 0000000000000000000000000000000000000000..86f99e6333f59388fe58d27bf3bd8728bd9c16a2
--- /dev/null
+++ b/src/td/td_interface_h.F90
@@ -0,0 +1,197 @@
+module td_interface_oct_m
+  use electron_space_oct_m
+  use grid_oct_m
+  use hamiltonian_elec_oct_m
+  use interaction_partner_oct_m
+  use ions_oct_m
+  use multicomm_oct_m
+  use namespace_oct_m
+  use output_low_oct_m
+  use space_oct_m
+  use states_elec_oct_m
+  use td_oct_m
+  use v_ks_oct_m
+  use, intrinsic :: iso_fortran_env
+
+  private
+  public ::                    &
+    td_run,                    &
+    td_run_init,               &
+    td_init,                   &
+    td_init_run,               &
+    td_end,                    &
+    td_end_run,                &
+    td_check_point,            &
+    td_dump,                   &
+    td_allocate_wavefunctions, &
+    td_init_gaugefield,        &
+    td_load_restart_from_gs,   &
+    td_load_restart_from_td,   &
+    td_init_with_wavefunctions,&
+    td_get_from_scratch,       &
+    td_set_from_scratch
+
+  ! Subroutine/Functions
+  interface
+    module subroutine td_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
+      type(td_t),               intent(inout) :: td
+      type(namespace_t),        intent(in)    :: namespace
+      type(multicomm_t),        intent(inout) :: mc
+      type(grid_t),             intent(inout) :: gr
+      type(ions_t),             intent(inout) :: ions
+      type(states_elec_t),      intent(inout) :: st
+      type(v_ks_t),             intent(inout) :: ks
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      type(partner_list_t),     intent(in)    :: ext_partners
+      type(output_t),           intent(inout) :: outp
+      type(electron_space_t),   intent(in)    :: space
+      logical,                  intent(inout) :: from_scratch
+    end subroutine td_run
+
+    module subroutine td_run_init()
+    end subroutine td_run_init
+
+    module subroutine td_init(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp)
+      type(td_t),               intent(inout) :: td
+      type(namespace_t),        intent(in)    :: namespace
+      class(space_t),           intent(in)    :: space
+      type(grid_t),             intent(in)    :: gr
+      type(ions_t),             intent(inout) :: ions
+      type(states_elec_t),      intent(in)    :: st
+      type(v_ks_t),             intent(in)    :: ks
+      type(hamiltonian_elec_t), intent(in)    :: hm
+      type(partner_list_t),     intent(in)    :: ext_partners
+      type(output_t),           intent(in)    :: outp
+    end subroutine td_init
+
+    module subroutine td_init_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
+      type(td_t),               intent(inout) :: td
+      type(namespace_t),        intent(in)    :: namespace
+      type(multicomm_t),        intent(inout) :: mc
+      type(grid_t),             intent(inout) :: gr
+      type(ions_t),             intent(inout) :: ions
+      type(states_elec_t),      intent(inout) :: st
+      type(v_ks_t),             intent(inout) :: ks
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      type(partner_list_t),     intent(in)    :: ext_partners
+      type(output_t),           intent(inout) :: outp
+      type(electron_space_t),   intent(in)    :: space
+      logical,                  intent(inout) :: from_scratch
+    end subroutine td_init_run
+
+    module subroutine td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
+      type(td_t),               intent(inout) :: td
+      type(namespace_t),        intent(in)    :: namespace
+      type(multicomm_t),        intent(inout) :: mc
+      type(grid_t),             intent(inout) :: gr
+      type(ions_t),             intent(inout) :: ions
+      type(states_elec_t),      intent(inout) :: st
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      class(space_t),           intent(in)    :: space
+    end subroutine td_allocate_wavefunctions
+
+    module subroutine td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
+      type(td_t),               intent(inout) :: td
+      type(namespace_t),        intent(in)    :: namespace
+      type(grid_t),             intent(inout) :: gr
+      type(states_elec_t),      intent(inout) :: st
+      type(v_ks_t),             intent(inout) :: ks
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      type(partner_list_t),     intent(in)    :: ext_partners
+      class(space_t),           intent(in)    :: space
+    end subroutine td_init_gaugefield
+
+    module subroutine td_end(td)
+      type(td_t),               intent(inout) :: td
+    end subroutine td_end
+
+    module subroutine td_end_run(td, st, hm)
+      type(td_t),               intent(inout) :: td
+      type(states_elec_t),      intent(inout) :: st
+      type(hamiltonian_elec_t), intent(inout) :: hm
+    end subroutine td_end_run
+
+    module subroutine td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
+      iter, scsteps, etime, stopping, from_scratch)
+      type(td_t),               intent(inout) :: td
+      type(namespace_t),        intent(in)    :: namespace
+      type(multicomm_t),        intent(in)    :: mc
+      type(grid_t),             intent(inout) :: gr
+      type(ions_t),             intent(inout) :: ions
+      type(states_elec_t),      intent(inout) :: st
+      type(v_ks_t),             intent(inout) :: ks
+      type(hamiltonian_elec_t), intent(inout) :: hm
+      type(partner_list_t),     intent(in)    :: ext_partners
+      type(output_t),           intent(in)    :: outp
+      type(electron_space_t),   intent(in)    :: space
+      integer,                  intent(in)    :: iter
+      integer,                  intent(in)    :: scsteps
+      real(real64),             intent(inout) :: etime
+      logical,                  intent(in)    :: stopping
+      logical,                  intent(inout) :: from_scratch
+    end subroutine td_check_point
+
+    module subroutine td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
+      type(td_t),                  intent(inout) :: td
+      type(namespace_t),           intent(in)    :: namespace
+      type(electron_space_t),      intent(in)    :: space
+      type(multicomm_t),           intent(in)    :: mc
+      type(grid_t),                intent(inout) :: gr
+      type(ions_t),                intent(inout) :: ions
+      type(partner_list_t),        intent(in)    :: ext_partners
+      type(states_elec_t), target, intent(inout) :: st
+      type(v_ks_t),                intent(inout) :: ks
+      type(hamiltonian_elec_t),    intent(inout) :: hm
+      type(output_t),              intent(inout) :: outp
+      logical,                     intent(in)    :: from_scratch
+    end subroutine td_init_with_wavefunctions
+
+    module subroutine td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
+      type(td_t),                  intent(inout) :: td
+      type(namespace_t),           intent(in)    :: namespace
+      class(space_t),              intent(in)    :: space
+      type(multicomm_t),           intent(in)    :: mc
+      type(grid_t),                intent(inout) :: gr
+      type(partner_list_t),        intent(in)    :: ext_partners
+      type(states_elec_t), target, intent(inout) :: st
+      type(v_ks_t),                intent(inout) :: ks
+      type(hamiltonian_elec_t),    intent(inout) :: hm
+    end subroutine td_load_restart_from_gs
+
+    module subroutine td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
+      type(td_t),               intent(in)  :: td
+      type(namespace_t),        intent(in)  :: namespace
+      class(space_t),           intent(in)  :: space
+      type(grid_t),             intent(in)  :: gr
+      type(states_elec_t),      intent(in)  :: st
+      type(hamiltonian_elec_t), intent(in)  :: hm
+      type(v_ks_t),             intent(in)  :: ks
+      type(partner_list_t),     intent(in)  :: ext_partners
+      integer,                  intent(in)  :: iter
+      integer,                  intent(out) :: ierr
+    end subroutine td_dump
+
+    module function td_get_from_scratch(td) result(res)
+      type(td_t), intent(in) :: td
+      logical :: res
+    end function td_get_from_scratch
+
+    module subroutine td_set_from_scratch(td, from_scratch)
+      type(td_t), intent(inout) :: td
+      logical,    intent(in)    :: from_scratch
+    end subroutine td_set_from_scratch
+
+    module subroutine td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
+      type(td_t),                  intent(inout) :: td
+      type(namespace_t),           intent(in)    :: namespace
+      class(space_t),              intent(in)    :: space
+      type(multicomm_t),           intent(in)    :: mc
+      type(grid_t),                intent(inout) :: gr
+      type(partner_list_t),        intent(in)    :: ext_partners
+      type(states_elec_t), target, intent(inout) :: st
+      type(v_ks_t),                intent(inout) :: ks
+      type(hamiltonian_elec_t),    intent(inout) :: hm
+      logical,                     intent(inout) :: from_scratch
+    end subroutine td_load_restart_from_td
+  end interface
+end module td_interface_oct_m