Refactor bonded interactions storage (#4350)

Fixes #4225
Fixes #4377
Fixes #4169

Description of changes:
- Bonded interactions are now stored as ScriptInterface objects, are immutable and can be removed
- It is no longer possible to overwrite a bond object by a different type of bond (avoids segmentation faults)

doc/sphinx/io.rst

@@ -103,8 +103,8 @@ Be aware of the following limitations:
 
 * Checkpointing only supports recursion on the head node. It is therefore
   impossible to checkpoint a :class:`espressomd.system.System` instance that
-  contains LB boundaries, constraints or auto-update accumulators, when the
-  simulation is running with 2 or more MPI nodes.
+  contains LB boundaries, constraints, bonded interactions or auto-update
+  accumulators, when the simulation is running with 2 or more MPI nodes.
 
 * The active actors, i.e., the content of ``system.actors``, are checkpointed. For lattice-Boltzmann fluids, this only includes the parameters such as the lattice constant (``agrid``). The actual flow field has to be saved separately with the lattice-Boltzmann specific methods
   :meth:`espressomd.lb.HydrodynamicInteraction.save_checkpoint`

samples/immersed_boundary/addBending.py

@@ -33,6 +33,6 @@ def AddBending(system, kb):
             id3 = int(line[2])
             id4 = int(line[3])
             tribend = espressomd.interactions.IBM_Tribend(
-                ind1=id1, ind2=id2, ind3=id3, ind4=id4, kb=kb, refShape="initial")
+                ind1=id1, ind2=id2, ind3=id3, ind4=id4, kb=kb, refShape="Initial")
             system.bonded_inter.add(tribend)
             system.part[id1].add_bond((tribend, id2, id3, id4))

src/core/Observable_stat.cpp

@@ -52,7 +52,7 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
 #else
   constexpr std::size_t n_vs = 0;
 #endif
-  auto const n_bonded = bonded_ia_params.size();
+  auto const n_bonded = bonded_ia_params.get_next_key();
   auto const n_non_bonded = max_non_bonded_pairs();
   constexpr std::size_t n_ext_fields = 1; // reduction over all fields
   constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions

src/core/bonded_interactions/CMakeLists.txt

@@ -2,15 +2,9 @@ target_sources(
   EspressoCore
   PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/angle_cosine.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/angle_cossquare.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/angle_harmonic.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/bonded_coulomb.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/bonded_coulomb_sr.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/bonded_interaction_data.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/bonded_tab.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/dihedral.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/fene.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/harmonic.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/quartic.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/rigid_bond.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/thermalized_bond.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/thermalized_bond_utils.cpp)

src/core/bonded_interactions/angle_cosine.hpp

