Breaking bonds

src/core/CMakeLists.txt

@@ -1,5 +1,6 @@
 set(EspressoCore_SRC
     accumulators.cpp
+    bond_breakage.cpp
     bond_error.cpp
     cells.cpp
     collision.cpp

src/core/bond_breakage.cpp

@@ -0,0 +1,232 @@
+/*
+ * Copyright (C) 2022 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * 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 "bond_breakage.hpp"
+#include "cells.hpp"
+#include "communication.hpp"
+#include "errorhandling.hpp"
+#include "particle_data.hpp"
+
+#include <utils/mpi/gather_buffer.hpp>
+
+#include <boost/functional/hash.hpp>
+#include <boost/variant.hpp>
+
+#include <cassert>
+#include <cstddef>
+#include <memory>
+#include <unordered_set>
+#include <vector>
+
+namespace BondBreakage {
+
+// Bond breakage specifications
+std::unordered_map<int, std::shared_ptr<BreakageSpec>> breakage_specs;
+
+// Delete Actions
+struct DeleteBond {
+  int particle_id;
+  int bond_partner_id;
+  int bond_type;
+  std::size_t hash_value() const {
+    std::size_t seed = 3875;
+    boost::hash_combine(seed, particle_id);
+    boost::hash_combine(seed, bond_partner_id);
+    boost::hash_combine(seed, bond_type);
+    return seed;
+  }
+  bool operator==(DeleteBond const &rhs) const {
+    return rhs.particle_id == particle_id and
+           rhs.bond_partner_id == bond_partner_id and
+           rhs.bond_type == bond_type;
+  }
+};
+
+struct DeleteAllBonds {
+  int particle_id_1;
+  int particle_id_2;
+  std::size_t hash_value() const {
+    std::size_t seed = 75;
+    boost::hash_combine(seed, particle_id_1);
+    boost::hash_combine(seed, particle_id_2);
+    return seed;
+  }
+  bool operator==(DeleteAllBonds const &rhs) const {
+    return rhs.particle_id_1 == particle_id_1 and
+           rhs.particle_id_2 == particle_id_2;
+  }
+};
+} // namespace BondBreakage
+
+// Hash support for std::unordered_set
+namespace boost {
+template <> struct hash<BondBreakage::DeleteBond> {
+  std::size_t operator()(BondBreakage::DeleteBond const &t) const noexcept {
+    return t.hash_value();
+  }
+};
+template <> struct hash<BondBreakage::DeleteAllBonds> {
+  std::size_t operator()(BondBreakage::DeleteAllBonds const &t) const noexcept {
+    return t.hash_value();
+  }
+};
+} // namespace boost
+
+namespace BondBreakage {
+// Variant holding any of the actions
+using Action = boost::variant<DeleteBond, DeleteAllBonds>;
+
+// Set of actions
+using ActionSet = std::unordered_set<Action>;
+
+// Broken bond record
+struct QueueEntry {
+  int particle_id;
+  int bond_partner_id;
+  int bond_type;
+
+  /// Serialization for synchronization across mpi ranks
+  friend class boost::serialization::access;
+  template <typename Archive>
+  void serialize(Archive &ar, const unsigned int version) {
+    ar &particle_id;
+    ar &bond_partner_id;
+    ar &bond_type;
+  }
+};
+
+/** @brief Queue to record bonds broken during a time step */
+using Queue = std::vector<QueueEntry>;
+Queue queue;
+
+/** @brief Retrieve breakage specification for the bond type */
+boost::optional<BreakageSpec> get_breakage_spec(int bond_type) {
+  if (breakage_specs.find(bond_type) != breakage_specs.end()) {
+    return {*(breakage_specs.at(bond_type))};
+  }
+  return {};
+}
+
+/** Add a particle+bond combination to the breakage queue */
+void queue_breakage(int particle_id, int bond_partner_id, int bond_type) {
+  queue.emplace_back(QueueEntry{particle_id, bond_partner_id, bond_type});
+}
+
+bool check_and_handle_breakage(int particle_id, int bond_partner_id,
+                               int bond_type, double distance) {
+  // Retrieve specification for this bond type
+  auto spec = get_breakage_spec(bond_type);
+  if (!spec)
+    return false; // No breakage rule for this bond type
+
+  // Is the bond length longer than the breakage length?
+  if (distance >= (*spec).breakage_length) {
+    queue_breakage(particle_id, bond_partner_id, bond_type);
+    return true;
+  }
+  return false;
+}
+
+void clear_queue() { queue.clear(); }
+
+/** @brief Gathers combined queue from all mpi ranks */
+Queue gather_global_queue(Queue const &local_queue) {
+  Queue res = local_queue;
+  Utils::Mpi::gather_buffer(res, comm_cart);
+  boost::mpi::broadcast(comm_cart, res, 0);
+  return res;
+}
+
+/** @brief Constructs the actions to take for a breakage queue entry */
+ActionSet actions_for_breakage(QueueEntry const &e) {
+  // Retrieve relevant breakage spec
+  auto const spec = get_breakage_spec(e.bond_type);
+  assert(spec);
+
+  // Handle different action types
+  if ((*spec).action_type == ActionType::DELETE_BOND)
+    return {DeleteBond{e.particle_id, e.bond_partner_id, e.bond_type}};
+#ifdef VIRTUAL_SITES_RELATIVE
+  if ((*spec).action_type == ActionType::REVERT_BIND_AT_POINT_OF_COLLISION) {
+    // We need to find the base particles for the two virtual sites
+    // between which the bond broke.
+    auto p1 = cell_structure.get_local_particle(e.particle_id);
+    auto p2 = cell_structure.get_local_particle(e.bond_partner_id);
+    if (!p1 || !p2)
+      return {}; // particles not on this mpi rank
+
+    if (!p1->p.is_virtual || !p2->p.is_virtual) {
+      runtimeErrorMsg() << "The REVERT_BIND_AT_POINT_OF_COLLISION bond "
+                           "breakage action has to be configured for the "
+                           "bond on the virtual site. Encountered a particle "
+                           "that is not virtual.";
+      return {};
+    }
+
+    return {
+        // Bond between virtual sites
+        DeleteBond{e.particle_id, e.bond_partner_id, e.bond_type},
+        // Bond between base particles. We do not know, on which of the two
+        // the bond is defined, since bonds are stored only on one partner
+        DeleteAllBonds{p1->p.vs_relative.to_particle_id,
+                       p2->p.vs_relative.to_particle_id},
+        DeleteAllBonds{p2->p.vs_relative.to_particle_id,
+                       p1->p.vs_relative.to_particle_id},
+    };
+  }
+#endif // VIRTUAL_SITES_RELATIVE
+  return {};
+}
+
+// Handler for the different delete events
+class execute : public boost::static_visitor<> {
+public:
+  void operator()(DeleteBond const &d) const {
+    auto p = cell_structure.get_local_particle(d.particle_id);
+    if (!p)
+      return;
+    local_remove_bond(*p, {d.bond_type, d.bond_partner_id});
+  }
+  void operator()(DeleteAllBonds const &d) const {
+    auto p = cell_structure.get_local_particle(d.particle_id_1);
+    if (!p)
+      return;
+    local_remove_pair_bonds_to(*p, d.particle_id_2);
+  }
+};
+
+void process_queue() {
+  if (breakage_specs.empty())
+    return;
+
+  auto global_queue = gather_global_queue(queue);
+
+  // Construct delete actions from breakage queue
+  ActionSet actions = {};
+  for (auto const &e : global_queue) {
+    // Convert to merge() once we are on C++17
+    auto to_add = actions_for_breakage(e);
+    actions.insert(to_add.begin(), to_add.end());
+  }
+
+  // Execute actions
+  for (auto const &a : actions) {
+    boost::apply_visitor(execute{}, a);
+  }
+}
+} // namespace BondBreakage

