diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index af41473123e..c6fd3ecc3d5 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -1,5 +1,6 @@ set(EspressoCore_SRC accumulators.cpp + bond_breakage.cpp bond_error.cpp cells.cpp collision.cpp diff --git a/src/core/bond_breakage.cpp b/src/core/bond_breakage.cpp new file mode 100644 index 00000000000..9549fa15417 --- /dev/null +++ b/src/core/bond_breakage.cpp @@ -0,0 +1,232 @@ +/* + * Copyright (C) 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 "bond_breakage.hpp" +#include "cells.hpp" +#include "communication.hpp" +#include "errorhandling.hpp" +#include "particle_data.hpp" + +#include + +#include +#include + +#include +#include +#include +#include +#include + +namespace BondBreakage { + +// Bond breakage specifications +std::unordered_map> breakage_specs; + +// Delete Actions +struct DeleteBond { + int particle_id; + int bond_partner_id; + int bond_type; + std::size_t hash_value() const { + std::size_t seed = 3875; + boost::hash_combine(seed, particle_id); + boost::hash_combine(seed, bond_partner_id); + boost::hash_combine(seed, bond_type); + return seed; + } + bool operator==(DeleteBond const &rhs) const { + return rhs.particle_id == particle_id and + rhs.bond_partner_id == bond_partner_id and + rhs.bond_type == bond_type; + } +}; + +struct DeleteAllBonds { + int particle_id_1; + int particle_id_2; + std::size_t hash_value() const { + std::size_t seed = 75; + boost::hash_combine(seed, particle_id_1); + boost::hash_combine(seed, particle_id_2); + return seed; + } + bool operator==(DeleteAllBonds const &rhs) const { + return rhs.particle_id_1 == particle_id_1 and + rhs.particle_id_2 == particle_id_2; + } +}; +} // namespace BondBreakage + +// Hash support for std::unordered_set +namespace boost { +template <> struct hash { + std::size_t operator()(BondBreakage::DeleteBond const &t) const noexcept { + return t.hash_value(); + } +}; +template <> struct hash { + std::size_t operator()(BondBreakage::DeleteAllBonds const &t) const noexcept { + return t.hash_value(); + } +}; +} // namespace boost + +namespace BondBreakage { +// Variant holding any of the actions +using Action = boost::variant; + +// Set of actions +using ActionSet = std::unordered_set; + +// Broken bond record +struct QueueEntry { + int particle_id; + int bond_partner_id; + int bond_type; + + /// Serialization for synchronization across mpi ranks + friend class boost::serialization::access; + template + void serialize(Archive &ar, const unsigned int version) { + ar &particle_id; + ar &bond_partner_id; + ar &bond_type; + } +}; + +/** @brief Queue to record bonds broken during a time step */ +using Queue = std::vector; +Queue queue; + +/** @brief Retrieve breakage specification for the bond type */ +boost::optional get_breakage_spec(int bond_type) { + if (breakage_specs.find(bond_type) != breakage_specs.end()) { + return {*(breakage_specs.at(bond_type))}; + } + return {}; +} + +/** Add a particle+bond combination to the breakage queue */ +void queue_breakage(int particle_id, int bond_partner_id, int bond_type) { + queue.emplace_back(QueueEntry{particle_id, bond_partner_id, bond_type}); +} + +bool check_and_handle_breakage(int particle_id, int bond_partner_id, + int bond_type, double distance) { + // Retrieve specification for this bond type + auto spec = get_breakage_spec(bond_type); + if (!spec) + return false; // No breakage rule for this bond type + + // Is the bond length longer than the breakage length? + if (distance >= (*spec).breakage_length) { + queue_breakage(particle_id, bond_partner_id, bond_type); + return true; + } + return false; +} + +void clear_queue() { queue.clear(); } + +/** @brief Gathers combined queue from all mpi ranks */ +Queue gather_global_queue(Queue const &local_queue) { + Queue res = local_queue; + Utils::Mpi::gather_buffer(res, comm_cart); + boost::mpi::broadcast(comm_cart, res, 0); + return res; +} + +/** @brief Constructs the actions to take for a breakage queue entry */ +ActionSet actions_for_breakage(QueueEntry const &e) { + // Retrieve relevant breakage spec + auto const spec = get_breakage_spec(e.bond_type); + assert(spec); + + // Handle different action types + if ((*spec).action_type == ActionType::DELETE_BOND) + return {DeleteBond{e.particle_id, e.bond_partner_id, e.bond_type}}; +#ifdef VIRTUAL_SITES_RELATIVE + if ((*spec).action_type == ActionType::REVERT_BIND_AT_POINT_OF_COLLISION) { + // We need to find the base particles for the two virtual sites + // between which the bond broke. + auto p1 = cell_structure.get_local_particle(e.particle_id); + auto p2 = cell_structure.get_local_particle(e.bond_partner_id); + if (!p1 || !p2) + return {}; // particles not on this mpi rank + + if (!p1->p.is_virtual || !p2->p.is_virtual) { + runtimeErrorMsg() << "The REVERT_BIND_AT_POINT_OF_COLLISION bond " + "breakage action has to be configured for the " + "bond on the virtual site. Encountered a particle " + "that is not virtual."; + return {}; + } + + return { + // Bond between virtual sites + DeleteBond{e.particle_id, e.bond_partner_id, e.bond_type}, + // Bond between base particles. We do not know, on which of the two + // the bond is defined, since bonds are stored only on one partner + DeleteAllBonds{p1->p.vs_relative.to_particle_id, + p2->p.vs_relative.to_particle_id}, + DeleteAllBonds{p2->p.vs_relative.to_particle_id, + p1->p.vs_relative.to_particle_id}, + }; + } +#endif // VIRTUAL_SITES_RELATIVE + return {}; +} + +// Handler for the different delete events +class execute : public boost::static_visitor<> { +public: + void operator()(DeleteBond const &d) const { + auto p = cell_structure.get_local_particle(d.particle_id); + if (!p) + return; + local_remove_bond(*p, {d.bond_type, d.bond_partner_id}); + } + void operator()(DeleteAllBonds const &d) const { + auto p = cell_structure.get_local_particle(d.particle_id_1); + if (!p) + return; + local_remove_pair_bonds_to(*p, d.particle_id_2); + } +}; + +void process_queue() { + if (breakage_specs.empty()) + return; + + auto global_queue = gather_global_queue(queue); + + // Construct delete actions from breakage queue + ActionSet actions = {}; + for (auto const &e : global_queue) { + // Convert to merge() once we are on C++17 + auto to_add = actions_for_breakage(e); + actions.insert(to_add.begin(), to_add.end()); + } + + // Execute actions + for (auto const &a : actions) { + boost::apply_visitor(execute{}, a); + } +} +} // namespace BondBreakage diff --git a/src/core/bond_breakage.hpp b/src/core/bond_breakage.hpp new file mode 100644 index 00000000000..5a6f0bcd689 --- /dev/null +++ b/src/core/bond_breakage.hpp @@ -0,0 +1,50 @@ +/* + * Copyright (C) 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 +#include + +namespace BondBreakage { + +enum class ActionType { + NONE = 0, + DELETE_BOND = 1, + REVERT_BIND_AT_POINT_OF_COLLISION = 2 +}; + +struct BreakageSpec { + double breakage_length; + ActionType action_type; +}; + +extern std::unordered_map> breakage_specs; + +/** @brief Check if the bond between the particles should break, if yes, queue + * it. + */ +bool check_and_handle_breakage(int particle_id, int bond_partner_id, + int bond_type, double distance); + +void clear_queue(); + +void process_queue(); + +} // namespace BondBreakage diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 526dc079d24..3847f17d123 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -26,6 +26,7 @@ #include "EspressoSystemInterface.hpp" +#include "bond_breakage.hpp" #include "cells.hpp" #include "collision.hpp" #include "comfixed_global.hpp" @@ -148,7 +149,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { #ifdef COLLISION_DETECTION prepare_local_collision_queue(); #endif - + BondBreakage::clear_queue(); auto particles = cell_structure.local_particles(); auto ghost_particles = cell_structure.ghost_particles(); #ifdef ELECTROSTATICS diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index c455c01debf..6e781a390da 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -28,6 +28,7 @@ #include "forces.hpp" +#include "bond_breakage.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" #include "bonded_interactions/thermalized_bond_kernel.hpp" #include "immersed_boundary/ibm_tribend.hpp" @@ -411,6 +412,15 @@ inline bool add_bonded_four_body_force(Bonded_IA_Parameters const &iaparams, inline bool add_bonded_force(Particle &p1, int bond_id, Utils::Span partners) { + + // Consider for bond breakage + if (partners.size() == 1) { + auto d = box_geo.get_mi_vector(p1.r.p, partners[0]->r.p).norm(); + if (BondBreakage::check_and_handle_breakage( + p1.p.identity, partners[0]->p.identity, bond_id, d)) + return false; + } + auto const &iaparams = *bonded_ia_params.at(bond_id); switch (number_of_partners(iaparams)) { diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 69950a960fd..581f03bd7bb 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -35,6 +35,7 @@ #include "ParticleRange.hpp" #include "accumulators.hpp" +#include "bond_breakage.hpp" #include "bonded_interactions/rigid_bond.hpp" #include "cells.hpp" #include "collision.hpp" @@ -325,6 +326,7 @@ int integrate(int n_steps, int reuse_forces) { #ifdef COLLISION_DETECTION handle_collisions(); #endif + BondBreakage::process_queue(); } integrated_steps++; diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 1444816d977..501f1c98296 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -199,6 +199,28 @@ struct RemoveBond { void serialize(Archive &ar, long int) { ar & bond; } }; +/** + * @brief Delete pair bonds to a specific partner + */ +struct RemovePairBondsTo { + int other_pid; + + void operator()(Particle &p) const { + using Bond = std::vector; + std::vector to_delete; + for (auto b: p.bonds()) { + if (b.partner_ids().size() == 1 and b.partner_ids()[0] == other_pid) + to_delete.push_back(Bond{b.bond_id(),other_pid}); + } + for (auto b: to_delete) { + RemoveBond{b}(p); + } + } + template + void serialize(Archive &ar, long int) { + ar & other_pid; + } +}; /** * @brief Delete all bonds. @@ -325,6 +347,14 @@ struct UpdateVisitor : public boost::static_visitor { }; } // namespace +void local_remove_bond(Particle &p, std::vector const &bond) { + RemoveBond{bond}(p); +} + +void local_remove_pair_bonds_to(Particle &p, int other_pid) { + RemovePairBondsTo{other_pid}(p); +} + void mpi_send_update_message_local(int node, int id) { if (node == comm_cart.rank()) { UpdateMessage msg{}; diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index 390ae8b7387..ce61404536e 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -299,10 +299,21 @@ void delete_particle_bond(int part, Utils::Span bond); */ void delete_particle_bonds(int part); +/** @brief Removs the specified bond from the particle + * @param bond The bond in the form + * {bond_id, partner_1, partner_2, ...} + */ +void local_remove_bond(Particle &p, std::vector const &bond); + +/** @brief Removes all pair bonds on the particle which have the specified + * particle id as partner + */ +void local_remove_pair_bonds_to(Particle &p, int other_pid); + /** Call only on the head node: Add bond to particle. * @param part identity of principal atom of the bond. - * @param bond field containing the bond type number and the - * identity of all bond partners (secondary atoms of the bond). + * @param bond field containing the bond type number and the identity + * of all bond partners (secondary atoms of the bond). */ void add_particle_bond(int part, Utils::Span bond); diff --git a/src/python/espressomd/bond_breakage.py b/src/python/espressomd/bond_breakage.py new file mode 100644 index 00000000000..ef7228a9ab4 --- /dev/null +++ b/src/python/espressomd/bond_breakage.py @@ -0,0 +1,30 @@ +# +# Copyright (C) 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 . +# +from .script_interface import script_interface_register, ScriptObjectMap, ScriptInterfaceHelper + + +@script_interface_register +class BreakageSpec(ScriptInterfaceHelper): + _so_name = "BondBreakage::BreakageSpec" + + +@script_interface_register +class BreakageSpecs(ScriptObjectMap): + _so_name = "BondBreakage::BreakageSpecs" + _key_type = int diff --git a/src/python/espressomd/script_interface.pyx b/src/python/espressomd/script_interface.pyx index 2d0a207240c..531639c8f1c 100644 --- a/src/python/espressomd/script_interface.pyx +++ b/src/python/espressomd/script_interface.pyx @@ -15,7 +15,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . import numpy as np -from .utils import to_char_pointer, to_str, handle_errors, array_locked +from .utils import to_char_pointer, to_str, handle_errors, array_locked, is_valid_type from .utils cimport Vector2d, Vector3d, Vector4d, make_array_locked cimport cpython.object @@ -385,6 +385,87 @@ class ScriptObjectRegistry(ScriptInterfaceHelper): return self.call_method("size") +class ScriptObjectMap(ScriptObjectRegistry): + + """ + Represents a script interface ObjectMap (dict)-like object + + Item acces via [key]. + """ + + _key_type = int + + def __init__(self, *args, **kwargs): + if args: + params, (_unpickle_so_class, (_so_name, bytestring)) = args + assert _so_name == self._so_name + self = _unpickle_so_class(_so_name, bytestring) + self.__setstate__(params) + else: + super().__init__(**kwargs) + + def remove(self, key): + """ + Removes the element with the given key + """ + # Validate key type + if not is_valid_type(key, self._key_type): + raise ValueError( + f"Key has to be of type {self._key_type.__name__}") + + self.call_method("erase", key=key) + + def clear(self): + """ + Remove all elements. + + """ + self.call_method("clear") + + def _get(self, key): + if not is_valid_type(key, self._key_type): + raise ValueError( + f"Key has to be of type {self._key_type.__name__}") + + return self.call_method("get", key=key) + + def __getitem__(self, key): + return self._get(key) + + def __setitem__(self, key, value): + self.call_method("insert", key=key, object=value) + + def __delitem__(self, key): + self.remove(key) + + def keys(self): + return self.call_method("keys") + + def __iter__(self): + for k in self.keys(): yield k + + def items(self): + for k in self.keys(): yield k, self[k] + + @classmethod + def _restore_object(cls, so_callback, so_callback_args, state): + so = so_callback(*so_callback_args) + so.__setstate__(state) + return so + + def __reduce__(self): + so_callback, (so_name, so_bytestring) = super().__reduce__() + return (ScriptObjectMap._restore_object, + (so_callback, (so_name, so_bytestring), self.__getstate__())) + + def __getstate__(self): + return dict(self.items()) + + def __setstate__(self, params): + for key, val in params.items(): + self[key] = val + + # Map from script object names to their corresponding python classes _python_class_by_so_name = {} diff --git a/src/python/espressomd/system.pyx b/src/python/espressomd/system.pyx index 995debf502d..9ad952757f8 100644 --- a/src/python/espressomd/system.pyx +++ b/src/python/espressomd/system.pyx @@ -44,6 +44,7 @@ if LB_BOUNDARIES or LB_BOUNDARIES_GPU: from .comfixed import ComFixed from .utils cimport check_type_or_throw_except from .utils import handle_errors, array_locked +from .bond_breakage import BreakageSpecs IF VIRTUAL_SITES: from .virtual_sites import ActiveVirtualSitesHandle, VirtualSitesOff @@ -143,6 +144,8 @@ cdef class System: """:class:`espressomd.actors.Actors`""" analysis """:class:`espressomd.analyze.Analysis`""" + bond_breakage + """:class:`espressomd.bond_breakage.BreakageSpecs`""" galilei """:class:`espressomd.galilei.GalileiTransform`""" integrator @@ -182,6 +185,7 @@ cdef class System: self.auto_update_accumulators = AutoUpdateAccumulators() self.bonded_inter = interactions.BondedInteractions() self.cell_system = CellSystem() + self.bond_breakage = BreakageSpecs() IF COLLISION_DETECTION == 1: self.collision_detection = CollisionDetection() self.comfixed = ComFixed() @@ -226,6 +230,7 @@ cdef class System: odict['lbboundaries'] = System.__getattribute__( self, "lbboundaries") odict['thermostat'] = System.__getattribute__(self, "thermostat") + odict['bond_breakage'] = System.__getattribute__(self, "bond_breakage") IF COLLISION_DETECTION: odict['collision_detection'] = System.__getattribute__( self, "collision_detection") diff --git a/src/script_interface/CMakeLists.txt b/src/script_interface/CMakeLists.txt index 6c657ae04ac..cc2a9e1373d 100644 --- a/src/script_interface/CMakeLists.txt +++ b/src/script_interface/CMakeLists.txt @@ -4,6 +4,7 @@ add_library( GlobalContext.cpp ContextManager.cpp) add_subdirectory(accumulators) +add_subdirectory(bond_breakage) add_subdirectory(collision_detection) add_subdirectory(constraints) add_subdirectory(cluster_analysis) diff --git a/src/script_interface/ObjectMap.hpp b/src/script_interface/ObjectMap.hpp index 76c32ad05cb..d6af7da67df 100644 --- a/src/script_interface/ObjectMap.hpp +++ b/src/script_interface/ObjectMap.hpp @@ -124,6 +124,11 @@ class ObjectMap : public BaseType { return none; } + if (method == "get") { + auto key = get_value(parameters.at("key")); + return Variant{m_elements.at(key)}; + } + if (method == "get_map") { std::unordered_map ret; @@ -133,6 +138,14 @@ class ObjectMap : public BaseType { return ret; } + if (method == "keys") { + std::vector res; + for (auto const &kv : m_elements) { + res.push_back(kv.first); + } + return res; + } + if (method == "clear") { clear(); return none; diff --git a/src/script_interface/bond_breakage/BreakageSpec.hpp b/src/script_interface/bond_breakage/BreakageSpec.hpp new file mode 100644 index 00000000000..d44f9900816 --- /dev/null +++ b/src/script_interface/bond_breakage/BreakageSpec.hpp @@ -0,0 +1,69 @@ +/* + * Copyright (C) 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 "bond_breakage.hpp" + +#include "script_interface/ScriptInterface.hpp" + +#include + +namespace ScriptInterface { +namespace BondBreakage { + +class BreakageSpec : public AutoParameters { +public: + BreakageSpec() : m_breakage_spec(new ::BondBreakage::BreakageSpec) { + add_parameters({ + {"breakage_length", m_breakage_spec->breakage_length}, + {"action_type", + [this](const Variant &v) { + m_breakage_spec->action_type = ::BondBreakage::ActionType{ + m_breakage_str_to_enum.at(boost::get(v))}; + }, + [this]() { + return Variant( + m_breakage_enum_to_str.at(m_breakage_spec->action_type)); + }}, + }); + } + + std::shared_ptr<::BondBreakage::BreakageSpec> breakage_spec() const { + return m_breakage_spec; + } + +private: + std::shared_ptr<::BondBreakage::BreakageSpec> m_breakage_spec; + std::unordered_map<::BondBreakage::ActionType, std::string> + m_breakage_enum_to_str = { + {::BondBreakage::ActionType::NONE, "none"}, + {::BondBreakage::ActionType::DELETE_BOND, "revert_center_bond"}, + {::BondBreakage::ActionType::REVERT_BIND_AT_POINT_OF_COLLISION, + "revert_vs_bond"}}; + std::unordered_map + m_breakage_str_to_enum = { + {"none", ::BondBreakage::ActionType::NONE}, + {"revert_center_bond", ::BondBreakage::ActionType::DELETE_BOND}, + {"revert_vs_bond", + ::BondBreakage::ActionType::REVERT_BIND_AT_POINT_OF_COLLISION}}; +}; + +} // namespace BondBreakage +} // namespace ScriptInterface diff --git a/src/script_interface/bond_breakage/BreakageSpecs.hpp b/src/script_interface/bond_breakage/BreakageSpecs.hpp new file mode 100644 index 00000000000..214ab7df2ac --- /dev/null +++ b/src/script_interface/bond_breakage/BreakageSpecs.hpp @@ -0,0 +1,64 @@ +/* + * Copyright (C) 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 "BreakageSpec.hpp" + +#include "core/bond_breakage.hpp" + +#include "script_interface/ObjectMap.hpp" +#include "script_interface/ScriptInterface.hpp" + +#include +#include +#include + +namespace ScriptInterface { +namespace BondBreakage { +class BreakageSpecs : public ObjectMap { + using container_type = std::unordered_map>; + +public: + using key_type = typename container_type::key_type; + using mapped_type = typename container_type::mapped_type; + + key_type insert_in_core(mapped_type const &) override { + if (context()->is_head_node()) { + throw std::runtime_error( + "Inserting breakage spec without a bond type is not permitted."); + } + return {}; + } + void insert_in_core(key_type const &key, + mapped_type const &obj_ptr) override { + auto core_spec = obj_ptr->breakage_spec(); + ::BondBreakage::breakage_specs.insert({key, core_spec}); + } + void erase_in_core(key_type const &key) override { + ::BondBreakage::breakage_specs.erase(key); + } + +private: + // disable serialization: bond breakage has its own pickling logic + std::string get_internal_state() const override { return {}; } + void set_internal_state(std::string const &state) override {} +}; +} // namespace BondBreakage +} // namespace ScriptInterface diff --git a/src/script_interface/bond_breakage/CMakeLists.txt b/src/script_interface/bond_breakage/CMakeLists.txt new file mode 100644 index 00000000000..43c571ff9a2 --- /dev/null +++ b/src/script_interface/bond_breakage/CMakeLists.txt @@ -0,0 +1,2 @@ +target_sources(ScriptInterface + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/initialize.cpp) diff --git a/src/script_interface/bond_breakage/initialize.cpp b/src/script_interface/bond_breakage/initialize.cpp new file mode 100644 index 00000000000..380665fdd13 --- /dev/null +++ b/src/script_interface/bond_breakage/initialize.cpp @@ -0,0 +1,31 @@ +/* + * Copyright (C) 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 "initialize.hpp" +#include "BreakageSpec.hpp" +#include "BreakageSpecs.hpp" + +namespace ScriptInterface { +namespace BondBreakage { +void initialize(Utils::Factory *f) { + f->register_new("BondBreakage::BreakageSpec"); + f->register_new("BondBreakage::BreakageSpecs"); +} +} // namespace BondBreakage +} // namespace ScriptInterface diff --git a/src/script_interface/bond_breakage/initialize.hpp b/src/script_interface/bond_breakage/initialize.hpp new file mode 100644 index 00000000000..43a06dd7ef6 --- /dev/null +++ b/src/script_interface/bond_breakage/initialize.hpp @@ -0,0 +1,31 @@ +/* + * Copyright (C) 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 "script_interface/ObjectHandle.hpp" + +#include + +namespace ScriptInterface { +namespace BondBreakage { +void initialize(Utils::Factory *f); + +} // namespace BondBreakage +} // namespace ScriptInterface diff --git a/src/script_interface/initialize.cpp b/src/script_interface/initialize.cpp index b6627905e84..2a55eb1e528 100644 --- a/src/script_interface/initialize.cpp +++ b/src/script_interface/initialize.cpp @@ -31,6 +31,7 @@ #include "ComFixed.hpp" #include "CylindricalTransformationParameters.hpp" #include "accumulators/initialize.hpp" +#include "bond_breakage/initialize.hpp" #include "collision_detection/initialize.hpp" #include "interactions/initialize.hpp" #include "lbboundaries/initialize.hpp" @@ -47,6 +48,7 @@ void initialize(Utils::Factory *f) { Writer::initialize(f); #endif Accumulators::initialize(f); + BondBreakage::initialize(f); Observables::initialize(f); ClusterAnalysis::initialize(f); Interactions::initialize(f); diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 31f3c5be59f..be1a68c21e1 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -85,6 +85,7 @@ foreach( save_checkpoint_${TEST_COMBINATION}_${TEST_BINARY}) endforeach(TEST_BINARY) endforeach(TEST_COMBINATION) +python_test(FILE bond_breakage.py MAX_NUM_PROC 4) python_test(FILE cellsystem.py MAX_NUM_PROC 4) python_test(FILE tune_skin.py MAX_NUM_PROC 1) python_test(FILE constraint_homogeneous_magnetic_field.py MAX_NUM_PROC 4) diff --git a/testsuite/python/bond_breakage.py b/testsuite/python/bond_breakage.py new file mode 100644 index 00000000000..ba538cb81c9 --- /dev/null +++ b/testsuite/python/bond_breakage.py @@ -0,0 +1,289 @@ +# +# Copyright (C) 2013-2019 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 . +# +import unittest as ut +import unittest_decorators as utx +import espressomd +import numpy as np +import espressomd.interactions +from espressomd.bond_breakage import BreakageSpec +from espressomd.interactions import HarmonicBond + + +class BondBreakageCommon: + system = espressomd.System(box_l=[10] * 3) + system.cell_system.skin = 0.4 + system.time_step = 0.01 + system.min_global_cut = 2 + + +@utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE") +class BondBreakage(BondBreakageCommon, ut.TestCase): + + @classmethod + def setUpClass(cls): + pos1 = cls.system.box_l / 2 - 0.5 + pos2 = cls.system.box_l / 2 + 0.5 + cls.p1 = cls.system.part.add(pos=pos1) + cls.p2 = cls.system.part.add(pos=pos2) + + cls.p1v = cls.system.part.add(pos=pos1) + cls.p1v.vs_auto_relate_to(cls.p1) + + cls.p2v = cls.system.part.add(pos=pos2) + cls.p2v.vs_auto_relate_to(cls.p2) + + cls.h1 = HarmonicBond(k=1, r_0=0) + cls.h2 = HarmonicBond(k=1, r_0=0) + cls.system.bonded_inter.add(cls.h1) + cls.system.bonded_inter.add(cls.h2) + + @classmethod + def tearDownClass(cls): + cls.system.part.clear() + cls.system.bonded_inter.clear() + + def tearDown(self): + self.system.bond_breakage.clear() + + def test_00_interface(self): + self.assertEqual(len(self.system.bond_breakage), 0) + + spec2 = BreakageSpec( + breakage_length=1.2, + action_type="revert_center_bond") + spec4 = BreakageSpec(breakage_length=0.2, action_type="revert_vs_bond") + self.system.bond_breakage[2] = spec2 + self.system.bond_breakage[4] = spec4 + self.assertEqual(self.system.bond_breakage[2], spec2) + self.assertEqual(self.system.bond_breakage[4], spec4) + self.assertEqual(len(self.system.bond_breakage), 2) + self.assertEqual(sorted(self.system.bond_breakage.keys()), [2, 4]) + self.assertEqual( + sorted(self.system.bond_breakage.items()), + [(2, spec2), (4, spec4)]) + + self.system.bond_breakage.clear() + self.assertEqual(len(self.system.bond_breakage), 0) + self.assertEqual(self.system.bond_breakage.keys(), []) + with self.assertRaisesRegex(ValueError, "Key has to be of type int"): + self.system.bond_breakage[self.h1] + with self.assertRaisesRegex(ValueError, "Key has to be of type int"): + self.system.bond_breakage.remove(self.h1) + with self.assertRaisesRegex(RuntimeError, "Inserting breakage spec without a bond type is not permitted"): + self.system.bond_breakage.call_method("insert", object=spec2) + + def test_ignore(self): + system = self.system + + # Particles closer than cutoff + system.bond_breakage[self.h1._bond_id] = BreakageSpec( + breakage_length=2, action_type="revert_center_bond") + + self.p1.bonds = ((self.h1, self.p2)) + system.integrator.run(1) + self.assertEqual(self.p1.bonds, ((self.h1, self.p2.id),)) + + self.p2.bonds = [(self.h1, self.p1)] + system.integrator.run(1) + self.assertEqual(self.p2.bonds, ((self.h1, self.p1.id),)) + + # Different bond type + system.bond_breakage[self.h1._bond_id] = BreakageSpec( + breakage_length=0.2, action_type="revert_center_bond") + self.p1.bonds = [(self.h2, self.p2)] + self.p2.bonds = [(self.h2, self.p1)] + system.integrator.run(1) + self.assertEqual(self.p1.bonds, ((self.h2, self.p2.id),)) + self.assertEqual(self.p2.bonds, ((self.h2, self.p1.id),)) + + def test_delete_bond(self): + system = self.system + + # Particles closer than cutoff + system.bond_breakage[self.h1._bond_id] = BreakageSpec( + breakage_length=0, action_type="revert_center_bond") + + self.p1.bonds = [(self.h1, self.p2)] + system.integrator.run(1) + self.assertEqual(self.p1.bonds, ()) + + self.p2.bonds = [(self.h1, self.p1)] + system.integrator.run(1) + self.assertEqual(self.p2.bonds, ()) + + def test_revert_bind_at_point_of_collision(self): + system = self.system + + # Particles closer than cutoff + system.bond_breakage[self.h1._bond_id] = BreakageSpec( + breakage_length=0.5, action_type="revert_vs_bond") + + self.p1.bonds = [(self.h2, self.p2)] + self.p1v.bonds = [(self.h1, self.p2v)] + system.integrator.run(1) + self.assertEqual(self.p1v.bonds, ()) + self.assertEqual(self.p1.bonds, ()) + + self.p2.bonds = [(self.h2, self.p1)] + self.p1v.bonds = [(self.h1, self.p2v)] + system.integrator.run(1) + self.assertEqual(self.p1.bonds, ()) + self.assertEqual(self.p1v.bonds, ()) + + def test_exceptions(self): + system = self.system + + # Particles closer than cutoff + system.bond_breakage[self.h2._bond_id] = BreakageSpec( + breakage_length=0.5, action_type="revert_vs_bond") + + self.p1.bonds = [(self.h2, self.p2)] + self.p1v.bonds = [(self.h1, self.p2v)] + with self.assertRaisesRegex(Exception, "The REVERT_BIND_AT_POINT_OF_COLLISION bond breakage action has to be configured for the bond on the virtual site"): + system.integrator.run(1) + + +@utx.skipIfMissingFeatures("LENNARD_JONES") +class NetworkBreakage(BondBreakageCommon, ut.TestCase): + + @classmethod + def setUpClass(cls): + cls.system.box_l = 3 * [20] + cls.system.min_global_cut = 0.6 + cls.system.time_step = 0.01 + cls.system.cell_system.skin = 0.4 + + def count_bonds(self, pairs): + bonds_count = 0 + for pair in pairs: + for bond in self.system.part.by_id(pair[0]).bonds: + if bond[1] == pair[1]: + bonds_count += 1 + for bond in self.system.part.by_id(pair[1]).bonds: + if bond[1] == pair[0]: + bonds_count += 1 + return bonds_count + + def setUp(self): + + box_vol = self.system.box_l[0]**3. + phi = 0.4 + + r = 1. + solid_vol = phi * box_vol + part_vol = 4 / 3 * np.pi * r**3 + part_num = int(solid_vol / part_vol) + + np.random.seed(seed=678) + for i in range(part_num): + pos = np.random.rand(3) * self.system.box_l[0] + self.system.part.add(pos=pos) + + self.system.non_bonded_inter[0, 0].lennard_jones.set_params( + sigma=1., epsilon=1., cutoff=2**(1 / 6), shift='auto') + self.system.integrator.set_steepest_descent(f_max=0, + gamma=0.1, + max_displacement=0.1) + self.system.integrator.run(100) + self.system.integrator.set_vv() + + for i in range(part_num): + self.system.part.by_id(i).fix = [1, 1, 1] + + self.system.thermostat.set_langevin(kT=0.0, gamma=1.0, seed=41) + + def tearDown(self): + self.system.part.clear() + self.system.bonded_inter.clear() + self.system.thermostat.turn_off() + + def test_center_bonds(self): + + harm = espressomd.interactions.HarmonicBond(k=1.0, r_0=0.0, r_cut=5) + self.system.bonded_inter.add(harm) + + 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.integrator.run(1) + + self.system.collision_detection.set_params(mode="off") + self.system.bond_breakage[harm._bond_id] = BreakageSpec( + breakage_length=crit, action_type="revert_center_bond") + self.system.integrator.run(1) + + bonds_dist = 0 + pairs = self.system.cell_system.get_pairs(crit, types=[0]) + for pair in pairs: + dist = self.system.distance( + self.system.part.by_id(pair[0]), + self.system.part.by_id(pair[1])) + if dist <= crit: + bonds_dist += 1 + + bonds_count = self.count_bonds(pairs) + np.testing.assert_equal(bonds_dist, bonds_count) + + @utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE") + def test_vs_bonds(self): + + harm = espressomd.interactions.HarmonicBond(k=1.0, r_0=0.0, r_cut=5) + virt = espressomd.interactions.Virtual() + self.system.bonded_inter.add(harm) + self.system.bonded_inter.add(virt) + + 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.integrator.run(1) + + self.system.collision_detection.set_params(mode="off") + self.system.bond_breakage[harm._bond_id] = BreakageSpec( + breakage_length=crit_vs, action_type="revert_vs_bond") + self.system.integrator.run(1) + + bonds_dist = 0 + pairs = self.system.cell_system.get_pairs( + 2**(1 / 6) * 2 / 3, types=[1]) + + for pair in pairs: + r1 = self.system.part.by_id(pair[0]).vs_relative[0] + r2 = self.system.part.by_id(pair[1]).vs_relative[0] + dist = self.system.distance( + self.system.part.by_id(r1), + self.system.part.by_id(r2)) + dist_vs = self.system.distance( + self.system.part.by_id(pair[0]), + self.system.part.by_id(pair[1])) + if dist_vs <= crit_vs: + if dist <= crit: + if dist > 0.0: + bonds_dist += 1 + + bonds_count = self.count_bonds(pairs) + + np.testing.assert_equal(bonds_dist, bonds_count) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 01db89bb88a..97862467dca 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -31,6 +31,7 @@ import espressomd.lbboundaries import espressomd.shapes import espressomd.constraints +import espressomd.bond_breakage modes = {x for mode in set("@TEST_COMBINATION@".upper().split('-')) for x in [mode, mode.split('.')[0]]} @@ -231,6 +232,9 @@ ibm_triel_bond = espressomd.interactions.IBM_Triel( ind1=p1.id, ind2=p2.id, ind3=p3.id, k1=1.1, k2=1.2, maxDist=1.6, elasticLaw="NeoHookean") +break_spec = espressomd.bond_breakage.BreakageSpec( + breakage_length=5., action_type="revert_center_bond") +system.bond_breakage[strong_harmonic_bond._bond_id] = break_spec checkpoint.register("system") checkpoint.register("acc_mean_variance") @@ -239,6 +243,7 @@ checkpoint.register("ibm_volcons_bond") checkpoint.register("ibm_tribend_bond") checkpoint.register("ibm_triel_bond") +checkpoint.register("break_spec") # calculate forces system.integrator.run(0) diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 9e36f088969..5eefa05cc53 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -341,6 +341,17 @@ def test_bonded_inter(self): bond_ids = system.bonded_inter.call_method('get_bond_ids') self.assertEqual(len(bond_ids), len(system.bonded_inter)) + def test_bond_breakage_specs(self): + # check the ObjectHandle was correctly initialized (including MPI) + spec_ids = list(system.bond_breakage.keys()) + self.assertEqual(len(spec_ids), 1) + cpt_spec = system.bond_breakage[spec_ids[0]] + self.assertAlmostEqual( + break_spec.breakage_length, + cpt_spec.breakage_length, + delta=1e-10) + self.assertEqual(break_spec.action_type, cpt_spec.action_type) + @ut.skipIf('THERM.LB' in modes, 'LB thermostat in modes') @utx.skipIfMissingFeatures(['ELECTROSTATICS', 'MASS', 'ROTATION']) def test_drude_helpers(self):