Encapsulate Lees-Edwards

src/core/integrate.cpp

@@ -92,54 +92,40 @@ volatile std::sig_atomic_t ctrl_C = 0;
 } // namespace
 
 namespace LeesEdwards {
-/** @brief Currently active Lees-Edwards protocol. */
-static std::shared_ptr<ActiveProtocol> protocol = nullptr;
-
-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, double sim_time) {
+void LeesEdwards::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),
-                                get_shear_velocity(sim_time, *protocol));
+    assert(m_protocol != nullptr);
+    box_geo.lees_edwards_update(get_pos_offset(sim_time, *m_protocol),
+                                get_shear_velocity(sim_time, *m_protocol));
   }
 }
 
-void set_protocol(std::shared_ptr<ActiveProtocol> new_protocol) {
-  auto &system = System::get_system();
+void LeesEdwards::set_protocol(std::shared_ptr<ActiveProtocol> protocol) {
+  auto &system = get_system();
   auto &cell_structure = *system.cell_structure;
   auto &box_geo = *system.box_geo;
   box_geo.set_type(BoxType::LEES_EDWARDS);
-  protocol = std::move(new_protocol);
-  LeesEdwards::update_box_params(box_geo, system.get_sim_time());
+  m_protocol = std::move(protocol);
+  update_box_params(box_geo, system.get_sim_time());
   system.propagation->recalc_forces = true;
   cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
 }
 
-void unset_protocol() {
-  auto &system = System::get_system();
+void LeesEdwards::unset_protocol() {
+  auto &system = get_system();
   auto &cell_structure = *system.cell_structure;
   auto &box_geo = *system.box_geo;
-  protocol = nullptr;
+  m_protocol = nullptr;
   box_geo.set_type(BoxType::CUBOID);
   system.propagation->recalc_forces = true;
   cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
 }
 
-template <class Kernel> void run_kernel(BoxGeometry const &box_geo) {
-  if (box_geo.type() == BoxType::LEES_EDWARDS) {
-    auto &system = System::get_system();
-    auto &cell_structure = *system.cell_structure;
-    auto const kernel = Kernel{box_geo};
-    auto const particles = cell_structure.local_particles();
-    std::for_each(particles.begin(), particles.end(),
-                  [&kernel](auto &p) { kernel(p); });
-  }
-}
 } // namespace LeesEdwards
 
 void Propagation::update_default_propagation() {
@@ -515,14 +501,19 @@ int System::System::integrate(int n_steps, int reuse_forces) {
       save_old_position(particles, cell_structure->ghost_particles());
 #endif
 
-    LeesEdwards::update_box_params(*box_geo, sim_time);
+    lees_edwards->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);
+    if (box_geo->type() == BoxType::LEES_EDWARDS) {
+      auto const kernel = LeesEdwards::Push{*box_geo};
+      for (auto &p : particles) {
+        kernel(p);
+      }
+    }
 
 #ifdef NPT
     if (propagation.integ_switch != INTEG_METHOD_NPT_ISO)
@@ -572,7 +563,12 @@ int System::System::integrate(int n_steps, int reuse_forces) {
     if (propagation.integ_switch == INTEG_METHOD_BD) {
       resort_particles_if_needed(*this);
     }
-    LeesEdwards::run_kernel<LeesEdwards::UpdateOffset>(*box_geo);
+    if (box_geo->type() == BoxType::LEES_EDWARDS) {
+      auto const kernel = LeesEdwards::UpdateOffset{*box_geo};
+      for (auto &p : particles) {
+        kernel(p);
+      }
+    }
 #ifdef BOND_CONSTRAINT
     // SHAKE velocity updates
     if (n_rigidbonds) {
@@ -647,7 +643,7 @@ int System::System::integrate(int n_steps, int reuse_forces) {
     }
 
   } // for-loop over integration steps
-  LeesEdwards::update_box_params(*box_geo, sim_time);
+  lees_edwards->update_box_params(*box_geo, sim_time);
 #ifdef CALIPER
   CALI_CXX_MARK_LOOP_END(integration_loop);
 #endif
@@ -734,5 +730,5 @@ int System::System::integrate_with_signal_handler(int n_steps, int reuse_forces,
 void System::System::set_sim_time(double value) {
   sim_time = value;
   propagation->recalc_forces = true;
-  LeesEdwards::update_box_params(*box_geo, sim_time);
+  lees_edwards->update_box_params(*box_geo, sim_time);
 }

src/core/lees_edwards/lees_edwards.hpp

@@ -16,17 +16,17 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef CORE_LEES_EDWARDS_LEES_EDWARDS_HPP
-#define CORE_LEES_EDWARDS_LEES_EDWARDS_HPP
+
+#pragma once
 
 #include "BoxGeometry.hpp"
 #include "Particle.hpp"
 #include "lees_edwards/protocols.hpp"
+#include "system/Leaf.hpp"
 
 #include <cmath>
 #include <memory>
 
-#include <iostream>
 namespace LeesEdwards {
 class UpdateOffset {
 protected:
@@ -90,14 +90,20 @@ inline Utils::Vector3d verlet_list_offset(BoxGeometry const &box,
   return {};
 }
 
-/** @brief Get currently active Lees-Edwards protocol. */
-std::weak_ptr<ActiveProtocol> get_protocol();
+class LeesEdwards : public System::Leaf<LeesEdwards> {
+  std::shared_ptr<ActiveProtocol> m_protocol;
+
+public:
+  /** @brief Get currently active Lees-Edwards protocol. */
+  auto get_protocol() const { return m_protocol; }
 
-/** @brief Set a new Lees-Edwards protocol. */
-void set_protocol(std::shared_ptr<ActiveProtocol> new_protocol);
+  /** @brief Set a new Lees-Edwards protocol. */
+  void set_protocol(std::shared_ptr<ActiveProtocol> protocol);
 
-/** @brief Delete the currently active Lees-Edwards protocol. */
-void unset_protocol();
+  /** @brief Delete the currently active Lees-Edwards protocol. */
+  void unset_protocol();
+
+  void update_box_params(BoxGeometry &box_geo, double sim_time);
+};
 
 } // namespace LeesEdwards
-#endif

src/core/lees_edwards/protocols.hpp

@@ -16,14 +16,13 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef CORE_LEES_EDWARDS_PROTOCOLS_HPP
-#define CORE_LEES_EDWARDS_PROTOCOLS_HPP
 
-#include <utils/Vector.hpp>
+#pragma once
 
-#include <boost/variant.hpp>
+#include <utils/Vector.hpp>
 
 #include <cmath>
+#include <variant>
 
 namespace LeesEdwards {
 
@@ -73,12 +72,12 @@ struct OscillatoryShear {
 };
 
 /** Type which holds the currently active protocol */
-using ActiveProtocol = boost::variant<Off, LinearShear, OscillatoryShear>;
+using ActiveProtocol = std::variant<Off, LinearShear, OscillatoryShear>;
 
-class PosOffsetGetter : public boost::static_visitor<double> {
+class PosOffsetGetter {
 public:
   PosOffsetGetter(double time) : m_time{time} {}
-  template <typename T> double operator()(const T &protocol) const {
+  template <typename T> double operator()(T const &protocol) const {
     return protocol.pos_offset(m_time);
   }
 
@@ -87,14 +86,14 @@ class PosOffsetGetter : public boost::static_visitor<double> {
 };
 
 inline double get_pos_offset(double time, ActiveProtocol const &protocol) {
-  return boost::apply_visitor(PosOffsetGetter(time), protocol);
+  return std::visit(PosOffsetGetter(time), protocol);
 }
 
 /** Visitor to get shear velocity from the Lees-Edwards protocol */
-class ShearVelocityGetter : public boost::static_visitor<double> {
+class ShearVelocityGetter {
 public:
   ShearVelocityGetter(double time) : m_time{time} {}
-  template <typename T> double operator()(const T &protocol) const {
+  template <typename T> double operator()(T const &protocol) const {
     return protocol.shear_velocity(m_time);
   }
 
@@ -104,9 +103,7 @@ class ShearVelocityGetter : public boost::static_visitor<double> {
 
 /** Calculation of current velocity */
 inline double get_shear_velocity(double time, ActiveProtocol const &protocol) {
-  return boost::apply_visitor(ShearVelocityGetter(time), protocol);
+  return std::visit(ShearVelocityGetter(time), protocol);
 }
 
 } // namespace LeesEdwards
-
-#endif

src/core/system/System.cpp

@@ -65,6 +65,7 @@ System::System(Private) {
   comfixed = std::make_shared<ComFixed>();
   galilei = std::make_shared<Galilei>();
   bond_breakage = std::make_shared<BondBreakage::BondBreakage>();
+  lees_edwards = std::make_shared<LeesEdwards::LeesEdwards>();
   reinit_thermo = true;
   time_step = -1.;
   sim_time = 0.;
@@ -75,6 +76,7 @@ System::System(Private) {
 void System::initialize() {
   auto handle = shared_from_this();
   cell_structure->bind_system(handle);
+  lees_edwards->bind_system(handle);
 #ifdef CUDA
   gpu.bind_system(handle);
   gpu.initialize();

src/core/system/System.hpp

@@ -47,6 +47,9 @@ class Observable_stat;
 namespace BondBreakage {
 class BondBreakage;
 }
+namespace LeesEdwards {
+class LeesEdwards;
+}
 
 namespace System {
 
@@ -254,6 +257,7 @@ class System : public std::enable_shared_from_this<System> {
   std::shared_ptr<ComFixed> comfixed;
   std::shared_ptr<Galilei> galilei;
   std::shared_ptr<BondBreakage::BondBreakage> bond_breakage;
+  std::shared_ptr<LeesEdwards::LeesEdwards> lees_edwards;
 
 protected:
   /** @brief Whether the thermostat has to be reinitialized before integration.

src/core/system/System.impl.hpp

@@ -30,6 +30,7 @@
 #include "galilei/ComFixed.hpp"
 #include "galilei/Galilei.hpp"
 #include "integrators/Propagation.hpp"
+#include "lees_edwards/lees_edwards.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
 #include "BoxGeometry.hpp"

src/core/virtual_sites.cpp

@@ -29,7 +29,6 @@
 #include "Particle.hpp"
 #include "communication.hpp"
 #include "errorhandling.hpp"
-#include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
 #include <utils/Vector.hpp>

src/python/espressomd/system.py

@@ -237,7 +237,7 @@ def _restore_object(cls, so_callback, so_callback_args, state):
 
     def __getstate__(self):
         checkpointable_properties = [
-            "non_bonded_inter", "bonded_inter", "lees_edwards",
+            "non_bonded_inter", "bonded_inter",
             "part", "auto_update_accumulators",
             "constraints",
         ]