Encapsulate electrostatics and magnetostatics

src/core/EspressoSystemStandAlone.cpp

@@ -23,7 +23,6 @@
 #include "EspressoSystemStandAlone.hpp"
 #include "cell_system/CellStructure.hpp"
 #include "cell_system/CellStructureType.hpp"
-#include "cells.hpp"
 #include "communication.hpp"
 #include "system/System.hpp"
 #include "system/System.impl.hpp"
@@ -55,8 +54,7 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) {
 }
 
 void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const {
-  m_instance->box_geo->set_length(box_l);
-  m_instance->on_boxl_change();
+  m_instance->set_box_l(box_l);
 }
 
 void EspressoSystemStandAlone::set_node_grid(

src/core/EspressoSystemStandAlone.hpp

@@ -19,6 +19,8 @@
 
 #pragma once
 
+#include <system/System.hpp>
+
 #include <utils/Vector.hpp>
 
 #include <memory>
@@ -37,6 +39,7 @@ class EspressoSystemStandAlone {
   void set_node_grid(Utils::Vector3i const &node_grid) const;
   void set_time_step(double time_step) const;
   void set_skin(double new_skin) const;
+  auto get_handle() { return m_instance; }
 
 private:
   bool head_node;

src/core/actor/registration.hpp

@@ -25,21 +25,41 @@
 #include <functional>
 #include <memory>
 #include <optional>
+#include <variant>
+
+namespace System {
+class System;
+}
 
 template <typename Variant, typename T, class F>
 void add_actor(boost::mpi::communicator const &comm,
+               std::shared_ptr<System::System> const &system,
                std::optional<Variant> &active_actor,
                std::shared_ptr<T> const &actor, F &&on_actor_change) {
   std::optional<Variant> other = actor;
+  auto const activate = [&system](auto &leaf) {
+    leaf->bind_system(system);
+    leaf->on_activation();
+  };
+  auto const deactivate = [&system](auto &leaf) {
+    leaf->detach_system(system);
+  };
   auto const cleanup_if_any_rank_failed = [&](bool failed) {
     if (boost::mpi::all_reduce(comm, failed, std::logical_or<>())) {
+      deactivate(actor);
       active_actor.swap(other);
+      if (active_actor) {
+        std::visit([&](auto &leaf) { activate(leaf); }, *active_actor);
+      }
       on_actor_change();
     }
   };
   try {
     active_actor.swap(other);
-    actor->on_activation();
+    if (other) {
+      std::visit([&](auto &leaf) { deactivate(leaf); }, *other);
+    }
+    activate(actor);
     on_actor_change();
     cleanup_if_any_rank_failed(false);
   } catch (...) {

src/core/electrostatics/actor.hpp

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 The ESPResSo project
+ * Copyright (C) 2022-2023 The ESPResSo project
  *
  * This file is part of ESPResSo.
  *
@@ -16,20 +16,24 @@
  * 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 ESPRESSO_SRC_CORE_ELECTROSTATICS_ACTOR_HPP
-#define ESPRESSO_SRC_CORE_ELECTROSTATICS_ACTOR_HPP
+
+#pragma once
 
 #include "config/config.hpp"
 
 #ifdef ELECTROSTATICS
 
-#include "electrostatics/solver.hpp"
+#include "system/Leaf.hpp"
 
 #include <stdexcept>
 
 namespace Coulomb {
 
-template <typename Class> class Actor {
+/** @brief Check if the system is charge-neutral. */
+void check_charge_neutrality(System::System const &system,
+                             double relative_tolerance);
+
+template <typename Class> class Actor : public System::Leaf<Actor<Class>> {
 public:
   static auto constexpr charge_neutrality_tolerance_default = 2e-12;
   /**
@@ -51,12 +55,12 @@ template <typename Class> class Actor {
 
   void sanity_checks_charge_neutrality() const {
     if (charge_neutrality_tolerance != -1.) {
-      check_charge_neutrality(charge_neutrality_tolerance);
+      check_charge_neutrality(System::Leaf<Actor<Class>>::get_system(),
+                              charge_neutrality_tolerance);
     }
   }
 };
 
 } // namespace Coulomb
 
 #endif // ELECTROSTATICS
-#endif

src/core/electrostatics/coulomb.cpp

@@ -325,11 +325,11 @@ static auto calc_charge_excess_ratio(std::vector<double> const &charges) {
   return std::abs(sum_kahan(q_sum)) / q_min;
 }
 
-void check_charge_neutrality(double relative_tolerance) {
+void check_charge_neutrality(System::System const &system,
+                             double relative_tolerance) {
   // collect non-zero particle charges from all nodes
-  auto &cell_structure = *System::get_system().cell_structure;
   std::vector<double> local_charges;
-  for (auto const &p : cell_structure.local_particles()) {
+  for (auto const &p : system.cell_structure->local_particles()) {
     local_charges.push_back(p.q());
   }
   std::vector<std::vector<double>> node_charges;

src/core/electrostatics/coulomb.hpp

@@ -16,8 +16,8 @@
  * 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 ESPRESSO_SRC_CORE_ELECTROSTATICS_COULOMB_HPP
-#define ESPRESSO_SRC_CORE_ELECTROSTATICS_COULOMB_HPP
+
+#pragma once
 
 #include "config/config.hpp"
 
@@ -109,4 +109,3 @@ template <> struct has_pressure<CoulombMMM1D> : std::false_type {};
 } // namespace traits
 } // namespace Coulomb
 #endif // ELECTROSTATICS
-#endif

src/core/electrostatics/coulomb_inline.hpp

@@ -17,16 +17,14 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#ifndef ESPRESSO_SRC_CORE_ELECTROSTATICS_COULOMB_INLINE_HPP
-#define ESPRESSO_SRC_CORE_ELECTROSTATICS_COULOMB_INLINE_HPP
+#pragma once
 
 #include "config/config.hpp"
 
 #include "electrostatics/coulomb.hpp"
 #include "electrostatics/solver.hpp"
 
 #include "Particle.hpp"
-#include "system/System.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/demangle.hpp>
@@ -86,11 +84,9 @@ struct ShortRangeForceCorrectionsKernel {
   result_type
   operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &ptr) const {
     auto const &actor = *ptr;
-    auto const &box_geo = *System::get_system().box_geo;
-    return kernel_type{
-        [&actor, &box_geo](Particle &p1, Particle &p2, double q1q2) {
-          actor.add_pair_force_corrections(p1, p2, q1q2, box_geo);
-        }};
+    return kernel_type{[&actor](Particle &p1, Particle &p2, double q1q2) {
+      actor.add_pair_force_corrections(p1, p2, q1q2);
+    }};
   }
 #endif // P3M
 };
@@ -137,16 +133,15 @@ struct ShortRangeEnergyKernel {
   result_type
   operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &ptr) const {
     auto const &actor = *ptr;
-    auto const &box_geo = *System::get_system().box_geo;
     auto const energy_kernel = std::visit(*this, actor.base_solver);
-    return kernel_type{[&actor, &box_geo, energy_kernel](
+    return kernel_type{[&actor, energy_kernel](
                            Particle const &p1, Particle const &p2, double q1q2,
                            Utils::Vector3d const &d, double dist) {
       auto energy = 0.;
       if (energy_kernel) {
         energy = (*energy_kernel)(p1, p2, q1q2, d, dist);
       }
-      return energy + actor.pair_energy_correction(p1, p2, q1q2, box_geo);
+      return energy + actor.pair_energy_correction(p1, p2, q1q2);
     }};
   }
 #endif // P3M
@@ -209,5 +204,3 @@ Solver::pair_energy_kernel() const {
 }
 
 } // namespace Coulomb
-
-#endif

src/core/electrostatics/debye_hueckel.hpp

@@ -24,8 +24,7 @@
  * Calculate the Debye-Hückel energy and force for a particle pair.
  */
 
-#ifndef ESPRESSO_SRC_CORE_ELECTROSTATICS_DEBYE_HUECKEL_HPP
-#define ESPRESSO_SRC_CORE_ELECTROSTATICS_DEBYE_HUECKEL_HPP
+#pragma once
 
 #include "config/config.hpp"
 
@@ -108,5 +107,3 @@ struct DebyeHueckel : public Coulomb::Actor<DebyeHueckel> {
 };
 
 #endif // ELECTROSTATICS
-
-#endif