diff --git a/src/core/PropagationMode.hpp b/src/core/PropagationMode.hpp
index ef8e029e81..ff87400f10 100644
--- a/src/core/PropagationMode.hpp
+++ b/src/core/PropagationMode.hpp
@@ -41,3 +41,12 @@ enum PropagationMode : int {
   ROT_STOKESIAN = 1 << 14,
 };
 } // namespace PropagationMode
+
+/** \name Integrator switches */
+enum IntegratorSwitch : int {
+  INTEG_METHOD_NPT_ISO = 0,
+  INTEG_METHOD_NVT = 1,
+  INTEG_METHOD_STEEPEST_DESCENT = 2,
+  INTEG_METHOD_BD = 3,
+  INTEG_METHOD_SD = 7,
+};
diff --git a/src/core/accumulators/Correlator.cpp b/src/core/accumulators/Correlator.cpp
index cbc1c453a0..fccadbd62c 100644
--- a/src/core/accumulators/Correlator.cpp
+++ b/src/core/accumulators/Correlator.cpp
@@ -18,8 +18,6 @@
  */
 #include "Correlator.hpp"
 
-#include "integrate.hpp"
-
 #include <utils/Vector.hpp>
 #include <utils/math/sqr.hpp>
 #include <utils/serialization/multi_array.hpp>
diff --git a/src/core/accumulators/Correlator.hpp b/src/core/accumulators/Correlator.hpp
index 4258871922..e60fa27d40 100644
--- a/src/core/accumulators/Correlator.hpp
+++ b/src/core/accumulators/Correlator.hpp
@@ -101,8 +101,8 @@
  */
 
 #include "AccumulatorBase.hpp"
-#include "integrate.hpp"
 #include "observables/Observable.hpp"
+#include "system/System.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -155,11 +155,12 @@ class Correlator : public AccumulatorBase {
              obs_ptr obs2, Utils::Vector3d correlation_args_ = {})
       : AccumulatorBase(delta_N), finalized(false), t(0),
         m_correlation_args(correlation_args_), m_tau_lin(tau_lin),
-        m_dt(delta_N * get_time_step()), m_tau_max(tau_max),
+        m_dt(::System::get_system().get_time_step()), m_tau_max(tau_max),
         compressA_name(std::move(compress1_)),
         compressB_name(std::move(compress2_)),
         corr_operation_name(std::move(corr_operation)), A_obs(std::move(obs1)),
         B_obs(std::move(obs2)) {
+    m_dt *= static_cast<double>(delta_N);
     initialize();
   }
 
diff --git a/src/core/bonded_interactions/thermalized_bond_utils.cpp b/src/core/bonded_interactions/thermalized_bond_utils.cpp
index edde9e90e0..1b87b98c06 100644
--- a/src/core/bonded_interactions/thermalized_bond_utils.cpp
+++ b/src/core/bonded_interactions/thermalized_bond_utils.cpp
@@ -20,7 +20,6 @@
  */
 
 #include "thermalized_bond_utils.hpp"
-#include "integrate.hpp"
 #include "thermalized_bond.hpp"
 
 #include "bonded_interactions/bonded_interaction_data.hpp"
diff --git a/src/core/cell_system/HybridDecomposition.hpp b/src/core/cell_system/HybridDecomposition.hpp
index f5f7220a92..c8993027b7 100644
--- a/src/core/cell_system/HybridDecomposition.hpp
+++ b/src/core/cell_system/HybridDecomposition.hpp
@@ -31,7 +31,6 @@
 #include "LocalBox.hpp"
 #include "Particle.hpp"
 #include "ghosts.hpp"
-#include "integrate.hpp"
 
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
diff --git a/src/core/cells.cpp b/src/core/cells.cpp
index 4c8676074b..a0f95342f8 100644
--- a/src/core/cells.cpp
+++ b/src/core/cells.cpp
@@ -35,7 +35,6 @@
 #include "Particle.hpp"
 #include "communication.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
 #include "particle_node.hpp"
 #include "system/System.hpp"
 
diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp
index da57ce7fe5..9ee236f798 100644
--- a/src/core/dpd.cpp
+++ b/src/core/dpd.cpp
@@ -30,7 +30,6 @@
 #include "BoxGeometry.hpp"
 #include "cell_system/CellStructure.hpp"
 #include "cells.hpp"
-#include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "random.hpp"
 #include "system/System.hpp"
diff --git a/src/core/ek/Solver.cpp b/src/core/ek/Solver.cpp
index a11aed8eb9..2de532a762 100644
--- a/src/core/ek/Solver.cpp
+++ b/src/core/ek/Solver.cpp
@@ -26,7 +26,6 @@
 #include "ek/EKNone.hpp"
 #include "ek/EKWalberla.hpp"
 
-#include "integrate.hpp"
 #include "system/System.hpp"
 
 #include <utils/Vector.hpp>
diff --git a/src/core/electrostatics/coulomb.cpp b/src/core/electrostatics/coulomb.cpp
index 870bd72ea9..5285093753 100644
--- a/src/core/electrostatics/coulomb.cpp
+++ b/src/core/electrostatics/coulomb.cpp
@@ -32,8 +32,6 @@
 #include "communication.hpp"
 #include "electrostatics/icc.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
-#include "npt.hpp"
 #include "system/System.hpp"
 
 #include <utils/Vector.hpp>
