Remove more global variables

doc/sphinx/analysis.rst

@@ -787,7 +787,7 @@ To obtain the cluster structure of a system, an instance of
 To to create a cluster structure with above criterion::
 
     import espressomd.cluster_analysis
-    cs = espressomd.cluster_analysis.ClusterStructure(distance_criterion=dc)
+    cs = espressomd.cluster_analysis.ClusterStructure(distance_criterion=dc, system=system)
 
 In most cases, the cluster analysis is carried out by calling the
 :any:`espressomd.cluster_analysis.ClusterStructure.run_for_all_pairs` method.

doc/tutorials/ferrofluid/ferrofluid_part1.ipynb

@@ -802,7 +802,9 @@
    "source": [
     "# SOLUTION CELL\n",
     "# Setup cluster analysis\n",
-    "cluster_structure = espressomd.cluster_analysis.ClusterStructure(pair_criterion=espressomd.pair_criteria.DistanceCriterion(cut_off=1.3 * LJ_SIGMA))"
+    "cluster_structure = espressomd.cluster_analysis.ClusterStructure(\n",
+    "    system=system,\n",
+    "    pair_criterion=espressomd.pair_criteria.DistanceCriterion(cut_off=1.3 * LJ_SIGMA))"
    ]
   },
   {

src/core/bonded_interactions/bonded_interaction_data.cpp

@@ -16,18 +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/>.
  */
+
 #include "bonded_interaction_data.hpp"
 #include "rigid_bond.hpp"
 #include "system/System.hpp"
 #include "thermalized_bond.hpp"
+#include "thermostat.hpp"
 
 #include <boost/variant.hpp>
 
 #include <algorithm>
-#include <cstddef>
-#include <vector>
-
-BondedInteractionsMap bonded_ia_params;
+#include <numeric>
 
 /** Visitor to get the bond cutoff from the bond parameter variant */
 class BondCutoff : public boost::static_visitor<double> {
@@ -37,20 +36,18 @@ class BondCutoff : public boost::static_visitor<double> {
   }
 };
 
-double maximal_cutoff_bonded() {
+double BondedInteractionsMap::maximal_cutoff() const {
   auto const max_cut_bonded = std::accumulate(
-      bonded_ia_params.begin(), bonded_ia_params.end(), BONDED_INACTIVE_CUTOFF,
-      [](auto max_cut, auto const &kv) {
+      begin(), end(), BONDED_INACTIVE_CUTOFF, [](auto max_cut, auto const &kv) {
         return std::max(max_cut,
                         boost::apply_visitor(BondCutoff(), *kv.second));
       });
 
   /* Check if there are dihedrals */
-  auto const any_dihedrals = std::any_of(
-      bonded_ia_params.begin(), bonded_ia_params.end(), [](auto const &kv) {
-        return (boost::get<DihedralBond>(&(*kv.second)) ||
-                boost::get<TabulatedDihedralBond>(&(*kv.second)));
-      });
+  auto const any_dihedrals = std::any_of(begin(), end(), [](auto const &kv) {
+    return (boost::get<DihedralBond>(&(*kv.second)) ||
+            boost::get<TabulatedDihedralBond>(&(*kv.second)));
+  });
 
   /* dihedrals: the central particle is indirectly connected to the fourth
    * particle via the third particle, so we have to double the cutoff */
@@ -72,9 +69,21 @@ void BondedInteractionsMap::on_ia_change() {
     }
 #endif
   }
-  if (System::is_system_set()) {
-    auto &system = System::get_system();
-    system.on_short_range_ia_change();
-    system.on_thermostat_param_change(); // thermalized bonds
+  if (auto system = m_system.lock()) {
+    system->on_short_range_ia_change();
+    system->on_thermostat_param_change(); // thermalized bonds
+  }
+}
+
+void BondedInteractionsMap::activate_bond(mapped_type const &ptr) {
+  auto &system = get_system();
+  if (auto bond = boost::get<ThermalizedBond>(ptr.get())) {
+    bond->set_thermostat_view(system.thermostat);
+  }
+}
+
+void BondedInteractionsMap::deactivate_bond(mapped_type const &ptr) {
+  if (auto bond = boost::get<ThermalizedBond>(ptr.get())) {
+    bond->unset_thermostat_view();
   }
 }

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -43,7 +43,9 @@
 #include "rigid_bond.hpp"
 #include "thermalized_bond.hpp"
 
+#include "BondList.hpp"
 #include "TabulatedPotential.hpp"
+#include "system/Leaf.hpp"
 
 #include <boost/serialization/access.hpp>
 #include <boost/serialization/variant.hpp>
@@ -52,6 +54,7 @@
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <optional>
 #include <stdexcept>
 #include <unordered_map>
 #include <vector>
@@ -100,7 +103,10 @@ using Bonded_IA_Parameters =
                    RigidBond, IBMTriel, IBMVolCons, IBMTribend,
                    OifGlobalForcesBond, OifLocalForcesBond, VirtualBond>;
 
-class BondedInteractionsMap {
+/**
+ * @brief container for bonded interactions.
+ */
+class BondedInteractionsMap : public System::Leaf<BondedInteractionsMap> {
   using container_type =
       std::unordered_map<int, std::shared_ptr<Bonded_IA_Parameters>>;
 
@@ -111,32 +117,43 @@ class BondedInteractionsMap {
   using iterator = typename container_type::iterator;
   using const_iterator = typename container_type::const_iterator;
 
-  explicit BondedInteractionsMap() = default;
+  BondedInteractionsMap() = default;
+  virtual ~BondedInteractionsMap() = default;
 
   iterator begin() { return m_params.begin(); }
   iterator end() { return m_params.end(); }
   const_iterator begin() const { return m_params.begin(); }
   const_iterator end() const { return m_params.end(); }
 
   void insert(key_type const &key, mapped_type const &ptr) {
+    if (m_params.contains(key)) {
+      deactivate_bond(m_params.at(key));
+    }
+    activate_bond(ptr);
     next_key = std::max(next_key, key + 1);
     m_params[key] = ptr;
     on_ia_change();
   }
   key_type insert(mapped_type const &ptr) {
+    activate_bond(ptr);
     auto const key = next_key++;
     m_params[key] = ptr;
     on_ia_change();
     return key;
   }
   auto erase(key_type const &key) {
+    if (m_params.contains(key)) {
+      deactivate_bond(m_params.at(key));
+    }
     auto &&obj = m_params.erase(key);
     on_ia_change();
     return obj;
   }
+  virtual void activate_bond(mapped_type const &ptr);
+  virtual void deactivate_bond(mapped_type const &ptr);
   mapped_type at(key_type const &key) const { return m_params.at(key); }
   auto count(key_type const &key) const { return m_params.count(key); }
-  bool contains(key_type const &key) const { return m_params.count(key); }
+  bool contains(key_type const &key) const { return m_params.contains(key); }
   bool empty() const { return m_params.empty(); }
   auto size() const { return m_params.size(); }
   auto get_next_key() const { return next_key; }
@@ -147,6 +164,69 @@ class BondedInteractionsMap {
 #ifdef BOND_CONSTRAINT
   auto get_n_rigid_bonds() const { return n_rigid_bonds; }
 #endif
+  std::optional<key_type> find_bond_id(mapped_type const &bond) const {
+    for (auto const &kv : m_params) {
+      if (kv.second == bond) {
+        return kv.first;
+      }
+    }
+    return std::nullopt;
+  }
+
+  /**
+   * @brief Calculate the maximal cutoff of bonded interactions, required to
+   * determine the cell size for communication.
+   *
+   * Bond angle and dihedral potentials do not contain a cutoff intrinsically.
+   * The cutoff for these potentials depends on the bond length potentials
+   * (it is assumed that particles participating in a bond angle or dihedral
+   * potential are bound to each other by some bond length potential). For bond
+   * angle potentials nothing has to be done. For dihedral potentials the cutoff
+   * is set to twice the maximal cutoff because the particle in which the bond
+   * is stored is only bonded to the first two partners, one of which has an
+   * additional bond to the third partner.
+   */
+  double maximal_cutoff() const;
+
+  /**
+   * @brief Checks both particles for a specific bond, even on ghost particles.
+   *
+   * @param p           particle to check for the bond
+   * @param p_partner   possible bond partner
+   * @tparam BondType   Bond type to check for. Must be of one of the types in
+   *                    @ref Bonded_IA_Parameters.
+   */
+  template <typename BondType>
+  bool pair_bond_exists_on(Particle const &p, Particle const &p_partner) const {
+    auto const &bonds = p.bonds();
+    return std::any_of(
+        bonds.begin(), bonds.end(),
+        [this, partner_id = p_partner.id()](BondView const &bond) {
+          auto const &bond_ptr = at(bond.bond_id());
+          return (boost::get<BondType>(bond_ptr.get()) != nullptr) and
+                 (bond.partner_ids()[0] == partner_id);
+        });
+  }
+
+  /**
+   * @brief Checks both particles for a specific bond, even on ghost particles.
+   *
+   * @param p1     particle on which the bond may be stored
+   * @param p2     particle on which the bond may be stored
+   * @tparam BondType Bond type to check for.
+   */
+  template <typename BondType>
+  bool pair_bond_exists_between(Particle const &p1, Particle const &p2) const {
+    if (&p1 == &p2)
+      return false;
+
+    // Check if particles have bonds and search for the bond of interest.
+    // Could be saved on either particle, so we need to check both.
+    return pair_bond_exists_on<BondType>(p1, p2) or
+           pair_bond_exists_on<BondType>(p2, p1);
+  }
+
+  void on_ia_change();
 
 private:
   container_type m_params = {};
@@ -155,27 +235,9 @@ class BondedInteractionsMap {
 #ifdef BOND_CONSTRAINT
   int n_rigid_bonds = 0;
 #endif
-  void on_ia_change();
 };
 
-/** Field containing the parameters of the bonded ia types */
-extern BondedInteractionsMap bonded_ia_params;
-
-/** Calculate the maximal cutoff of bonded interactions, required to
- *  determine the cell size for communication.
- *
- *  Bond angle and dihedral potentials do not contain a cutoff intrinsically.
- *  The cutoff for these potentials depends on the bond length potentials
- *  (it is assumed that particles participating in a bond angle or dihedral
- *  potential are bound to each other by some bond length potential). For bond
- *  angle potentials nothing has to be done. For dihedral potentials the cutoff
- *  is set to twice the maximal cutoff because the particle in which the bond
- *  is stored is only bonded to the first two partners, one of which has an
- *  additional bond to the third partner.
- */
-double maximal_cutoff_bonded();
-
-/** Return the number of bonded partners for the specified bond */
+/** @brief Get the number of bonded partners for the specified bond. */
 inline int number_of_partners(Bonded_IA_Parameters const &iaparams) {
   return boost::apply_visitor(BondNumPartners(), iaparams);
 }

src/core/bonded_interactions/thermalized_bond.hpp

@@ -29,10 +29,16 @@
 
 #include <utils/Vector.hpp>
 
+#include <cassert>
 #include <cmath>
+#include <memory>
 #include <optional>
 #include <tuple>
 
+namespace Thermostat {
+class Thermostat;
+}
+
 /** Parameters for Thermalized bond */
 struct ThermalizedBond {
   double temp_com;
@@ -70,10 +76,20 @@ struct ThermalizedBond {
     pref2_dist = std::sqrt(24. * gamma_distance / time_step * temp_distance);
   }
 
+  void set_thermostat_view(
+      std::weak_ptr<Thermostat::Thermostat const> const &thermostat) {
+    m_thermostat = thermostat;
+  }
+
+  void unset_thermostat_view() { m_thermostat.reset(); }
+
   std::optional<std::tuple<Utils::Vector3d, Utils::Vector3d>>
   forces(Particle const &p1, Particle const &p2,
          Utils::Vector3d const &dx) const;
 
+private:
+  std::weak_ptr<Thermostat::Thermostat const> m_thermostat;
+
 private:
   friend boost::serialization::access;
   template <typename Archive>

src/core/bonded_interactions/thermalized_bond_kernel.hpp

@@ -25,7 +25,6 @@
 
 #include "Particle.hpp"
 #include "random.hpp"
-#include "system/System.hpp"
 #include "thermostat.hpp"
 
 #include <utils/Vector.hpp>
@@ -54,8 +53,9 @@ ThermalizedBond::forces(Particle const &p1, Particle const &p2,
   auto const sqrt_mass_red = sqrt(p1.mass() * p2.mass() / mass_tot);
   auto const com_vel = mass_tot_inv * (p1.mass() * p1.v() + p2.mass() * p2.v());
   auto const dist_vel = p2.v() - p1.v();
-  auto const &thermalized_bond =
-      *::System::get_system().thermostat->thermalized_bond;
+  auto const thermostat_view = m_thermostat.lock();
+  assert(thermostat_view);
+  auto const &thermalized_bond = *thermostat_view->thermalized_bond;
 
   Utils::Vector3d force1{};
   Utils::Vector3d force2{};