diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index a4553a0a13..9ae4054353 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -19,123 +19,168 @@ automatically during the simulation, every time two particles collide. This is useful for simulations of chemical reactions and irreversible adhesion processes. Both, sliding and non-sliding contacts can be created. -The collision detection is controlled via the system -:attr:`~espressomd.system.System.collision_detection` attribute, +The collision detection is controlled via the +:attr:`system.collision_detection +` attribute, which is an instance of the class :class:`~espressomd.collision_detection.CollisionDetection`. -Several modes are available for different types of binding. - -* ``"bind_centers"``: adds a pair-bond between two particles at their first collision. - By making the bonded interaction *stiff* enough, the particles can be held together - after the collision. Note that the particles can still slide on each others' surface, - as the pair bond is not directional. This mode is set up as follows:: - - import espressomd - import espressomd.interactions - - system = espressomd.System(box_l=[1, 1, 1]) - bond_centers = espressomd.interactions.HarmonicBond(k=1000, r_0=1.5) - system.bonded_inter.add(bond_centers) - system.collision_detection.set_params(mode="bind_centers", distance=1.5, - bond_centers=bond_centers) - - The parameters are as follows: - - * ``distance`` is the distance between two particles at which the binding is triggered. - This cutoff distance, ``1.5`` in the example above, is typically chosen slightly larger - than the particle diameter. It is also a good choice for the equilibrium length of the bond. - * ``bond_centers`` is the bonded interaction to be created between the particles - (an instance of :class:`~espressomd.interactions.HarmonicBond` in the example above). - No guarantees are made regarding which of the two colliding particles gets the bond. - Once there is a bond of this type on any of the colliding particles, - no further binding occurs for this pair of particles. - -* ``"bind_at_point_of_collision"``: this mode prevents sliding of the colliding particles at the contact. - This is achieved by creating two virtual sites at the point of collision. - They are rigidly connected to each of the colliding particles. - A bond is then created between the virtual sites, or an angular bond between - the two colliding particles and the virtual particles. In the latter case, - the virtual particles are the centers of the angle potentials - (particle 2 in the description of the angle potential, see :ref:`Bond-angle interactions`). - Due to the rigid connection between each of the - particles in the collision and its respective virtual site, a sliding - at the contact point is no longer possible. See the documentation on - :ref:`Rigid arrangements of particles` for details. In addition to the bond between the virtual - sites, the bond between the colliding particles is also created, i.e., - the ``"bind_at_point_of_collision"`` mode implicitly includes the ``"bind_centers"`` mode. - You can either use a real bonded interaction to prevent wobbling around - the point of contact or you can use :class:`espressomd.interactions.Virtual` which acts as a marker, only. - The method is setup as follows:: - - system.collision_detection.set_params( - mode="bind_at_point_of_collision", - distance=0.1, - bond_centers=harmonic_bond1, - bond_vs=harmonic_bond2, - part_type_vs=1, - vs_placement=0) - - The parameters ``distance`` and ``bond_centers`` have the same meaning - as in the ``"bind_centers"`` mode. The remaining parameters are as follows: - - * ``bond_vs`` is the bond to be added between the two virtual sites created on collision. - This is either a pair-bond with an equilibrium length matching the distance between - the virtual sites, or an angle bond fully stretched in its equilibrium configuration. - * ``part_type_vs`` is the particle type assigned to the virtual sites created on collision. - In nearly all cases, no non-bonded interactions should be defined for this particle type. - * ``vs_placement`` controls, where on the line connecting the centers of the colliding - particles, the virtual sites are placed. A value of 0 means that the virtual sites are - placed at the same position as the colliding particles on which they are based. - A value of 0.5 will result in the virtual sites being placed at the mid-point between - the two colliding particles. A value of 1 will result the virtual site associated - to the first colliding particle to be placed at the position of the second colliding - particle. In most cases, 0.5, is a good choice. Then, the bond connecting the virtual - sites should have an equilibrium length of zero. - -* ``"glue_to_surface"``: This mode is used to irreversibly attach small particles - to the surface of a big particle. It is asymmetric in that several small particles - can be bound to a big particle but not vice versa. The small particles can change type - after collision to make them *inert*. On collision, a single virtual site is placed - and related to the big particle. Then, a bond (``bond_centers``) connects the big - and the small particle. A second bond (``bond_vs``) connects the virtual site and - the small particle. Further required parameters are: - - * ``part_type_to_attach_vs_to``: Type of the particle to which the virtual site is attached, i.e., the *big* particle. - * ``part_type_to_be_glued``: Type of the particle bound to the virtual site (the *small* particle). - * ``part_type_after_glueing``: The type assigned to the particle bound to the virtual site (*small* particle) after the collision. - * ``part_type_vs``: Particle type assigned to the virtual site created during the collision. - * ``distance_glued_particle_to_vs``: Distance of the virtual site to the particle being bound to it (*small* particle). - - Note: When the type of a particle is changed on collision, this makes the - particle inert with regards to further collision. Should a particle of - type ``part_type_to_be_glued`` collide with two particles in a single - time step, no guarantees are made with regards to which partner is selected. - In particular, there is no guarantee that the choice is unbiased. - - The method is used as follows:: - - system.collision_detection.set_params( - mode="glue_to_surface", - distance=0.1, - distance_glued_particle_to_vs=0.02, - bond_centers=harmonic_bond1, - bond_vs=harmonic_bond2, - part_type_vs=1, - part_type_to_attach_vs_to=2, - part_type_to_be_glued=3, - part_type_after_glueing=4) - - -The following limitations currently apply for the collision detection: - -* No distinction is currently made between different particle types for the ``"bind_centers"`` method. - -* The ``"bind_at_point_of_collision"`` and ``"glue_to_surface"`` approaches require - the feature ``VIRTUAL_SITES_RELATIVE`` to be activated in :file:`myconfig.hpp`. - -* The ``"bind_at_point_of_collision"`` approach cannot handle collisions - between virtual sites +Several protocols are available for different types of dynamic binding. +The currently active collision mode can be removed by assigning ``None`` +or :class:`~espressomd.collision_detection.Off` to the +:attr:`system.collision_detection.protocol +` attribute. + +.. _Bind centers: + +Bind centers +~~~~~~~~~~~~ + +Add a pair-bond between two particles at their first collision. +By making the bonded interaction *stiff* enough, the particles can be held together +after the collision. Note that the particles can still slide on each others' surface, +as the pair bond is not directional. This protocol affects all particle types. +This protocol is set up with :class:`~espressomd.collision_detection.BindCenters` as follows:: + + import espressomd + import espressomd.interactions + import espressomd.collision_detection + system = espressomd.System(box_l=[1, 1, 1]) + bond_centers = espressomd.interactions.HarmonicBond(k=1000, r_0=0.1) + system.bonded_inter.add(bond_centers) + system.collision_detection.protocol = espressomd.collision_detection.BindCenters( + distance=0.1, bond_centers=bond_centers) + +The parameters are as follows: + +* ``distance`` is the distance between two particles at which the binding is triggered. + This cutoff distance, ``0.1`` in the example above, is typically chosen slightly larger + than the particle diameter. It is also a good choice for the equilibrium length of the bond. +* ``bond_centers`` is the bonded interaction to be created between the particles + (an instance of :class:`~espressomd.interactions.HarmonicBond` in the example above). + No guarantees are made regarding which of the two colliding particles gets the bond. + Once there is a bond of this type on any of the colliding particles, + no further binding occurs for this pair of particles. + +.. note:: + + The following features are required: + ``COLLISION_DETECTION``. + +.. _Bind at point of collision: + +Bind at point of collision +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Add two pair-bonds between two particles and two automatically generated virtual sites at their first collision. + +This protocol prevents sliding of the colliding particles at the contact point. +This is achieved by creating two virtual sites at the point of collision. +They are rigidly connected to each of the colliding particles. +Then, either a pair bond is added between the virtual sites, or an angular bond +is added between the two colliding particles and the virtual particles. +In the latter case, the virtual particles are the centers of the angle potentials +(particle 2 in the description of the angle potential, see :ref:`Bond-angle interactions`). +Due to the rigid connection between each of the colliding particles and their +respective virtual sites, sliding at the contact point is no longer possible. +See :ref:`Rigid arrangements of particles` for details. This protocol affects all particle types. + +In addition to the bond between the virtual sites, a bond between the colliding +particles is also created. You can either use a real bonded interaction to prevent wobbling +around the point of contact or you can use :class:`espressomd.interactions.Virtual` +which acts as a marker, only. + +This protocol is set up with :class:`~espressomd.collision_detection.BindAtPointOfCollision` as follows:: + + import espressomd + import espressomd.interactions + import espressomd.collision_detection + system = espressomd.System(box_l=[1, 1, 1]) + bond_centers = espressomd.interactions.HarmonicBond(k=1000, r_0=0.1) + bond_vs = espressomd.interactions.HarmonicBond(k=10000, r_0=0.02) + system.bonded_inter.add(bond_centers) + system.bonded_inter.add(bond_vs) + system.collision_detection.protocol = espressomd.collision_detection.BindAtPointOfCollision( + distance=0.1, + bond_centers=bond_centers, + bond_vs=bond_vs, + part_type_vs=1, + vs_placement=0.5) + +The parameters ``distance`` and ``bond_centers`` have the same meaning +as in the :ref:`Bind centers` protocol. The remaining parameters are as follows: + +* ``bond_vs`` is the bond to be added between the two virtual sites created on collision. + This is either a pair-bond with an equilibrium length matching the distance between + the virtual sites, or an angle bond fully stretched in its equilibrium configuration. +* ``part_type_vs`` is the particle type assigned to the virtual sites created on collision. + In nearly all cases, no non-bonded interactions should be defined for this particle type. +* ``vs_placement`` controls where the virtual sites are placed on the line connecting + the colliding particles. A value of 0 means that the virtual sites are + placed at the same position as the colliding particles on which they are based. + A value of 0.5 will result in the virtual sites being placed at the mid-point between + the two colliding particles. A value of 1 will result the virtual site associated + to the first colliding particle to be placed at the position of the second colliding + particle. In most cases, 0.5, is a good choice. Then, the bond connecting the virtual + sites should have an equilibrium length of zero. + +.. note:: + + The following features are required: + ``COLLISION_DETECTION``, ``VIRTUAL_SITES_RELATIVE``. + +.. _Glue to surface: + +Glue to surface +~~~~~~~~~~~~~~~ + +Attach small particles to the surface of a large particle. +The bond can be made irreversible. + +Several small particles can be bound to a large particle but not vice versa. +The small particles can change type after collision to become *inert*. + +This protocol is set up with :class:`~espressomd.collision_detection.GlueToSurface` as follows:: + + import espressomd + import espressomd.interactions + import espressomd.collision_detection + system = espressomd.System(box_l=[1, 1, 1]) + bond_centers = espressomd.interactions.HarmonicBond(k=1000, r_0=0.1) + bond_vs = espressomd.interactions.HarmonicBond(k=10000, r_0=0.02) + system.bonded_inter.add(bond_centers) + system.bonded_inter.add(bond_vs) + system.collision_detection.protocol = espressomd.collision_detection.GlueToSurface( + distance=0.1, + distance_glued_particle_to_vs=0.02, + bond_centers=bond_centers, + bond_vs=bond_vs, + part_type_vs=1, + part_type_to_attach_vs_to=2, + part_type_to_be_glued=3, + part_type_after_glueing=4) + +On collision, a single virtual site is placed and related to the large particle. +Then a bond (``bond_centers``) connects the large and the small particle. +A second bond (``bond_vs``) connects the virtual site and the small particle. +Further required parameters are: + +* ``part_type_to_attach_vs_to``: Type of the particle to which the virtual site is attached, i.e., the *large* particle. +* ``part_type_to_be_glued``: Type of the particle bound to the virtual site (the *small* particle). +* ``part_type_after_glueing``: The type assigned to the particle bound to the virtual site (*small* particle) after the collision. +* ``part_type_vs``: Particle type assigned to the virtual site created during the collision. +* ``distance_glued_particle_to_vs``: Distance of the virtual site to the particle being bound to it (*small* particle), as a fraction of the pair distance. + +Note: When the type of a particle is changed on collision, this makes the +particle inert with regards to further collisions. Should a particle of +type ``part_type_to_be_glued`` collide with two particles in a single +time step, no guarantees are made with regards to which partner is selected. +In particular, there is no guarantee that the choice is unbiased. + +.. note:: + + The following features are required: + ``COLLISION_DETECTION``, ``VIRTUAL_SITES_RELATIVE``. .. _Deleting bonds when particles are pulled apart: @@ -209,10 +254,10 @@ features can be combined to model reversible bonds. Two combinations are possible: * ``"delete_bond"`` mode for breakable bonds together with - ``"bind_centers"`` mode for collision detection: + the :ref:`Bind centers` protocol of collision detection: used to create or delete a bond between two real particles * ``"revert_bind_at_point_of_collision"`` mode for breakable bonds together - with ``"bind_at_point_of_collision"`` mode for collision detection: + with the :ref:`Bind at point of collision` protocol of collision detection: used to create or delete virtual sites (the implicitly created bond between the real particles isn't affected) diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index ce44167449..d97b70a983 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -23,7 +23,6 @@ add_library( espresso_core SHARED bond_error.cpp cells.cpp - collision.cpp communication.cpp dpd.cpp energy.cpp @@ -102,6 +101,7 @@ add_subdirectory(bond_breakage) add_subdirectory(bonded_interactions) add_subdirectory(cell_system) add_subdirectory(cluster_analysis) +add_subdirectory(collision_detection) add_subdirectory(constraints) add_subdirectory(ek) add_subdirectory(electrostatics) diff --git a/src/core/collision.cpp b/src/core/collision.cpp deleted file mode 100644 index 0070d90100..0000000000 --- a/src/core/collision.cpp +++ /dev/null @@ -1,465 +0,0 @@ -/* - * Copyright (C) 2011-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 . - */ - -#include "config/config.hpp" - -#ifdef COLLISION_DETECTION - -#include "collision.hpp" - -#include "BoxGeometry.hpp" -#include "Particle.hpp" -#include "bonded_interactions/bonded_interaction_data.hpp" -#include "cell_system/Cell.hpp" -#include "cell_system/CellStructure.hpp" -#include "communication.hpp" -#include "errorhandling.hpp" -#include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "system/System.hpp" -#include "virtual_sites.hpp" - -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace boost { -namespace serialization { -template -void serialize(Archive &ar, CollisionPair &c, const unsigned int) { - ar & c.pp1; - ar & c.pp2; -} -} // namespace serialization -} // namespace boost - -namespace { -Particle &get_part(CellStructure &cell_structure, int id) { - auto const p = cell_structure.get_local_particle(id); - - if (not p) { - throw std::runtime_error("Could not handle collision because particle " + - std::to_string(id) + " was not found."); - } - - return *p; -} -} // namespace - -void CollisionDetection::initialize() { - // If mode is OFF, no further checks - if (mode == CollisionModeType::OFF) { - return; - } - // Validate distance - if (mode != CollisionModeType::OFF) { - if (distance <= 0.) { - throw std::domain_error("Parameter 'distance' must be > 0"); - } - - // Cache square of cutoff - distance_sq = Utils::sqr(distance); - } - -#ifndef VIRTUAL_SITES_RELATIVE - // The collision modes involving virtual sites also require the creation of a - // bond between the colliding - // If we don't have virtual sites, virtual site binding isn't possible. - if ((mode == CollisionModeType::BIND_VS) or - (mode == CollisionModeType::GLUE_TO_SURF)) { - throw std::runtime_error("collision modes based on virtual sites require " - "the VIRTUAL_SITES_RELATIVE feature"); - } -#endif - -#ifdef VIRTUAL_SITES - // Check vs placement parameter - if (mode == CollisionModeType::BIND_VS) { - if (vs_placement < 0. or vs_placement > 1.) { - throw std::domain_error( - "Parameter 'vs_placement' must be between 0 and 1"); - } - } -#endif - - auto &system = get_system(); - auto &bonded_ias = *system.bonded_ias; - auto &nonbonded_ias = *system.nonbonded_ias; - - // Check if bond exists - assert(mode != CollisionModeType::BIND_CENTERS or - bonded_ias.contains(bond_centers)); - assert(mode != CollisionModeType::BIND_VS or bonded_ias.contains(bond_vs)); - - // If the bond type to bind particle centers is not a pair bond... - // Check that the bonds have the right number of partners - if ((mode == CollisionModeType::BIND_CENTERS) and - (number_of_partners(*bonded_ias.at(bond_centers)) != 1)) { - throw std::runtime_error("The bond type to be used for binding particle " - "centers needs to be a pair bond"); - } - - // The bond between the virtual sites can be pair or triple - if ((mode == CollisionModeType::BIND_VS) and - (number_of_partners(*bonded_ias.at(bond_vs)) != 1 and - number_of_partners(*bonded_ias.at(bond_vs)) != 2)) { - throw std::runtime_error("The bond type to be used for binding virtual " - "sites needs to be a pair bond"); - } - - // Create particle types - - if (mode == CollisionModeType::BIND_VS) { - if (vs_particle_type < 0) { - throw std::domain_error("Collision detection particle type for virtual " - "sites needs to be >=0"); - } - nonbonded_ias.make_particle_type_exist(vs_particle_type); - } - - if (mode == CollisionModeType::GLUE_TO_SURF) { - if (vs_particle_type < 0) { - throw std::domain_error("Collision detection particle type for virtual " - "sites needs to be >=0"); - } - nonbonded_ias.make_particle_type_exist(vs_particle_type); - - if (part_type_to_be_glued < 0) { - throw std::domain_error("Collision detection particle type to be glued " - "needs to be >=0"); - } - nonbonded_ias.make_particle_type_exist(part_type_to_be_glued); - - if (part_type_to_attach_vs_to < 0) { - throw std::domain_error("Collision detection particle type to attach " - "the virtual site to needs to be >=0"); - } - nonbonded_ias.make_particle_type_exist(part_type_to_attach_vs_to); - - if (part_type_after_glueing < 0) { - throw std::domain_error("Collision detection particle type after gluing " - "needs to be >=0"); - } - nonbonded_ias.make_particle_type_exist(part_type_after_glueing); - } - - system.on_short_range_ia_change(); -} - -/** @brief Calculate position of vs for GLUE_TO_SURFACE mode. - * Returns id of particle to bind vs to. - */ -static auto const &glue_to_surface_calc_vs_pos( - Particle const &p1, Particle const &p2, BoxGeometry const &box_geo, - CollisionDetection const &collision_params, Utils::Vector3d &pos) { - double ratio; - auto const vec21 = box_geo.get_mi_vector(p1.pos(), p2.pos()); - auto const dist = vec21.norm(); - - // Find out, which is the particle to be glued. - if ((p1.type() == collision_params.part_type_to_be_glued) and - (p2.type() == collision_params.part_type_to_attach_vs_to)) { - ratio = 1. - collision_params.dist_glued_part_to_vs / dist; - } else if ((p2.type() == collision_params.part_type_to_be_glued) and - (p1.type() == collision_params.part_type_to_attach_vs_to)) { - ratio = collision_params.dist_glued_part_to_vs / dist; - } else { - throw std::runtime_error("This should never be thrown. Bug."); - } - pos = p2.pos() + vec21 * ratio; - if (p1.type() == collision_params.part_type_to_attach_vs_to) - return p1; - - return p2; -} - -static void bind_at_point_of_collision_calc_vs_pos( - Particle const &p1, Particle const &p2, BoxGeometry const &box_geo, - CollisionDetection const &collision_params, Utils::Vector3d &pos1, - Utils::Vector3d &pos2) { - auto const vec21 = box_geo.get_mi_vector(p1.pos(), p2.pos()); - pos1 = p1.pos() - vec21 * collision_params.vs_placement; - pos2 = p1.pos() - vec21 * (1. - collision_params.vs_placement); -} - -#ifdef VIRTUAL_SITES_RELATIVE -static void place_vs_and_relate_to_particle( - CellStructure &cell_structure, BoxGeometry const &box_geo, - CollisionDetection const &collision_params, double const min_global_cut, - int const current_vs_pid, Utils::Vector3d const &pos, int const relate_to) { - Particle new_part; - new_part.id() = current_vs_pid; - new_part.pos() = pos; - auto p_vs = cell_structure.add_particle(std::move(new_part)); - vs_relate_to(*p_vs, get_part(cell_structure, relate_to), box_geo, - min_global_cut); - p_vs->type() = collision_params.vs_particle_type; -} - -static void bind_at_poc_create_bond_between_vs( - CellStructure &cell_structure, BondedInteractionsMap const &bonded_ias, - CollisionDetection const &collision_params, int const current_vs_pid, - CollisionPair const &c) { - switch (number_of_partners(*bonded_ias.at(collision_params.bond_vs))) { - case 1: { - // Create bond between the virtual particles - const int bondG[] = {current_vs_pid - 2}; - // Only add bond if vs was created on this node - if (auto p = cell_structure.get_local_particle(current_vs_pid - 1)) - p->bonds().insert({collision_params.bond_vs, bondG}); - break; - } - case 2: { - // Create 1st bond between the virtual particles - const int bondG[] = {c.pp1, c.pp2}; - // Only add bond if vs was created on this node - if (auto p = cell_structure.get_local_particle(current_vs_pid - 1)) - p->bonds().insert({collision_params.bond_vs, bondG}); - if (auto p = cell_structure.get_local_particle(current_vs_pid - 2)) - p->bonds().insert({collision_params.bond_vs, bondG}); - break; - } - } -} - -static void glue_to_surface_bind_part_to_vs( - Particle const *const p1, Particle const *const p2, - int const vs_pid_plus_one, CollisionDetection const &collision_params, - CellStructure &cell_structure) { - // Create bond between the virtual particles - const int bondG[] = {vs_pid_plus_one - 1}; - - if (p1->type() == collision_params.part_type_after_glueing) { - get_part(cell_structure, p1->id()) - .bonds() - .insert({collision_params.bond_vs, bondG}); - } else { - get_part(cell_structure, p2->id()) - .bonds() - .insert({collision_params.bond_vs, bondG}); - } -} - -#endif // VIRTUAL_SITES_RELATIVE - -static auto gather_collision_queue(std::vector const &local) { - auto global = local; - if (comm_cart.size() > 1) { - Utils::Mpi::gather_buffer(global, ::comm_cart); - boost::mpi::broadcast(::comm_cart, global, 0); - } - return global; -} - -void CollisionDetection::handle_collisions(CellStructure &cell_structure) { - // 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 - // time step - auto const bind_centers = mode != CollisionModeType::OFF and - mode != CollisionModeType::GLUE_TO_SURF; - if (bind_centers) { - for (auto &c : local_collision_queue) { - // put the bond to the non-ghost particle; at least one partner always is - if (cell_structure.get_local_particle(c.pp1)->is_ghost()) { - std::swap(c.pp1, c.pp2); - } - - const int bondG[] = {c.pp2}; - - get_part(cell_structure, c.pp1).bonds().insert({bond_centers, bondG}); - } - } - -#ifdef VIRTUAL_SITES_RELATIVE - auto &system = get_system(); - auto const &box_geo = *system.box_geo; - auto const min_global_cut = system.get_min_global_cut(); - if ((mode == CollisionModeType::BIND_VS) || - (mode == CollisionModeType::GLUE_TO_SURF)) { - // Gather the global collision queue, because only one node has a collision - // across node boundaries in its queue. - // The other node might still have to change particle properties on its - // non-ghost particle - auto global_collision_queue = gather_collision_queue(local_collision_queue); - - // Sync max_seen_part - auto const global_max_seen_particle = boost::mpi::all_reduce( - ::comm_cart, cell_structure.get_max_local_particle_id(), - boost::mpi::maximum()); - - int current_vs_pid = global_max_seen_particle + 1; - - // Iterate over global collision queue - for (auto &c : global_collision_queue) { - - // Get particle pointers - Particle *p1 = cell_structure.get_local_particle(c.pp1); - Particle *p2 = cell_structure.get_local_particle(c.pp2); - - // Only nodes take part in particle creation and binding - // that see both particles - - // If we cannot access both particles, both are ghosts, - // or one is ghost and one is not accessible - // we only increase the counter for the ext id to use based on the - // number of particles created by other nodes - if (((!p1 or p1->is_ghost()) and (!p2 or p2->is_ghost())) or !p1 or !p2) { - // Increase local counters - if (mode == CollisionModeType::BIND_VS) { - current_vs_pid++; - } - // For glue to surface, we have only one vs - current_vs_pid++; - if (mode == CollisionModeType::GLUE_TO_SURF) { - if (p1 and p1->type() == part_type_to_be_glued) { - p1->type() = part_type_after_glueing; - } - if (p2 and p2->type() == part_type_to_be_glued) { - p2->type() = part_type_after_glueing; - } - } // mode glue to surface - - } else { // We consider the pair because one particle - // is local to the node and the other is local or ghost - // If we are in the two vs mode - // Virtual site related to first particle in the collision - if (mode == CollisionModeType::BIND_VS) { - Utils::Vector3d pos1, pos2; - - // Enable rotation on the particles to which vs will be attached - p1->set_can_rotate_all_axes(); - p2->set_can_rotate_all_axes(); - - // Positions of the virtual sites - bind_at_point_of_collision_calc_vs_pos(*p1, *p2, box_geo, *this, pos1, - pos2); - - auto handle_particle = [&](Particle *p, Utils::Vector3d const &pos) { - if (not p->is_ghost()) { - place_vs_and_relate_to_particle(cell_structure, box_geo, *this, - min_global_cut, current_vs_pid, - pos, p->id()); - // Particle storage locations may have changed due to - // added particle - p1 = cell_structure.get_local_particle(c.pp1); - p2 = cell_structure.get_local_particle(c.pp2); - } - }; - - // place virtual sites on the node where the base particle is not a - // ghost - handle_particle(p1, pos1); - // Increment counter - current_vs_pid++; - - handle_particle(p2, pos2); - // Increment counter - current_vs_pid++; - // Create bonds between the vs. - - bind_at_poc_create_bond_between_vs(cell_structure, *system.bonded_ias, - *this, current_vs_pid, c); - } // mode VS - - if (mode == CollisionModeType::GLUE_TO_SURF) { - // If particles are made inert by a type change on collision: - // We skip the pair if one of the particles has already reacted - // but we still increase the particle counters, as other nodes - // can not always know whether or not a vs is placed - if (part_type_after_glueing != part_type_to_be_glued) { - if ((p1->type() == part_type_after_glueing) || - (p2->type() == part_type_after_glueing)) { - current_vs_pid++; - continue; - } - } - - Utils::Vector3d pos; - Particle const &attach_vs_to = - glue_to_surface_calc_vs_pos(*p1, *p2, box_geo, *this, pos); - - // Add a bond between the centers of the colliding particles - // The bond is placed on the node that has p1 - if (!p1->is_ghost()) { - const int bondG[] = {c.pp2}; - get_part(cell_structure, c.pp1) - .bonds() - .insert({bond_centers, bondG}); - } - - // Change type of particle being attached, to make it inert - if (p1->type() == part_type_to_be_glued) { - p1->type() = part_type_after_glueing; - } - if (p2->type() == part_type_to_be_glued) { - p2->type() = part_type_after_glueing; - } - - // Vs placement happens on the node that has p1 - if (!attach_vs_to.is_ghost()) { - place_vs_and_relate_to_particle(cell_structure, box_geo, *this, - min_global_cut, current_vs_pid, pos, - attach_vs_to.id()); - // Particle storage locations may have changed due to - // added particle - p1 = cell_structure.get_local_particle(c.pp1); - p2 = cell_structure.get_local_particle(c.pp2); - current_vs_pid++; - } else { // Just update the books - current_vs_pid++; - } - glue_to_surface_bind_part_to_vs(p1, p2, current_vs_pid, *this, - cell_structure); - } - } // we considered the pair - } // Loop over all collisions in the queue -#ifdef ADDITIONAL_CHECKS - if (!Utils::Mpi::all_compare(comm_cart, current_vs_pid)) { - throw std::runtime_error("Nodes disagree about current_vs_pid"); - } -#endif - - // If any node had a collision, all nodes need to resort - if (not global_collision_queue.empty()) { - cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); - cell_structure.update_ghosts_and_resort_particle( - Cells::DATA_PART_PROPERTIES | Cells::DATA_PART_BONDS); - } - system.update_used_propagations(); - } // are we in one of the vs_based methods -#endif // defined VIRTUAL_SITES_RELATIVE - - clear_queue(); -} - -#endif // COLLISION_DETECTION diff --git a/src/core/collision.hpp b/src/core/collision.hpp deleted file mode 100644 index 0a5545c16e..0000000000 --- a/src/core/collision.hpp +++ /dev/null @@ -1,151 +0,0 @@ -/* - * Copyright (C) 2011-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 . - */ - -#pragma once - -#include "config/config.hpp" - -#include "BondList.hpp" -#include "Particle.hpp" -#include "cell_system/CellStructure.hpp" -#include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "system/Leaf.hpp" - -/** @brief Protocols for collision handling. */ -enum class CollisionModeType : int { - /** @brief Deactivate collision detection. */ - OFF = 0, - /** @brief Create bond between centers of colliding particles. */ - BIND_CENTERS = 1, - /** - * @brief Create a bond between the centers of the colliding particles, - * plus two virtual sites at the point of collision and bind them - * together. This prevents the particles from sliding against each - * other. Requires VIRTUAL_SITES_RELATIVE. - */ - BIND_VS = 2, - /** @brief Glue a particle to a specific spot on another particle. */ - GLUE_TO_SURF = 3, -}; - -/// Data type holding the info about a single collision -struct CollisionPair { - int pp1; // 1st particle id - int pp2; // 2nd particle id -}; - -class CollisionDetection : public System::Leaf { -public: - CollisionDetection() - : mode(CollisionModeType::OFF), distance(0.), distance_sq(0.), - bond_centers(-1), bond_vs(-1) {} - - /// collision protocol - CollisionModeType mode; - /// distance at which particles are bound - double distance; - // Square of distance at which particle are bound - double distance_sq; - - /// bond type used between centers of colliding particles - int bond_centers; - /// bond type used between virtual sites - int bond_vs; - /// particle type for virtual sites created on collision - int vs_particle_type; - - /// For mode "glue to surface": The distance from the particle which is to be - /// glued to the new virtual site - double dist_glued_part_to_vs; - /// For mode "glue to surface": The particle type being glued - int part_type_to_be_glued; - /// For mode "glue to surface": The particle type to which the virtual site is - /// attached - int part_type_to_attach_vs_to; - /// Particle type to which the newly glued particle is converted - int part_type_after_glueing; - /** Placement of virtual sites for MODE_VS. - * 0=on same particle as related to, - * 1=on collision partner, - * 0.5=in the middle between - */ - double vs_placement; - - /** @brief Validates parameters and creates particle types if needed. */ - void initialize(); - - auto cutoff() const { - if (mode != CollisionModeType::OFF) { - return distance; - } - return INACTIVE_CUTOFF; - } - - /// Handle queued collisions - void handle_collisions(CellStructure &cell_structure); - - void clear_queue() { local_collision_queue.clear(); } - - /** @brief Detect (and queue) a collision between the given particles. */ - void detect_collision(Particle const &p1, Particle const &p2, - double const dist_sq) { - if (dist_sq > distance_sq) - return; - - // If we are in the glue to surface mode, check that the particles - // are of the right type - if (mode == CollisionModeType::GLUE_TO_SURF) - if (!glue_to_surface_criterion(p1, p2)) - return; - - // Ignore virtual particles - if (p1.is_virtual() or p2.is_virtual()) - return; - - // Check, if there's already a bond between the particles - if (pair_bond_exists_on(p1.bonds(), p2.id(), bond_centers)) - return; - - if (pair_bond_exists_on(p2.bonds(), p1.id(), bond_centers)) - return; - - /* If we're still here, there is no previous bond between the particles, - we have a new collision */ - - // do not create bond between ghost particles - if (p1.is_ghost() and p2.is_ghost()) { - return; - } - local_collision_queue.push_back({p1.id(), p2.id()}); - } - - // private: - /// During force calculation, colliding particles are recorded in the queue. - /// The queue is processed after force calculation, when it is safe to add - /// particles. - std::vector local_collision_queue; - - /** @brief Check additional criteria for the glue_to_surface collision mode */ - bool glue_to_surface_criterion(Particle const &p1, Particle const &p2) const { - return (((p1.type() == part_type_to_be_glued) and - (p2.type() == part_type_to_attach_vs_to)) or - ((p2.type() == part_type_to_be_glued) and - (p1.type() == part_type_to_attach_vs_to))); - } -}; diff --git a/src/core/collision_detection/ActiveProtocol.hpp b/src/core/collision_detection/ActiveProtocol.hpp new file mode 100644 index 0000000000..6d23c74a76 --- /dev/null +++ b/src/core/collision_detection/ActiveProtocol.hpp @@ -0,0 +1,40 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "BindAtPointOfCollision.hpp" +#include "BindCenters.hpp" +#include "GlueToSurface.hpp" +#include "Off.hpp" + +#include + +namespace CollisionDetection { +using ActiveProtocol = std::variant; +} // namespace CollisionDetection +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/BindAtPointOfCollision.cpp b/src/core/collision_detection/BindAtPointOfCollision.cpp new file mode 100644 index 0000000000..da9bea7c55 --- /dev/null +++ b/src/core/collision_detection/BindAtPointOfCollision.cpp @@ -0,0 +1,193 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#include + +#ifdef COLLISION_DETECTION +#ifdef VIRTUAL_SITES_RELATIVE + +#include "BindAtPointOfCollision.hpp" +#include "CollisionPair.hpp" +#include "utils.hpp" + +#include "BoxGeometry.hpp" +#include "Particle.hpp" +#include "bonded_interactions/bonded_interaction_data.hpp" +#include "cell_system/Cell.hpp" +#include "cell_system/CellStructure.hpp" +#include "communication.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "system/System.hpp" +#include "virtual_sites.hpp" + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace CollisionDetection { + +void BindAtPointOfCollision::initialize(System::System &system) { + // Validate distance + if (distance <= 0.) { + throw std::domain_error("Parameter 'distance' must be > 0"); + } + // Cache square of cutoff + distance_sq = Utils::sqr(distance); + + // Check vs placement parameter + if (vs_placement < 0. or vs_placement > 1.) { + throw std::domain_error("Parameter 'vs_placement' must be between 0 and 1"); + } + + // Check if bond exists + assert(system.bonded_ias->contains(bond_vs)); + // The bond between the virtual sites can be pair or triple + auto const n_partners = number_of_partners(*system.bonded_ias->at(bond_vs)); + if (n_partners != 1 and n_partners != 2) { + throw std::runtime_error("The bond type to be used for binding virtual " + "sites needs to be a pair bond"); + } + + // Create particle types + if (part_type_vs < 0) { + throw std::domain_error("Collision detection particle type for virtual " + "sites needs to be >=0"); + } + system.nonbonded_ias->make_particle_type_exist(part_type_vs); +} + +void BindAtPointOfCollision::handle_collisions( + System::System &system, std::vector &local_collision_queue) { + auto &cell_structure = *system.cell_structure; + auto const min_global_cut = system.get_min_global_cut(); + auto const &box_geo = *system.box_geo; + + add_bind_centers(local_collision_queue, cell_structure, bond_centers); + + // Gather the global collision queue, because only one node has a collision + // across node boundaries in its queue. + // The other node might still have to change particle properties on its + // non-ghost particle + auto global_collision_queue = gather_collision_queue(local_collision_queue); + + // Synchornize max_seen_part + auto const global_max_seen_particle = boost::mpi::all_reduce( + ::comm_cart, cell_structure.get_max_local_particle_id(), + boost::mpi::maximum()); + + int current_vs_pid = global_max_seen_particle + 1; + + // Iterate over global collision queue + for (auto &c : global_collision_queue) { + + // Get particle pointers + Particle *p1 = cell_structure.get_local_particle(c.first); + Particle *p2 = cell_structure.get_local_particle(c.second); + + // Only nodes take part in particle creation and binding + // that see both particles + + // If we cannot access both particles, both are ghosts, + // or one is ghost and one is not accessible + // we only increase the counter for the ext id to use based on the + // number of particles created by other nodes + if (((!p1 or p1->is_ghost()) and (!p2 or p2->is_ghost())) or !p1 or !p2) { + // Increase local counters + current_vs_pid += 2; + continue; + } + // We consider the pair because one particle is local to the node and + // the other is local or ghost. + + // Enable rotation on the particles to which vs will be attached + p1->set_can_rotate_all_axes(); + p2->set_can_rotate_all_axes(); + + // Positions of the virtual sites + auto const vec21 = box_geo.get_mi_vector(p1->pos(), p2->pos()); + auto const pos1 = p1->pos() - vec21 * vs_placement; + auto const pos2 = p1->pos() - vec21 * (1. - vs_placement); + + auto handle_particle = [&](Particle *p, Utils::Vector3d const &pos) { + if (not p->is_ghost()) { + place_vs_and_relate_to_particle(cell_structure, box_geo, part_type_vs, + min_global_cut, current_vs_pid, pos, + p->id()); + // Particle storage locations may have changed due to + // added particle + p1 = cell_structure.get_local_particle(c.first); + p2 = cell_structure.get_local_particle(c.second); + } + }; + + // place virtual sites on the node where the base particle is not a ghost + handle_particle(p1, pos1); + // Increment counter + current_vs_pid++; + + handle_particle(p2, pos2); + // Increment counter + current_vs_pid++; + + // Create bonds + auto const n_partners = number_of_partners(*system.bonded_ias->at(bond_vs)); + if (n_partners == 1) { + // Create bond between the virtual particles + const int bondG[] = {current_vs_pid - 2}; + // Only add bond if vs was created on this node + if (auto p = cell_structure.get_local_particle(current_vs_pid - 1)) + p->bonds().insert({bond_vs, bondG}); + } + if (n_partners == 2) { + // Create 1st bond between the virtual particles + const int bondG[] = {c.first, c.second}; + // Only add bond if vs was created on this node + if (auto p = cell_structure.get_local_particle(current_vs_pid - 1)) + p->bonds().insert({bond_vs, bondG}); + if (auto p = cell_structure.get_local_particle(current_vs_pid - 2)) + p->bonds().insert({bond_vs, bondG}); + } + } // Loop over all collisions in the queue + +#ifdef ADDITIONAL_CHECKS + assert(Utils::Mpi::all_compare(::comm_cart, current_vs_pid) && + "Nodes disagree about current_vs_pid"); +#endif + + // If any node had a collision, all nodes need to resort + if (not global_collision_queue.empty()) { + cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); + cell_structure.update_ghosts_and_resort_particle( + Cells::DATA_PART_PROPERTIES | Cells::DATA_PART_BONDS); + system.update_used_propagations(); + } +} + +} // namespace CollisionDetection + +#endif // VIRTUAL_SITES_RELATIVE +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/BindAtPointOfCollision.hpp b/src/core/collision_detection/BindAtPointOfCollision.hpp new file mode 100644 index 0000000000..f6c4a8007f --- /dev/null +++ b/src/core/collision_detection/BindAtPointOfCollision.hpp @@ -0,0 +1,81 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION +#ifdef VIRTUAL_SITES_RELATIVE + +#include "CollisionPair.hpp" +#include "common.hpp" + +#include "Particle.hpp" +#include "system/System.hpp" + +#include + +namespace CollisionDetection { + +class BindAtPointOfCollision { +public: + /// Distance at which particle are bound + double distance; + /// Square of distance at which particle are bound + double distance_sq; + /// bond type used between centers of colliding particles + int bond_centers; + /// bond type used between virtual sites + int bond_vs; + /** + * @brief Barycenter of the virtual site. + * 0=on same particle as related to, + * 1=on collision partner, + * 0.5=in the middle between the two particles + */ + double vs_placement; + /// particle type for virtual sites created on collision + int part_type_vs; + + BindAtPointOfCollision(double distance, int bond_centers, int bond_vs, + double vs_placement, int part_type_vs) + : distance{distance}, distance_sq{distance * distance}, + bond_centers{bond_centers}, bond_vs{bond_vs}, + vs_placement{vs_placement}, part_type_vs{part_type_vs} {} + + void initialize(System::System &system); + + auto cutoff() const { return distance; } + + void handle_collisions(System::System &system, + std::vector &local_collision_queue); + + bool detect_collision(Particle const &p1, Particle const &p2, + double const dist_sq) const { + if (dist_sq > distance_sq) + return false; + + return detect_collision_common(p1, p2, bond_centers); + } +}; + +} // namespace CollisionDetection +#endif // VIRTUAL_SITES_RELATIVE +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/BindCenters.cpp b/src/core/collision_detection/BindCenters.cpp new file mode 100644 index 0000000000..d648895e76 --- /dev/null +++ b/src/core/collision_detection/BindCenters.cpp @@ -0,0 +1,66 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#include + +#ifdef COLLISION_DETECTION + +#include "BindCenters.hpp" +#include "CollisionPair.hpp" +#include "utils.hpp" + +#include "Particle.hpp" +#include "bonded_interactions/bonded_interaction_data.hpp" +#include "cell_system/CellStructure.hpp" +#include "system/System.hpp" + +#include + +#include +#include +#include +#include + +namespace CollisionDetection { + +void BindCenters::initialize(System::System &system) { + // Validate distance + if (distance <= 0.) { + throw std::domain_error("Parameter 'distance' must be > 0"); + } + // Cache square of cutoff + distance_sq = Utils::sqr(distance); + // Check if bond exists + assert(system.bonded_ias->contains(bond_centers)); + // If the bond type to bind particle centers is not a pair bond... + // Check that the bonds have the right number of partners + if (number_of_partners(*system.bonded_ias->at(bond_centers)) != 1) { + throw std::runtime_error("The bond type to be used for binding particle " + "centers needs to be a pair bond"); + } +} + +void BindCenters::handle_collisions( + System::System &system, std::vector &local_collision_queue) { + add_bind_centers(local_collision_queue, *system.cell_structure, bond_centers); +} + +} // namespace CollisionDetection + +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/BindCenters.hpp b/src/core/collision_detection/BindCenters.hpp new file mode 100644 index 0000000000..e2f58b2fc9 --- /dev/null +++ b/src/core/collision_detection/BindCenters.hpp @@ -0,0 +1,66 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "CollisionPair.hpp" +#include "common.hpp" + +#include "Particle.hpp" +#include "system/System.hpp" + +#include + +namespace CollisionDetection { + +class BindCenters { +public: + /// Distance at which particle are bound + double distance; + /// Square of distance at which particle are bound + double distance_sq; + /// bond type used between centers of colliding particles + int bond_centers; + + BindCenters(double distance, int bond_centers) + : distance{distance}, distance_sq{distance * distance}, + bond_centers{bond_centers} {} + + void initialize(System::System &system); + + auto cutoff() const { return distance; } + + void handle_collisions(System::System &system, + std::vector &local_collision_queue); + + bool detect_collision(Particle const &p1, Particle const &p2, + double const dist_sq) const { + if (dist_sq > distance_sq) + return false; + + return detect_collision_common(p1, p2, bond_centers); + } +}; + +} // namespace CollisionDetection +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/CMakeLists.txt b/src/core/collision_detection/CMakeLists.txt new file mode 100644 index 0000000000..57f5001766 --- /dev/null +++ b/src/core/collision_detection/CMakeLists.txt @@ -0,0 +1,25 @@ +# +# Copyright (C) 2024 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 . +# + +target_sources( + espresso_core + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/CollisionDetection.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/BindCenters.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/BindAtPointOfCollision.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GlueToSurface.cpp) diff --git a/src/core/collision_detection/CollisionDetection.cpp b/src/core/collision_detection/CollisionDetection.cpp new file mode 100644 index 0000000000..9a1dc15f45 --- /dev/null +++ b/src/core/collision_detection/CollisionDetection.cpp @@ -0,0 +1,69 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#include + +#ifdef COLLISION_DETECTION + +#include "CollisionDetection.hpp" + +#include "system/System.hpp" + +#include +#include +#include + +namespace CollisionDetection { + +void CollisionDetection::handle_collisions() { + auto &system = get_system(); + if (m_protocol) { + std::visit( + [&system, this](auto &protocol) { + protocol.handle_collisions(system, local_collision_queue); + }, + *m_protocol); + } + clear_queue(); +} + +void CollisionDetection::initialize() { + auto &system = get_system(); + if (m_protocol) { + std::visit([&system](auto &protocol) { protocol.initialize(system); }, + *m_protocol); + } + system.on_short_range_ia_change(); +} + +void CollisionDetection::set_protocol( + std::shared_ptr protocol) { + m_protocol = std::move(protocol); + initialize(); +} + +void CollisionDetection::unset_protocol() { + m_protocol = nullptr; + auto &system = get_system(); + system.on_short_range_ia_change(); +} + +} // namespace CollisionDetection + +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/CollisionDetection.hpp b/src/core/collision_detection/CollisionDetection.hpp new file mode 100644 index 0000000000..2683a99cd4 --- /dev/null +++ b/src/core/collision_detection/CollisionDetection.hpp @@ -0,0 +1,93 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "ActiveProtocol.hpp" +#include "CollisionPair.hpp" + +#include "Particle.hpp" +#include "system/Leaf.hpp" + +#include +#include +#include +#include + +namespace CollisionDetection { + +class CollisionDetection : public System::Leaf { + std::shared_ptr m_protocol; + +public: + CollisionDetection() = default; + /** @brief Get currently active collision protocol. */ + auto get_protocol() const { return m_protocol; } + /** @brief Set a new collision protocol. */ + void set_protocol(std::shared_ptr protocol); + /** @brief Delete the currently active collision protocol. */ + void unset_protocol(); + /** @brief Validate parameters and create particle types if needed. */ + void initialize(); + + auto is_off() const { + return m_protocol == nullptr or std::holds_alternative(*m_protocol); + } + + auto cutoff() const { + if (m_protocol == nullptr) { + return INACTIVE_CUTOFF; + } + return std::visit([](auto const &protocol) { return protocol.cutoff(); }, + *m_protocol); + } + + /// Handle queued collisions + void handle_collisions(); + + void clear_queue() { local_collision_queue.clear(); } + + /** @brief Detect (and queue) a collision between the given particles. */ + void detect_collision(Particle const &p1, Particle const &p2, + double const dist_sq) { + if (m_protocol) { + bool collision_detected = std::visit( + [&p1, &p2, dist_sq](auto const &protocol) { + return protocol.detect_collision(p1, p2, dist_sq); + }, + *m_protocol); + + if (collision_detected) { + local_collision_queue.emplace_back(p1.id(), p2.id()); + } + } + } + /// During force calculation, colliding particles are recorded in the queue. + /// The queue is processed after force calculation, when it is safe to add + /// particles. + std::vector local_collision_queue; +}; + +} // namespace CollisionDetection + +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/CollisionPair.hpp b/src/core/collision_detection/CollisionPair.hpp new file mode 100644 index 0000000000..d943336105 --- /dev/null +++ b/src/core/collision_detection/CollisionPair.hpp @@ -0,0 +1,31 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include + +namespace CollisionDetection { +using CollisionPair = std::pair; +} // namespace CollisionDetection +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/GlueToSurface.cpp b/src/core/collision_detection/GlueToSurface.cpp new file mode 100644 index 0000000000..b4cad55c02 --- /dev/null +++ b/src/core/collision_detection/GlueToSurface.cpp @@ -0,0 +1,214 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#include + +#ifdef COLLISION_DETECTION +#ifdef VIRTUAL_SITES_RELATIVE + +#include "CollisionPair.hpp" +#include "GlueToSurface.hpp" +#include "utils.hpp" + +#include "BoxGeometry.hpp" +#include "Particle.hpp" +#include "bonded_interactions/bonded_interaction_data.hpp" +#include "cell_system/Cell.hpp" +#include "cell_system/CellStructure.hpp" +#include "communication.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "system/System.hpp" +#include "virtual_sites.hpp" + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace CollisionDetection { + +void GlueToSurface::initialize(System::System &system) { + // Validate distance + if (distance <= 0.) { + throw std::domain_error("Parameter 'distance' must be > 0"); + } + // Cache square of cutoff + distance_sq = Utils::sqr(distance); + + if (part_type_vs < 0) { + throw std::domain_error("Collision detection particle type for virtual " + "sites needs to be >=0"); + } + system.nonbonded_ias->make_particle_type_exist(part_type_vs); + + if (part_type_to_be_glued < 0) { + throw std::domain_error("Collision detection particle type to be glued " + "needs to be >=0"); + } + system.nonbonded_ias->make_particle_type_exist(part_type_to_be_glued); + + if (part_type_to_attach_vs_to < 0) { + throw std::domain_error("Collision detection particle type to attach " + "the virtual site to needs to be >=0"); + } + system.nonbonded_ias->make_particle_type_exist(part_type_to_attach_vs_to); + + if (part_type_after_glueing < 0) { + throw std::domain_error("Collision detection particle type after gluing " + "needs to be >=0"); + } + system.nonbonded_ias->make_particle_type_exist(part_type_after_glueing); +} + +void GlueToSurface::handle_collisions( + System::System &system, std::vector &local_collision_queue) { + auto &cell_structure = *system.cell_structure; + auto const min_global_cut = system.get_min_global_cut(); + 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 + // time step + + // Gather the global collision queue, because only one node has a collision + // across node boundaries in its queue. + // The other node might still have to change particle properties on its + // non-ghost particle + auto global_collision_queue = gather_collision_queue(local_collision_queue); + + // Sync max_seen_part + auto const global_max_seen_particle = boost::mpi::all_reduce( + ::comm_cart, cell_structure.get_max_local_particle_id(), + boost::mpi::maximum()); + + int current_vs_pid = global_max_seen_particle + 1; + + // Iterate over global collision queue + for (auto &c : global_collision_queue) { + + // Get particle pointers + Particle *p1 = cell_structure.get_local_particle(c.first); + Particle *p2 = cell_structure.get_local_particle(c.second); + + // Only nodes take part in particle creation and binding + // that see both particles + + // If we cannot access both particles, both are ghosts, + // or one is ghost and one is not accessible + // we only increase the counter for the ext id to use based on the + // number of particles created by other nodes + if (((!p1 or p1->is_ghost()) and (!p2 or p2->is_ghost())) or !p1 or !p2) { + // For glue to surface, we have only one vs + current_vs_pid++; + + if (p1 and p1->type() == part_type_to_be_glued) { + p1->type() = part_type_after_glueing; + } + if (p2 and p2->type() == part_type_to_be_glued) { + p2->type() = part_type_after_glueing; + } + continue; + } + // If particles are made inert by a type change on collision: + // We skip the pair if one of the particles has already reacted + // but we still increase the particle counters, as other nodes + // can not always know whether or not a vs is placed + if (part_type_after_glueing != part_type_to_be_glued) { + if ((p1->type() == part_type_after_glueing) or + (p2->type() == part_type_after_glueing)) { + current_vs_pid++; + continue; + } + } + + double ratio = -1.; + auto const vec21 = box_geo.get_mi_vector(p1->pos(), p2->pos()); + auto const dist = vec21.norm(); + + // Find out, which is the particle to be glued. + if ((p1->type() == part_type_to_be_glued) and + (p2->type() == part_type_to_attach_vs_to)) { + ratio = 1. - dist_glued_part_to_vs / dist; + } else if ((p2->type() == part_type_to_be_glued) and + (p1->type() == part_type_to_attach_vs_to)) { + ratio = dist_glued_part_to_vs / dist; + } + assert(ratio != -1.); + auto const pos = p2->pos() + vec21 * ratio; + auto const &attach_vs_to = + (p1->type() == part_type_to_attach_vs_to) ? *p1 : *p2; + + // Add a bond between the centers of the colliding particles + // The bond is placed on the node that has p1 + if (!p1->is_ghost()) { + const int bondG[] = {c.second}; + get_part(cell_structure, c.first).bonds().insert({bond_centers, bondG}); + } + + // Change type of particle being attached, to make it inert + if (p1->type() == part_type_to_be_glued) { + p1->type() = part_type_after_glueing; + } + if (p2->type() == part_type_to_be_glued) { + p2->type() = part_type_after_glueing; + } + + if (attach_vs_to.is_ghost()) { + current_vs_pid++; + } else { + // VS placement happens on the node that has p1 + place_vs_and_relate_to_particle(cell_structure, box_geo, part_type_vs, + min_global_cut, current_vs_pid, pos, + attach_vs_to.id()); + // Particle storage locations may have changed due to added particle + p1 = cell_structure.get_local_particle(c.first); + p2 = cell_structure.get_local_particle(c.second); + current_vs_pid++; + } + // Create bond between the virtual particles + auto const p = (p1->type() == part_type_after_glueing) ? p1 : p2; + int const bondG[] = {current_vs_pid - 1}; + get_part(cell_structure, p->id()).bonds().insert({bond_vs, bondG}); + } // Loop over all collisions in the queue + +#ifdef ADDITIONAL_CHECKS + assert(Utils::Mpi::all_compare(::comm_cart, current_vs_pid) && + "Nodes disagree about current_vs_pid"); +#endif + + // If any node had a collision, all nodes need to resort + if (not global_collision_queue.empty()) { + cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); + cell_structure.update_ghosts_and_resort_particle( + Cells::DATA_PART_PROPERTIES | Cells::DATA_PART_BONDS); + system.update_used_propagations(); + } +} + +} // namespace CollisionDetection + +#endif // VIRTUAL_SITES_RELATIVE +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/GlueToSurface.hpp b/src/core/collision_detection/GlueToSurface.hpp new file mode 100644 index 0000000000..999f15ffaf --- /dev/null +++ b/src/core/collision_detection/GlueToSurface.hpp @@ -0,0 +1,99 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION +#ifdef VIRTUAL_SITES_RELATIVE + +#include "CollisionPair.hpp" +#include "common.hpp" + +#include "Particle.hpp" +#include "system/System.hpp" + +#include + +namespace CollisionDetection { + +class GlueToSurface { +public: + /// Distance at which particle are bound + double distance; + /// Square of distance at which particle are bound + double distance_sq; + /// bond type used between centers of colliding particles + int bond_centers; + /// bond type used between virtual sites + int bond_vs; + /// particle type for virtual sites created on collision + int part_type_vs; + /// Distance from the particle which is to be glued to the new virtual site + double dist_glued_part_to_vs; + /// The particle type being glued (the small particle) + int part_type_to_be_glued; + /// The particle type to which virtuals site are attached (the large particle) + int part_type_to_attach_vs_to; + /// Particle type to which the newly glued particle is converted + int part_type_after_glueing; + + GlueToSurface(double distance, int bond_centers, int bond_vs, + int part_type_vs, double dist_glued_part_to_vs, + int part_type_to_be_glued, int part_type_to_attach_vs_to, + int part_type_after_glueing) + : distance{distance}, distance_sq{distance * distance}, + bond_centers{bond_centers}, bond_vs{bond_vs}, + part_type_vs{part_type_vs}, + dist_glued_part_to_vs{dist_glued_part_to_vs}, + part_type_to_be_glued{part_type_to_be_glued}, + part_type_to_attach_vs_to{part_type_to_attach_vs_to}, + part_type_after_glueing{part_type_after_glueing} {} + + void initialize(System::System &system); + + auto cutoff() const { return distance; } + + /** @brief Check additional criteria for the glue_to_surface collision mode */ + bool glue_to_surface_criterion(Particle const &p1, Particle const &p2) const { + return ((p1.type() == part_type_to_be_glued and + p2.type() == part_type_to_attach_vs_to) or + (p2.type() == part_type_to_be_glued and + p1.type() == part_type_to_attach_vs_to)); + } + + void handle_collisions(System::System &system, + std::vector &local_collision_queue); + + bool detect_collision(Particle const &p1, Particle const &p2, + double const dist_sq) const { + if (dist_sq > distance_sq) + return false; + + if (!glue_to_surface_criterion(p1, p2)) + return false; + + return detect_collision_common(p1, p2, bond_centers); + } +}; + +} // namespace CollisionDetection +#endif // VIRTUAL_SITES_RELATIVE +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/Off.hpp b/src/core/collision_detection/Off.hpp new file mode 100644 index 0000000000..21f6e65c3f --- /dev/null +++ b/src/core/collision_detection/Off.hpp @@ -0,0 +1,53 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "CollisionPair.hpp" + +#include "Particle.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "system/System.hpp" + +#include + +namespace CollisionDetection { + +class Off { +public: + Off() = default; + + auto cutoff() const { return INACTIVE_CUTOFF; } + + void initialize(System::System &) {} + + void handle_collisions(System::System &, std::vector &) {} + + bool detect_collision(Particle const &, Particle const &, + double const) const { + return false; + } +}; + +} // namespace CollisionDetection +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/common.hpp b/src/core/collision_detection/common.hpp new file mode 100644 index 0000000000..921477cb6f --- /dev/null +++ b/src/core/collision_detection/common.hpp @@ -0,0 +1,50 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "BondList.hpp" +#include "Particle.hpp" + +namespace CollisionDetection { + +inline auto detect_collision_common(Particle const &p1, Particle const &p2, + int const bond_centers) { + // Ignore virtual particles + if (p1.is_virtual() or p2.is_virtual()) + return false; + + // Check, if there's already a bond between the particles + if (pair_bond_exists_on(p1.bonds(), p2.id(), bond_centers) or + pair_bond_exists_on(p2.bonds(), p1.id(), bond_centers)) + return false; + + /* If we're still here, there is no previous bond between the particles, + we have a new collision */ + + // do not create bond between ghost particles + return not(p1.is_ghost() and p2.is_ghost()); +} + +} // namespace CollisionDetection +#endif // COLLISION_DETECTION diff --git a/src/core/collision_detection/utils.hpp b/src/core/collision_detection/utils.hpp new file mode 100644 index 0000000000..6c9577a140 --- /dev/null +++ b/src/core/collision_detection/utils.hpp @@ -0,0 +1,99 @@ +/* + * Copyright (C) 2011-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "CollisionPair.hpp" + +#include "BoxGeometry.hpp" +#include "Particle.hpp" +#include "cell_system/CellStructure.hpp" +#include "communication.hpp" +#include "virtual_sites.hpp" + +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace CollisionDetection { + +inline auto &get_part(CellStructure &cell_structure, int id) { + auto const p = cell_structure.get_local_particle(id); + + if (not p) { + throw std::runtime_error("Could not handle collision because particle " + + std::to_string(id) + " was not found."); + } + + return *p; +} + +#ifdef VIRTUAL_SITES_RELATIVE +inline void place_vs_and_relate_to_particle( + CellStructure &cell_structure, BoxGeometry const &box_geo, + int const part_type_vs, double const min_global_cut, + int const current_vs_pid, Utils::Vector3d const &pos, int const relate_to) { + Particle new_part; + new_part.id() = current_vs_pid; + new_part.pos() = pos; + auto p_vs = cell_structure.add_particle(std::move(new_part)); + vs_relate_to(*p_vs, get_part(cell_structure, relate_to), box_geo, + min_global_cut); + p_vs->type() = part_type_vs; +} +#endif // VIRTUAL_SITES_RELATIVE + +inline auto gather_collision_queue(std::vector const &local) { + auto global = local; + if (::comm_cart.size() > 1) { + Utils::Mpi::gather_buffer(global, ::comm_cart); + boost::mpi::broadcast(::comm_cart, global, 0); + } + return global; +} + +inline void add_bind_centers(std::vector &collision_queue, + CellStructure &cell_structure, int bond_centers) { + for (auto &c : collision_queue) { + // Ensure that the bond is associated with the non-ghost particle + if (cell_structure.get_local_particle(c.first)->is_ghost()) { + std::swap(c.first, c.second); + } + + const int bondG[] = {c.second}; + + // Insert the bond for the non-ghost particle + get_part(cell_structure, c.first).bonds().insert({bond_centers, bondG}); + } +} + +} // namespace CollisionDetection + +#endif // COLLISION_DETECTION diff --git a/src/core/forces.cpp b/src/core/forces.cpp index b9881e166c..d5a8c0e8ba 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -31,7 +31,7 @@ #include "bond_breakage/bond_breakage.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" -#include "collision.hpp" +#include "collision_detection/CollisionDetection.hpp" #include "communication.hpp" #include "constraints/Constraints.hpp" #include "electrostatics/icc.hpp" @@ -197,7 +197,7 @@ void System::System::calculate_forces() { coulomb_kernel_ptr, dipoles_kernel_ptr, elc_kernel_ptr); #ifdef COLLISION_DETECTION - if (collision_detection.mode != CollisionModeType::OFF) { + if (not collision_detection.is_off()) { collision_detection.detect_collision(p1, p2, d.dist2); } #endif diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 08b514efb1..34e2abe9f2 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -42,7 +42,7 @@ #include "bonded_interactions/bonded_interaction_data.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" -#include "collision.hpp" +#include "collision_detection/CollisionDetection.hpp" #include "communication.hpp" #include "errorhandling.hpp" #include "forces.hpp" @@ -659,7 +659,7 @@ int System::System::integrate(int n_steps, int reuse_forces) { #endif #ifdef COLLISION_DETECTION - collision_detection->handle_collisions(*cell_structure); + collision_detection->handle_collisions(); #endif bond_breakage->process_queue(*this); } diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp index 676b824155..5243f82fc6 100644 --- a/src/core/system/System.cpp +++ b/src/core/system/System.cpp @@ -31,7 +31,7 @@ #include "cell_system/CellStructure.hpp" #include "cell_system/CellStructureType.hpp" #include "cell_system/HybridDecomposition.hpp" -#include "collision.hpp" +#include "collision_detection/CollisionDetection.hpp" #include "communication.hpp" #include "electrostatics/icc.hpp" #include "errorhandling.hpp" @@ -75,7 +75,8 @@ System::System(Private) { oif_global = std::make_shared(); immersed_boundaries = std::make_shared(); #ifdef COLLISION_DETECTION - collision_detection = std::make_shared(); + collision_detection = + std::make_shared(); #endif bond_breakage = std::make_shared(); lees_edwards = std::make_shared(); @@ -503,7 +504,7 @@ unsigned System::get_global_ghost_flags() const { } #ifdef COLLISION_DETECTION - if (collision_detection->mode != CollisionModeType::OFF) { + if (not collision_detection->is_off()) { data_parts |= Cells::DATA_PART_BONDS; } #endif diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index bc576045c5..55be810194 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -50,7 +50,9 @@ class Galilei; class Observable_stat; class OifGlobal; class ImmersedBoundaries; +namespace CollisionDetection { class CollisionDetection; +} namespace BondBreakage { class BondBreakage; } @@ -284,7 +286,7 @@ class System : public std::enable_shared_from_this { std::shared_ptr oif_global; std::shared_ptr immersed_boundaries; #ifdef COLLISION_DETECTION - std::shared_ptr collision_detection; + std::shared_ptr collision_detection; #endif std::shared_ptr bond_breakage; std::shared_ptr lees_edwards; diff --git a/src/core/system/System.impl.hpp b/src/core/system/System.impl.hpp index 24a3f8c9c0..b3762721a2 100644 --- a/src/core/system/System.impl.hpp +++ b/src/core/system/System.impl.hpp @@ -29,7 +29,7 @@ #include "bond_breakage/bond_breakage.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" #include "cell_system/CellStructure.hpp" -#include "collision.hpp" +#include "collision_detection/CollisionDetection.hpp" #include "constraints/Constraints.hpp" #include "galilei/ComFixed.hpp" #include "galilei/Galilei.hpp" diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 58a9bba30e..829aedbfc6 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -37,6 +37,7 @@ namespace utf = boost::unit_test; #include "bonded_interactions/harmonic.hpp" #include "cell_system/CellStructure.hpp" #include "cell_system/CellStructureType.hpp" +#include "collision_detection/utils.hpp" #include "communication.hpp" #include "cuda/utils.hpp" #include "electrostatics/coulomb.hpp" @@ -559,6 +560,10 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { BOOST_CHECK_THROW(throw BondUnknownTypeError(), std::exception); BOOST_CHECK_THROW(throw BondInvalidSizeError(2), std::exception); BOOST_CHECK_EQUAL(BondInvalidSizeError(2).size, 2); +#ifdef COLLISION_DETECTION + BOOST_CHECK_THROW(CollisionDetection::get_part(*system.cell_structure, 777), + std::runtime_error); +#endif } // check exceptions diff --git a/src/python/espressomd/collision_detection.py b/src/python/espressomd/collision_detection.py index ac90777f50..cce17f48a5 100644 --- a/src/python/espressomd/collision_detection.py +++ b/src/python/espressomd/collision_detection.py @@ -1,5 +1,5 @@ # -# Copyright (C) 2010-2022 The ESPResSo project +# Copyright (C) 2010-2024 The ESPResSo project # # This file is part of ESPResSo. # @@ -32,79 +32,127 @@ class CollisionDetection(ScriptInterfaceHelper): the :attr:`~espressomd.system.System.collision_detection` attribute of the system class to access the collision detection. - Use method :meth:`~espressomd.collision_detection.CollisionDetection.set_params` - to change the parameters of the collision detection. + Attributes + ---------- + protocol : + Protocol instance, or ``None`` if not set. """ _so_name = "CollisionDetection::CollisionDetection" _so_features = ("COLLISION_DETECTION",) - # Do not allow setting of individual attributes - def __setattr__(self, *args, **kwargs): - raise Exception( - "Please set all parameters at once via collision_detection.set_params()") - def set_params(self, **kwargs): - """ - Set the parameters for the collision detection +@script_interface_register +class Off(ScriptInterfaceHelper): - See :ref:`Creating bonds when particles collide` for detailed instructions. + """ + Disable collision detection. + """ + _so_name = "CollisionDetection::Off" + _so_features = ("COLLISION_DETECTION",) - Parameters - ---------- - mode : :obj:`str`, {"off", "bind_centers", "bind_at_point_of_collision", "glue_to_surface"} - Collision detection mode - distance : :obj:`float` - Distance below which a pair of particles is considered in the - collision detection +@script_interface_register +class BindCenters(ScriptInterfaceHelper): - bond_centers : :obj:`espressomd.interactions.BondedInteraction` - Bond to add between the colliding particles + """ + Add a pair bond between two particles upon collision. + Particles can still slide around their contact point. - bond_vs : :obj:`espressomd.interactions.BondedInteraction` - Bond to add between virtual sites (for modes using virtual sites) + Parameters + ---------- + distance : :obj:`float` + Distance below which two particles are considered to have collided. - part_type_vs : :obj:`int` - Particle type of the virtual sites being created on collision - (virtual sites based modes) + bond_centers : :obj:`espressomd.interactions.BondedInteraction` + Bond to add between the colliding particles. - part_type_to_be_glued : :obj:`int` - particle type for ``"glue_to_surface"`` mode. See user guide. + """ + _so_name = "CollisionDetection::BindCenters" + _so_features = ("COLLISION_DETECTION",) - part_type_to_attach_vs_to : :obj:`int` - particle type for ``"glue_to_surface"`` mode. See user guide. - part_type_after_glueing : :obj:`int` - particle type for ``"glue_to_surface"`` mode. See user guide. +@script_interface_register +class BindAtPointOfCollision(ScriptInterfaceHelper): + + """ + Add a pair bond between two particles upon collision, and two extra + pair bonds or angle bonds with two automatically-generated virtual sites. + This protocol prevents sliding of the particles at the contact point. - distance_glued_particle_to_vs : :obj:`float` - Distance for ``"glue_to_surface"`` mode. See user guide. + Parameters + ---------- + distance : :obj:`float` + Distance below which two particles are considered to have collided. - """ + bond_centers : :obj:`espressomd.interactions.BondedInteraction` + Bond to add between the colliding particles. - if "mode" not in kwargs: - raise ValueError( - "Collision mode must be specified via the 'mode' argument") - self.call_method("set_params", **kwargs) + bond_vs : :obj:`espressomd.interactions.BondedInteraction` + Bond to add between virtual sites. - def get_parameter(self, name): - value = super().get_parameter(name) - if name in ["bond_centers", "bond_vs"]: - if value == -1: # Not defined - value = None - else: - value = self.call_method("get_bond_by_id", bond_id=value) - return value + part_type_vs : :obj:`int` + Particle type of the virtual sites created on collision. - def get_params(self): - """Returns the parameters of the collision detection as dict. + vs_placement : :obj:`float` + Barycenter of the virtual sites. A value of 0 means that the virtual sites + are placed at the same position as the colliding particles on which they are based. + A value of 0.5 will result in the virtual sites being placed at the mid-point between + the two colliding particles. A value of 1 will result the virtual site associated to + the first colliding particle to be placed at the position of the second colliding particle. + In most cases, 0.5, is a good choice. Then, the bond connecting the virtual sites + should have an equilibrium length of zero. - """ - params = {} - mode = super().get_parameter("mode") - for name in self.call_method("params_for_mode", mode=mode): - params[name] = self.get_parameter(name) - return params + """ + _so_name = "CollisionDetection::BindAtPointOfCollision" + _so_features = ("COLLISION_DETECTION", "VIRTUAL_SITES_RELATIVE") + + +@script_interface_register +class GlueToSurface(ScriptInterfaceHelper): + + """ + Attach small particles to the surface of a larger particle. + + It is asymmetric: several small particles can be bound to a large particle + but not vice versa. It can be made irreversible: the small particles can + change type after collision to become *inert*. + + On collision, a single virtual site is placed and related to the large + particle. Then a bond (``bond_centers``) connects the large and the small + particle. A second bond (``bond_vs``) connects the virtual site and + the small particle. + + Parameters + ---------- + distance : :obj:`float` + Distance below which two particles are considered to have collided. + + bond_centers : :obj:`espressomd.interactions.BondedInteraction` + Bond to add between the colliding particles. + + bond_vs : :obj:`espressomd.interactions.BondedInteraction` + Bond to add between virtual sites. + + part_type_vs : :obj:`int` + Particle type of the virtual sites created on collision. + + part_type_to_attach_vs_to : :obj:`int` + Type of the large particle. + + part_type_to_be_glued : :obj:`int` + Type of the small particle. + + part_type_after_glueing : :obj:`int` + Type of the small particle after collision. If different from + ``part_type_to_be_glued``, the bond is irreversible. + + distance_glued_particle_to_vs : :obj:`float` + Distance of the virtual site to the small particle, + as a fraction of the pair distance. + + """ + _so_name = "CollisionDetection::GlueToSurface" + _so_features = ("COLLISION_DETECTION", "VIRTUAL_SITES_RELATIVE") diff --git a/src/script_interface/collision_detection/BindAtPointOfCollision.hpp b/src/script_interface/collision_detection/BindAtPointOfCollision.hpp new file mode 100644 index 0000000000..aaf586fb55 --- /dev/null +++ b/src/script_interface/collision_detection/BindAtPointOfCollision.hpp @@ -0,0 +1,92 @@ +/* + * Copyright (C) 2021-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION +#ifdef VIRTUAL_SITES_RELATIVE + +#include "Protocol.hpp" + +#include "core/collision_detection/BindAtPointOfCollision.hpp" + +#include "script_interface/ScriptInterface.hpp" +#include "script_interface/auto_parameters/AutoParameters.hpp" + +#include +#include + +namespace ScriptInterface::CollisionDetection { + +class BindAtPointOfCollision : public Protocol { + using CoreClass = ::CollisionDetection::BindAtPointOfCollision; + +public: + BindAtPointOfCollision() { + add_parameters({ + {"bond_centers", AutoParameter::read_only, + [this]() { + return get_bond_variant_by_id( + std::get(*m_protocol).bond_centers); + }}, + {"distance", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).distance; }}, + {"bond_vs", AutoParameter::read_only, + [this]() { + return get_bond_variant_by_id( + std::get(*m_protocol).bond_vs); + }}, + {"vs_placement", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).vs_placement; }}, + {"part_type_vs", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).part_type_vs; }}, + {"_bond_centers", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).bond_centers; }}, + {"_bond_vs", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).bond_vs; }}, + }); + } + std::shared_ptr<::CollisionDetection::ActiveProtocol> protocol() override { + return m_protocol; + } + +private: + std::shared_ptr<::CollisionDetection::ActiveProtocol> m_protocol; + +protected: + void do_initialize(VariantMap const ¶ms) override { + auto const bond_center = find_bond_id(params.contains("_bond_centers") + ? params.at("_bond_centers") + : params.at("bond_centers")); + auto const bond_vs = + find_bond_id(params.contains("_bond_vs") ? params.at("_bond_vs") + : params.at("bond_vs")); + m_protocol = std::make_shared<::CollisionDetection::ActiveProtocol>( + CoreClass(get_value(params, "distance"), bond_center, bond_vs, + get_value(params, "vs_placement"), + get_value(params, "part_type_vs"))); + } +}; + +} // namespace ScriptInterface::CollisionDetection + +#endif // VIRTUAL_SITES_RELATIVE +#endif // COLLISION_DETECTION diff --git a/src/script_interface/collision_detection/BindCenters.hpp b/src/script_interface/collision_detection/BindCenters.hpp new file mode 100644 index 0000000000..87926d13ee --- /dev/null +++ b/src/script_interface/collision_detection/BindCenters.hpp @@ -0,0 +1,73 @@ +/* + * Copyright (C) 2021-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "Protocol.hpp" + +#include "core/collision_detection/BindCenters.hpp" + +#include "script_interface/ScriptInterface.hpp" +#include "script_interface/auto_parameters/AutoParameters.hpp" + +#include +#include + +namespace ScriptInterface::CollisionDetection { + +class BindCenters : public Protocol { + using CoreClass = ::CollisionDetection::BindCenters; + +public: + BindCenters() { + add_parameters( + {{"bond_centers", AutoParameter::read_only, + [this]() { + return get_bond_variant_by_id( + std::get(*m_protocol).bond_centers); + }}, + {"_bond_centers", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).bond_centers; }}, + {"distance", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).distance; }}}); + } + std::shared_ptr<::CollisionDetection::ActiveProtocol> protocol() override { + return m_protocol; + } + +private: + std::shared_ptr<::CollisionDetection::ActiveProtocol> m_protocol; + +protected: + void do_initialize(VariantMap const ¶ms) override { + auto const bond_center = find_bond_id(params.contains("_bond_centers") + ? params.at("_bond_centers") + : params.at("bond_centers")); + m_protocol = std::make_shared<::CollisionDetection::ActiveProtocol>( + CoreClass(get_value(params, "distance"), bond_center)); + } +}; + +} // namespace ScriptInterface::CollisionDetection + +#endif // COLLISION_DETECTION diff --git a/src/script_interface/collision_detection/CollisionDetection.hpp b/src/script_interface/collision_detection/CollisionDetection.hpp index bb1727fd0d..7003f3f09c 100644 --- a/src/script_interface/collision_detection/CollisionDetection.hpp +++ b/src/script_interface/collision_detection/CollisionDetection.hpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2010-2022 The ESPResSo project + * Copyright (C) 2010-2024 The ESPResSo project * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 * Max-Planck-Institute for Polymer Research, Theory Group * @@ -21,10 +21,12 @@ #pragma once -#include "config/config.hpp" +#include #ifdef COLLISION_DETECTION +#include "Protocol.hpp" + #include "script_interface/ScriptInterface.hpp" #include "script_interface/interactions/BondedInteraction.hpp" #include "script_interface/interactions/BondedInteractions.hpp" @@ -32,239 +34,74 @@ #include "script_interface/system/System.hpp" #include "core/bonded_interactions/bonded_interaction_data.hpp" -#include "core/collision.hpp" +#include "core/collision_detection/CollisionDetection.hpp" -#include #include -#include -#include -#include -#include -#include -#include -#include -namespace ScriptInterface { -namespace CollisionDetection { +namespace ScriptInterface::CollisionDetection { class CollisionDetection : public AutoParameters { - std::shared_ptr<::CollisionDetection> m_handle; + std::shared_ptr<::CollisionDetection::CollisionDetection> m_handle; std::unique_ptr m_params; - std::weak_ptr m_bonded_ias; - - std::unordered_map const cd_mode_to_name = { - {CollisionModeType::OFF, "off"}, - {CollisionModeType::BIND_CENTERS, "bind_centers"}, - {CollisionModeType::BIND_VS, "bind_at_point_of_collision"}, - {CollisionModeType::GLUE_TO_SURF, "glue_to_surface"}, - }; - std::unordered_map cd_name_to_mode; - std::unordered_map> const cd_mode_to_parameters = { - {CollisionModeType::OFF, {"mode"}}, - {CollisionModeType::BIND_CENTERS, {"mode", "bond_centers", "distance"}}, - {CollisionModeType::BIND_VS, - {"mode", "bond_centers", "bond_vs", "part_type_vs", "distance", - "vs_placement"}}, - {CollisionModeType::GLUE_TO_SURF, - {"mode", "bond_centers", "bond_vs", "part_type_vs", - "part_type_to_be_glued", "part_type_to_attach_vs_to", - "part_type_after_glueing", "distance", - "distance_glued_particle_to_vs"}}, - }; - - auto find_bond_id(Variant const &v) const { - auto &system = get_system(); - if (is_type(v)) { - auto const bond_id = get_value(v); - std::optional retval = std::nullopt; - if (system.bonded_ias->contains(bond_id)) { - retval = bond_id; - } - return retval; - } - auto obj = get_value>(v); - return system.bonded_ias->find_bond_id(obj->bonded_ia()); - } + std::shared_ptr m_protocol; + std::weak_ptr<::BondedInteractionsMap> m_bonded_ias; + std::weak_ptr m_so_bonded_ias; public: CollisionDetection() { - for (auto const &kv : cd_mode_to_name) { - cd_name_to_mode[kv.second] = kv.first; - } - add_parameters({ - {"mode", - [this](Variant const &v) { - auto const name = get_value(v); - check_mode_name(name); - m_handle->mode = cd_name_to_mode.at(name); - }, - [this]() { return cd_mode_to_name.at(m_handle->mode); }}, - {"bond_centers", - [this](Variant const &v) { - auto const bond_id = find_bond_id(v); - if (not bond_id) { - throw std::invalid_argument("Bond in parameter 'bond_centers' was " - "not added to the system"); - } - m_handle->bond_centers = bond_id.value(); - }, - [this]() { return m_handle->bond_centers; }}, - {"bond_vs", - [this](Variant const &v) { - auto const bond_id = find_bond_id(v); - if (not bond_id) { - throw std::invalid_argument( - "Bond in parameter 'bond_vs' was not added to the system"); - } - m_handle->bond_vs = bond_id.value(); - }, - [this]() { return m_handle->bond_vs; }}, - {"distance", - [this](Variant const &v) { - m_handle->distance = get_value(v); - }, - [this]() { return m_handle->distance; }}, - {"distance_glued_particle_to_vs", - [this](Variant const &v) { - m_handle->dist_glued_part_to_vs = get_value(v); - }, - [this]() { return m_handle->dist_glued_part_to_vs; }}, - {"vs_placement", - [this](Variant const &v) { - m_handle->vs_placement = get_value(v); - }, - [this]() { return m_handle->vs_placement; }}, - {"part_type_vs", - [this](Variant const &v) { - m_handle->vs_particle_type = get_value(v); - }, - [this]() { return m_handle->vs_particle_type; }}, - {"part_type_to_be_glued", - [this](Variant const &v) { - m_handle->part_type_to_be_glued = get_value(v); - }, - [this]() { return m_handle->part_type_to_be_glued; }}, - {"part_type_to_attach_vs_to", - [this](Variant const &v) { - m_handle->part_type_to_attach_vs_to = get_value(v); - }, - [this]() { return m_handle->part_type_to_attach_vs_to; }}, - {"part_type_after_glueing", - [this](Variant const &v) { - m_handle->part_type_after_glueing = get_value(v); - }, - [this]() { return m_handle->part_type_after_glueing; }}, - }); + add_parameters( + {{"protocol", + [this](Variant const &value) { + if (is_none(value)) { + m_protocol = nullptr; + m_handle->unset_protocol(); + return; + } + auto const m_protocol_backup = m_protocol; + try { + context()->parallel_try_catch([&]() { + m_protocol = get_value>(value); + m_protocol->bind_bonded_ias(m_bonded_ias, m_so_bonded_ias); + m_handle->set_protocol(m_protocol->protocol()); + }); + } catch (...) { + m_protocol = m_protocol_backup; + if (m_protocol) { + m_handle->set_protocol(m_protocol->protocol()); + } else { + m_protocol = nullptr; + m_handle->unset_protocol(); + } + throw; + } + }, + [this]() { + if (m_protocol) + return make_variant(m_protocol); + return make_variant(none); + }}}); } void do_construct(VariantMap const ¶ms) override { m_params = std::make_unique(params); - if (params.empty()) { - (*m_params)["mode"] = std::string("off"); - } else { - // Assume we are reloading from a checkpoint file. - // This else branch can be removed once issue #4483 is fixed. - m_params = std::make_unique(); - auto const name = get_value(params, "mode"); - check_mode_name(name); - for (auto const &name : - cd_mode_to_parameters.at(cd_name_to_mode.at(name))) { - (*m_params)[name] = params.at(name); - } - (*m_params)["mode"] = params.at("mode"); - } - } - - Variant do_call_method(const std::string &name, - const VariantMap ¶ms) override { - if (name == "set_params") { - context()->parallel_try_catch([this, ¶ms]() { - auto const backup = std::make_shared<::CollisionDetection>(*m_handle); - auto &system = get_system(); - try { - // check provided parameters - check_input_parameters(params); - // set parameters - for (auto const &kv : params) { - do_set_parameter(get_value(kv.first), kv.second); - } - // sanitize parameters and calculate derived parameters - m_handle->initialize(); - return none; - } catch (...) { - // restore original parameters and re-throw exception - m_handle = system.collision_detection = backup; - m_handle->initialize(); - throw; - } - }); - } - if (name == "params_for_mode") { - auto const name = get_value(params, "mode"); - check_mode_name(name); - auto const mode = cd_name_to_mode.at(name); - return make_vector_of_variants(cd_mode_to_parameters.at(mode)); - } - if (name == "get_bond_by_id") { - if (not context()->is_head_node()) { - return {}; - } - return m_bonded_ias.lock()->call_method("get_bond", params); - } - return none; - } - - void attach(std::weak_ptr bonded_ias) { - m_bonded_ias = bonded_ias; - } - -private: - void check_mode_name(std::string const &name) const { - if (not cd_name_to_mode.contains(name)) { - throw std::invalid_argument("Unknown collision mode '" + name + "'"); - } } void on_bind_system(::System::System &system) override { m_handle = system.collision_detection; m_handle->bind_system(m_system.lock()); - if (not m_params->empty()) { - do_call_method("set_params", *m_params); - } - m_params.reset(); + m_bonded_ias = system.bonded_ias; } - void check_input_parameters(VariantMap const ¶ms) const { - auto const name = get_value(params, "mode"); - check_mode_name(name); - auto const mode = cd_name_to_mode.at(name); - auto const expected_params = cd_mode_to_parameters.at(mode); - auto const expected_param_names = - std::set{expected_params.begin(), expected_params.end()}; - std::set input_parameter_names = {}; - for (auto const &kv : params) { - auto const ¶m_name = kv.first; - if (not expected_param_names.contains(param_name)) { - // throw if parameter is unknown - std::ignore = get_parameter(param_name); - // throw if parameter is known but doesn't match the selected mode - throw std::runtime_error("Parameter '" + param_name + "' is not " + - "required for mode '" + name + "'"); - } - input_parameter_names.insert(param_name); - } - for (auto const ¶m_name : expected_param_names) { - if (not input_parameter_names.contains(param_name)) { - throw std::runtime_error("Parameter '" + param_name + "' is " + - "required for mode '" + name + "'"); - } + void attach(std::weak_ptr bonded_ias) { + m_so_bonded_ias = bonded_ias; + if (m_params) { + AutoParameters::do_construct(*m_params); + m_params.reset(); } } }; -} /* namespace CollisionDetection */ -} /* namespace ScriptInterface */ +} // namespace ScriptInterface::CollisionDetection #endif // COLLISION_DETECTION diff --git a/src/script_interface/collision_detection/GlueToSurface.hpp b/src/script_interface/collision_detection/GlueToSurface.hpp new file mode 100644 index 0000000000..b506667fcf --- /dev/null +++ b/src/script_interface/collision_detection/GlueToSurface.hpp @@ -0,0 +1,109 @@ +/* + * Copyright (C) 2021-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION +#ifdef VIRTUAL_SITES_RELATIVE + +#include "Protocol.hpp" + +#include "core/collision_detection/GlueToSurface.hpp" + +#include "script_interface/ScriptInterface.hpp" +#include "script_interface/auto_parameters/AutoParameters.hpp" + +#include +#include + +namespace ScriptInterface::CollisionDetection { + +class GlueToSurface : public Protocol { + using CoreClass = ::CollisionDetection::GlueToSurface; + +public: + GlueToSurface() { + add_parameters({ + {"bond_centers", AutoParameter::read_only, + [this]() { + return get_bond_variant_by_id( + std::get(*m_protocol).bond_centers); + }}, + {"distance", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).distance; }}, + {"bond_vs", AutoParameter::read_only, + [this]() { + return get_bond_variant_by_id( + std::get(*m_protocol).bond_vs); + }}, + {"part_type_vs", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).part_type_vs; }}, + {"distance_glued_particle_to_vs", AutoParameter::read_only, + [this]() { + return std::get(*m_protocol).dist_glued_part_to_vs; + }}, + {"part_type_to_be_glued", AutoParameter::read_only, + [this]() { + return std::get(*m_protocol).part_type_to_be_glued; + }}, + {"part_type_to_attach_vs_to", AutoParameter::read_only, + [this]() { + return std::get(*m_protocol).part_type_to_attach_vs_to; + }}, + {"part_type_after_glueing", AutoParameter::read_only, + [this]() { + return std::get(*m_protocol).part_type_after_glueing; + }}, + {"_bond_centers", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).bond_centers; }}, + {"_bond_vs", AutoParameter::read_only, + [this]() { return std::get(*m_protocol).bond_vs; }}, + }); + } + std::shared_ptr<::CollisionDetection::ActiveProtocol> protocol() override { + return m_protocol; + } + +private: + std::shared_ptr<::CollisionDetection::ActiveProtocol> m_protocol; + +protected: + void do_initialize(VariantMap const ¶ms) override { + auto const bond_center = find_bond_id(params.contains("_bond_centers") + ? params.at("_bond_centers") + : params.at("bond_centers")); + auto const bond_vs = + find_bond_id(params.contains("_bond_vs") ? params.at("_bond_vs") + : params.at("bond_vs")); + m_protocol = std::make_shared<::CollisionDetection::ActiveProtocol>( + CoreClass(get_value(params, "distance"), bond_center, bond_vs, + get_value(params, "part_type_vs"), + get_value(params, "distance_glued_particle_to_vs"), + get_value(params, "part_type_to_be_glued"), + get_value(params, "part_type_to_attach_vs_to"), + get_value(params, "part_type_after_glueing"))); + } +}; + +} // namespace ScriptInterface::CollisionDetection + +#endif // COLLISION_DETECTION +#endif // VIRTUAL_SITES_RELATIVE diff --git a/src/script_interface/collision_detection/Off.hpp b/src/script_interface/collision_detection/Off.hpp new file mode 100644 index 0000000000..7ee0d97f2e --- /dev/null +++ b/src/script_interface/collision_detection/Off.hpp @@ -0,0 +1,55 @@ +/* + * Copyright (C) 2021-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "Protocol.hpp" + +#include "core/collision_detection/Off.hpp" + +#include "script_interface/ScriptInterface.hpp" +#include "script_interface/auto_parameters/AutoParameters.hpp" + +#include + +namespace ScriptInterface::CollisionDetection { + +class Off : public Protocol { +public: + Off() + : m_protocol{std::make_shared<::CollisionDetection::ActiveProtocol>( + ::CollisionDetection::Off())} {} + + std::shared_ptr<::CollisionDetection::ActiveProtocol> protocol() override { + return m_protocol; + } + +protected: + void do_initialize(VariantMap const &) override {} + +private: + std::shared_ptr<::CollisionDetection::ActiveProtocol> m_protocol; +}; + +} // namespace ScriptInterface::CollisionDetection +#endif // COLLISION_DETECTION diff --git a/src/script_interface/collision_detection/Protocol.hpp b/src/script_interface/collision_detection/Protocol.hpp new file mode 100644 index 0000000000..817d38cd9b --- /dev/null +++ b/src/script_interface/collision_detection/Protocol.hpp @@ -0,0 +1,94 @@ +/* + * Copyright (C) 2021-2024 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 . + */ + +#pragma once + +#include + +#ifdef COLLISION_DETECTION + +#include "script_interface/interactions/BondedInteractions.hpp" + +#include "core/bonded_interactions/bonded_interaction_data.hpp" +#include "core/collision_detection/ActiveProtocol.hpp" + +#include "script_interface/ScriptInterface.hpp" +#include "script_interface/auto_parameters/AutoParameters.hpp" + +#include +#include +#include +#include +#include + +namespace ScriptInterface::CollisionDetection { + +class Protocol : public AutoParameters { + std::weak_ptr m_bonded_ias; + std::weak_ptr m_so_bonded_ias; + std::unique_ptr m_params; + +public: + virtual std::shared_ptr<::CollisionDetection::ActiveProtocol> protocol() = 0; + void bind_bonded_ias( + std::weak_ptr bonded_ias, + std::weak_ptr so_bonded_ias) { + m_bonded_ias = std::move(bonded_ias); + m_so_bonded_ias = std::move(so_bonded_ias); + if (m_params) { + do_initialize(*m_params); + m_params.reset(); + } + } + + void do_construct(VariantMap const ¶ms) override { + m_params = std::make_unique(params); + } + +protected: + auto find_bond_id(Variant const &v) const { + auto bonded_ias = m_bonded_ias.lock(); + assert(bonded_ias && "This protocol is not bound to a system"); + std::optional retval = std::nullopt; + if (is_type(v)) { + auto const bond_id = get_value(v); + if (bonded_ias->contains(bond_id)) { + retval = bond_id; + } + } else { + auto obj = get_value>(v); + retval = bonded_ias->find_bond_id(obj->bonded_ia()); + } + if (not retval) { + throw std::invalid_argument("Bond in parameter list was " + "not added to the system"); + } + return *retval; + } + auto get_bond_variant_by_id(int bond_id) { + auto so_bonded_ias = m_so_bonded_ias.lock(); + assert(so_bonded_ias != nullptr); + return so_bonded_ias->do_call_method("get", {{"key", bond_id}}); + } + virtual void do_initialize(VariantMap const ¶ms) = 0; +}; + +} // namespace ScriptInterface::CollisionDetection + +#endif // COLLISION_DETECTION diff --git a/src/script_interface/collision_detection/initialize.cpp b/src/script_interface/collision_detection/initialize.cpp index daa06ed44a..79fdd06184 100644 --- a/src/script_interface/collision_detection/initialize.cpp +++ b/src/script_interface/collision_detection/initialize.cpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2015-2022 The ESPResSo project + * Copyright (C) 2015-2024 The ESPResSo project * * This file is part of ESPResSo. * @@ -16,20 +16,28 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#include "config/config.hpp" -#include "initialize.hpp" +#include +#include "BindAtPointOfCollision.hpp" +#include "BindCenters.hpp" #include "CollisionDetection.hpp" +#include "GlueToSurface.hpp" +#include "Off.hpp" +#include "initialize.hpp" -namespace ScriptInterface { -namespace CollisionDetection { - +namespace ScriptInterface::CollisionDetection { void initialize(Utils::Factory *om) { #ifdef COLLISION_DETECTION om->register_new( "CollisionDetection::CollisionDetection"); -#endif + om->register_new("CollisionDetection::Off"); + om->register_new("CollisionDetection::BindCenters"); +#ifdef VIRTUAL_SITES_RELATIVE + om->register_new( + "CollisionDetection::BindAtPointOfCollision"); + om->register_new("CollisionDetection::GlueToSurface"); +#endif // VIRTUAL_SITES_RELATIVE +#endif // COLLISION_DETECTION } -} /* namespace CollisionDetection */ -} /* namespace ScriptInterface */ +} // namespace ScriptInterface::CollisionDetection diff --git a/src/script_interface/collision_detection/initialize.hpp b/src/script_interface/collision_detection/initialize.hpp index 5eb93c727b..4663833bc1 100644 --- a/src/script_interface/collision_detection/initialize.hpp +++ b/src/script_interface/collision_detection/initialize.hpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2015-2022 The ESPResSo project + * Copyright (C) 2015-2024 The ESPResSo project * * This file is part of ESPResSo. * @@ -17,19 +17,12 @@ * along with this program. If not, see . */ -#ifndef SCRIPT_INTERFACE_COLLISIONDETECTION_INITIALIZE_HPP -#define SCRIPT_INTERFACE_COLLISIONDETECTION_INITIALIZE_HPP +#pragma once #include #include -namespace ScriptInterface { -namespace CollisionDetection { - +namespace ScriptInterface::CollisionDetection { void initialize(Utils::Factory *om); - -} /* namespace CollisionDetection */ -} /* namespace ScriptInterface */ - -#endif +} // namespace ScriptInterface::CollisionDetection diff --git a/src/script_interface/interactions/BondedInteractions.hpp b/src/script_interface/interactions/BondedInteractions.hpp index ac9f978875..9bb8ebc73d 100644 --- a/src/script_interface/interactions/BondedInteractions.hpp +++ b/src/script_interface/interactions/BondedInteractions.hpp @@ -92,7 +92,7 @@ class BondedInteractions : public BondedInteractionsBase_t { m_bonds.erase(key); } -protected: +public: Variant do_call_method(std::string const &name, VariantMap const ¶ms) override { if (name == "get_size") { diff --git a/testsuite/python/bond_breakage.py b/testsuite/python/bond_breakage.py index c0565641d8..090f93f96e 100644 --- a/testsuite/python/bond_breakage.py +++ b/testsuite/python/bond_breakage.py @@ -265,6 +265,7 @@ def tearDown(self): self.system.part.clear() self.system.bonded_inter.clear() self.system.thermostat.turn_off() + self.system.min_global_cut = 0.6 @utx.skipIfMissingFeatures(["COLLISION_DETECTION"]) def test_center_bonds(self): @@ -274,16 +275,17 @@ def test_center_bonds(self): crit = 2**(1 / 6) * 2. - self.system.collision_detection.set_params(mode="bind_centers", - distance=2**(1 / 6) * 2.2, bond_centers=harm) + self.system.collision_detection.protocol = espressomd.collision_detection.BindCenters( + distance=2**(1 / 6) * 2.2, bond_centers=harm) self.system.integrator.run(1) - self.system.collision_detection.set_params(mode="off") + self.system.collision_detection.protocol = espressomd.collision_detection.Off() self.system.bond_breakage[harm] = BreakageSpec( breakage_length=crit, action_type="delete_bond") self.system.integrator.run(1) bonds_dist = 0 + self.system.min_global_cut = crit pairs = self.system.cell_system.get_pairs(crit, types=[0]) for pair in pairs: dist = self.system.distance( @@ -307,12 +309,11 @@ def test_vs_bonds(self): crit = 2**(1 / 6) * 1.5 crit_vs = 2**(1 / 6) * 1 / 3 * 1.2 - self.system.collision_detection.set_params(mode="bind_at_point_of_collision", - distance=crit, bond_centers=virt, bond_vs=harm, - part_type_vs=1, vs_placement=1 / 3) + self.system.collision_detection.protocol = espressomd.collision_detection.BindAtPointOfCollision( + distance=crit, bond_centers=virt, bond_vs=harm, part_type_vs=1, vs_placement=1 / 3) self.system.integrator.run(1) - self.system.collision_detection.set_params(mode="off") + self.system.collision_detection.protocol = espressomd.collision_detection.Off() self.system.bond_breakage[harm] = BreakageSpec( breakage_length=crit_vs, action_type="revert_bind_at_point_of_collision") self.system.integrator.run(1) diff --git a/testsuite/python/collision_detection.py b/testsuite/python/collision_detection.py index 1f0a2ae241..fa68dc47e4 100644 --- a/testsuite/python/collision_detection.py +++ b/testsuite/python/collision_detection.py @@ -64,7 +64,7 @@ def get_state_set_state_consistency(self): def test_bind_centers(self): system = self.system # Check that it leaves particles alone, when off - system.collision_detection.set_params(mode="off") + system.collision_detection.protocol = espressomd.collision_detection.Off() system.part.clear() p0 = system.part.add(pos=(0, 0, 0), id=0) @@ -74,10 +74,9 @@ def test_bind_centers(self): self.assertEqual(p0.bonds, ()) self.assertEqual(p1.bonds, ()) self.assertEqual(p2.bonds, ()) - # Check that it cannot be activated - system.collision_detection.set_params( - mode="bind_centers", distance=0.11, bond_centers=self.bond_center) + system.collision_detection.protocol = espressomd.collision_detection.BindCenters( + distance=0.11, bond_centers=self.bond_center) self.get_state_set_state_consistency() system.integrator.run(1, recalc_forces=True) bond0 = ((system.bonded_inter[0], 1),) @@ -93,9 +92,10 @@ def test_bind_centers(self): self.assertEqual(p2.bonds, ()) # Check turning it off - system.collision_detection.set_params(mode="off") + system.collision_detection.protocol = espressomd.collision_detection.Off() self.get_state_set_state_consistency() - self.assertEqual(system.collision_detection.mode, "off") + self.assertIsInstance( + self.system.collision_detection.protocol, espressomd.collision_detection.Off) def run_test_bind_at_point_of_collision_for_pos(self, *positions): system = self.system @@ -113,9 +113,8 @@ def run_test_bind_at_point_of_collision_for_pos(self, *positions): # 2 non-virtual + 2 virtual + one that doesn't take part expected_np = 4 * len(positions) + 1 - system.collision_detection.set_params( - mode="bind_at_point_of_collision", bond_centers=self.bond_center, - bond_vs=self.bond_vs, part_type_vs=1, vs_placement=0.4, distance=0.11) + system.collision_detection.protocol = espressomd.collision_detection.BindAtPointOfCollision( + bond_centers=self.bond_center, bond_vs=self.bond_vs, part_type_vs=1, vs_placement=0.4, distance=0.11) self.get_state_set_state_consistency() system.integrator.run(1, recalc_forces=True) self.verify_state_after_bind_at_poc(expected_np) @@ -132,7 +131,7 @@ def run_test_bind_at_point_of_collision_for_pos(self, *positions): self.verify_state_after_bind_at_poc(expected_np) def verify_state_after_bind_at_poc(self, expected_np): - if self.system.collision_detection.bond_vs == self.bond_angle_vs: + if self.system.collision_detection.protocol.bond_vs == self.bond_angle_vs: self.verify_state_after_bind_at_poc_triplet(expected_np) else: self.verify_state_after_bind_at_poc_pair(expected_np) @@ -227,7 +226,7 @@ def verify_bind_at_poc(self, p1, p2, vs1, vs2): system = self.system # Check for presence of vs # Check for bond between vs - if self.system.collision_detection.bond_vs == self.bond_angle_vs: + if self.system.collision_detection.protocol.bond_vs == self.bond_angle_vs: bond_p1 = ((self.bond_pair, p2.id), (self.bond_center, p2.id),) bond_p2 = ((self.bond_pair, p1.id), (self.bond_center, p1.id),) self.assertTrue(p1.bonds == bond_p1 or p2.bonds == bond_p2) @@ -265,7 +264,7 @@ def verify_bind_at_poc(self, p1, p2, vs1, vs2): else: dist_centers = p1.pos - p2.pos expected_pos = system.part.by_id(rel_to).pos_folded + \ - system.collision_detection.vs_placement * dist_centers + system.collision_detection.protocol.vs_placement * dist_centers dist = expected_pos - p.pos_folded dist -= np.round(dist / system.box_l) * system.box_l self.assertLess(np.linalg.norm(dist), 1E-12) @@ -311,8 +310,8 @@ def test_bind_at_point_of_collision_triplet(self): # 2 non-virtual + 2 virtual + one that doesn't take part expected_np = 4 * len(positions) + 1 - system.collision_detection.set_params( - mode="bind_at_point_of_collision", bond_centers=self.bond_center, + system.collision_detection.protocol = espressomd.collision_detection.BindAtPointOfCollision( + bond_centers=self.bond_center, bond_vs=self.bond_angle_vs, part_type_vs=1, vs_placement=0.4, distance=0.11) self.get_state_set_state_consistency() system.integrator.run(1, recalc_forces=True) @@ -355,13 +354,8 @@ def test_bind_at_point_of_collision_random(self): system.integrator.run(10) # Collision detection - system.collision_detection.set_params( - mode="bind_at_point_of_collision", - distance=0.11, - bond_centers=self.bond_center, - bond_vs=self.bond_vs, - part_type_vs=1, - vs_placement=0.4) + system.collision_detection.protocol = espressomd.collision_detection.BindAtPointOfCollision( + bond_centers=self.bond_center, bond_vs=self.bond_vs, part_type_vs=1, vs_placement=0.4, distance=0.11) self.get_state_set_state_consistency() # Integrate lj liquid @@ -436,8 +430,8 @@ def run_test_glue_to_surface_for_pos(self, *positions): # 2 non-virtual + 1 virtual + one that doesn't take part expected_np = 3 * len(positions) + 1 - system.collision_detection.set_params( - mode="glue_to_surface", distance=0.11, + system.collision_detection.protocol = espressomd.collision_detection.GlueToSurface( + distance=0.11, distance_glued_particle_to_vs=0.02, bond_centers=self.bond_center, bond_vs=self.bond_vs, part_type_vs=self.part_type_vs, part_type_to_attach_vs_to=self.part_type_to_attach_vs_to, @@ -595,8 +589,8 @@ def test_glue_to_surface_random(self): system.integrator.run(10) # Collision detection - system.collision_detection.set_params( - mode="glue_to_surface", distance=0.11, + system.collision_detection.protocol = espressomd.collision_detection.GlueToSurface( + distance=0.11, distance_glued_particle_to_vs=0.02, bond_centers=self.bond_center, bond_vs=self.bond_vs, part_type_vs=self.part_type_vs, part_type_to_attach_vs_to=self.part_type_to_attach_vs_to, diff --git a/testsuite/python/collision_detection_interface.py b/testsuite/python/collision_detection_interface.py index 3e2e5dbe7a..47005002d9 100644 --- a/testsuite/python/collision_detection_interface.py +++ b/testsuite/python/collision_detection_interface.py @@ -53,35 +53,22 @@ class CollisionDetection(ut.TestCase): "part_type_to_be_glued": 2, "part_type_after_glueing": 3 }, } + classes = {"bind_centers": espressomd.collision_detection.BindCenters, + "bind_at_point_of_collision": espressomd.collision_detection.BindAtPointOfCollision, + "glue_to_surface": espressomd.collision_detection.GlueToSurface} def tearDown(self): - self.system.collision_detection.set_params(mode="off") + self.system.collision_detection.protocol = None def test_00_interface_and_defaults(self): - # Is it off by default - self.assertEqual(self.system.collision_detection.mode, "off") - - # Make sure params cannot be set individually - with self.assertRaises(Exception): - self.system.collision_detection.mode = "bind_centers" - - # Verify exception throwing for unknown collision modes - with self.assertRaisesRegex(ValueError, "Unknown collision mode 'unknown'"): - self.system.collision_detection.set_params(mode="unknown") - with self.assertRaisesRegex(ValueError, "Collision mode must be specified via the 'mode' argument"): - self.system.collision_detection.set_params() - + self.assertIsNone(self.system.collision_detection.protocol) self.assertIsNone(self.system.collision_detection.call_method("none")) - # That should work - self.system.collision_detection.set_params(mode="off") - self.assertEqual(self.system.collision_detection.mode, "off") - def check_stored_parameters(self, mode, **kwargs): """ Check if collision detection stored parameters match input values. """ - parameters = self.system.collision_detection.get_params() + parameters = self.system.collision_detection.protocol.get_params() parameters_ref = self.valid_coldet_params[mode].copy() parameters_ref.update(kwargs) for key, value_ref in parameters_ref.items(): @@ -96,7 +83,14 @@ def set_coldet(self, mode, **invalid_params): """ params = self.valid_coldet_params.get(mode).copy() params.update(invalid_params) - self.system.collision_detection.set_params(mode=mode, **params) + self.system.collision_detection.protocol = self.classes[mode](**params) + + def test_off(self): + self.system.collision_detection.protocol = None + with self.assertRaisesRegex(ValueError, "Parameter 'distance' must be > 0"): + self.set_coldet("bind_centers", distance=-2.) + # check if original parameters have been preserved + self.assertIsNone(self.system.collision_detection.protocol) def test_bind_centers(self): self.set_coldet("bind_centers", distance=0.5) @@ -104,31 +98,29 @@ def test_bind_centers(self): self.set_coldet("bind_centers", distance=-2.) with self.assertRaisesRegex(ValueError, "Parameter 'distance' must be > 0"): self.set_coldet("bind_centers", distance=0.) - with self.assertRaisesRegex(ValueError, "Bond in parameter 'bond_centers' was not added to the system"): + with self.assertRaisesRegex(ValueError, "Bond in parameter list was not added to the system"): bond = espressomd.interactions.HarmonicBond(k=1., r_0=0.1) self.set_coldet("bind_centers", bond_centers=bond) with self.assertRaisesRegex(RuntimeError, "The bond type to be used for binding particle centers needs to be a pair bond"): self.set_coldet("bind_centers", bond_centers=self.bond_angle) - with self.assertRaisesRegex(RuntimeError, "Unknown parameter 'unknown'"): - self.set_coldet("bind_centers", unknown=1) - with self.assertRaisesRegex(RuntimeError, "Parameter 'part_type_vs' is not required for mode 'bind_centers'"): - self.set_coldet("bind_centers", part_type_vs=1) - with self.assertRaisesRegex(RuntimeError, "Parameter 'distance' is required for mode 'bind_centers'"): - self.system.collision_detection.set_params( - mode="bind_centers", bond_centers=self.bond_harmonic) - with self.assertRaisesRegex(Exception, "Please set all parameters at once via collision_detection.set_params"): - self.system.collision_detection.mode = "bind_at_point_of_collision" + with self.assertRaisesRegex(RuntimeError, "Parameter 'distance' is missing"): + self.system.collision_detection.protocol = espressomd.collision_detection.BindCenters( + bond_centers=self.bond_harmonic) # check if original parameters have been preserved self.check_stored_parameters("bind_centers", distance=0.5) @utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE") def test_bind_at_point_of_collision(self): self.set_coldet("bind_at_point_of_collision", distance=0.5) + with self.assertRaisesRegex(ValueError, "Parameter 'distance' must be > 0"): + self.set_coldet("bind_at_point_of_collision", distance=-2.) + with self.assertRaisesRegex(ValueError, "Parameter 'distance' must be > 0"): + self.set_coldet("bind_at_point_of_collision", distance=0.) with self.assertRaisesRegex(ValueError, "Parameter 'vs_placement' must be between 0 and 1"): self.set_coldet("bind_at_point_of_collision", vs_placement=-0.01) with self.assertRaisesRegex(ValueError, "Parameter 'vs_placement' must be between 0 and 1"): self.set_coldet("bind_at_point_of_collision", vs_placement=1.01) - with self.assertRaisesRegex(ValueError, "Bond in parameter 'bond_vs' was not added to the system"): + with self.assertRaisesRegex(ValueError, "Bond in parameter list 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 bond"): @@ -143,12 +135,16 @@ def test_bind_at_point_of_collision(self): @utx.skipIfMissingFeatures("VIRTUAL_SITES") def test_bind_at_point_of_collision_norotation(self): if not espressomd.has_features("VIRTUAL_SITES_RELATIVE"): - with self.assertRaisesRegex(RuntimeError, "require the VIRTUAL_SITES_RELATIVE feature"): + with self.assertRaisesRegex(RuntimeError, "Missing features VIRTUAL_SITES_RELATIVE"): self.set_coldet("bind_at_point_of_collision") @utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE") def test_glue_to_surface(self): self.set_coldet("glue_to_surface", distance=0.5) + with self.assertRaisesRegex(ValueError, "Parameter 'distance' must be > 0"): + self.set_coldet("glue_to_surface", distance=-2.) + with self.assertRaisesRegex(ValueError, "Parameter 'distance' must be > 0"): + self.set_coldet("glue_to_surface", distance=0.) with self.assertRaisesRegex(ValueError, "type for virtual sites needs to be >=0"): self.set_coldet("glue_to_surface", part_type_vs=-1) with self.assertRaisesRegex(ValueError, "type to be glued needs to be >=0"): diff --git a/testsuite/python/lees_edwards.py b/testsuite/python/lees_edwards.py index 5c6fc10400..f2259dafb8 100644 --- a/testsuite/python/lees_edwards.py +++ b/testsuite/python/lees_edwards.py @@ -79,7 +79,7 @@ def tearDown(self): system.bonded_inter.clear() system.lees_edwards.protocol = None if espressomd.has_features("COLLISION_DETECTION"): - system.collision_detection.set_params(mode="off") + system.collision_detection.protocol = espressomd.collision_detection.Off() def test_00_is_none_by_default(self): @@ -617,8 +617,8 @@ def test_le_colldet(self): virt = espressomd.interactions.Virtual() system.bonded_inter.add(virt) - system.collision_detection.set_params( - mode="bind_centers", distance=1., bond_centers=harm) + system.collision_detection.protocol = espressomd.collision_detection.BindCenters( + distance=1., bond_centers=harm) # After two integration steps we should not have a bond, # as the collision detection uses the distant calculation @@ -635,15 +635,15 @@ def test_le_colldet(self): np.testing.assert_array_equal(len(bond_list), 1) system.part.clear() - system.collision_detection.set_params(mode="off") + system.collision_detection.protocol = espressomd.collision_detection.Off() system.time = 0 system.lees_edwards.protocol = espressomd.lees_edwards.LinearShear( shear_velocity=-1.0, initial_pos_offset=0.0) - system.collision_detection.set_params( - mode="bind_at_point_of_collision", distance=1., bond_centers=virt, - bond_vs=harm, part_type_vs=31, vs_placement=1 / 3) + system.collision_detection.protocol = espressomd.collision_detection.BindAtPointOfCollision( + distance=1., bond_centers=virt, bond_vs=harm, part_type_vs=31, + vs_placement=1. / 3.) col_part1 = system.part.add( pos=(2.5, 4.5, 2.5), type=30, fix=[True, True, True]) diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index c40e54caba..29e8d387de 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -277,6 +277,8 @@ # bonded interactions harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0.0, k=1.0) system.bonded_inter.add(harmonic_bond) +strong_harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0.0, k=5e5) +system.bonded_inter.add(strong_harmonic_bond) p2.add_bond((harmonic_bond, p1)) if 'THERM.LB' in modes or 'THERM.LANGEVIN' in modes: # create Drude particles @@ -298,8 +300,6 @@ thermalized_bond=therm_bond1, p_core=p2, type_drude=10, alpha=1., mass_drude=0.6, coulomb_prefactor=0.8, thole_damping=2.) checkpoint.register("dh") -strong_harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0.0, k=5e5) -system.bonded_inter.add(strong_harmonic_bond) p4.add_bond((strong_harmonic_bond, p3)) ibm_volcons_bond = espressomd.interactions.IBM_VolCons(softID=15, kappaV=0.01) ibm_tribend_bond = espressomd.interactions.IBM_Tribend( @@ -325,8 +325,14 @@ checkpoint.register("particle_force0") checkpoint.register("particle_force1") if espressomd.has_features("COLLISION_DETECTION"): - system.collision_detection.set_params( - mode="bind_centers", distance=0.11, bond_centers=harmonic_bond) + if espressomd.has_features("VIRTUAL_SITES_RELATIVE"): + protocol = espressomd.collision_detection.BindAtPointOfCollision( + distance=0.12, bond_centers=harmonic_bond, + bond_vs=strong_harmonic_bond, part_type_vs=2, vs_placement=1. / 3.) + else: + protocol = espressomd.collision_detection.BindCenters( + distance=0.11, bond_centers=harmonic_bond) + system.collision_detection.protocol = protocol particle_propagation0 = p1.propagation particle_propagation1 = p2.propagation diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 1594902557..8ac84fda12 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -682,15 +682,21 @@ def test_bonded_inter(self): bond_ids = system.bonded_inter.call_method('get_bond_ids') self.assertEqual(len(bond_ids), len(system.bonded_inter)) # check bonded interactions - partcl_1 = system.part.by_id(1) + p1 = system.part.by_id(1) + p4 = system.part.by_id(4) reference = {'r_0': 0.0, 'k': 1.0, 'r_cut': 0.0} - self.assertEqual(partcl_1.bonds[0][0].params, reference) - self.assertEqual(system.bonded_inter[0].params, reference) + self.assertEqual(p1.bonds[0][0].params, reference) + self.assertAlmostEqual( + system.bonded_inter[0].params, reference, delta=1e-6) + reference = {'r_0': 0.0, 'k': 5e5, 'r_cut': 0.0} + self.assertEqual(p4.bonds[0][0].params, reference) + self.assertAlmostEqual( + system.bonded_inter[1].params, reference, delta=1e-1) # all thermalized bonds should be identical if has_drude: reference = therm_params - self.assertEqual(partcl_1.bonds[1][0].params, reference) - self.assertEqual(system.bonded_inter[1].params, reference) + self.assertEqual(p1.bonds[1][0].params, reference) + self.assertEqual(system.bonded_inter[2].params, reference) self.assertEqual(therm_bond2.params, reference) # immersed boundary bonds self.assertEqual( @@ -921,10 +927,20 @@ def test_comfixed(self): @utx.skipIfMissingFeatures('COLLISION_DETECTION') def test_collision_detection(self): - coldet = system.collision_detection - self.assertEqual(coldet.mode, "bind_centers") - self.assertAlmostEqual(coldet.distance, 0.11, delta=1E-9) - self.assertEqual(coldet.bond_centers, system.bonded_inter[0]) + protocol = system.collision_detection.protocol + if espressomd.has_features("VIRTUAL_SITES_RELATIVE"): + self.assertIsInstance( + protocol, espressomd.collision_detection.BindAtPointOfCollision) + self.assertAlmostEqual(protocol.distance, 0.12, delta=1E-9) + self.assertEqual(protocol.bond_centers, system.bonded_inter[0]) + self.assertEqual(protocol.bond_vs, system.bonded_inter[1]) + self.assertEqual(protocol.part_type_vs, 2) + self.assertAlmostEqual(protocol.vs_placement, 1. / 3., delta=1e-6) + else: + self.assertIsInstance( + protocol, espressomd.collision_detection.BindCenters) + self.assertAlmostEqual(protocol.distance, 0.11, delta=1E-9) + self.assertEqual(protocol.bond_centers, system.bonded_inter[0]) @utx.skipIfMissingFeatures('EXCLUSIONS') def test_exclusions(self):