src/core/bond_breakage.hpp

@@ -0,0 +1,50 @@
+/*
+ * Copyright (C) 2022 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <memory>
+#include <unordered_map>
+
+namespace BondBreakage {
+
+enum class ActionType {
+  NONE = 0,
+  DELETE_BOND = 1,
+  REVERT_BIND_AT_POINT_OF_COLLISION = 2
+};
+
+struct BreakageSpec {
+  double breakage_length;
+  ActionType action_type;
+};
+
+extern std::unordered_map<int, std::shared_ptr<BreakageSpec>> breakage_specs;
+
+/** @brief Check if the bond between the particles should break, if yes, queue
+ *  it.
+ */
+bool check_and_handle_breakage(int particle_id, int bond_partner_id,
+                               int bond_type, double distance);
+
+void clear_queue();
+
+void process_queue();
+
+} // namespace BondBreakage

src/core/forces.cpp

@@ -26,6 +26,7 @@
 
 #include "EspressoSystemInterface.hpp"
 
+#include "bond_breakage.hpp"
 #include "cells.hpp"
 #include "collision.hpp"
 #include "comfixed_global.hpp"
@@ -148,7 +149,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) {
 #ifdef COLLISION_DETECTION
   prepare_local_collision_queue();
 #endif
-
+  BondBreakage::clear_queue();
   auto particles = cell_structure.local_particles();
   auto ghost_particles = cell_structure.ghost_particles();
 #ifdef ELECTROSTATICS

src/core/forces_inline.hpp

@@ -28,6 +28,7 @@
 
 #include "forces.hpp"
 
+#include "bond_breakage.hpp"
 #include "bonded_interactions/bonded_interaction_data.hpp"
 #include "bonded_interactions/thermalized_bond_kernel.hpp"
 #include "immersed_boundary/ibm_tribend.hpp"
@@ -411,6 +412,15 @@ inline bool add_bonded_four_body_force(Bonded_IA_Parameters const &iaparams,
 
 inline bool add_bonded_force(Particle &p1, int bond_id,
                              Utils::Span<Particle *> partners) {
+
+  // Consider for bond breakage
+  if (partners.size() == 1) {
+    auto d = box_geo.get_mi_vector(p1.r.p, partners[0]->r.p).norm();
+    if (BondBreakage::check_and_handle_breakage(
+            p1.p.identity, partners[0]->p.identity, bond_id, d))
+      return false;
+  }
+
   auto const &iaparams = *bonded_ia_params.at(bond_id);
 
   switch (number_of_partners(iaparams)) {

src/core/integrate.cpp

@@ -35,6 +35,7 @@
 
 #include "ParticleRange.hpp"
 #include "accumulators.hpp"
+#include "bond_breakage.hpp"
 #include "bonded_interactions/rigid_bond.hpp"
 #include "cells.hpp"
 #include "collision.hpp"
@@ -325,6 +326,7 @@ int integrate(int n_steps, int reuse_forces) {
 #ifdef COLLISION_DETECTION
       handle_collisions();
 #endif
+      BondBreakage::process_queue();
     }
 
     integrated_steps++;

src/core/particle_data.cpp

@@ -199,6 +199,28 @@ struct RemoveBond {
   void serialize(Archive &ar, long int) { ar & bond; }
 };
 
+/**
+ * @brief Delete pair bonds to a specific partner
+ */
+struct RemovePairBondsTo {
+   int other_pid;
+
+    void operator()(Particle &p) const {
+      using Bond = std::vector<int>;
+      std::vector<Bond> to_delete;
+      for (auto b: p.bonds()) {
+         if (b.partner_ids().size() == 1 and b.partner_ids()[0] == other_pid)
+           to_delete.push_back(Bond{b.bond_id(),other_pid});
+      }
+      for (auto b: to_delete) {
+        RemoveBond{b}(p);
+      }
+    }
+    template<class Archive>
+            void serialize(Archive &ar, long int) {
+        ar & other_pid;
+    }
+};
 
 /**
  * @brief Delete all bonds.
@@ -325,6 +347,14 @@ struct UpdateVisitor : public boost::static_visitor<void> {
 };
 } // namespace
 
+void local_remove_bond(Particle &p, std::vector<int> const &bond) {
+  RemoveBond{bond}(p);
+}
+
+void local_remove_pair_bonds_to(Particle &p, int other_pid) {
+  RemovePairBondsTo{other_pid}(p);
+}
+
 void mpi_send_update_message_local(int node, int id) {
   if (node == comm_cart.rank()) {
     UpdateMessage msg{};