ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
bond_breakage.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
22
24#include "communication.hpp"
25#include "errorhandling.hpp"
26#include "system/System.hpp"
27
30
31#include <boost/mpi.hpp>
32#include <boost/serialization/access.hpp>
33#include <boost/variant.hpp>
34
35#include <cassert>
36#include <memory>
37#include <unordered_set>
38#include <utility>
39#include <vector>
40
41namespace BondBreakage {
42
43// Variant holding any of the actions
44using Action = boost::variant<DeleteBond, DeleteAngleBond, DeleteAllBonds>;
45
46// Set of actions
47using ActionSet = std::unordered_set<Action>;
48
49/** Add a particle+bond combination to the breakage queue */
50void BondBreakage::queue_breakage(int particle_id,
51 BondPartners const &bond_partners,
52 int bond_type) {
53 m_queue.emplace_back(QueueEntry{particle_id, bond_partners, bond_type});
54}
55
56/** @brief Gathers combined queue from all mpi ranks */
57static auto gather_global_queue(Queue const &local_queue) {
58 Queue res = local_queue;
59 if (comm_cart.size() > 1) {
61 boost::mpi::broadcast(comm_cart, res, 0);
62 }
63 return res;
64}
65
66/** @brief Constructs the actions to take for a breakage queue entry */
67static ActionSet actions_for_breakage(CellStructure const &cell_structure,
68 QueueEntry const &e,
69 BreakageSpec const &spec) {
70 auto is_angle_bond = [](auto const &bond_partners) {
71 return bond_partners[1];
72 }; // optional for second partner engaged
73
74 // Handle different action types
76 if (is_angle_bond(e.bond_partners)) {
78 {{*(e.bond_partners[0]), *(e.bond_partners[1])}},
79 e.bond_type}};
80 }
81 return {DeleteBond{e.particle_id, *(e.bond_partners[0]), e.bond_type}};
82 }
83#ifdef VIRTUAL_SITES_RELATIVE
84 // revert bind at point of collision for pair bonds
86 not is_angle_bond(e.bond_partners)) {
87 // We need to find the base particles for the two virtual sites
88 // between which the bond broke.
89 auto p1 = cell_structure.get_local_particle(e.particle_id);
90 auto p2 = cell_structure.get_local_particle(*(e.bond_partners[0]));
91 if (p1 and p2) {
92 if (not p1->is_virtual() or not p2->is_virtual()) {
93 runtimeErrorMsg() << "The REVERT_BIND_AT_POINT_OF_COLLISION bond "
94 "breakage action has to be configured for the "
95 "bond on the virtual site. Encountered a particle "
96 "that is not virtual.";
97 return {};
98 }
99
100 return {
101 // Bond between virtual sites
103 // Bond between base particles. We do not know, on which of these
104 // the bond is defined, since bonds are stored only on one partner
105 DeleteAllBonds{p1->vs_relative().to_particle_id,
106 p2->vs_relative().to_particle_id},
107 DeleteAllBonds{p2->vs_relative().to_particle_id,
108 p1->vs_relative().to_particle_id},
109 };
110 }
111 } else {
112 // revert bind at point of collision for angle bonds
113 auto vs = cell_structure.get_local_particle(e.particle_id);
114 auto p1 = cell_structure.get_local_particle(*(e.bond_partners[0]));
115 auto p2 = cell_structure.get_local_particle(*(e.bond_partners[1]));
116 if (p1 and p2) {
117 if (not vs->is_virtual()) {
118 runtimeErrorMsg() << "The REVERT_BIND_AT_POINT_OF_COLLISION bond "
119 "breakage action has to be configured for the "
120 "bond on the virtual site. Encountered a particle "
121 "that is not virtual.";
122 return {};
123 }
124
125 return {// Angle bond on the virtual site
126 DeleteAngleBond{e.particle_id, {p1->id(), p2->id()}, e.bond_type},
127 // Bond between base particles. We do not know, on which of these
128 // the bond is defined, since bonds are stored only on one partner
129 DeleteAllBonds{p1->id(), p2->id()},
130 DeleteAllBonds{p2->id(), p1->id()}};
131 }
132 }
133#endif // VIRTUAL_SITES_RELATIVE
134 return {};
135}
136
137/**
138 * @brief Delete specific bond.
139 */
140static void remove_bond(Particle &p, BondView const &view) {
141 auto &bond_list = p.bonds();
142 auto it = std::find(bond_list.begin(), bond_list.end(), view);
143 if (it != bond_list.end()) {
144 bond_list.erase(it);
145 }
146}
147
148/**
149 * @brief Delete pair bonds to a specific partner
150 */
151static void remove_pair_bonds_to(Particle &p, int other_pid) {
152 std::vector<std::pair<int, int>> to_delete;
153 for (auto b : p.bonds()) {
154 if (b.partner_ids().size() == 1 and b.partner_ids()[0] == other_pid)
155 to_delete.emplace_back(b.bond_id(), other_pid);
156 }
157 for (auto const &b : to_delete) {
158 remove_bond(p, BondView(b.first, {&b.second, 1}));
159 }
160}
161
162// Handler for the different delete events
163class execute : public boost::static_visitor<> {
164 CellStructure &cell_structure;
165
166public:
167 explicit execute(CellStructure &cell_structure)
168 : cell_structure{cell_structure} {}
169
170 void operator()(DeleteBond const &d) const {
171 if (auto p = cell_structure.get_local_particle(d.particle_id)) {
172 remove_bond(*p, BondView(d.bond_type, {&d.bond_partner_id, 1}));
173 }
174 }
175 void operator()(DeleteAngleBond const &d) const {
176 if (auto p = cell_structure.get_local_particle(d.particle_id)) {
177 remove_bond(*p, BondView(d.bond_type, {&d.bond_partner_id[0], 2}));
178 }
179 }
180 void operator()(DeleteAllBonds const &d) const {
181 if (auto p = cell_structure.get_local_particle(d.particle_id_1)) {
183 }
184 }
185};
186
187void BondBreakage::process_queue_impl(System::System &system) {
188 auto global_queue = gather_global_queue(m_queue);
189 auto &cell_structure = *system.cell_structure;
190
191 // Construct delete actions from breakage queue
192 ActionSet actions = {};
193 for (auto const &e : global_queue) {
194 // Retrieve relevant breakage spec
195 assert(breakage_specs.count(e.bond_type) != 0);
196 auto const &spec = breakage_specs.at(e.bond_type);
197 actions.merge(actions_for_breakage(cell_structure, e, *spec));
198 }
199
200 // Execute actions
201 for (auto const &a : actions) {
202 boost::apply_visitor(execute(cell_structure), a);
203 system.on_particle_change();
204 }
205}
206
207} // namespace BondBreakage
execute(CellStructure &cell_structure)
void operator()(DeleteAngleBond const &d) const
void operator()(DeleteBond const &d) const
void operator()(DeleteAllBonds const &d) const
Immutable view on a bond.
Definition BondList.hpp:44
Main system class.
void on_particle_change()
Called every time a particle property changes.
std::shared_ptr< CellStructure > cell_structure
boost::mpi::communicator comm_cart
The communicator.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
static void remove_pair_bonds_to(Particle &p, int other_pid)
Delete pair bonds to a specific partner.
static void remove_bond(Particle &p, BondView const &view)
Delete specific bond.
static ActionSet actions_for_breakage(CellStructure const &cell_structure, QueueEntry const &e, BreakageSpec const &spec)
Constructs the actions to take for a breakage queue entry.
static auto gather_global_queue(Queue const &local_queue)
Gathers combined queue from all mpi ranks.
boost::variant< DeleteBond, DeleteAngleBond, DeleteAllBonds > Action
std::array< std::optional< int >, 2 > BondPartners
Stores one or two bond parnters for pair/angle bonds.
std::vector< QueueEntry > Queue
Record bonds broken during a time step.
std::unordered_set< Action > ActionSet
void gather_buffer(std::vector< T, Allocator > &buffer, boost::mpi::communicator const &comm, int root=0)
Gather buffer with different size on each node.
Describes a cell structure / cell system.
Particle * get_local_particle(int id)
Get a local particle by id.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & bonds() const
Definition Particle.hpp:428