Implement manual trigger of bond breakage

Co-authored-by: Rudolf Weeber <weeber@icp.uni-stuttgart.de>
Co-authored-by: Jean-Noël Grad <jgrad@icp.uni-stuttgart.de>

src/core/bond_breakage/bond_breakage.cpp

@@ -204,4 +204,42 @@ void BondBreakage::process_queue_impl(System::System &system) {
   }
 }
 
+static bool bond_handler(BondBreakage &bond_breakage, Particle &p,
+                         std::span<Particle *> partners, int bond_id,
+                         BoxGeometry const &box_geo) {
+  auto retval = false;
+  if (partners.size() == 1u) { // pair bonds
+    auto d = box_geo.get_mi_vector(p.pos(), partners[0]->pos()).norm();
+    retval = bond_breakage.check_and_handle_breakage(
+        p.id(), {{partners[0]->id(), std::nullopt}}, bond_id, d);
+  } else if (partners.size() == 2u) { // angle bond
+    auto d =
+        box_geo.get_mi_vector(partners[0]->pos(), partners[1]->pos()).norm();
+    retval = bond_breakage.check_and_handle_breakage(
+        p.id(), {{partners[0]->id(), partners[1]->id()}}, bond_id, d);
+  }
+  return retval;
+}
+
+void BondBreakage::execute_bond_breakage(System::System &system) {
+  system.cell_structure->update_ghosts_and_resort_particle(
+      system.get_global_ghost_flags());
+
+  // Clear the bond breakage queue
+  clear_queue();
+
+  // Create the bond kernel function (the bond handler)
+  auto bond_kernel = [&](Particle &p, int bond_id,
+                         std::span<Particle *> partners) {
+    bond_handler(*this, p, partners, bond_id, *system.box_geo);
+    return false;
+  };
+
+  // Use the CellStructure::bond_loop to process bonds
+  system.cell_structure->bond_loop(bond_kernel);
+
+  // Process the bond breakage queue
+  process_queue(system);
+}
+
 } // namespace BondBreakage

src/core/bond_breakage/bond_breakage.hpp

@@ -96,6 +96,8 @@ class BondBreakage {
 
   void clear_queue() { m_queue.clear(); }
 
+  void execute_bond_breakage(System::System &system);
+
   void process_queue(System::System &system) {
     if (not breakage_specs.empty()) {
       process_queue_impl(system);

src/python/espressomd/bond_breakage.py

@@ -42,6 +42,10 @@ class BreakageSpec(ScriptInterfaceHelper):
 class BreakageSpecs(ScriptObjectMap):
     _so_name = "BondBreakage::BreakageSpecs"
 
+    def execute(self):
+        """Execute the bond breakage on the current state of the system"""
+        self.call_method("execute")
+
     def _get_key(self, key):
         """Convert a bond object to a bond id."""
         if isinstance(key, interactions.BondedInteraction):

src/script_interface/bond_breakage/BreakageSpecs.hpp

@@ -55,6 +55,15 @@ class BreakageSpecs : public BreakageSpecsBase_t {
     m_bond_breakage = std::make_shared<::BondBreakage::BondBreakage>();
     restore_from_checkpoint(params);
   }
+  Variant do_call_method(std::string const &name,
+                         VariantMap const &parameters) override {
+    if (name == "execute") {
+      context()->parallel_try_catch(
+          [this]() { m_bond_breakage->execute_bond_breakage(get_system()); });
+      return {};
+    }
+    return Base::do_call_method(name, parameters);
+  }
 
 private:
   void on_bind_system(::System::System &system) override {

testsuite/python/bond_breakage.py

@@ -22,7 +22,7 @@
 import numpy as np
 import espressomd.interactions
 from espressomd.bond_breakage import BreakageSpec
-from espressomd.interactions import HarmonicBond, AngleHarmonic
+from espressomd.interactions import HarmonicBond, AngleHarmonic, Dihedral
 
 
 class BondBreakageCommon:
@@ -33,16 +33,18 @@ class BondBreakageCommon:
 
 
 @utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE")
-class BondBreakage(BondBreakageCommon, ut.TestCase):
+class BondBreakageTest(BondBreakageCommon, ut.TestCase):
 
     @classmethod
     def setUpClass(cls):
         pos1 = cls.system.box_l / 2 - 0.5
         pos2 = cls.system.box_l / 2 + 0.5
         pos3 = cls.system.box_l / 2 - 1.5
+        pos4 = cls.system.box_l / 2 - 1.5 + [0.5, 0., 0.]
         cls.p1 = cls.system.part.add(pos=pos1)
         cls.p2 = cls.system.part.add(pos=pos2)
         cls.p3 = cls.system.part.add(pos=pos3)
+        cls.p4 = cls.system.part.add(pos=pos4)
 
         cls.p1v = cls.system.part.add(pos=pos1)
         cls.p1v.vs_auto_relate_to(cls.p1)
@@ -53,9 +55,11 @@ def setUpClass(cls):
         cls.h1 = HarmonicBond(k=1, r_0=0)
         cls.h2 = HarmonicBond(k=1, r_0=0)
         cls.angle = AngleHarmonic(bend=1, phi0=np.pi)
+        cls.dihe = Dihedral(bend=1, mult=1, phase=np.pi)
         cls.system.bonded_inter.add(cls.h1)
         cls.system.bonded_inter.add(cls.h2)
         cls.system.bonded_inter.add(cls.angle)
+        cls.system.bonded_inter.add(cls.dihe)
 
     @classmethod
     def tearDownClass(cls):
@@ -84,6 +88,7 @@ def test_00_interface(self):
         self.system.bond_breakage.clear()
         self.assertEqual(len(self.system.bond_breakage), 0)
         self.assertEqual(self.system.bond_breakage.keys(), [])
+        self.assertIsNone(self.system.bond_breakage.call_method("unknown"))
         with self.assertRaisesRegex(ValueError, "Key has to be of type 'int'"):
             self.system.bond_breakage[None]
         with self.assertRaisesRegex(ValueError, "Key has to be of type 'int'"):
@@ -134,6 +139,10 @@ def test_delete_bond(self):
         system.integrator.run(1)
         self.assertEqual(self.p2.bonds, ())
 
+        self.p1.bonds = [(self.h1, self.p2)]
+        system.bond_breakage.execute()
+        self.assertEqual(self.p1.bonds, ())
+
     def test_delete_angle_bond(self):
         system = self.system
 
@@ -143,13 +152,23 @@ def test_delete_angle_bond(self):
         self.p1.bonds = ((self.angle, self.p2, self.p3),)
         bonds = self.p1.bonds
         system.integrator.run(1)
-        # should still be there. Not bfeyond breakage disst
+        # should still be there. Not beyond breakage dist
         self.assertEqual(self.p1.bonds, bonds)
         self.system.bond_breakage.clear()
         system.bond_breakage[self.angle] = BreakageSpec(
             breakage_length=1, action_type="delete_bond")
         system.integrator.run(1)
         self.assertEqual(self.p1.bonds, ())
+        # manual trigger
+        self.p1.bonds = ((self.angle, self.p2, self.p3),)
+        system.bond_breakage.execute()
+        self.assertEqual(self.p1.bonds, ())
+        # manual trigger for dihedral angle bond (no-op)
+        self.p2.bonds = ((self.dihe, self.p1, self.p3, self.p4),)
+        bonds = self.p2.bonds
+        system.bond_breakage.execute()
+        self.assertEqual(self.p2.bonds, bonds)
+        self.p2.bonds = ()
 
     def test_revert_bind_at_point_of_collision_pair(self):
         system = self.system