@@ -49,7 +49,6 @@ struct AngleCosineBond {
 
   static constexpr int num = 2;
 
-  AngleCosineBond() = default;
   AngleCosineBond(double bend, double phi0);
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>

src/core/bonded_interactions/angle_cossquare.hpp

@@ -46,7 +46,6 @@ struct AngleCossquareBond {
 
   static constexpr int num = 2;
 
-  AngleCossquareBond() = default;
   AngleCossquareBond(double bend, double phi0);
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>

src/core/bonded_interactions/angle_harmonic.hpp

@@ -45,8 +45,10 @@ struct AngleHarmonicBond {
 
   static constexpr int num = 2;
 
-  AngleHarmonicBond() = default;
-  AngleHarmonicBond(double bend, double phi0);
+  AngleHarmonicBond(double bend, double phi0) {
+    this->bend = bend;
+    this->phi0 = phi0;
+  }
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
   forces(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left,

src/core/bonded_interactions/bonded_coulomb.hpp

@@ -23,8 +23,6 @@
 /** \file
  *  Routines to calculate the bonded Coulomb potential between
  *  particle pairs.
- *
- *  Implementation in \ref bonded_coulomb.cpp
  */
 
 #include "config.hpp"
@@ -44,8 +42,7 @@ struct BondedCoulomb {
 
   static constexpr int num = 1;
 
-  BondedCoulomb() = default;
-  BondedCoulomb(double prefactor);
+  BondedCoulomb(double prefactor) { this->prefactor = prefactor; }
 
   boost::optional<Utils::Vector3d> force(double q1q2,
                                          Utils::Vector3d const &dx) const;

src/core/bonded_interactions/bonded_coulomb_sr.hpp

@@ -24,8 +24,6 @@
  *  Routines to calculate the short-range part of the bonded Coulomb potential
  *  between particle pairs. Can be used to subtract certain intramolecular
  *  interactions in combination with Thole damping.
- *
- *  Implementation in \ref bonded_coulomb_sr.cpp.
  */
 
 #include "config.hpp"
@@ -48,8 +46,7 @@ struct BondedCoulombSR {
 
   static constexpr int num = 1;
 
-  BondedCoulombSR() = default;
-  BondedCoulombSR(double q1q2);
+  BondedCoulombSR(double q1q2) { this->q1q2 = q1q2; }
 
   boost::optional<Utils::Vector3d> force(Utils::Vector3d const &dx) const;
   boost::optional<double> energy(Particle const &p1, Particle const &p2,

src/core/bonded_interactions/bonded_interaction_data.cpp

@@ -28,7 +28,7 @@
 #include <cstddef>
 #include <vector>
 
-std::vector<Bonded_IA_Parameters> bonded_ia_params;
+BondedInteractionsMap bonded_ia_params;
 
 /** Visitor to get the bond cutoff from the bond parameter variant */
 class BondCutoff : public boost::static_visitor<double> {
@@ -41,28 +41,19 @@ class BondCutoff : public boost::static_visitor<double> {
 double maximal_cutoff_bonded() {
   auto const max_cut_bonded = boost::accumulate(
       bonded_ia_params, BONDED_INACTIVE_CUTOFF,
-      [](auto max_cut, Bonded_IA_Parameters const &bond) {
-        return std::max(max_cut, boost::apply_visitor(BondCutoff(), bond));
+      [](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 &bond) {
-        return (boost::get<DihedralBond>(&bond) ||
-                boost::get<TabulatedDihedralBond>(&bond));
+      bonded_ia_params.begin(), bonded_ia_params.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 */
   return (any_dihedrals) ? 2 * max_cut_bonded : max_cut_bonded;
 }
-
-void make_bond_type_exist(int type) {
-  std::size_t ns = static_cast<std::size_t>(type) + 1;
-  auto const old_size = bonded_ia_params.size();
-  if (ns <= old_size) {
-    return;
-  }
-  /* else allocate new memory */
-  bonded_ia_params.resize(ns, NoneBond());
-}

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -49,9 +49,11 @@
 #include <boost/serialization/variant.hpp>
 #include <boost/variant.hpp>
 
+#include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <stdexcept>
+#include <unordered_map>
 #include <vector>
 
 /* Special cutoff value for a disabled bond.
@@ -100,15 +102,48 @@ using Bonded_IA_Parameters =
                    RigidBond, IBMTriel, IBMVolCons, IBMTribend,
                    OifGlobalForcesBond, OifLocalForcesBond, VirtualBond>;
 
-/** Field containing the parameters of the bonded ia types */
-extern std::vector<Bonded_IA_Parameters> bonded_ia_params;
+class BondedInteractionsMap {
+  using container_type =
+      std::unordered_map<int, std::shared_ptr<Bonded_IA_Parameters>>;
 
-/** Makes sure that \ref bonded_ia_params is large enough to cover the
- *  parameters for the bonded interaction type.
- *  Attention: 1: There is no initialization done here.
- *  2: Use only in connection with creating new or overwriting old bond types
- */
-void make_bond_type_exist(int type);
+public:
+  using key_type = typename container_type::key_type;
+  using mapped_type = typename container_type::mapped_type;
+  using value_type = typename container_type::value_type;
+  using iterator = typename container_type::iterator;
+  using const_iterator = typename container_type::const_iterator;
+
+  explicit 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) {
+    next_key = std::max(next_key, key + 1);
+    m_params[key] = ptr;
+  }
+  key_type insert(mapped_type const &ptr) {
+    auto const key = next_key++;
+    m_params[key] = ptr;
+    return key;
+  }
+  auto erase(key_type const &key) { return m_params.erase(key); }
+  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 empty() const { return m_params.empty(); }
+  auto size() const { return m_params.size(); }
+  auto get_next_key() const { return next_key; }
+
+private:
+  container_type m_params = {};
+  key_type next_key = static_cast<key_type>(0);
+};
+
+/** 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.

src/core/bonded_interactions/bonded_interaction_utils.hpp

@@ -40,7 +40,8 @@ inline bool pair_bond_enum_exists_on(Particle const &p,
                                      Particle const &p_partner) {
   return boost::algorithm::any_of(
       p.bonds(), [partner_id = p_partner.identity()](BondView const &bond) {
-        return (boost::get<BondType>(&bonded_ia_params[bond.bond_id()])) and
+        auto const &bond_ptr = bonded_ia_params.at(bond.bond_id());
+        return (boost::get<BondType>(bond_ptr.get()) != nullptr) and
                (bond.partner_ids()[0] == partner_id);
       });
 }
@@ -65,19 +66,4 @@ inline bool pair_bond_enum_exists_between(Particle const &p1,
          pair_bond_enum_exists_on<BondType>(p2, p1);
 }
 
-/** Sets bond parameters.
- *  @param bond_id  ID of the bond to be set. If the bond ID doesn't exist yet,
- *                  it will be created.
- *  @param iaparams Bond to be set.
- *  @tparam BondType One of the types in @ref Bonded_IA_Parameters.
- */
-template <typename BondType>
-void set_bonded_ia_params(int bond_id, BondType const &iaparams) {
-  make_bond_type_exist(bond_id);
-  bonded_ia_params[bond_id] = iaparams;
-
-  /* broadcast interaction parameters */
-  mpi_bcast_ia_params(bond_id, -1);
-}
-
 #endif