Skip to content

Commit

Permalink
Remove the bind_three_particles collision method (#4823)
Browse files Browse the repository at this point in the history 
Step one of #4483 

Description of changes:
- API change: the bind_three_particles feature was removed
  • Loading branch information
@kodiakhq
kodiakhq[bot] authored Jul 4, 2024
2 parents e4d63b2 + a4a3e6b commit 98eca27
Show file tree
Hide file tree
Showing 7 changed files with 8 additions and 302 deletions.
37 changes: 0 additions & 37 deletions doc/sphinx/advanced_methods.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand Down
164 changes: 3 additions & 161 deletions src/core/collision.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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<int>(
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<int, 2> 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,
Expand Down Expand Up @@ -414,74 +327,8 @@ std::vector<CollisionPair> 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<CollisionPair> 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
Expand All @@ -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)) {
Expand Down Expand Up @@ -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();
}

Expand Down
11 changes: 1 addition & 10 deletions src/core/collision.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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,
Expand Down
15 changes: 3 additions & 12 deletions src/python/espressomd/collision_detection.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`
Expand Down Expand Up @@ -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
Expand All @@ -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:
Expand Down
7 changes: 0 additions & 7 deletions src/script_interface/collision_detection/CollisionDetection.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ class CollisionDetection : public AutoParameters<CollisionDetection> {
{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<std::string, CollisionModeType> cd_name_to_mode;
std::unordered_map<CollisionModeType,
Expand All @@ -62,9 +61,6 @@ class CollisionDetection : public AutoParameters<CollisionDetection> {
"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:
Expand All @@ -83,9 +79,6 @@ class CollisionDetection : public AutoParameters<CollisionDetection> {

{"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",
Expand Down
Loading

0 comments on commit 98eca27

Please sign in to comment.