From a4a3e6bf6e2db6f95a9fb246724b392e7ba097e6 Mon Sep 17 00:00:00 2001 From: slottsborgen Date: Fri, 24 Nov 2023 10:56:30 +0100 Subject: [PATCH] Remove the bind_three_particles collision method --- doc/sphinx/advanced_methods.rst | 37 ---- src/core/collision.cpp | 164 +----------------- src/core/collision.hpp | 11 +- src/python/espressomd/collision_detection.py | 15 +- .../CollisionDetection.hpp | 7 - testsuite/python/collision_detection.py | 55 ------ .../python/collision_detection_interface.py | 21 +-- 7 files changed, 8 insertions(+), 302 deletions(-) diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index d5874a9a8c8..a4553a0a138 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -126,43 +126,6 @@ Several modes are available for different types of binding. part_type_to_be_glued=3, part_type_after_glueing=4) -* ``"bind_three_particles"`` allows for the creation of agglomerates which maintain - their shape similarly to those create by the mode ``"bind_at_point_of_collision"``. - The present approach works without virtual sites. Instead, for each two-particle - collision, the surrounding is searched for a third particle. If one is found, - angular bonds are placed to maintain the local shape. - If all three particles are within the cutoff distance, an angle bond is added - on each of the three particles in addition - to the distance based bonds between the particle centers. - If two particles are within the cutoff of a central particle (e.g., chain of three particles) - an angle bond is placed on the central particle. - The angular bonds being added are determined from the angle between the particles. - This method does not depend on the particles' rotational - degrees of freedom being integrated. Virtual sites are not required. - The method, along with the corresponding bonds are setup as follows:: - - first_angle_bond_id = 0 - n_angle_bonds = 181 # 0 to 180 degrees in one degree steps - for i in range(0, n_angle_bonds, 1): - bond_id = first_angle_bond_id + i - system.bonded_inter[bond_id] = espressomd.interactions.AngleHarmonic( - bend=1., phi0=float(i) / float(n_angle_bonds - 1) * np.pi) - - bond_centers = espressomd.interactions.HarmonicBond(k=1., r_0=0.1, r_cut=0.5) - system.bonded_inter.add(bond_centers) - - system.collision_detection.set_params( - mode="bind_three_particles", - bond_centers=bond_centers, - bond_three_particles=first_angle_bond_id, - three_particle_binding_angle_resolution=n_angle_bonds, - distance=0.1) - - Important: The bonds for the angles are mapped via their numerical bond ids. - In this example, ids from 0 to 180 are used. All other bonds required for - the simulation need to be added to the system after those bonds. In particular, - this applies to the bonded interaction passed via ``bond_centers`` - The following limitations currently apply for the collision detection: diff --git a/src/core/collision.cpp b/src/core/collision.cpp index 73602544cc4..62b4d6b3a4f 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -167,26 +167,7 @@ void Collision_parameters::initialize() { (get_bond_num_partners(collision_params.bond_vs) != 1 and get_bond_num_partners(collision_params.bond_vs) != 2)) { throw std::runtime_error("The bond type to be used for binding virtual " - "sites needs to be a pair or three-particle bond"); - } - - if (collision_params.mode == CollisionModeType::BIND_THREE_PARTICLES) { - if (collision_params.bond_three_particles + - collision_params.three_particle_angle_resolution > - bonded_ia_params.size()) { - throw std::runtime_error( - "Insufficient bonds defined for three particle binding"); - } - - for (int i = collision_params.bond_three_particles; - i < collision_params.bond_three_particles + - collision_params.three_particle_angle_resolution; - i++) { - if (get_bond_num_partners(i) != 2) { - throw std::runtime_error( - "The bonds for three particle binding need to be angle bonds."); - } - } + "sites needs to be a pair bond"); } // Create particle types @@ -275,74 +256,6 @@ static void bind_at_point_of_collision_calc_vs_pos(Particle const &p1, pos2 = p1.pos() - vec21 * (1. - collision_params.vs_placement); } -// Considers three particles for three_particle_binding and performs -// the binding if the criteria are met -static void coldet_do_three_particle_bond(Particle &p, Particle const &p1, - Particle const &p2, - BoxGeometry const &box_geo) { - // If p1 and p2 are not closer or equal to the cutoff distance, skip - // p1: - if (box_geo.get_mi_vector(p.pos(), p1.pos()).norm() > - collision_params.distance) - return; - // p2: - if (box_geo.get_mi_vector(p.pos(), p2.pos()).norm() > - collision_params.distance) - return; - - // Check, if there already is a three-particle bond centered on p - // with p1 and p2 as partners. If so, skip this triplet. - // Note that the bond partners can appear in any order. - - /* Check if a bond is a bond placed by the collision detection */ - auto is_collision_bond = [](BondView const &bond) { - return (bond.bond_id() >= collision_params.bond_three_particles) and - (bond.bond_id() <= - collision_params.bond_three_particles + - collision_params.three_particle_angle_resolution); - }; - /* Check if the bond is between the particles we are currently considering */ - auto has_same_partners = [id1 = p1.id(), - id2 = p2.id()](BondView const &bond) { - auto const partner_ids = bond.partner_ids(); - - return ((partner_ids[0] == id1) and (partner_ids[1] == id2)) or - ((partner_ids[0] == id2) and (partner_ids[1] == id1)); - }; - - auto const &bonds = p.bonds(); - if (std::any_of(bonds.begin(), bonds.end(), [=](auto const &bond) { - return is_collision_bond(bond) and has_same_partners(bond); - })) { - return; - } - - // If we are still here, we need to create angular bond - // First, find the angle between the particle p, p1 and p2 - - /* vector from p to p1 */ - auto const vec1 = box_geo.get_mi_vector(p.pos(), p1.pos()).normalize(); - /* vector from p to p2 */ - auto const vec2 = box_geo.get_mi_vector(p.pos(), p2.pos()).normalize(); - - auto const cosine = std::clamp(vec1 * vec2, -TINY_COS_VALUE, TINY_COS_VALUE); - - // Bond angle - auto const phi = acos(cosine); - - // We find the bond id by dividing the range from 0 to pi in - // three_particle_angle_resolution steps and by adding the id - // of the bond for zero degrees. - auto const bond_id = static_cast( - floor(0.5 + phi / std::numbers::pi * - (collision_params.three_particle_angle_resolution - 1)) + - collision_params.bond_three_particles); - - // Create the bond - const std::array bondT = {{p1.id(), p2.id()}}; - p.bonds().insert({bond_id, bondT}); -} - #ifdef VIRTUAL_SITES_RELATIVE static void place_vs_and_relate_to_particle(CellStructure &cell_structure, BoxGeometry const &box_geo, @@ -414,74 +327,8 @@ std::vector gather_global_collision_queue() { return res; } -static void three_particle_binding_do_search(Cell *basecell, Particle &p1, - Particle &p2, - BoxGeometry const &box_geo) { - auto handle_cell = [&p1, &p2, &box_geo](Cell *c) { - for (auto &p : c->particles()) { - // Skip collided particles themselves - if ((p.id() == p1.id()) or (p.id() == p2.id())) { - continue; - } - - // We need all cyclical permutations, here (bond is placed on 1st - // particle, order of bond partners does not matter, so we don't need - // non-cyclic permutations). - // @ref coldet_do_three_particle_bond checks the bonding criterion and if - // the involved particles are not already bonded before it binds them. - if (!p.is_ghost()) { - coldet_do_three_particle_bond(p, p1, p2, box_geo); - } - - if (!p1.is_ghost()) { - coldet_do_three_particle_bond(p1, p, p2, box_geo); - } - - if (!p2.is_ghost()) { - coldet_do_three_particle_bond(p2, p, p1, box_geo); - } - } - }; - - /* Search the base cell ... */ - handle_cell(basecell); - - /* ... and all the neighbors. */ - for (auto &n : basecell->neighbors().all()) { - handle_cell(n); - } -} - -// Goes through the collision queue and for each pair in it -// looks for a third particle by using the domain decomposition -// cell system. If found, it performs three particle binding -static void three_particle_binding_domain_decomposition( - CellStructure &cell_structure, BoxGeometry const &box_geo, - std::vector const &gathered_queue) { - - for (auto &c : gathered_queue) { - // If we have both particles, at least as ghosts, Get the corresponding cell - // indices - if (cell_structure.get_local_particle(c.pp1) and - cell_structure.get_local_particle(c.pp2)) { - Particle &p1 = *cell_structure.get_local_particle(c.pp1); - Particle &p2 = *cell_structure.get_local_particle(c.pp2); - auto cell1 = cell_structure.find_current_cell(p1); - auto cell2 = cell_structure.find_current_cell(p2); - - if (cell1) - three_particle_binding_do_search(cell1, p1, p2, box_geo); - if (cell2 and cell1 != cell2) - three_particle_binding_do_search(cell2, p1, p2, box_geo); - - } // If local particles exist - } // Loop over total collisions -} - // Handle the collisions stored in the queue void handle_collisions(CellStructure &cell_structure) { - auto &system = System::get_system(); - auto const &box_geo = *system.box_geo; // Note that the glue to surface mode adds bonds between the centers // but does so later in the process. This is needed to guarantee that // a particle can only be glued once, even if queued twice in a single @@ -503,6 +350,8 @@ void handle_collisions(CellStructure &cell_structure) { // Virtual sites based collision schemes #ifdef VIRTUAL_SITES_RELATIVE + auto &system = System::get_system(); + auto const &box_geo = *system.box_geo; auto const min_global_cut = system.get_min_global_cut(); if ((collision_params.mode == CollisionModeType::BIND_VS) || (collision_params.mode == CollisionModeType::GLUE_TO_SURF)) { @@ -660,13 +509,6 @@ void handle_collisions(CellStructure &cell_structure) { } // are we in one of the vs_based methods #endif // defined VIRTUAL_SITES_RELATIVE - // three-particle-binding part - if (collision_params.mode == CollisionModeType::BIND_THREE_PARTICLES) { - auto gathered_queue = gather_global_collision_queue(); - three_particle_binding_domain_decomposition(cell_structure, box_geo, - gathered_queue); - } // if TPB - local_collision_queue.clear(); } diff --git a/src/core/collision.hpp b/src/core/collision.hpp index 339e47f81a8..29c0bfe9bd1 100644 --- a/src/core/collision.hpp +++ b/src/core/collision.hpp @@ -41,15 +41,13 @@ enum class CollisionModeType : int { BIND_VS = 2, /** @brief Glue a particle to a specific spot on another particle. */ GLUE_TO_SURF = 3, - /** @brief Three particle binding mode. */ - BIND_THREE_PARTICLES = 4 }; class Collision_parameters { public: Collision_parameters() : mode(CollisionModeType::OFF), distance(0.), distance2(0.), - bond_centers(-1), bond_vs(-1), bond_three_particles(-1) {} + bond_centers(-1), bond_vs(-1) {} /// collision protocol CollisionModeType mode; @@ -75,13 +73,6 @@ class Collision_parameters { int part_type_to_attach_vs_to; /// Particle type to which the newly glued particle is converted int part_type_after_glueing; - /// First bond type (for zero degrees) used for the three-particle bond - /// (angle potential) - int bond_three_particles; - /// Number of angle bonds to use (angular resolution) - /// different angle bonds with different equilibrium angles - /// Are expected to have ids immediately following to bond_three_particles - int three_particle_angle_resolution; /** Placement of virtual sites for MODE_VS. * 0=on same particle as related to, * 1=on collision partner, diff --git a/src/python/espressomd/collision_detection.py b/src/python/espressomd/collision_detection.py index d7470e5b150..24dde5de809 100644 --- a/src/python/espressomd/collision_detection.py +++ b/src/python/espressomd/collision_detection.py @@ -63,7 +63,7 @@ def set_params(self, **kwargs): Parameters ---------- - mode : :obj:`str`, {"off", "bind_centers", "bind_at_point_of_collision", "bind_three_particles", "glue_to_surface"} + mode : :obj:`str`, {"off", "bind_centers", "bind_at_point_of_collision", "glue_to_surface"} Collision detection mode distance : :obj:`float` @@ -92,22 +92,13 @@ def set_params(self, **kwargs): distance_glued_particle_to_vs : :obj:`float` Distance for ``"glue_to_surface"`` mode. See user guide. - bond_three_particles : :obj:`espressomd.interactions.BondedInteraction` - First angular bond for the ``"bind_three_particles"`` mode. See - user guide - - three_particle_binding_angle_resolution : :obj:`int` - Resolution for the angular bonds (mode ``"bind_three_particles"``). - Resolution+1 bonds are needed to accommodate the case of 180 degrees - angles - """ if "mode" not in kwargs: raise ValueError( "Collision mode must be specified via the 'mode' argument") # Convert bonds to bond ids - for name in ["bond_centers", "bond_vs", "bond_three_particle_binding"]: + for name in ["bond_centers", "bond_vs"]: if name in kwargs: if isinstance(kwargs[name], BondedInteraction): kwargs[name] = kwargs[name]._bond_id @@ -117,7 +108,7 @@ def get_parameter(self, name): """Gets a single parameter from the collision detection.""" value = super().get_parameter(name) - if name in ["bond_centers", "bond_vs", "bond_three_particle_binding"]: + if name in ["bond_centers", "bond_vs"]: if value == -1: # Not defined value = None else: diff --git a/src/script_interface/collision_detection/CollisionDetection.hpp b/src/script_interface/collision_detection/CollisionDetection.hpp index 2ab4e93752b..b523c3d96f9 100644 --- a/src/script_interface/collision_detection/CollisionDetection.hpp +++ b/src/script_interface/collision_detection/CollisionDetection.hpp @@ -47,7 +47,6 @@ class CollisionDetection : public AutoParameters { {CollisionModeType::BIND_CENTERS, "bind_centers"}, {CollisionModeType::BIND_VS, "bind_at_point_of_collision"}, {CollisionModeType::GLUE_TO_SURF, "glue_to_surface"}, - {CollisionModeType::BIND_THREE_PARTICLES, "bind_three_particles"}, }; std::unordered_map cd_name_to_mode; std::unordered_map { "part_type_to_be_glued", "part_type_to_attach_vs_to", "part_type_after_glueing", "distance", "distance_glued_particle_to_vs"}}, - {CollisionModeType::BIND_THREE_PARTICLES, - {"mode", "bond_centers", "distance", "bond_three_particles", - "three_particle_binding_angle_resolution"}}, }; public: @@ -83,9 +79,6 @@ class CollisionDetection : public AutoParameters { {"bond_centers", collision_params.bond_centers}, {"bond_vs", collision_params.bond_vs}, - {"bond_three_particles", collision_params.bond_three_particles}, - {"three_particle_binding_angle_resolution", - collision_params.three_particle_angle_resolution}, {"distance", collision_params.distance}, {"distance_glued_particle_to_vs", diff --git a/testsuite/python/collision_detection.py b/testsuite/python/collision_detection.py index 1f67d676f4c..1f0a2ae2411 100644 --- a/testsuite/python/collision_detection.py +++ b/testsuite/python/collision_detection.py @@ -684,61 +684,6 @@ def test_glue_to_surface_random(self): # Tidy system.non_bonded_inter[0, 0].lennard_jones.deactivate() - def test_bind_three_particles(self): - system = self.system - # Setup particles - system.part.clear() - dx = np.array((1, 0, 0)) - dy = np.array((0, 1, 0)) - a = np.array((0.499, 0.499, 0.499)) - b = a + 0.1 * dx - c = a + 0.03 * dx + 0.03 * dy - d = a + 0.03 * dx - 0.03 * dy - e = a - 0.1 * dx - - system.part.add(id=0, pos=a) - system.part.add(id=1, pos=b) - system.part.add(id=2, pos=c) - system.part.add(id=3, pos=d) - system.part.add(id=4, pos=e) - - # Setup bonds - res = 181 - for i in range(res): - system.bonded_inter[i + 4] = espressomd.interactions.AngleHarmonic( - bend=1, phi0=float(i) / (res - 1) * np.pi) - cutoff = 0.11 - system.collision_detection.set_params( - mode="bind_three_particles", bond_centers=self.bond_center, - bond_three_particles=4, three_particle_binding_angle_resolution=res, distance=cutoff) - self.get_state_set_state_consistency() - - system.time_step = 1E-6 - system.integrator.run(1, recalc_forces=True) - self.verify_triangle_binding(cutoff, system.bonded_inter[4], res) - # Make sure no extra bonds appear - system.integrator.run(1, recalc_forces=True) - self.verify_triangle_binding(cutoff, system.bonded_inter[4], res) - - # Place the particles in two steps and make sure, the bonds are the - # same - system.part.clear() - system.part.add(id=0, pos=a) - system.part.add(id=2, pos=c) - system.part.add(id=3, pos=d) - system.integrator.run(1, recalc_forces=True) - - system.part.add(id=4, pos=e) - system.part.add(id=1, pos=b) - system.cell_system.set_regular_decomposition() - system.integrator.run(1, recalc_forces=True) - self.verify_triangle_binding(cutoff, system.bonded_inter[4], res) - system.cell_system.set_n_square() - system.part.all().bonds = () - system.integrator.run(1, recalc_forces=True) - self.verify_triangle_binding(cutoff, system.bonded_inter[4], res) - system.time_step = self.time_step - def verify_triangle_binding(self, distance, first_bond, angle_res): system = self.system # Gather pairs diff --git a/testsuite/python/collision_detection_interface.py b/testsuite/python/collision_detection_interface.py index c035a413bc2..25df82dad93 100644 --- a/testsuite/python/collision_detection_interface.py +++ b/testsuite/python/collision_detection_interface.py @@ -46,11 +46,6 @@ class CollisionDetection(ut.TestCase): "bond_vs": bond_harmonic, "bond_centers": bond_harmonic, "part_type_vs": 1, "distance": 0.1, "vs_placement": 0.1 }, - "bind_three_particles": { - "bond_centers": bond_harmonic, "bond_three_particles": 0, - "three_particle_binding_angle_resolution": bond_angle_resolution, - "distance": 0.1 - }, "glue_to_surface": { "distance": 0.1, "distance_glued_particle_to_vs": 0.02, "bond_centers": bond_harmonic, "bond_vs": bond_harmonic, @@ -136,7 +131,7 @@ def test_bind_at_point_of_collision(self): with self.assertRaisesRegex(RuntimeError, "Bond in parameter 'bond_vs' was not added to the system"): bond = espressomd.interactions.HarmonicBond(k=1., r_0=0.1) self.set_coldet("bind_at_point_of_collision", bond_vs=bond) - with self.assertRaisesRegex(RuntimeError, "bond type to be used for binding virtual sites needs to be a pair or three-particle bond"): + with self.assertRaisesRegex(RuntimeError, "bond type to be used for binding virtual sites needs to be a pair bond"): self.set_coldet( "bind_at_point_of_collision", bond_vs=self.bond_dihe) with self.assertRaisesRegex(ValueError, "type for virtual sites needs to be >=0"): @@ -151,20 +146,6 @@ def test_bind_at_point_of_collision_norotation(self): with self.assertRaisesRegex(RuntimeError, "require the VIRTUAL_SITES_RELATIVE feature"): self.set_coldet("bind_at_point_of_collision") - def test_bind_three_particles(self): - self.set_coldet("bind_three_particles", distance=0.5) - with self.assertRaisesRegex(RuntimeError, "Insufficient bonds defined for three particle binding"): - self.set_coldet( - "bind_three_particles", - three_particle_binding_angle_resolution=self.bond_angle_resolution + - 1000) - with self.assertRaisesRegex(RuntimeError, "The bonds for three particle binding need to be angle bonds"): - self.set_coldet( - "bind_three_particles", - three_particle_binding_angle_resolution=self.bond_angle_resolution + 1) - # check if original parameters have been preserved - self.check_stored_parameters("bind_three_particles", distance=0.5) - @utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE") def test_glue_to_surface(self): self.set_coldet("glue_to_surface", distance=0.5)