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{};