Encapsulate propagator

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,
+};

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>

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();
   }
 

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"

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>

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"
 

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"

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>

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

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

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,14 +406,16 @@ 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;
 
   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.;

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) {

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,

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"

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) {