@@ -114,8 +112,7 @@ struct LongRangePressure {
 
 #ifdef P3M
   auto operator()(std::shared_ptr<CoulombP3M> const &actor) const {
-    actor->charge_assign(m_particles);
-    return actor->p3m_calc_kspace_pressure_tensor();
+    return actor->long_range_pressure(m_particles);
   }
 #endif // P3M
 
@@ -215,24 +212,10 @@ struct LongRangeForce {
 
 #ifdef P3M
   void operator()(std::shared_ptr<CoulombP3M> const &actor) const {
-    actor->charge_assign(m_particles);
-#ifdef NPT
-    if (integ_switch == INTEG_METHOD_NPT_ISO) {
-      auto const energy = actor->long_range_kernel(true, true, m_particles);
-      npt_add_virial_contribution(energy);
-    } else
-#endif // NPT
-      actor->add_long_range_forces(m_particles);
+    actor->add_long_range_forces(m_particles);
   }
 #ifdef CUDA
   void operator()(std::shared_ptr<CoulombP3MGPU> const &actor) const {
-#ifdef NPT
-    if (integ_switch == INTEG_METHOD_NPT_ISO) {
-      actor->charge_assign(m_particles);
-      auto const energy = actor->long_range_energy(m_particles);
-      npt_add_virial_contribution(energy);
-    }
-#endif // NPT
     actor->add_long_range_forces(m_particles);
   }
 #endif // CUDA
@@ -266,7 +249,6 @@ struct LongRangeEnergy {
 
 #ifdef P3M
   auto operator()(std::shared_ptr<CoulombP3M> const &actor) const {
-    actor->charge_assign(m_particles);
     return actor->long_range_energy(m_particles);
   }
   auto
diff --git a/src/core/electrostatics/icc.cpp b/src/core/electrostatics/icc.cpp
index dd7348163a..975e76d4f4 100644
--- a/src/core/electrostatics/icc.cpp
+++ b/src/core/electrostatics/icc.cpp
@@ -34,13 +34,14 @@
 
 #include "Particle.hpp"
 #include "ParticleRange.hpp"
+#include "PropagationMode.hpp"
 #include "actor/visitors.hpp"
 #include "cell_system/CellStructure.hpp"
 #include "communication.hpp"
 #include "electrostatics/coulomb.hpp"
 #include "electrostatics/coulomb_inline.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
+#include "integrators/Propagation.hpp"
 #include "system/System.hpp"
 
 #include <utils/constants.hpp>
@@ -269,7 +270,7 @@ struct SanityChecksICC {
 void ICCStar::sanity_check() const {
   sanity_checks_active_solver();
 #ifdef NPT
-  if (integ_switch == INTEG_METHOD_NPT_ISO) {
+  if (get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO) {
     throw std::runtime_error("ICC does not work in the NPT ensemble");
   }
 #endif
diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp
index e62d5ff772..1a78b36d10 100644
--- a/src/core/electrostatics/p3m.cpp
+++ b/src/core/electrostatics/p3m.cpp
@@ -43,12 +43,14 @@
 #include "Particle.hpp"
 #include "ParticlePropertyIterator.hpp"
 #include "ParticleRange.hpp"
+#include "PropagationMode.hpp"
 #include "actor/visitors.hpp"
 #include "cell_system/CellStructure.hpp"
 #include "cell_system/CellStructureType.hpp"
 #include "communication.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
+#include "integrators/Propagation.hpp"
+#include "npt.hpp"
 #include "system/System.hpp"
 #include "tuning.hpp"
 
@@ -404,7 +406,8 @@ static auto calc_dipole_moment(boost::mpi::communicator const &comm,
  *  in @cite essmann95a. The part \f$\Pi_{\textrm{corr}, \alpha, \beta}\f$
  *  eq. (2.8) is not present here since M is the empty set in our simulations.
  */
-Utils::Vector9d CoulombP3M::p3m_calc_kspace_pressure_tensor() {
+Utils::Vector9d
+CoulombP3M::long_range_pressure(ParticleRange const &particles) {
   using namespace detail::FFT_indexing;
 
   auto const &box_geo = *get_system().box_geo;
@@ -412,6 +415,7 @@ Utils::Vector9d CoulombP3M::p3m_calc_kspace_pressure_tensor() {
   Utils::Vector9d node_k_space_pressure_tensor{};
 
   if (p3m.sum_q2 > 0.) {
+    charge_assign(particles);
     p3m.sm.gather_grid(p3m.rs_mesh.data(), comm_cart, p3m.local_mesh.dim);
     fft_perform_forw(p3m.rs_mesh.data(), p3m.fft, comm_cart);
 
@@ -464,12 +468,25 @@ Utils::Vector9d CoulombP3M::p3m_calc_kspace_pressure_tensor() {
 
 double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
                                      ParticleRange const &particles) {
-  auto const &box_geo = *get_system().box_geo;
+  auto const &system = get_system();
+  auto const &box_geo = *system.box_geo;
+#ifdef NPT
+  auto const npt_flag =
+      force_flag and (system.propagation->integ_switch == INTEG_METHOD_NPT_ISO);
+#else
+  auto constexpr npt_flag = false;
+#endif
 
-  /* Gather information for FFT grid inside the nodes domain (inner local mesh)
-   * and perform forward 3D FFT (Charge Assignment Mesh). */
-  p3m.sm.gather_grid(p3m.rs_mesh.data(), comm_cart, p3m.local_mesh.dim);
-  fft_perform_forw(p3m.rs_mesh.data(), p3m.fft, comm_cart);
+  if (p3m.sum_q2 > 0.) {
+    if (not has_actor_of_type<ElectrostaticLayerCorrection>(
+            system.coulomb.impl->solver)) {
+      charge_assign(particles);
+    }
+    /* Gather information for FFT grid inside the nodes domain (inner local
+     * mesh) and perform forward 3D FFT (Charge Assignment Mesh). */
+    p3m.sm.gather_grid(p3m.rs_mesh.data(), comm_cart, p3m.local_mesh.dim);
+    fft_perform_forw(p3m.rs_mesh.data(), p3m.fft, comm_cart);
+  }
 
   auto p_q_range = ParticlePropertyRange::charge_range(particles);
   auto p_force_range = ParticlePropertyRange::force_range(particles);
@@ -549,7 +566,7 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
   }
 
   /* === k-space energy calculation  === */
-  if (energy_flag) {
+  if (energy_flag or npt_flag) {
     auto node_energy = 0.;
     for (int i = 0; i < p3m.fft.plan[3].new_size; i++) {
       // Use the energy optimized influence function for energy!
@@ -575,7 +592,14 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
         energy += pref * box_dipole.value().norm2();
       }
     }
-    return prefactor * energy;
+    energy *= prefactor;
+    if (npt_flag) {
+      npt_add_virial_contribution(energy);
+    }
+    if (not energy_flag) {
+      energy = 0.;
+    }
+    return energy;
   }
 
   return 0.;
diff --git a/src/core/electrostatics/p3m.hpp b/src/core/electrostatics/p3m.hpp
index 82506a9359..776eb57cb3 100644
--- a/src/core/electrostatics/p3m.hpp
+++ b/src/core/electrostatics/p3m.hpp
@@ -218,7 +218,7 @@ struct CoulombP3M : public Coulomb::Actor<CoulombP3M> {
   }
 
   /** Compute the k-space part of the pressure tensor */
-  Utils::Vector9d p3m_calc_kspace_pressure_tensor();
+  Utils::Vector9d long_range_pressure(ParticleRange const &particles);
 
   /** Compute the k-space part of energies. */
   double long_range_energy(ParticleRange const &particles) {
diff --git a/src/core/electrostatics/p3m_gpu.cpp b/src/core/electrostatics/p3m_gpu.cpp
index 425adff720..f8cf8abfe8 100644
--- a/src/core/electrostatics/p3m_gpu.cpp
+++ b/src/core/electrostatics/p3m_gpu.cpp
@@ -29,6 +29,9 @@
 #include "electrostatics/coulomb.hpp"
 
 #include "ParticleRange.hpp"
+#include "PropagationMode.hpp"
+#include "integrators/Propagation.hpp"
+#include "npt.hpp"
 #include "system/GpuParticleData.hpp"
 #include "system/System.hpp"
 
@@ -46,7 +49,15 @@ static auto get_n_part_safe(GpuParticleData const &gpu) {
   return static_cast<unsigned int>(n_part);
 }
 
-void CoulombP3MGPU::add_long_range_forces(ParticleRange const &) {
+void CoulombP3MGPU::add_long_range_forces(ParticleRange const &particles) {
+#ifdef NPT
+  if (get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO) {
+    auto const energy = long_range_energy(particles);
+    npt_add_virial_contribution(energy);
+  }
+#else
+  static_cast<void>(particles);
+#endif
   if (this_node == 0) {
     auto &gpu = get_system().gpu;
     p3m_gpu_add_farfield_force(*m_gpu_data, gpu, prefactor,
diff --git a/src/core/energy.cpp b/src/core/energy.cpp
index a62df7bdb9..ae66fd23e5 100644
--- a/src/core/energy.cpp
+++ b/src/core/energy.cpp
@@ -25,7 +25,6 @@
 #include "constraints.hpp"
 #include "energy_inline.hpp"
 #include "global_ghost_flags.hpp"
-#include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "short_range_loop.hpp"
 #include "system/System.hpp"
diff --git a/src/core/forces.cpp b/src/core/forces.cpp
index c7a10f0a4e..1bc21bffcb 100644
--- a/src/core/forces.cpp
+++ b/src/core/forces.cpp
@@ -38,7 +38,7 @@
 #include "forces_inline.hpp"
 #include "galilei/ComFixed.hpp"
 #include "immersed_boundaries.hpp"
-#include "integrate.hpp"
+#include "integrators/Propagation.hpp"
 #include "lb/particle_coupling.hpp"
 #include "magnetostatics/dipoles.hpp"
 #include "nonbonded_interactions/VerletCriterion.hpp"
@@ -228,7 +228,10 @@ void System::System::calculate_forces(double kT) {
 #endif // CUDA
 
 #ifdef VIRTUAL_SITES_RELATIVE
-  vs_relative_back_transfer_forces_and_torques(*cell_structure);
+  if (propagation->used_propagations &
+      (PropagationMode::TRANS_VS_RELATIVE | PropagationMode::ROT_VS_RELATIVE)) {
+    vs_relative_back_transfer_forces_and_torques(*cell_structure);
+  }
 #endif
 
   // Communication step: ghost forces
@@ -241,7 +244,7 @@ void System::System::calculate_forces(double kT) {
   force_capping(particles, force_cap);
 
   // mark that forces are now up-to-date
-  recalc_forces = false;
+  propagation->recalc_forces = false;
 }
 
 void calc_long_range_forces(const ParticleRange &particles) {
diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp
index bd9bc3309e..dcb5e81a11 100644
--- a/src/core/integrate.cpp
+++ b/src/core/integrate.cpp
@@ -27,6 +27,7 @@
  */
 
 #include "integrate.hpp"
+#include "integrators/Propagation.hpp"
 #include "integrators/brownian_inline.hpp"
 #include "integrators/steepest_descent.hpp"
 #include "integrators/stokesian_dynamics_inline.hpp"
@@ -86,18 +87,6 @@
 #endif
 #endif
 
-int integ_switch = INTEG_METHOD_NVT;
-static int default_propagation = PropagationMode::NONE;
-static int used_propagations = PropagationMode::NONE;
-
-/** Actual simulation time. */
-static double sim_time = 0.0;
-
-bool recalc_forces = true;
-
-static int lb_skipped_md_steps = 0;
-static int ek_skipped_md_steps = 0;
-
 namespace {
 volatile std::sig_atomic_t ctrl_C = 0;
 } // namespace
@@ -112,7 +101,7 @@ std::weak_ptr<ActiveProtocol> get_protocol() { return protocol; }
  * @brief Update the Lees-Edwards parameters of the box geometry
  * for the current simulation time.
  */
-static void update_box_params(BoxGeometry &box_geo) {
+static void update_box_params(BoxGeometry &box_geo, double sim_time) {
   if (box_geo.type() == BoxType::LEES_EDWARDS) {
     assert(protocol != nullptr);
     box_geo.lees_edwards_update(get_pos_offset(sim_time, *protocol),
@@ -126,8 +115,8 @@ void set_protocol(std::shared_ptr<ActiveProtocol> new_protocol) {
   auto &box_geo = *system.box_geo;
   box_geo.set_type(BoxType::LEES_EDWARDS);
   protocol = std::move(new_protocol);
-  LeesEdwards::update_box_params(box_geo);
-  ::recalc_forces = true;
+  LeesEdwards::update_box_params(box_geo, system.get_sim_time());
+  system.propagation->recalc_forces = true;
   cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
 }
 
@@ -137,7 +126,7 @@ void unset_protocol() {
   auto &box_geo = *system.box_geo;
   protocol = nullptr;
   box_geo.set_type(BoxType::CUBOID);
-  ::recalc_forces = true;
+  system.propagation->recalc_forces = true;
   cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
 }
 
@@ -153,7 +142,7 @@ template <class Kernel> void run_kernel(BoxGeometry const &box_geo) {
 }
 } // namespace LeesEdwards
 
-static void update_default_propagation() {
+void Propagation::update_default_propagation() {
   switch (integ_switch) {
   case INTEG_METHOD_STEEPEST_DESCENT:
     default_propagation = PropagationMode::NONE;
@@ -205,35 +194,36 @@ static void update_default_propagation() {
 }
 
 void System::System::update_used_propagations() {
-  ::used_propagations = PropagationMode::NONE;
+  int used_propagations = PropagationMode::NONE;
   for (auto &p : cell_structure->local_particles()) {
-    ::used_propagations |= p.propagation();
+    used_propagations |= p.propagation();
   }
-  if (::used_propagations & PropagationMode::SYSTEM_DEFAULT) {
-    ::used_propagations |= default_propagation;
+  if (used_propagations & PropagationMode::SYSTEM_DEFAULT) {
+    used_propagations |= propagation->default_propagation;
   }
-  ::used_propagations = boost::mpi::all_reduce(::comm_cart, ::used_propagations,
-                                               std::bit_or<int>());
+  used_propagations = boost::mpi::all_reduce(::comm_cart, used_propagations,
+                                             std::bit_or<int>());
+  propagation->used_propagations = used_propagations;
 }
 
 void System::System::integrator_sanity_checks() const {
   if (time_step <= 0.) {
     runtimeErrorMsg() << "time_step not set";
   }
-  if (integ_switch == INTEG_METHOD_STEEPEST_DESCENT) {
+  if (propagation->integ_switch == INTEG_METHOD_STEEPEST_DESCENT) {
     if (thermo_switch != THERMO_OFF) {
       runtimeErrorMsg()
           << "The steepest descent integrator is incompatible with thermostats";
     }
   }
-  if (integ_switch == INTEG_METHOD_NVT) {
+  if (propagation->integ_switch == INTEG_METHOD_NVT) {
     if (thermo_switch & (THERMO_NPT_ISO | THERMO_BROWNIAN | THERMO_SD)) {
       runtimeErrorMsg() << "The VV integrator is incompatible with the "
                            "currently active combination of thermostats";
     }
   }
 #ifdef NPT
-  if (::used_propagations & PropagationMode::TRANS_LANGEVIN_NPT) {
+  if (propagation->used_propagations & PropagationMode::TRANS_LANGEVIN_NPT) {
     if (thermo_switch != THERMO_OFF and thermo_switch != THERMO_NPT_ISO) {
       runtimeErrorMsg() << "The NpT integrator requires the NpT thermostat";
     }
@@ -242,12 +232,12 @@ void System::System::integrator_sanity_checks() const {
     }
   }
 #endif
-  if (::used_propagations & PropagationMode::TRANS_BROWNIAN) {
+  if (propagation->used_propagations & PropagationMode::TRANS_BROWNIAN) {
     if (thermo_switch != THERMO_BROWNIAN) {
       runtimeErrorMsg() << "The BD integrator requires the BD thermostat";
     }
   }
-  if (::used_propagations & PropagationMode::TRANS_STOKESIAN) {
+  if (propagation->used_propagations & PropagationMode::TRANS_STOKESIAN) {
 #ifdef STOKESIAN_DYNAMICS
     if (thermo_switch != THERMO_OFF && thermo_switch != THERMO_SD) {
       runtimeErrorMsg() << "The SD integrator requires the SD thermostat";
@@ -320,18 +310,13 @@ static void resort_particles_if_needed(System::System &system) {
   }
 }
 
-static bool should_propagate_with(Particle const &p, int mode) {
-  return p.propagation() & mode or
-         ((default_propagation & mode) and
-          (p.propagation() & PropagationMode::SYSTEM_DEFAULT));
-}
-
 void System::System::thermostats_force_init(double kT) {
+  auto const &propagation = *this->propagation;
   for (auto &p : cell_structure->local_particles()) {
-    if (should_propagate_with(p, PropagationMode::TRANS_LANGEVIN))
+    if (propagation.should_propagate_with(p, PropagationMode::TRANS_LANGEVIN))
       p.force() += friction_thermo_langevin(langevin, p, time_step, kT);
 #ifdef ROTATION
-    if (should_propagate_with(p, PropagationMode::ROT_LANGEVIN))
+    if (propagation.should_propagate_with(p, PropagationMode::ROT_LANGEVIN))
       p.torque() += convert_vector_body_to_space(
           p, friction_thermo_langevin_rotation(langevin, p, time_step, kT));
 #endif
@@ -341,11 +326,11 @@ void System::System::thermostats_force_init(double kT) {
 /** @brief Calls the hook for propagation kernels before the force calculation
  *  @return whether or not to stop the integration loop early.
  */
-static bool integrator_step_1(ParticleRange const &particles, double kT,
+static bool integrator_step_1(ParticleRange const &particles,
+                              Propagation const &propagation, double kT,
                               double time_step) {
-  using namespace PropagationMode;
   // steepest decent
-  if (integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
+  if (propagation.integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
     return steepest_descent_step(particles);
 
   for (auto &p : particles) {
@@ -354,103 +339,100 @@ static bool integrator_step_1(ParticleRange const &particles, double kT,
     if (p.is_virtual())
       continue;
 #endif
-    if (should_propagate_with(p, TRANS_LB_MOMENTUM_EXCHANGE))
+    if (propagation.should_propagate_with(
+            p, PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE))
       velocity_verlet_propagator_1(p, time_step);
-    if (should_propagate_with(p, TRANS_NEWTON))
+    if (propagation.should_propagate_with(p, PropagationMode::TRANS_NEWTON))
       velocity_verlet_propagator_1(p, time_step);
 #ifdef ROTATION
-    if (should_propagate_with(p, ROT_EULER))
+    if (propagation.should_propagate_with(p, PropagationMode::ROT_EULER))
       velocity_verlet_rotator_1(p, time_step);
 #endif
-    if (should_propagate_with(p, TRANS_LANGEVIN))
+    if (propagation.should_propagate_with(p, PropagationMode::TRANS_LANGEVIN))
       velocity_verlet_propagator_1(p, time_step);
 #ifdef ROTATION
-    if (should_propagate_with(p, ROT_LANGEVIN))
+    if (propagation.should_propagate_with(p, PropagationMode::ROT_LANGEVIN))
       velocity_verlet_rotator_1(p, time_step);
 #endif
-    if (should_propagate_with(p, TRANS_BROWNIAN))
+    if (propagation.should_propagate_with(p, PropagationMode::TRANS_BROWNIAN))
       brownian_dynamics_propagator(brownian, p, time_step, kT);
 #ifdef ROTATION
-    if (should_propagate_with(p, ROT_BROWNIAN))
+    if (propagation.should_propagate_with(p, PropagationMode::ROT_BROWNIAN))
       brownian_dynamics_rotator(brownian, p, time_step, kT);
 #endif
   }
 
 #ifdef NPT
-  if (::used_propagations & TRANS_LANGEVIN_NPT) {
-    if (default_propagation & TRANS_LANGEVIN_NPT)
-      velocity_verlet_npt_step_1(
-          particles.filter(PropagationPredicateNPT(default_propagation)),
-          time_step);
+  if ((propagation.used_propagations & PropagationMode::TRANS_LANGEVIN_NPT) and
+      (propagation.default_propagation & PropagationMode::TRANS_LANGEVIN_NPT)) {
+    auto pred = PropagationPredicateNPT(propagation.default_propagation);
+    velocity_verlet_npt_step_1(particles.filter(pred), time_step);
   }
 #endif
 
 #ifdef STOKESIAN_DYNAMICS
-  if (::used_propagations & TRANS_STOKESIAN) {
-    if (default_propagation & TRANS_STOKESIAN) {
-      stokesian_dynamics_step_1(
-          particles.filter(PropagationPredicateStokesian(default_propagation)),
-          time_step);
-    }
+  if ((propagation.used_propagations & PropagationMode::TRANS_STOKESIAN) and
+      (propagation.default_propagation & PropagationMode::TRANS_STOKESIAN)) {
+    auto pred = PropagationPredicateStokesian(propagation.default_propagation);
+    stokesian_dynamics_step_1(particles.filter(pred), time_step);
   }
 #endif // STOKESIAN_DYNAMICS
 
-  sim_time += time_step;
   return false;
 }
 
-static void integrator_step_2(ParticleRange const &particles, double kT,
+static void integrator_step_2(ParticleRange const &particles,
+                              Propagation const &propagation, double kT,
                               double time_step) {
-  using namespace PropagationMode;
-  if (integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
+  if (propagation.integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
     return;
+
   for (auto &p : particles) {
 #ifdef VIRTUAL_SITES
     // virtual sites are updated later in the integration loop
     if (p.is_virtual())
       continue;
 #endif
-    if (should_propagate_with(p, TRANS_LB_MOMENTUM_EXCHANGE))
+    if (propagation.should_propagate_with(
+            p, PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE))
       velocity_verlet_propagator_2(p, time_step);
-    if (should_propagate_with(p, TRANS_NEWTON))
+    if (propagation.should_propagate_with(p, PropagationMode::TRANS_NEWTON))
       velocity_verlet_propagator_2(p, time_step);
 #ifdef ROTATION
-    if (should_propagate_with(p, ROT_EULER))
+    if (propagation.should_propagate_with(p, PropagationMode::ROT_EULER))
       velocity_verlet_rotator_2(p, time_step);
 #endif
-    if (should_propagate_with(p, TRANS_LANGEVIN))
+    if (propagation.should_propagate_with(p, PropagationMode::TRANS_LANGEVIN))
       velocity_verlet_propagator_2(p, time_step);
 #ifdef ROTATION
-    if (should_propagate_with(p, ROT_LANGEVIN))
+    if (propagation.should_propagate_with(p, PropagationMode::ROT_LANGEVIN))
       velocity_verlet_rotator_2(p, time_step);
 #endif
   }
 
 #ifdef NPT
-  if (::used_propagations & TRANS_LANGEVIN_NPT) {
-    if (default_propagation & TRANS_LANGEVIN_NPT)
-      velocity_verlet_npt_step_2(
-          particles.filter(PropagationPredicateNPT(default_propagation)),
-          time_step);
+  if ((propagation.used_propagations & PropagationMode::TRANS_LANGEVIN_NPT) and
+      (propagation.default_propagation & PropagationMode::TRANS_LANGEVIN_NPT)) {
+    auto pred = PropagationPredicateNPT(propagation.default_propagation);
+    velocity_verlet_npt_step_2(particles.filter(pred), time_step);
   }
 #endif
 }
 
-namespace System {
-
-int System::integrate(int n_steps, int reuse_forces) {
+int System::System::integrate(int n_steps, int reuse_forces) {
 #ifdef CALIPER
   CALI_CXX_MARK_FUNCTION;
 #endif
+  auto &propagation = *this->propagation;
 #ifdef VIRTUAL_SITES_RELATIVE
-  auto const has_vs_rel = []() {
-    return ::used_propagations & (PropagationMode::ROT_VS_RELATIVE |
-                                  PropagationMode::TRANS_VS_RELATIVE);
+  auto const has_vs_rel = [&propagation]() {
+    return propagation.used_propagations & (PropagationMode::ROT_VS_RELATIVE |
+                                            PropagationMode::TRANS_VS_RELATIVE);
   };
 #endif
 
   // Prepare particle structure and run sanity checks of all active algorithms
-  update_default_propagation();
+  propagation.update_default_propagation();
   update_used_propagations();
   on_integration_start();
 
@@ -460,7 +442,8 @@ int System::integrate(int n_steps, int reuse_forces) {
 
   // Additional preparations for the first integration step
   if (reuse_forces == INTEG_REUSE_FORCES_NEVER or
-      (recalc_forces and reuse_forces != INTEG_REUSE_FORCES_ALWAYS)) {
+      ((reuse_forces != INTEG_REUSE_FORCES_ALWAYS) and
+       propagation.recalc_forces)) {
 #ifdef CALIPER
     CALI_MARK_BEGIN("Initial Force Calculation");
 #endif
@@ -477,7 +460,7 @@ int System::integrate(int n_steps, int reuse_forces) {
 
     calculate_forces(::temperature);
 
-    if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
+    if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
 #ifdef ROTATION
       convert_initial_torques(cell_structure->local_particles());
 #endif
@@ -504,7 +487,7 @@ int System::integrate(int n_steps, int reuse_forces) {
 
   auto lb_active = false;
   auto ek_active = false;
-  if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
+  if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
     lb_active = lb.is_solver_set();
     ek_active = ek.is_ready_for_propagation();
   }
@@ -532,15 +515,17 @@ int System::integrate(int n_steps, int reuse_forces) {
       save_old_position(particles, cell_structure->ghost_particles());
 #endif
 
-    LeesEdwards::update_box_params(*box_geo);
-    bool early_exit = integrator_step_1(particles, temperature, time_step);
+    LeesEdwards::update_box_params(*box_geo, sim_time);
+    bool early_exit =
+        integrator_step_1(particles, propagation, ::temperature, time_step);
     if (early_exit)
       break;
 
+    sim_time += time_step;
     LeesEdwards::run_kernel<LeesEdwards::Push>(*box_geo);
 
 #ifdef NPT
-    if (integ_switch != INTEG_METHOD_NPT_ISO)
+    if (propagation.integ_switch != INTEG_METHOD_NPT_ISO)
 #endif
     {
       resort_particles_if_needed(*this);
@@ -559,7 +544,7 @@ int System::integrate(int n_steps, int reuse_forces) {
 #ifdef VIRTUAL_SITES_RELATIVE
     if (has_vs_rel()) {
 #ifdef NPT
-      if (integ_switch == INTEG_METHOD_NPT_ISO) {
+      if (propagation.integ_switch == INTEG_METHOD_NPT_ISO) {
         cell_structure->update_ghosts_and_resort_particle(
             Cells::DATA_PART_PROPERTIES);
       }
@@ -579,12 +564,12 @@ int System::integrate(int n_steps, int reuse_forces) {
     calculate_forces(::temperature);
 
 #ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
-    if (::used_propagations & PropagationMode::TRANS_LB_TRACER) {
+    if (propagation.used_propagations & PropagationMode::TRANS_LB_TRACER) {
       lb_tracers_add_particle_force_to_fluid(*cell_structure, time_step);
     }
 #endif
-    integrator_step_2(particles, temperature, time_step);
-    if (integ_switch == INTEG_METHOD_BD) {
+    integrator_step_2(particles, propagation, ::temperature, time_step);
+    if (propagation.integ_switch == INTEG_METHOD_BD) {
       resort_particles_if_needed(*this);
     }
     LeesEdwards::run_kernel<LeesEdwards::UpdateOffset>(*box_geo);
@@ -596,7 +581,7 @@ int System::integrate(int n_steps, int reuse_forces) {
 #endif
 
     // propagate one-step functionalities
-    if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
+    if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
       if (lb_active and ek_active) {
         // assume that they are coupled, which is not necessarily true
         auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
@@ -607,36 +592,37 @@ int System::integrate(int n_steps, int reuse_forces) {
               << "LB and EK are active but with different time steps.";
         }
 
-        assert(lb_skipped_md_steps == ek_skipped_md_steps);
+        assert(propagation.lb_skipped_md_steps ==
+               propagation.ek_skipped_md_steps);
 
-        lb_skipped_md_steps += 1;
-        ek_skipped_md_steps += 1;
-        if (lb_skipped_md_steps >= md_steps_per_lb_step) {
-          lb_skipped_md_steps = 0;
-          ek_skipped_md_steps = 0;
+        propagation.lb_skipped_md_steps += 1;
+        propagation.ek_skipped_md_steps += 1;
+        if (propagation.lb_skipped_md_steps >= md_steps_per_lb_step) {
+          propagation.lb_skipped_md_steps = 0;
+          propagation.ek_skipped_md_steps = 0;
           lb.propagate();
           ek.propagate();
         }
         lb_lbcoupling_propagate();
       } else if (lb_active) {
         auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
-        lb_skipped_md_steps += 1;
-        if (lb_skipped_md_steps >= md_steps_per_lb_step) {
-          lb_skipped_md_steps = 0;
+        propagation.lb_skipped_md_steps += 1;
+        if (propagation.lb_skipped_md_steps >= md_steps_per_lb_step) {
+          propagation.lb_skipped_md_steps = 0;
           lb.propagate();
         }
         lb_lbcoupling_propagate();
       } else if (ek_active) {
         auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
-        ek_skipped_md_steps += 1;
-        if (ek_skipped_md_steps >= md_steps_per_ek_step) {
-          ek_skipped_md_steps = 0;
+        propagation.ek_skipped_md_steps += 1;
+        if (propagation.ek_skipped_md_steps >= md_steps_per_ek_step) {
+          propagation.ek_skipped_md_steps = 0;
           ek.propagate();
         }
       }
 
 #ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
-      if (::used_propagations & PropagationMode::TRANS_LB_TRACER) {
+      if (propagation.used_propagations & PropagationMode::TRANS_LB_TRACER) {
         lb_tracers_propagate(*cell_structure, time_step);
       }
 #endif
@@ -661,7 +647,7 @@ int System::integrate(int n_steps, int reuse_forces) {
     }
 
   } // for-loop over integration steps
-  LeesEdwards::update_box_params(*box_geo);
+  LeesEdwards::update_box_params(*box_geo, sim_time);
 #ifdef CALIPER
   CALI_CXX_MARK_LOOP_END(integration_loop);
 #endif
@@ -680,7 +666,7 @@ int System::integrate(int n_steps, int reuse_forces) {
   cell_structure->update_verlet_stats(n_steps, n_verlet_updates);
 
 #ifdef NPT
-  if (integ_switch == INTEG_METHOD_NPT_ISO) {
+  if (propagation.integ_switch == INTEG_METHOD_NPT_ISO) {
     synchronize_npt_state();
   }
 #endif
@@ -694,8 +680,8 @@ int System::integrate(int n_steps, int reuse_forces) {
   return integrated_steps;
 }
 
-int System::integrate_with_signal_handler(int n_steps, int reuse_forces,
-                                          bool update_accumulators) {
+int System::System::integrate_with_signal_handler(int n_steps, int reuse_forces,
+                                                  bool update_accumulators) {
   assert(n_steps >= 0);
 
   // Override the signal handler so that the integrator obeys Ctrl+C
@@ -745,21 +731,8 @@ int System::integrate_with_signal_handler(int n_steps, int reuse_forces,
   return 0;
 }
 
-} // namespace System
-
-double get_time_step() { return System::get_system().get_time_step(); }
-
-double get_sim_time() { return sim_time; }
-
-void set_time(double value) {
-  ::sim_time = value;
-  ::recalc_forces = true;
-  auto &system = System::get_system();
-  auto &box_geo = *system.box_geo;
-  LeesEdwards::update_box_params(box_geo);
-}
-
-void set_integ_switch(int value) {
-  ::integ_switch = value;
-  ::recalc_forces = true;
+void System::System::set_sim_time(double value) {
+  sim_time = value;
+  propagation->recalc_forces = true;
+  LeesEdwards::update_box_params(*box_geo, sim_time);
 }
diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp
index 53252b0e88..5ef1282524 100644
--- a/src/core/integrate.hpp
+++ b/src/core/integrate.hpp
@@ -29,21 +29,14 @@
 
 #include "config/config.hpp"
 
+#include "PropagationMode.hpp"
+
 #ifdef WALBERLA
 #include <string>
 
 #include <utils/Vector.hpp>
 #endif
 
-/** \name Integrator switches */
-/**@{*/
-#define INTEG_METHOD_NPT_ISO 0
-#define INTEG_METHOD_NVT 1
-#define INTEG_METHOD_STEEPEST_DESCENT 2
-#define INTEG_METHOD_BD 3
-#define INTEG_METHOD_SD 7
-/**@}*/
-
 /** \name Integrator error codes */
 /**@{*/
 #define INTEG_ERROR_RUNTIME -1
@@ -54,18 +47,12 @@
 /**@{*/
 /// recalculate forces unconditionally (mostly used for timing)
 #define INTEG_REUSE_FORCES_NEVER -1
-/// recalculate forces if @ref recalc_forces is set
+/// recalculate forces only if @ref Propagation::recalc_forces is set
 #define INTEG_REUSE_FORCES_CONDITIONALLY 0
 /// do not recalculate forces (mostly when reading checkpoints with forces)
 #define INTEG_REUSE_FORCES_ALWAYS 1
 /**@}*/
 
-/** Switch determining which integrator to use. */
-extern int integ_switch;
-
-/** If true, the forces will be recalculated before the next integration. */
-extern bool recalc_forces;
-
 #ifdef WALBERLA
 void walberla_tau_sanity_checks(std::string method, double tau,
                                 double time_step);
@@ -76,14 +63,3 @@ void walberla_agrid_sanity_checks(std::string method,
                                   Utils::Vector3d const &lattice_right,
                                   double agrid);
 #endif // WALBERLA
-
-/** Get time step */
-double get_time_step();
-
-/** Get simulation time */
-double get_sim_time();
-
-/** @brief Set the simulation time. */
-void set_time(double value);
-
-void set_integ_switch(int value);
diff --git a/src/core/integrators/Propagation.hpp b/src/core/integrators/Propagation.hpp
new file mode 100644
index 0000000000..f810c5f1ba
--- /dev/null
+++ b/src/core/integrators/Propagation.hpp
@@ -0,0 +1,49 @@
+/*
+ * Copyright (C) 2010-2023 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "PropagationMode.hpp"
+
+class Propagation {
+public:
+  int integ_switch = INTEG_METHOD_NVT;
+  int used_propagations = PropagationMode::NONE;
+  int default_propagation = PropagationMode::NONE;
+  int lb_skipped_md_steps = 0;
+  int ek_skipped_md_steps = 0;
+  /** If true, forces will be recalculated before the next integration. */
+  bool recalc_forces = true;
+
+  void update_default_propagation();
+
+  template <typename Particle>
+  bool should_propagate_with(Particle const &p, int mode) const {
+    return (p.propagation() & mode) or
+           ((default_propagation & mode) and
+            (p.propagation() & PropagationMode::SYSTEM_DEFAULT));
+  }
+
+  void set_integ_switch(int value) {
+    integ_switch = value;
+    recalc_forces = true;
+  }
+};
diff --git a/src/core/integrators/velocity_verlet_inline.hpp b/src/core/integrators/velocity_verlet_inline.hpp
index 67458932b1..909dd1e2d0 100644
--- a/src/core/integrators/velocity_verlet_inline.hpp
+++ b/src/core/integrators/velocity_verlet_inline.hpp
@@ -24,7 +24,6 @@
 #include "Particle.hpp"
 #include "ParticleRange.hpp"
 #include "cell_system/CellStructure.hpp"
-#include "integrate.hpp"
 #include "rotation.hpp"
 
 /** Propagate the velocities and positions. Integration steps before force
diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp
index 083815d595..328bf31859 100644
--- a/src/core/integrators/velocity_verlet_npt.cpp
+++ b/src/core/integrators/velocity_verlet_npt.cpp
@@ -28,7 +28,6 @@
 #include "cell_system/CellStructure.hpp"
 #include "communication.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
 #include "npt.hpp"
 #include "system/System.hpp"
 #include "thermostat.hpp"
diff --git a/src/core/magnetostatics/dipoles.cpp b/src/core/magnetostatics/dipoles.cpp
index 18f6eedcde..4f60e1fa8b 100644
--- a/src/core/magnetostatics/dipoles.cpp
+++ b/src/core/magnetostatics/dipoles.cpp
@@ -31,8 +31,6 @@
 #include "actor/visitors.hpp"
 #include "communication.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
-#include "npt.hpp"
 #include "system/System.hpp"
 
 #include <utils/Vector.hpp>
@@ -40,7 +38,6 @@
 #include <utils/demangle.hpp>
 
 #include <cassert>
-#include <cstdio>
 #include <optional>
 #include <stdexcept>
 
@@ -123,15 +120,7 @@ struct LongRangeForce {
 
 #ifdef DP3M
   void operator()(std::shared_ptr<DipolarP3M> const &actor) const {
-    actor->dipole_assign(m_particles);
-#ifdef NPT
-    if (integ_switch == INTEG_METHOD_NPT_ISO) {
-      auto const energy = actor->kernel(true, true, m_particles);
-      npt_add_virial_contribution(energy);
-      fprintf(stderr, "dipolar_P3M at this moment is added to p_vir[0]\n");
-    } else
-#endif // NPT
-      actor->kernel(true, false, m_particles);
+    actor->add_long_range_forces(m_particles);
   }
 #endif // DP3M
   void operator()(std::shared_ptr<DipolarLayerCorrection> const &actor) const {
@@ -165,8 +154,7 @@ struct LongRangeEnergy {
 
 #ifdef DP3M
   double operator()(std::shared_ptr<DipolarP3M> const &actor) const {
-    actor->dipole_assign(m_particles);
-    return actor->kernel(false, true, m_particles);
+    return actor->long_range_energy(m_particles);
   }
 #endif // DP3M
   double
diff --git a/src/core/magnetostatics/dp3m.cpp b/src/core/magnetostatics/dp3m.cpp
index 7fb65ec0be..d262e857b0 100644
--- a/src/core/magnetostatics/dp3m.cpp
+++ b/src/core/magnetostatics/dp3m.cpp
@@ -43,11 +43,12 @@
 #include "LocalBox.hpp"
 #include "Particle.hpp"
 #include "ParticleRange.hpp"
+#include "PropagationMode.hpp"
 #include "cell_system/CellStructure.hpp"
 #include "cell_system/CellStructureType.hpp"
 #include "communication.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
+#include "integrators/Propagation.hpp"
 #include "npt.hpp"
 #include "system/System.hpp"
 #include "tuning.hpp"
@@ -64,6 +65,7 @@
 
 #include <algorithm>
 #include <array>
+#include <cstdio>
 #include <functional>
 #include <optional>
 #include <sstream>
@@ -251,17 +253,22 @@ template <int cao> struct AssignForces {
 };
 } // namespace
 
-double DipolarP3M::kernel(bool force_flag, bool energy_flag,
-                          ParticleRange const &particles) {
-  int i, ind, j[3];
+double DipolarP3M::long_range_kernel(bool force_flag, bool energy_flag,
+                                     ParticleRange const &particles) {
   /* k-space energy */
-  double k_space_energy_dip = 0.;
-  double tmp0, tmp1;
-
-  auto const &box_geo = *get_system().box_geo;
+  double energy = 0.;
+  auto const &system = get_system();
+  auto const &box_geo = *system.box_geo;
   auto const dipole_prefac = prefactor / Utils::int_pow<3>(dp3m.params.mesh[0]);
-
-  if (dp3m.sum_mu2 > 0) {
+#ifdef NPT
+  auto const npt_flag =
+      force_flag and (system.propagation->integ_switch == INTEG_METHOD_NPT_ISO);
+#else
+  auto constexpr npt_flag = false;
+#endif
+
+  if (dp3m.sum_mu2 > 0.) {
+    dipole_assign(particles);
     /* Gather information for FFT grid inside the nodes domain (inner local
      * mesh) and perform forward 3D FFT (Charge Assignment Mesh). */
     std::array<double *, 3> meshes = {{dp3m.rs_mesh_dip[0].data(),
@@ -279,20 +286,21 @@ double DipolarP3M::kernel(bool force_flag, bool energy_flag,
   }
 
   /* === k-space energy calculation  === */
-  if (energy_flag) {
+  if (energy_flag or npt_flag) {
     /*********************
        Dipolar energy
     **********************/
-    if (dp3m.sum_mu2 > 0) {
+    if (dp3m.sum_mu2 > 0.) {
       /* i*k differentiation for dipolar gradients:
        * |(\Fourier{\vect{mu}}(k)\cdot \vect{k})|^2 */
-      ind = 0;
-      i = 0;
-      double node_k_space_energy_dip = 0.0;
+      int ind = 0;
+      int i = 0;
+      int j[3];
+      double node_energy = 0.0;
       for (j[0] = 0; j[0] < dp3m.fft.plan[3].new_mesh[0]; j[0]++) {
         for (j[1] = 0; j[1] < dp3m.fft.plan[3].new_mesh[1]; j[1]++) {
           for (j[2] = 0; j[2] < dp3m.fft.plan[3].new_mesh[2]; j[2]++) {
-            node_k_space_energy_dip +=
+            node_energy +=
                 dp3m.g_energy[i] *
                 (Utils::sqr(
                      dp3m.rs_mesh_dip[0][ind] *
@@ -313,23 +321,19 @@ double DipolarP3M::kernel(bool force_flag, bool energy_flag,
           }
         }
       }
-      node_k_space_energy_dip *=
-          dipole_prefac * Utils::pi() * box_geo.length_inv()[0];
-      boost::mpi::reduce(comm_cart, node_k_space_energy_dip, k_space_energy_dip,
-                         std::plus<>(), 0);
+      node_energy *= dipole_prefac * Utils::pi() * box_geo.length_inv()[0];
+      boost::mpi::reduce(comm_cart, node_energy, energy, std::plus<>(), 0);
 
       if (dp3m.energy_correction == 0.)
         calc_energy_correction();
 
       if (this_node == 0) {
         /* self energy correction */
-        k_space_energy_dip -= prefactor * dp3m.sum_mu2 * 2. *
-                              Utils::int_pow<3>(dp3m.params.alpha) *
-                              Utils::sqrt_pi_i() / 3.;
+        energy -= prefactor * dp3m.sum_mu2 * Utils::sqrt_pi_i() * (2. / 3.) *
+                  Utils::int_pow<3>(dp3m.params.alpha);
 
         /* dipolar energy correction due to systematic Madelung-self effects */
-        auto const volume = box_geo.volume();
-        k_space_energy_dip += prefactor * dp3m.energy_correction / volume;
+        energy += prefactor * dp3m.energy_correction / box_geo.volume();
       }
     }
   } // if (energy_flag)
@@ -339,11 +343,13 @@ double DipolarP3M::kernel(bool force_flag, bool energy_flag,
     /****************************
      * DIPOLAR TORQUES (k-space)
      ****************************/
-    if (dp3m.sum_mu2 > 0) {
+    if (dp3m.sum_mu2 > 0.) {
       auto const two_pi_L_i = 2. * Utils::pi() * box_geo.length_inv()[0];
       /* fill in ks_mesh array for torque calculation */
-      ind = 0;
-      i = 0;
+      int ind = 0;
+      int i = 0;
+      int j[3];
+      double tmp0, tmp1;
 
       for (j[0] = 0; j[0] < dp3m.fft.plan[3].new_mesh[0]; j[0]++) { // j[0]=n_y
         for (j[1] = 0; j[1] < dp3m.fft.plan[3].new_mesh[1];
@@ -497,12 +503,20 @@ double DipolarP3M::kernel(bool force_flag, bool energy_flag,
 
   if (dp3m.params.epsilon != P3M_EPSILON_METALLIC) {
     auto const surface_term =
-        calc_surface_term(force_flag, energy_flag, particles);
-    if (this_node == 0)
-      k_space_energy_dip += surface_term;
+        calc_surface_term(force_flag, energy_flag or npt_flag, particles);
+    if (this_node == 0) {
+      energy += surface_term;
+    }
+  }
+  if (npt_flag) {
+    npt_add_virial_contribution(energy);
+    fprintf(stderr, "dipolar_P3M at this moment is added to p_vir[0]\n");
+  }
+  if (not energy_flag) {
+    energy = 0.;
   }
 
-  return k_space_energy_dip;
+  return energy;
 }
 
 double DipolarP3M::calc_surface_term(bool force_flag, bool energy_flag,
diff --git a/src/core/magnetostatics/dp3m.hpp b/src/core/magnetostatics/dp3m.hpp
index fedfb2c22a..85bfcdeb63 100644
--- a/src/core/magnetostatics/dp3m.hpp
+++ b/src/core/magnetostatics/dp3m.hpp
@@ -169,9 +169,19 @@ struct DipolarP3M : public Dipoles::Actor<DipolarP3M> {
   void tune();
   bool is_tuned() const { return m_is_tuned; }
 
+  /** Compute the k-space part of energies. */
+  double long_range_energy(ParticleRange const &particles) {
+    return long_range_kernel(false, true, particles);
+  }
+
+  /** Compute the k-space part of forces. */
+  void add_long_range_forces(ParticleRange const &particles) {
+    long_range_kernel(true, false, particles);
+  }
+
   /** Compute the k-space part of forces and energies. */
-  double kernel(bool force_flag, bool energy_flag,
-                ParticleRange const &particles);
+  double long_range_kernel(bool force_flag, bool energy_flag,
+                           ParticleRange const &particles);
 
   /** Calculate real-space contribution of p3m dipolar pair forces and torques.
    *  If NPT is compiled in, update the NpT virial.
diff --git a/src/core/npt.cpp b/src/core/npt.cpp
index 0ea9016b6a..d75e4bbe06 100644
--- a/src/core/npt.cpp
+++ b/src/core/npt.cpp
@@ -20,11 +20,12 @@
 
 #ifdef NPT
 
+#include "PropagationMode.hpp"
 #include "communication.hpp"
 #include "config/config.hpp"
 #include "electrostatics/coulomb.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
+#include "integrators/Propagation.hpp"
 #include "magnetostatics/dipoles.hpp"
 #include "system/System.hpp"
 
@@ -108,21 +109,19 @@ Utils::Vector<bool, 3> NptIsoParameters::get_direction() const {
           static_cast<bool>(geometry & ::nptgeom_dir[2])};
 }
 
-void npt_ensemble_init(const BoxGeometry &box) {
-  if (integ_switch == INTEG_METHOD_NPT_ISO) {
-    /* prepare NpT-integration */
-    nptiso.inv_piston = 1. / nptiso.piston;
-    nptiso.volume = pow(box.length()[nptiso.non_const_dim], nptiso.dimension);
-    if (recalc_forces) {
-      nptiso.p_inst = 0.0;
-      nptiso.p_vir = Utils::Vector3d{};
-      nptiso.p_vel = Utils::Vector3d{};
-    }
+void npt_ensemble_init(Utils::Vector3d const &box_l, bool recalc_forces) {
+  nptiso.inv_piston = 1. / nptiso.piston;
+  nptiso.volume = std::pow(box_l[nptiso.non_const_dim], nptiso.dimension);
+  if (recalc_forces) {
+    nptiso.p_inst = 0.0;
+    nptiso.p_vir = Utils::Vector3d{};
+    nptiso.p_vel = Utils::Vector3d{};
   }
 }
 
 void integrator_npt_sanity_checks() {
-  if (integ_switch == INTEG_METHOD_NPT_ISO) {
+  if (::System::get_system().propagation->used_propagations &
+      PropagationMode::TRANS_LANGEVIN_NPT) {
     try {
       nptiso.coulomb_dipole_sanity_checks();
     } catch (std::runtime_error const &err) {
@@ -133,19 +132,23 @@ void integrator_npt_sanity_checks() {
 
 /** reset virial part of instantaneous pressure */
 void npt_reset_instantaneous_virials() {
-  if (integ_switch == INTEG_METHOD_NPT_ISO)
+  if (::System::get_system().propagation->used_propagations &
+      PropagationMode::TRANS_LANGEVIN_NPT) {
     nptiso.p_vir = Utils::Vector3d{};
+  }
 }
 
 void npt_add_virial_contribution(double energy) {
-  if (integ_switch == INTEG_METHOD_NPT_ISO) {
+  if (::System::get_system().propagation->used_propagations &
+      PropagationMode::TRANS_LANGEVIN_NPT) {
     nptiso.p_vir[0] += energy;
   }
 }
 
 void npt_add_virial_contribution(const Utils::Vector3d &force,
                                  const Utils::Vector3d &d) {
-  if (integ_switch == INTEG_METHOD_NPT_ISO) {
+  if (::System::get_system().propagation->used_propagations &
+      PropagationMode::TRANS_LANGEVIN_NPT) {
     nptiso.p_vir += hadamard_product(force, d);
   }
 }
diff --git a/src/core/npt.hpp b/src/core/npt.hpp
index 658526c033..defec034b7 100644
--- a/src/core/npt.hpp
+++ b/src/core/npt.hpp
@@ -22,13 +22,11 @@
  *  Exports for the NpT code.
  */
 
-#ifndef NPT_H
-#define NPT_H
+#pragma once
 
 #include "config/config.hpp"
 
 #ifdef NPT
-#include "BoxGeometry.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -84,7 +82,7 @@ extern NptIsoParameters nptiso;
 /** @brief Synchronizes NpT state such as instantaneous and average pressure
  */
 void synchronize_npt_state();
-void npt_ensemble_init(const BoxGeometry &box);
+void npt_ensemble_init(Utils::Vector3d const &box_l, bool recalc_forces);
 void integrator_npt_sanity_checks();
 void npt_reset_instantaneous_virials();
 void npt_add_virial_contribution(double energy);
@@ -92,4 +90,3 @@ void npt_add_virial_contribution(const Utils::Vector3d &force,
                                  const Utils::Vector3d &d);
 
 #endif // NPT
-#endif
diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp
index f65c4ef1be..9ed35a4753 100644
--- a/src/core/system/System.cpp
+++ b/src/core/system/System.cpp
@@ -24,6 +24,7 @@
 
 #include "BoxGeometry.hpp"
 #include "LocalBox.hpp"
+#include "PropagationMode.hpp"
 #include "bonded_interactions/bonded_interaction_data.hpp"
 #include "cell_system/CellStructure.hpp"
 #include "cell_system/CellStructureType.hpp"
@@ -35,7 +36,6 @@
 #include "errorhandling.hpp"
 #include "global_ghost_flags.hpp"
 #include "immersed_boundaries.hpp"
-#include "integrate.hpp"
 #include "npt.hpp"
 #include "particle_node.hpp"
 #include "thermostat.hpp"
@@ -60,12 +60,14 @@ System::System(Private) {
   box_geo = std::make_shared<BoxGeometry>();
   local_geo = std::make_shared<LocalBox>();
   cell_structure = std::make_shared<CellStructure>(*box_geo);
+  propagation = std::make_shared<Propagation>();
   nonbonded_ias = std::make_shared<InteractionsNonBonded>();
   comfixed = std::make_shared<ComFixed>();
   galilei = std::make_shared<Galilei>();
   bond_breakage = std::make_shared<BondBreakage::BondBreakage>();
   reinit_thermo = true;
   time_step = -1.;
+  sim_time = 0.;
   force_cap = 0.;
   min_global_cut = INACTIVE_CUTOFF;
 }
@@ -106,7 +108,7 @@ void System::set_time_step(double value) {
 
 void System::set_force_cap(double value) {
   force_cap = value;
-  ::recalc_forces = true;
+  propagation->recalc_forces = true;
 }
 
 void System::set_min_global_cut(double value) {
@@ -175,7 +177,7 @@ void System::on_periodicity_change() {
 #endif
 
 #ifdef STOKESIAN_DYNAMICS
-  if (::integ_switch == INTEG_METHOD_SD) {
+  if (propagation->integ_switch == INTEG_METHOD_SD) {
     if (box_geo->periodic(0u) or box_geo->periodic(1u) or box_geo->periodic(2u))
       runtimeErrorMsg() << "Stokesian Dynamics requires periodicity "
                         << "(False, False, False)\n";
@@ -222,13 +224,13 @@ void System::on_timestep_change() {
 
 void System::on_short_range_ia_change() {
   rebuild_cell_structure();
-  ::recalc_forces = true;
+  propagation->recalc_forces = true;
 }
 
 void System::on_non_bonded_ia_change() {
   nonbonded_ias->recalc_maximal_cutoffs();
   rebuild_cell_structure();
-  ::recalc_forces = true;
+  propagation->recalc_forces = true;
 }
 
 void System::on_coulomb_change() {
@@ -245,13 +247,15 @@ void System::on_dipoles_change() {
   on_short_range_ia_change();
 }
 
-void System::on_constraint_change() { ::recalc_forces = true; }
+void System::on_constraint_change() { propagation->recalc_forces = true; }
 
-void System::on_lb_boundary_conditions_change() { ::recalc_forces = true; }
+void System::on_lb_boundary_conditions_change() {
+  propagation->recalc_forces = true;
+}
 
 void System::on_particle_local_change() {
   cell_structure->update_ghosts_and_resort_particle(global_ghost_flags());
-  ::recalc_forces = true;
+  propagation->recalc_forces = true;
 }
 
 void System::on_particle_change() {
@@ -266,7 +270,7 @@ void System::on_particle_change() {
 #ifdef DIPOLES
   dipoles.on_particle_change();
 #endif
-  ::recalc_forces = true;
+  propagation->recalc_forces = true;
 
   /* the particle information is no longer valid */
   invalidate_fetch_cache();
@@ -380,11 +384,13 @@ void System::on_integration_start() {
   if (reinit_thermo) {
     thermo_init(time_step);
     reinit_thermo = false;
-    ::recalc_forces = true;
+    propagation->recalc_forces = true;
   }
 
 #ifdef NPT
-  npt_ensemble_init(*box_geo);
+  if (propagation->integ_switch == INTEG_METHOD_NPT_ISO) {
+    npt_ensemble_init(box_geo->length(), propagation->recalc_forces);
+  }
 #endif
 
   invalidate_fetch_cache();
diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp
index 32003dbbd7..41772a84f2 100644
--- a/src/core/system/System.hpp
+++ b/src/core/system/System.hpp
@@ -39,6 +39,7 @@
 class BoxGeometry;
 class LocalBox;
 struct CellStructure;
+class Propagation;
 class InteractionsNonBonded;
 class ComFixed;
 class Galilei;
@@ -76,6 +77,12 @@ class System : public std::enable_shared_from_this<System> {
   /** @brief Set @ref time_step. */
   void set_time_step(double value);
 
+  /** @brief Get @ref sim_time. */
+  auto get_sim_time() const { return sim_time; }
+
+  /** @brief Set @ref sim_time. */
+  void set_sim_time(double value);
+
   /** @brief Get @ref force_cap. */
   auto get_force_cap() const { return force_cap; }
 
@@ -151,7 +158,7 @@ class System : public std::enable_shared_from_this<System> {
    *  It is up to the propagation kernels to increment the simulation time.
    *
    *  This function propagates the system according to the choice of integrator
-   *  stored in @ref integ_switch. The general structure is:
+   *  stored in @ref Propagation::integ_switch. The general structure is:
    *  - if reuse_forces is zero, recalculate the forces based on the current
    *    state of the system
    *  - Loop over the number of simulation steps:
@@ -242,6 +249,7 @@ class System : public std::enable_shared_from_this<System> {
   std::shared_ptr<BoxGeometry> box_geo;
   std::shared_ptr<LocalBox> local_geo;
   std::shared_ptr<CellStructure> cell_structure;
+  std::shared_ptr<Propagation> propagation;
   std::shared_ptr<InteractionsNonBonded> nonbonded_ias;
   std::shared_ptr<ComFixed> comfixed;
   std::shared_ptr<Galilei> galilei;
@@ -253,6 +261,8 @@ class System : public std::enable_shared_from_this<System> {
   bool reinit_thermo;
   /** @brief Molecular dynamics integrator time step. */
   double time_step;
+  /** @brief Molecular dynamics integrator simulation time. */
+  double sim_time;
   /** @brief Molecular dynamics integrator force capping. */
   double force_cap;
   /**
diff --git a/src/core/system/System.impl.hpp b/src/core/system/System.impl.hpp
index fdb950d2b4..3a4962573e 100644
--- a/src/core/system/System.impl.hpp
+++ b/src/core/system/System.impl.hpp
@@ -29,6 +29,7 @@
 #include "cell_system/CellStructure.hpp"
 #include "galilei/ComFixed.hpp"
 #include "galilei/Galilei.hpp"
+#include "integrators/Propagation.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
 #include "BoxGeometry.hpp"
diff --git a/src/core/thermostat.cpp b/src/core/thermostat.cpp
index cf8eab195c..cc44101ac3 100644
--- a/src/core/thermostat.cpp
+++ b/src/core/thermostat.cpp
@@ -29,7 +29,6 @@
 #include "communication.hpp"
 #include "dpd.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
 #include "npt.hpp"
 #include "system/System.hpp"
 #include "thermostat.hpp"
diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
index db3f81fc75..ee604d1f2b 100644
--- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
+++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
@@ -29,6 +29,7 @@ namespace utf = boost::unit_test;
 
 #include "Observable_stat.hpp"
 #include "Particle.hpp"
+#include "PropagationMode.hpp"
 #include "accumulators/TimeSeries.hpp"
 #include "actor/registration.hpp"
 #include "bonded_interactions/bonded_interaction_utils.hpp"
@@ -42,6 +43,7 @@ namespace utf = boost::unit_test;
 #include "electrostatics/p3m.hpp"
 #include "galilei/Galilei.hpp"
 #include "integrate.hpp"
+#include "integrators/Propagation.hpp"
 #include "magnetostatics/dipoles.hpp"
 #include "nonbonded_interactions/lj.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
@@ -340,7 +342,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
   // check integration
   {
     // set up velocity-Verlet integrator
-    set_integ_switch(INTEG_METHOD_NVT);
+    espresso::system->propagation->set_integ_switch(INTEG_METHOD_NVT);
 
     // reset system
     remove_translational_motion(system);
diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp
index 143147720f..04d2ca1f4d 100644
--- a/src/core/unit_tests/Verlet_list_test.cpp
+++ b/src/core/unit_tests/Verlet_list_test.cpp
@@ -36,14 +36,17 @@ namespace bdata = boost::unit_test::data;
 #include "particle_management.hpp"
 
 #include "Particle.hpp"
+#include "PropagationMode.hpp"
 #include "cell_system/CellStructureType.hpp"
 #include "cells.hpp"
 #include "communication.hpp"
 #include "integrate.hpp"
+#include "integrators/Propagation.hpp"
 #include "integrators/steepest_descent.hpp"
 #include "nonbonded_interactions/lj.hpp"
 #include "npt.hpp"
 #include "particle_node.hpp"
+#include "system/System.hpp"
 #include "thermostat.hpp"
 
 #include <utils/Vector.hpp>
@@ -85,7 +88,8 @@ struct : public IntegratorHelper {
   void set_integrator() const override {
     mpi_set_thermo_switch_local(THERMO_OFF);
     register_integrator(SteepestDescentParameters(0., 0.01, 100.));
-    set_integ_switch(INTEG_METHOD_STEEPEST_DESCENT);
+    espresso::system->propagation->set_integ_switch(
+        INTEG_METHOD_STEEPEST_DESCENT);
   }
   void set_particle_properties(int pid) const override {
     set_particle_property(pid, &Particle::ext_force,
@@ -98,7 +102,7 @@ struct : public IntegratorHelper {
 struct : public IntegratorHelper {
   void set_integrator() const override {
     mpi_set_thermo_switch_local(THERMO_OFF);
-    set_integ_switch(INTEG_METHOD_NVT);
+    espresso::system->propagation->set_integ_switch(INTEG_METHOD_NVT);
   }
   void set_particle_properties(int pid) const override {
     set_particle_v(pid, {20., 0., 0.});
@@ -110,7 +114,7 @@ struct : public IntegratorHelper {
 struct : public IntegratorHelper {
   void set_integrator() const override {
     ::nptiso = NptIsoParameters(1., 1e9, {true, true, true}, true);
-    set_integ_switch(INTEG_METHOD_NPT_ISO);
+    espresso::system->propagation->set_integ_switch(INTEG_METHOD_NPT_ISO);
     mpi_set_temperature_local(1.);
     mpi_npt_iso_set_rng_seed(0);
     mpi_set_thermo_switch_local(thermo_switch | THERMO_NPT_ISO);
diff --git a/src/core/virtual_sites/relative.cpp b/src/core/virtual_sites/relative.cpp
index 379744b4bf..2a8887f7cd 100644
--- a/src/core/virtual_sites/relative.cpp
+++ b/src/core/virtual_sites/relative.cpp
@@ -27,7 +27,6 @@
 #include "cells.hpp"
 #include "errorhandling.hpp"
 #include "forces.hpp"
-#include "integrate.hpp"
 #include "lees_edwards/lees_edwards.hpp"
 #include "rotation.hpp"
 
diff --git a/src/script_interface/cell_system/CellSystem.cpp b/src/script_interface/cell_system/CellSystem.cpp
index 5a88f0d161..42388dc9bc 100644
--- a/src/script_interface/cell_system/CellSystem.cpp
+++ b/src/script_interface/cell_system/CellSystem.cpp
@@ -28,7 +28,6 @@
 #include "core/cell_system/RegularDecomposition.hpp"
 #include "core/cells.hpp"
 #include "core/communication.hpp"
-#include "core/integrate.hpp"
 #include "core/nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "core/particle_node.hpp"
 #include "core/system/System.hpp"
diff --git a/src/script_interface/h5md/h5md.cpp b/src/script_interface/h5md/h5md.cpp
index 2bace3a2e8..e3218ce870 100644
--- a/src/script_interface/h5md/h5md.cpp
+++ b/src/script_interface/h5md/h5md.cpp
@@ -26,7 +26,6 @@
 #include "h5md.hpp"
 
 #include "cell_system/CellStructure.hpp"
-#include "core/integrate.hpp"
 #include "core/system/System.hpp"
 
 #include <cmath>
@@ -38,12 +37,12 @@ namespace Writer {
 Variant H5md::do_call_method(const std::string &name,
                              const VariantMap &parameters) {
   if (name == "write") {
-    auto const &system = System::get_system();
-    auto &cell_structure = *system.cell_structure;
-    m_h5md->write(
-        cell_structure.local_particles(), get_sim_time(),
-        static_cast<int>(std::round(get_sim_time() / get_time_step())),
-        *system.box_geo);
+    auto const &system = ::System::get_system();
+    auto const particles = system.cell_structure->local_particles();
+    auto const sim_time = system.get_sim_time();
+    auto const time_step = system.get_time_step();
+    auto const n_steps = static_cast<int>(std::round(sim_time / time_step));
+    m_h5md->write(particles, sim_time, n_steps, *system.box_geo);
   } else if (name == "flush") {
     m_h5md->flush();
   } else if (name == "close") {
diff --git a/src/script_interface/integrators/BrownianDynamics.cpp b/src/script_interface/integrators/BrownianDynamics.cpp
index ecc7e7699c..e4af62a7b2 100644
--- a/src/script_interface/integrators/BrownianDynamics.cpp
+++ b/src/script_interface/integrators/BrownianDynamics.cpp
@@ -21,12 +21,15 @@
 
 #include "script_interface/ScriptInterface.hpp"
 
-#include "core/integrate.hpp"
+#include "core/PropagationMode.hpp"
+#include "core/integrators/Propagation.hpp"
 
 namespace ScriptInterface {
 namespace Integrators {
 
-void BrownianDynamics::activate() { set_integ_switch(INTEG_METHOD_BD); }
+void ScriptInterface::Integrators::BrownianDynamics::activate() {
+  get_system().propagation->set_integ_switch(INTEG_METHOD_BD);
+}
 
 } // namespace Integrators
 } // namespace ScriptInterface
diff --git a/src/script_interface/integrators/IntegratorHandle.cpp b/src/script_interface/integrators/IntegratorHandle.cpp
index 59dacfe483..3b08badf6e 100644
--- a/src/script_interface/integrators/IntegratorHandle.cpp
+++ b/src/script_interface/integrators/IntegratorHandle.cpp
@@ -27,7 +27,8 @@
 #include "VelocityVerlet.hpp"
 #include "VelocityVerletIsoNPT.hpp"
 
-#include "core/integrate.hpp"
+#include "core/PropagationMode.hpp"
+#include "core/integrators/Propagation.hpp"
 #include "core/system/System.hpp"
 
 #include <memory>
@@ -44,8 +45,11 @@ IntegratorHandle::IntegratorHandle() {
              [&]() { get_system().set_time_step(get_value<double>(v)); });
        },
        [this]() { return get_system().get_time_step(); }},
-      {"time", [](Variant const &v) { set_time(get_value<double>(v)); },
-       []() { return get_sim_time(); }},
+      {"time",
+       [&](Variant const &v) {
+         get_system().set_sim_time(get_value<double>(v));
+       },
+       [this]() { return get_system().get_sim_time(); }},
       {"force_cap",
        [this](Variant const &v) {
          get_system().set_force_cap(get_value<double>(v));
@@ -62,7 +66,7 @@ IntegratorHandle::IntegratorHandle() {
          m_instance->activate();
        },
        [this]() {
-         switch (::integ_switch) {
+         switch (get_system().propagation->integ_switch) {
          case INTEG_METHOD_STEEPEST_DESCENT:
            return Variant{
                std::dynamic_pointer_cast<SteepestDescent>(m_instance)};
@@ -93,7 +97,8 @@ void IntegratorHandle::on_bind_system(::System::System &system) {
   auto const &params = *m_params;
   for (auto const &key : get_parameter_insertion_order()) {
     if (params.count(key) != 0ul) {
-      if (not(key == "time_step" and ::integ_switch == INTEG_METHOD_NVT and
+      if (not(key == "time_step" and
+              system.propagation->integ_switch == INTEG_METHOD_NVT and
               system.get_time_step() == -1. and
               is_type<double>(params.at(key)) and
               get_value<double>(is_type<double>(params.at(key))) == -1.)) {
diff --git a/src/script_interface/integrators/SteepestDescent.cpp b/src/script_interface/integrators/SteepestDescent.cpp
index c3a4b1bb02..5b8f7ddce1 100644
--- a/src/script_interface/integrators/SteepestDescent.cpp
+++ b/src/script_interface/integrators/SteepestDescent.cpp
@@ -21,7 +21,9 @@
 
 #include "script_interface/ScriptInterface.hpp"
 
+#include "core/PropagationMode.hpp"
 #include "core/integrate.hpp"
+#include "core/integrators/Propagation.hpp"
 #include "core/integrators/steepest_descent.hpp"
 
 #include <utils/Vector.hpp>
@@ -69,7 +71,7 @@ void SteepestDescent::do_construct(VariantMap const &params) {
 
 void SteepestDescent::activate() {
   register_integrator(get_instance());
-  set_integ_switch(INTEG_METHOD_STEEPEST_DESCENT);
+  get_system().propagation->set_integ_switch(INTEG_METHOD_STEEPEST_DESCENT);
 }
 
 } // namespace Integrators
diff --git a/src/script_interface/integrators/StokesianDynamics.cpp b/src/script_interface/integrators/StokesianDynamics.cpp
index 54e7ab8c8f..40d0c7b784 100644
--- a/src/script_interface/integrators/StokesianDynamics.cpp
+++ b/src/script_interface/integrators/StokesianDynamics.cpp
@@ -25,7 +25,8 @@
 
 #include "script_interface/ScriptInterface.hpp"
 
-#include "core/integrate.hpp"
+#include "core/PropagationMode.hpp"
+#include "core/integrators/Propagation.hpp"
 #include "core/stokesian_dynamics/sd_interface.hpp"
 
 #include <memory>
@@ -92,7 +93,7 @@ void StokesianDynamics::do_construct(VariantMap const &params) {
 void StokesianDynamics::activate() {
   context()->parallel_try_catch([&]() {
     register_integrator(get_instance());
-    set_integ_switch(INTEG_METHOD_SD);
+    get_system().propagation->set_integ_switch(INTEG_METHOD_SD);
   });
 }
 
diff --git a/src/script_interface/integrators/VelocityVerlet.cpp b/src/script_interface/integrators/VelocityVerlet.cpp
index 4d7c7e32ea..325ea42052 100644
--- a/src/script_interface/integrators/VelocityVerlet.cpp
+++ b/src/script_interface/integrators/VelocityVerlet.cpp
@@ -21,12 +21,15 @@
 
 #include "script_interface/ScriptInterface.hpp"
 
-#include "core/integrate.hpp"
+#include "core/PropagationMode.hpp"
+#include "core/integrators/Propagation.hpp"
 
 namespace ScriptInterface {
 namespace Integrators {
 
-void VelocityVerlet::activate() { set_integ_switch(INTEG_METHOD_NVT); }
+void VelocityVerlet::activate() {
+  get_system().propagation->set_integ_switch(INTEG_METHOD_NVT);
+}
 
 } // namespace Integrators
 } // namespace ScriptInterface
diff --git a/src/script_interface/integrators/VelocityVerletIsoNPT.cpp b/src/script_interface/integrators/VelocityVerletIsoNPT.cpp
index 8464f75247..52c43b9fae 100644
--- a/src/script_interface/integrators/VelocityVerletIsoNPT.cpp
+++ b/src/script_interface/integrators/VelocityVerletIsoNPT.cpp
@@ -25,7 +25,8 @@
 
 #include "script_interface/ScriptInterface.hpp"
 
-#include "core/integrate.hpp"
+#include "core/PropagationMode.hpp"
+#include "core/integrators/Propagation.hpp"
 #include "core/npt.hpp"
 
 #include <utils/Vector.hpp>
@@ -64,7 +65,7 @@ void VelocityVerletIsoNPT::do_construct(VariantMap const &params) {
 
 void VelocityVerletIsoNPT::activate() {
   ::nptiso = get_instance();
-  set_integ_switch(INTEG_METHOD_NPT_ISO);
+  get_system().propagation->set_integ_switch(INTEG_METHOD_NPT_ISO);
   get_system().on_thermostat_param_change();
 }
 
diff --git a/src/script_interface/walberla/LBFluid.cpp b/src/script_interface/walberla/LBFluid.cpp
index 744ac5a8e0..e800962bbe 100644
--- a/src/script_interface/walberla/LBFluid.cpp
+++ b/src/script_interface/walberla/LBFluid.cpp
@@ -25,7 +25,6 @@
 #include "WalberlaCheckpoint.hpp"
 
 #include "core/BoxGeometry.hpp"
-#include "core/integrate.hpp"
 #include "core/lb/LBWalberla.hpp"
 #include "core/lees_edwards/lees_edwards.hpp"
 #include "core/lees_edwards/protocols.hpp"
@@ -180,11 +179,13 @@ void LBFluid::do_construct(VariantMap const &params) {
       auto lees_edwards_object = std::make_unique<LeesEdwardsPack>(
           le_bc.shear_direction, le_bc.shear_plane_normal,
           [this, le_protocol]() {
-            return get_pos_offset(get_sim_time(), *le_protocol) /
+            auto const &system = ::System::get_system();
+            return get_pos_offset(system.get_sim_time(), *le_protocol) /
                    m_lb_params->get_agrid();
           },
           [this, le_protocol]() {
-            return get_shear_velocity(get_sim_time(), *le_protocol) *
+            auto const &system = ::System::get_system();
+            return get_shear_velocity(system.get_sim_time(), *le_protocol) *
                    (m_lb_params->get_tau() / m_lb_params->get_agrid());
           });
       m_instance->set_collision_model(std::move(lees_edwards_object));