diff --git a/src/core/bond_breakage/bond_breakage.cpp b/src/core/bond_breakage/bond_breakage.cpp index 296eec5282..9dcaf7837b 100644 --- a/src/core/bond_breakage/bond_breakage.cpp +++ b/src/core/bond_breakage/bond_breakage.cpp @@ -79,7 +79,6 @@ static ActionSet actions_for_breakage(CellStructure const &cell_structure, return bond_partners[1]; }; // optional for second partner engaged - // Handle different action types if (spec.action_type == ActionType::DELETE_BOND) { if (is_angle_bond(e.bond_partners)) { return {DeleteAngleBond{e.particle_id, @@ -88,54 +87,59 @@ static ActionSet actions_for_breakage(CellStructure const &cell_structure, } return {DeleteBond{e.particle_id, *(e.bond_partners[0]), e.bond_type}}; } + #ifdef ESPRESSO_VIRTUAL_SITES_RELATIVE // revert bind at point of collision for pair bonds - if (spec.action_type == ActionType::REVERT_BIND_AT_POINT_OF_COLLISION and - not is_angle_bond(e.bond_partners)) { - // 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_partners[0])); - if (p1 and p2) { - if (not p1->is_virtual() or not p2->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 {}; - } + if (spec.action_type == ActionType::REVERT_BIND_AT_POINT_OF_COLLISION) { + if (not is_angle_bond(e.bond_partners)) { + // 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_partners[0])); + if (p1 and p2) { + if (not p1->is_virtual() or not p2->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_partners[0]), e.bond_type}, - // Bond between base particles. We do not know, on which of these - // the bond is defined, since bonds are stored only on one partner - DeleteAllBonds{p1->vs_relative().to_particle_id, - p2->vs_relative().to_particle_id}, - DeleteAllBonds{p2->vs_relative().to_particle_id, - p1->vs_relative().to_particle_id}, - }; - } - } else { - // revert bind at point of collision for angle bonds - auto vs = cell_structure.get_local_particle(e.particle_id); - auto p1 = cell_structure.get_local_particle(*(e.bond_partners[0])); - auto p2 = cell_structure.get_local_particle(*(e.bond_partners[1])); - if (p1 and p2) { - if (not vs->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_partners[0]), e.bond_type}, + // Bond between base particles. We do not know, on which of these + // the bond is defined, since bonds are stored only on one partner + DeleteAllBonds{p1->vs_relative().to_particle_id, + p2->vs_relative().to_particle_id}, + DeleteAllBonds{p2->vs_relative().to_particle_id, + p1->vs_relative().to_particle_id}, + }; } + } else { + // revert bind at point of collision for angle bonds + auto vs = cell_structure.get_local_particle(e.particle_id); + auto p1 = cell_structure.get_local_particle(*(e.bond_partners[0])); + auto p2 = cell_structure.get_local_particle(*(e.bond_partners[1])); + if (p1 and p2) { + if (not vs->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 {// Angle bond on the virtual site - DeleteAngleBond{e.particle_id, {p1->id(), p2->id()}, e.bond_type}, - // Bond between base particles. We do not know, on which of these - // the bond is defined, since bonds are stored only on one partner - DeleteAllBonds{p1->id(), p2->id()}, - DeleteAllBonds{p2->id(), p1->id()}}; + return { + // Angle bond on the virtual site + DeleteAngleBond{e.particle_id, {p1->id(), p2->id()}, e.bond_type}, + // Bond between base particles. We do not know, on which of these + // the bond is defined, since bonds are stored only on one partner + DeleteAllBonds{p1->id(), p2->id()}, + DeleteAllBonds{p2->id(), p1->id()}}; + } } } #endif // ESPRESSO_VIRTUAL_SITES_RELATIVE diff --git a/testsuite/python/bond_breakage.py b/testsuite/python/bond_breakage.py index a071431331..e0e57491de 100644 --- a/testsuite/python/bond_breakage.py +++ b/testsuite/python/bond_breakage.py @@ -253,18 +253,15 @@ def count_bonds(self, pairs): return bonds_count def setUp(self): - - box_vol = self.system.box_l[0]**3. phi = 0.4 - r = 1. - solid_vol = phi * box_vol + solid_vol = phi * self.system.volume() part_vol = 4 / 3 * np.pi * r**3 part_num = int(solid_vol / part_vol) np.random.seed(seed=678) for i in range(part_num): - pos = np.random.rand(3) * self.system.box_l[0] + pos = np.random.rand(3) * self.system.box_l self.system.part.add(pos=pos) self.system.non_bonded_inter[0, 0].lennard_jones.set_params( @@ -360,5 +357,42 @@ def test_vs_bonds(self): np.testing.assert_equal(bonds_dist, bonds_count) +class BreakageAPI(BondBreakageCommon, ut.TestCase): + + def setUp(self): + self.system.box_l = 3 * [20] + self.system.min_global_cut = 0.1 + self.system.time_step = 0.01 + self.system.cell_system.skin = 0.4 + + def tearDown(self): + self.system.part.clear() + self.system.bonded_inter.clear() + + def test_bond_deletion(self): + p1 = self.system.part.add(pos=[0., 0., 0.], v=[0., 0., 0.]) + p2 = self.system.part.add(pos=[0., 0., 1.], v=[0., 0., 1.]) + bond = espressomd.interactions.HarmonicBond(k=1., r_0=1., r_cut=1.2) + self.system.bonded_inter.add(bond) + p1.bonds = ((bond, p2)) + self.system.bond_breakage[bond] = BreakageSpec( + breakage_length=1.1, action_type="delete_bond") + self.system.integrator.run(100) + assert np.linalg.norm(p2.pos - p1.pos) > bond.r_cut + self.assertEqual(len(p1.bonds), 0) + + def test_none(self): + p1 = self.system.part.add(pos=[0., 0., 0.], v=[0., 0., 0.]) + p2 = self.system.part.add(pos=[0., 0., 1.], v=[0., 0., 1.]) + bond = espressomd.interactions.HarmonicBond(k=1., r_0=1., r_cut=1.2) + self.system.bonded_inter.add(bond) + p1.bonds = ((bond, p2)) + self.system.bond_breakage[bond] = BreakageSpec( + breakage_length=1.1, action_type="none") + self.system.integrator.run(100) + assert np.linalg.norm(p2.pos - p1.pos) > bond.r_cut + self.assertEqual(len(p1.bonds), 1) + + if __name__ == "__main__": ut.main()