Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Breaking bonds #4456

Merged
merged 5 commits into from
Feb 24, 2022
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Implement bond breakage
RudolfWeeber authored and jngrad committed Feb 24, 2022
commit af9055a6984226d97fb3afce78b2ff08961be877
1 change: 1 addition & 0 deletions src/core/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
set(EspressoCore_SRC
accumulators.cpp
bond_breakage.cpp
bond_error.cpp
cells.cpp
collision.cpp
232 changes: 232 additions & 0 deletions src/core/bond_breakage.cpp
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/
#include "bond_breakage.hpp"
#include "cells.hpp"
#include "communication.hpp"
#include "errorhandling.hpp"
#include "particle_data.hpp"

#include <utils/mpi/gather_buffer.hpp>

#include <boost/functional/hash.hpp>
#include <boost/variant.hpp>

#include <cassert>
#include <cstddef>
#include <memory>
#include <unordered_set>
#include <vector>

namespace BondBreakage {

// Bond breakage specifications
std::unordered_map<int, std::shared_ptr<BreakageSpec>> 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<BondBreakage::DeleteBond> {
std::size_t operator()(BondBreakage::DeleteBond const &t) const noexcept {
return t.hash_value();
}
};
template <> struct hash<BondBreakage::DeleteAllBonds> {
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<DeleteBond, DeleteAllBonds>;

// Set of actions
using ActionSet = std::unordered_set<Action>;

// 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 <typename Archive>
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<QueueEntry>;
Queue queue;

/** @brief Retrieve breakage specification for the bond type */
boost::optional<BreakageSpec> 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
50 changes: 50 additions & 0 deletions src/core/bond_breakage.hpp
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

#pragma once

#include <memory>
#include <unordered_map>

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<int, std::shared_ptr<BreakageSpec>> 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
3 changes: 2 additions & 1 deletion src/core/forces.cpp
Original file line number Diff line number Diff line change
@@ -27,6 +27,7 @@

#include "EspressoSystemInterface.hpp"

#include "bond_breakage.hpp"
#include "cells.hpp"
#include "collision.hpp"
#include "comfixed_global.hpp"
@@ -149,7 +150,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
10 changes: 10 additions & 0 deletions src/core/forces_inline.hpp
Original file line number Diff line number Diff line change
@@ -25,6 +25,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"
@@ -403,6 +404,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<Particle *> 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)) {
2 changes: 2 additions & 0 deletions src/core/integrate.cpp
Original file line number Diff line number Diff line change
@@ -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++;
30 changes: 30 additions & 0 deletions src/core/particle_data.cpp
Original file line number Diff line number Diff line change
@@ -194,6 +194,28 @@ struct RemoveBond {
}
};

/**
* @brief Delete pair bonds to a specific partner
*/
struct RemovePairBondsTo {
int other_pid;

void operator()(Particle &p) const {
using Bond = std::vector<int>;
std::vector<Bond> 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<class Archive>
void serialize(Archive &ar, long int) {
ar & other_pid;
}
};

/**
* @brief Delete all bonds.
@@ -323,6 +345,14 @@ struct UpdateVisitor : public boost::static_visitor<void> {
};
} // namespace

void local_remove_bond(Particle &p, std::vector<int> 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{};
15 changes: 13 additions & 2 deletions src/core/particle_data.hpp
Original file line number Diff line number Diff line change
@@ -299,10 +299,21 @@ void delete_particle_bond(int part, Utils::Span<const int> bond);
*/
void delete_particle_bonds(int part);

/** @brief Removs the specified bond from the particle
* @param bond The bond in the form
* <tt>{bond_id, partner_1, partner_2, ...}</tt>
*/
void local_remove_bond(Particle &p, std::vector<int> 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<const int> bond);

30 changes: 30 additions & 0 deletions src/python/espressomd/bond_breakage.py
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
#
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
83 changes: 82 additions & 1 deletion src/python/espressomd/script_interface.pyx
Original file line number Diff line number Diff line change
@@ -15,7 +15,7 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import numpy as np
from .utils import to_char_pointer, to_str, handle_errors
from .utils import to_char_pointer, to_str, handle_errors, is_valid_type
from .utils cimport Vector3d, make_array_locked
cimport cpython.object

@@ -377,6 +377,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 = {}

5 changes: 5 additions & 0 deletions src/python/espressomd/system.pyx
Original file line number Diff line number Diff line change
@@ -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")
1 change: 1 addition & 0 deletions src/script_interface/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
13 changes: 13 additions & 0 deletions src/script_interface/ObjectMap.hpp
Original file line number Diff line number Diff line change
@@ -124,6 +124,11 @@ class ObjectMap : public BaseType {
return none;
}

if (method == "get") {
auto key = get_value<KeyType>(parameters.at("key"));
return Variant{m_elements.at(key)};
}

if (method == "get_map") {
std::unordered_map<KeyType, Variant> ret;

@@ -133,6 +138,14 @@ class ObjectMap : public BaseType {
return ret;
}

if (method == "keys") {
std::vector<Variant> res;
for (auto const &kv : m_elements) {
res.push_back(kv.first);
}
return res;
}

if (method == "clear") {
clear();
return none;
54 changes: 54 additions & 0 deletions src/script_interface/bond_breakage/BreakageSpec.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
/*
* 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 <http://www.gnu.org/licenses/>.
*/

#pragma once

#include "bond_breakage.hpp"

#include "script_interface/ScriptInterface.hpp"

#include <memory>

namespace ScriptInterface {
namespace BondBreakage {

class BreakageSpec : public AutoParameters<BreakageSpec> {
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(boost::get<int>(v));
},
[this]() { return Variant(int(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;
};

} // namespace BondBreakage
} // namespace ScriptInterface
64 changes: 64 additions & 0 deletions src/script_interface/bond_breakage/BreakageSpecs.hpp
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

#pragma once

#include "BreakageSpec.hpp"

#include "core/bond_breakage.hpp"

#include "script_interface/ObjectMap.hpp"
#include "script_interface/ScriptInterface.hpp"

#include <memory>
#include <stdexcept>
#include <unordered_map>

namespace ScriptInterface {
namespace BondBreakage {
class BreakageSpecs : public ObjectMap<BreakageSpec> {
using container_type = std::unordered_map<int, std::shared_ptr<BreakageSpec>>;

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
2 changes: 2 additions & 0 deletions src/script_interface/bond_breakage/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
target_sources(ScriptInterface
PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/initialize.cpp)
31 changes: 31 additions & 0 deletions src/script_interface/bond_breakage/initialize.cpp
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

#include "initialize.hpp"
#include "BreakageSpec.hpp"
#include "BreakageSpecs.hpp"

namespace ScriptInterface {
namespace BondBreakage {
void initialize(Utils::Factory<ObjectHandle> *f) {
f->register_new<BreakageSpec>("BondBreakage::BreakageSpec");
f->register_new<BreakageSpecs>("BondBreakage::BreakageSpecs");
}
} // namespace BondBreakage
} // namespace ScriptInterface
31 changes: 31 additions & 0 deletions src/script_interface/bond_breakage/initialize.hpp
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

#pragma once

#include "script_interface/ObjectHandle.hpp"

#include <utils/Factory.hpp>

namespace ScriptInterface {
namespace BondBreakage {
void initialize(Utils::Factory<ObjectHandle> *f);

} // namespace BondBreakage
} // namespace ScriptInterface
2 changes: 2 additions & 0 deletions src/script_interface/initialize.cpp
Original file line number Diff line number Diff line change
@@ -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"
@@ -46,6 +47,7 @@ void initialize(Utils::Factory<ObjectHandle> *f) {
Writer::initialize(f);
#endif
Accumulators::initialize(f);
BondBreakage::initialize(f);
Observables::initialize(f);
ClusterAnalysis::initialize(f);
Interactions::initialize(f);
1 change: 1 addition & 0 deletions testsuite/python/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -83,6 +83,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)
152 changes: 152 additions & 0 deletions testsuite/python/bond_breakage.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
#
# 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 <http://www.gnu.org/licenses/>.
#
import unittest as ut
import unittest_decorators as utx
import espressomd
from espressomd.bond_breakage import BreakageSpec
from espressomd.interactions import HarmonicBond


@utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE")
class BondBreakage(ut.TestCase):

@classmethod
def setUpClass(cls):
system = cls.system = espressomd.System(box_l=[10] * 3)
system.cell_system.skin = 0.4
system.time_step = 0.01
system.min_global_cut = 2

pos1 = system.box_l / 2 - 0.5
pos2 = system.box_l / 2 + 0.5
cls.p1 = system.part.add(pos=pos1)
cls.p2 = system.part.add(pos=pos2)

cls.p1v = system.part.add(pos=pos1)
cls.p1v.vs_auto_relate_to(cls.p1)

cls.p2v = 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)
system.bonded_inter.add(cls.h1)
system.bonded_inter.add(cls.h2)

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=1)
spec4 = BreakageSpec(breakage_length=0.2, action_type=2)
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=1)

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=1)
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=1)

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=2)

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=2)

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)


if __name__ == "__main__":
ut.main()
5 changes: 5 additions & 0 deletions testsuite/python/save_checkpoint.py
Original file line number Diff line number Diff line change
@@ -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=2)
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)
11 changes: 11 additions & 0 deletions testsuite/python/test_checkpoint.py
Original file line number Diff line number Diff line change
@@ -339,6 +339,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):