ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
core/reaction_methods/ReactionAlgorithm.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-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
20#include "config/config.hpp"
21
22#include "reaction_methods/ReactionAlgorithm.hpp"
23
24#include "BoxGeometry.hpp"
25#include "Observable_stat.hpp"
27#include "cells.hpp"
28#include "particle_node.hpp"
29#include "system/System.hpp"
30
31#include <utils/Vector.hpp>
32#include <utils/contains.hpp>
33
34#include <boost/mpi/collectives/all_reduce.hpp>
35#include <boost/mpi/collectives/broadcast.hpp>
36#include <boost/serialization/serialization.hpp>
37
38#include <algorithm>
39#include <cassert>
40#include <cmath>
41#include <functional>
42#include <iterator>
43#include <limits>
44#include <map>
45#include <numbers>
46#include <optional>
47#include <stdexcept>
48#include <string>
49#include <tuple>
50#include <utility>
51#include <vector>
52
53namespace boost::serialization {
54template <typename Archive, typename... Types>
55void serialize(Archive &ar, std::tuple<Types...> &tuple, const unsigned int) {
56 std::apply([&](auto &...item) { ((ar & item), ...); }, tuple);
57}
58} // namespace boost::serialization
59
60namespace ReactionMethods {
61
62/**
63 * Adds a reaction to the reaction system
64 */
66 std::shared_ptr<SingleReaction> const &new_reaction) {
67
68 // make ESPResSo count the particle numbers which take part in the reactions
69 for (int reactant_type : new_reaction->reactant_types)
70 init_type_map(reactant_type);
71 for (int product_type : new_reaction->product_types)
72 init_type_map(product_type);
73
75
76 reactions.push_back(new_reaction);
77}
78
79/**
80 * @details This method tries to keep the cell system overhead to a minimum.
81 * Event callbacks are only called once after all particles are updated,
82 * except for particle deletion (the cell structure is still reinitialized
83 * after each deletion).
84 */
86 auto &system = System::get_system();
87 auto const &old_state = get_old_system_state();
88 auto const &box_geo = *system.box_geo;
89 // restore the properties of changed and hidden particles
90 for (auto const &state : {old_state.changed, old_state.hidden}) {
91 for (auto const &[p_id, p_type] : state) {
93 if (auto p = get_local_particle(p_id)) {
94 p->type() = p_type;
95#ifdef ELECTROSTATICS
96 p->q() = charges_of_types.at(p_type);
97#endif
98 }
99 }
100 }
101 // delete created particles
102 for (auto const p_id : old_state.created) {
103 delete_particle(p_id);
104 }
105 // restore original positions and velocities
106 for (auto const &[p_id, pos, vel] : old_state.moved) {
107 if (auto p = get_local_particle(p_id)) {
108 p->v() = vel;
109 auto folded_pos = pos;
110 auto image_box = Utils::Vector3i{};
111 box_geo.fold_position(folded_pos, image_box);
112 p->pos() = folded_pos;
113 p->image_box() = image_box;
114 }
115 }
116 if (not old_state.moved.empty()) {
117 auto const &system = System::get_system();
118 auto &cell_structure = *system.cell_structure;
119 cell_structure.set_resort_particles(Cells::RESORT_GLOBAL);
120 }
121 system.on_particle_change();
123}
124
125/**
126 * Automatically sets the volume which is used by the reaction ensemble to the
127 * volume of a cuboid box.
128 */
130 auto const &box_geo = *System::get_system().box_geo;
131 volume = box_geo.volume();
132}
133
134/**
135 * Checks whether all particles exist for the provided reaction.
136 */
138 SingleReaction const &current_reaction) const {
139 bool enough_particles = true;
140 for (int i = 0; i < current_reaction.reactant_types.size(); i++) {
141 int current_number =
143 if (current_number < current_reaction.reactant_coefficients[i]) {
144 enough_particles = false;
145 break;
146 }
147 }
148 return enough_particles;
149}
150
152 ParticleChanges &bookkeeping) {
153 // create or hide particles of types with corresponding types in reaction
154 auto const n_product_types = reaction.product_types.size();
155 auto const n_reactant_types = reaction.reactant_types.size();
156 auto const get_random_p_id_of_type = [this](int type) {
157 auto const random_index = i_random(number_of_particles_with_type(type));
158 return get_random_p_id(type, random_index);
159 };
160 auto only_local_changes = true;
161 for (int i = 0; i < std::min(n_product_types, n_reactant_types); i++) {
162 auto const n_product_coef = reaction.product_coefficients[i];
163 auto const n_reactant_coef = reaction.reactant_coefficients[i];
164 // change std::min(reactant_coefficients(i),product_coefficients(i)) many
165 // particles of reactant_types(i) to product_types(i)
166 auto const old_type = reaction.reactant_types[i];
167 auto const new_type = reaction.product_types[i];
168#ifdef ELECTROSTATICS
169 if (charges_of_types.at(new_type) != charges_of_types.at(old_type)) {
170 only_local_changes = false;
171 }
172#endif
173 for (int j = 0; j < std::min(n_product_coef, n_reactant_coef); j++) {
174 auto const p_id = get_random_p_id_of_type(old_type);
175 on_particle_type_change(p_id, old_type, new_type);
176 if (auto p = get_local_particle(p_id)) {
177 p->type() = new_type;
178#ifdef ELECTROSTATICS
179 p->q() = charges_of_types.at(new_type);
180#endif
181 }
182 bookkeeping.changed.emplace_back(p_id, old_type);
183 }
184 // create product_coefficients(i)-reactant_coefficients(i) many product
185 // particles iff product_coefficients(i)-reactant_coefficients(i)>0,
186 // iff product_coefficients(i)-reactant_coefficients(i)<0, hide this number
187 // of reactant particles
188 auto const delta_n = n_product_coef - n_reactant_coef;
189 if (delta_n > 0) {
190 auto const type = reaction.product_types[i];
191 for (int j = 0; j < delta_n; j++) {
192 auto const p_id = create_particle(type);
193 check_exclusion_range(p_id, type);
194 bookkeeping.created.emplace_back(p_id);
195 }
196 only_local_changes = false;
197 } else if (delta_n < 0) {
198 auto const type = reaction.reactant_types[i];
199 for (int j = 0; j < -delta_n; j++) {
200 auto const p_id = get_random_p_id_of_type(type);
201 bookkeeping.hidden.emplace_back(p_id, type);
202 check_exclusion_range(p_id, type);
203 hide_particle(p_id, type);
204 }
205 only_local_changes = false;
206 }
207 }
208 // create or hide particles of types with noncorresponding replacement types
209 for (auto i = std::min(n_product_types, n_reactant_types);
210 i < std::max(n_product_types, n_reactant_types); i++) {
211 if (n_product_types < n_reactant_types) {
212 // hide superfluous reactant_types particles
213 auto const type = reaction.reactant_types[i];
214 for (int j = 0; j < reaction.reactant_coefficients[i]; j++) {
215 auto const p_id = get_random_p_id_of_type(type);
216 bookkeeping.hidden.emplace_back(p_id, type);
217 check_exclusion_range(p_id, type);
218 hide_particle(p_id, type);
219 }
220 } else {
221 // create additional product_types particles
222 auto const type = reaction.product_types[i];
223 for (int j = 0; j < reaction.product_coefficients[i]; j++) {
224 auto const p_id = create_particle(type);
225 check_exclusion_range(p_id, type);
226 bookkeeping.created.emplace_back(p_id);
227 }
228 }
229 }
230 // determine which fine-grained event to trigger
231 if (n_product_types == n_reactant_types and only_local_changes) {
233 } else {
235 }
236}
237
238std::unordered_map<int, int>
239ReactionAlgorithm::get_particle_numbers(SingleReaction const &reaction) const {
240 std::unordered_map<int, int> particle_numbers;
241 // reactants
242 for (int type : reaction.reactant_types) {
243 particle_numbers[type] = ::number_of_particles_with_type(type);
244 }
245 // products
246 for (int type : reaction.product_types) {
247 particle_numbers[type] = ::number_of_particles_with_type(type);
248 }
249 return particle_numbers;
250}
251
252std::optional<double>
254 auto &reaction = *reactions[reaction_id];
255 reaction.tried_moves++;
257 if (!all_reactant_particles_exist(reaction)) {
258 // make sure that no incomplete reaction is performed -> only need to
259 // consider rollback of complete reactions
260 return {};
261 }
262 auto &bookkeeping = make_new_system_state();
263 bookkeeping.reaction_id = reaction_id;
264 bookkeeping.old_particle_numbers = get_particle_numbers(reaction);
265 make_reaction_attempt(reaction, bookkeeping);
266 auto E_pot_new = std::numeric_limits<double>::max();
268 E_pot_new = calculate_potential_energy();
269 }
270 return {E_pot_new};
271}
272
274 double bf,
275 double E_pot_old,
276 double E_pot_new) {
277 auto constexpr exp_min = -708.4; // for IEEE-compatible double
278 auto const exponent = -(E_pot_new - E_pot_old) / kT;
279 auto const exponential = (exponent < exp_min) ? 0. : std::exp(exponent);
280 auto &reaction = *reactions[reaction_id];
281 reaction.accumulator_potential_energy_difference_exponential(
282 std::vector<double>{exponential});
283 if (get_random_uniform_number() >= bf) {
284 // reject trial move: restore previous state, energy is unchanged
286 return E_pot_old;
287 }
288 // accept trial move: delete hidden particles and return new system energy
289 for (auto const &[p_id, p_type] : get_old_system_state().hidden) {
290 delete_particle(p_id);
291 }
292 reaction.accepted_moves++;
294 return E_pot_new;
295}
296
297/**
298 * Hides a particle from short ranged interactions and from the electrostatic
299 * interaction. Additional hiding from interactions would need to be implemented
300 * here.
301 *
302 * Removing the particle charge and changing its type to a non existing one
303 * deactivates all interactions with other particles, as if the particle was
304 * inexistent (currently only type-based interactions are switched off, as well
305 * as the electrostatic interaction).
306 * This function does not break bonds for simple reactions, as long as there
307 * are no reactions like 2A -->B where one of the reacting A particles occurs
308 * in the polymer (think of bond breakages if the monomer in the polymer gets
309 * deleted in the reaction). This constraint is not of fundamental reason, but
310 * there would be a need for a rule for such "collision" reactions (a reaction
311 * like the one above).
312 */
313void ReactionAlgorithm::hide_particle(int p_id, int p_type) const {
315 if (auto p = get_local_particle(p_id)) {
316 p->type() = non_interacting_type;
317#ifdef ELECTROSTATICS
318 p->q() = 0.;
319#endif
320 }
321}
322
323/**
324 * Check if the inserted particle is too close to neighboring particles.
325 */
326void ReactionAlgorithm::check_exclusion_range(int p_id, int p_type) {
327
328 /* Check the exclusion radius of the inserted particle */
329 if (exclusion_radius_per_type.count(p_type) != 0) {
330 if (exclusion_radius_per_type[p_type] == 0.) {
331 return;
332 }
333 }
334
335 auto p1_ptr = get_real_particle(p_id);
336
337 std::vector<int> particle_ids;
339 auto all_ids = get_particle_ids_parallel();
340 /* remove the inserted particle id */
341 std::erase(all_ids, p_id);
342 particle_ids = all_ids;
343 } else {
344 auto &system = System::get_system();
345 system.on_observable_calc();
346 auto const local_ids =
347 get_short_range_neighbors(system, p_id, m_max_exclusion_range);
348 assert(p1_ptr == nullptr or !!local_ids);
349 if (local_ids) {
350 particle_ids = std::move(*local_ids);
351 }
352 }
353
354 if (p1_ptr != nullptr) {
355 auto &p1 = *p1_ptr;
356 auto const &system = System::get_system();
357 auto const &box_geo = *system.box_geo;
358 auto &cell_structure = *system.cell_structure;
359
360 /* Check if the inserted particle within the exclusion radius of any other
361 * particle */
362 for (auto const p2_id : particle_ids) {
363 if (auto const p2_ptr = cell_structure.get_local_particle(p2_id)) {
364 auto const &p2 = *p2_ptr;
365 double excluded_distance;
366 if (exclusion_radius_per_type.count(p_type) == 0 ||
367 exclusion_radius_per_type.count(p2.type()) == 0) {
368 excluded_distance = exclusion_range;
369 } else if (exclusion_radius_per_type[p2.type()] == 0.) {
370 continue;
371 } else {
372 excluded_distance = exclusion_radius_per_type[p_type] +
373 exclusion_radius_per_type[p2.type()];
374 }
375
376 auto const d_min = box_geo.get_mi_vector(p2.pos(), p1.pos()).norm();
377
378 if (d_min < excluded_distance) {
380 break;
381 }
382 }
383 }
384 if (m_comm.rank() != 0) {
386 }
387 } else if (m_comm.rank() == 0) {
388 m_comm.recv(boost::mpi::any_source, 1,
390 }
391 boost::mpi::broadcast(m_comm, particle_inside_exclusion_range_touched, 0);
392}
393
394/**
395 * Deletes the particle with the given p_id and stores the id if the deletion
396 * created a hole in the particle id range. This method is intended to only
397 * delete unbonded particles since bonds are coupled to ids. This is used to
398 * avoid the id range becoming excessively huge.
399 */
401 if (p_id < 0) {
402 throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
403 }
404 auto const old_max_seen_id = ::get_maximal_particle_id();
405 if (p_id == old_max_seen_id) {
406 // last particle, just delete
407 remove_particle(p_id);
408 // remove all saved empty p_ids which are greater than the max_seen_particle
409 // this is needed in order to avoid the creation of holes
410 for (auto p_id_iter = m_empty_p_ids_smaller_than_max_seen_particle.begin();
412 if ((*p_id_iter) >= old_max_seen_id)
414 p_id_iter); // update iterator after container was modified
415 else
416 ++p_id_iter;
417 }
418 } else if (p_id <= old_max_seen_id) {
419 remove_particle(p_id);
421 } else {
422 throw std::runtime_error(
423 "Particle id is greater than the max seen particle id");
424 }
425}
426
427void ReactionAlgorithm::set_cyl_constraint(double center_x, double center_y,
428 double radius) {
429 auto const &box_geo = *System::get_system().box_geo;
430 if (center_x < 0. or center_x > box_geo.length()[0])
431 throw std::domain_error("center_x is outside the box");
432 if (center_y < 0. or center_y > box_geo.length()[1])
433 throw std::domain_error("center_y is outside the box");
434 if (radius < 0.)
435 throw std::domain_error("radius is invalid");
436 m_cyl_x = center_x;
437 m_cyl_y = center_y;
438 m_cyl_radius = radius;
439 m_reaction_constraint = ReactionConstraint::CYL_Z;
440}
441
443 double slab_end_z) {
444 auto const &box_geo = *System::get_system().box_geo;
445 if (slab_start_z < 0. or slab_start_z > box_geo.length()[2])
446 throw std::domain_error("slab_start_z is outside the box");
447 if (slab_end_z < 0. or slab_end_z > box_geo.length()[2])
448 throw std::domain_error("slab_end_z is outside the box");
449 if (slab_end_z < slab_start_z)
450 throw std::domain_error("slab_end_z must be >= slab_start_z");
451 m_slab_start_z = slab_start_z;
452 m_slab_end_z = slab_end_z;
453 m_reaction_constraint = ReactionConstraint::SLAB_Z;
454}
455
456/**
457 * Writes a random position inside the central box into the provided array.
458 */
460 auto const &box_geo = *System::get_system().box_geo;
461 Utils::Vector3d out_pos{};
462
463 if (m_reaction_constraint == ReactionConstraint::CYL_Z) {
464 // see http://mathworld.wolfram.com/DiskPointPicking.html
465 // for uniform disk point picking in cylinder
466 auto const random_radius =
467 m_cyl_radius * std::sqrt(m_uniform_real_distribution(m_generator));
468 auto const random_phi =
469 2. * std::numbers::pi * m_uniform_real_distribution(m_generator);
470 out_pos[0] = m_cyl_x + random_radius * cos(random_phi);
471 out_pos[1] = m_cyl_y + random_radius * sin(random_phi);
472 out_pos[2] = box_geo.length()[2] * m_uniform_real_distribution(m_generator);
473 } else if (m_reaction_constraint == ReactionConstraint::SLAB_Z) {
474 out_pos[0] = box_geo.length()[0] * m_uniform_real_distribution(m_generator);
475 out_pos[1] = box_geo.length()[1] * m_uniform_real_distribution(m_generator);
476 out_pos[2] = m_slab_start_z + (m_slab_end_z - m_slab_start_z) *
477 m_uniform_real_distribution(m_generator);
478 } else {
479 assert(m_reaction_constraint == ReactionConstraint::NONE);
480 out_pos[0] = box_geo.length()[0] * m_uniform_real_distribution(m_generator);
481 out_pos[1] = box_geo.length()[1] * m_uniform_real_distribution(m_generator);
482 out_pos[2] = box_geo.length()[2] * m_uniform_real_distribution(m_generator);
483 }
484 return out_pos;
485}
486
487/**
488 * Creates a particle at the end of the observed particle id range.
489 */
490int ReactionAlgorithm::create_particle(int p_type) {
491 int p_id;
493 auto p_id_iter = std::min_element(
496 p_id = *p_id_iter;
498 } else {
499 p_id = ::get_maximal_particle_id() + 1;
500 }
501
502 // create random velocity vector according to Maxwell-Boltzmann distribution
503 auto pos = get_random_position_in_box();
504 auto vel = get_random_velocity_vector();
505
506 ::make_new_particle(p_id, pos);
507 if (auto p = get_local_particle(p_id)) {
508 p->v() = std::sqrt(kT / p->mass()) * vel;
509 p->type() = p_type;
510#ifdef ELECTROSTATICS
511 p->q() = charges_of_types.at(p_type);
512#endif
513 }
515 return p_id;
516}
517
518void ReactionAlgorithm::displacement_mc_move(int type, int n_particles) {
519 auto &bookkeeping = make_new_system_state();
520 // draw particle ids at random without replacement
521 int p_id = -1;
522 std::vector<int> drawn_pids{p_id};
523 for (int i = 0; i < n_particles; i++) {
524 // draw a new particle id
525 while (Utils::contains(drawn_pids, p_id)) {
526 auto const random_index = i_random(number_of_particles_with_type(type));
527 p_id = get_random_p_id(type, random_index);
528 }
529 drawn_pids.emplace_back(p_id);
530 // write new position and new velocity
531 typename decltype(ParticleChanges::moved)::value_type old_state;
532 auto const new_pos = get_random_position_in_box();
533 auto vel = get_random_velocity_vector();
534 if (auto p = get_real_particle(p_id)) {
535 old_state = {p_id, p->pos(), p->v()};
536 p->v() = std::sqrt(kT / p->mass()) * vel;
537 if (m_comm.rank() != 0) {
538 m_comm.send(0, 42, old_state);
539 }
540 } else if (m_comm.rank() == 0) {
541 m_comm.recv(boost::mpi::any_source, 42, old_state);
542 }
543 boost::mpi::broadcast(m_comm, old_state, 0);
544 bookkeeping.moved.emplace_back(old_state);
545 ::set_particle_pos(p_id, new_pos);
546
547 check_exclusion_range(p_id, type);
549 break;
550 }
551 }
552}
553
555 int n_particles) {
556
557 if (type < 0) {
558 throw std::domain_error("Parameter 'type_mc' must be >= 0");
559 }
560 if (n_particles < 0) {
561 throw std::domain_error(
562 "Parameter 'particle_number_to_be_changed' must be >= 0");
563 }
564
565 if (n_particles == 0) {
566 // reject
567 return false;
568 }
569
572
573 auto const n_particles_of_type = ::number_of_particles_with_type(type);
574 if (n_particles > n_particles_of_type) {
575 // reject
576 return false;
577 }
578
579 auto const E_pot_old = calculate_potential_energy();
580 displacement_mc_move(type, n_particles);
581 auto const E_pot_new = (particle_inside_exclusion_range_touched)
582 ? std::numeric_limits<double>::max()
584 auto constexpr exp_min = -708.4; // for IEEE-compatible double
585 auto const exponent = -(E_pot_new - E_pot_old) / kT;
586 auto const exponential = (exponent < exp_min) ? 0. : std::exp(exponent);
587
588 // Metropolis algorithm since proposal density is symmetric
589 auto const bf = std::min(1., exponential);
590
591 // // correct for enhanced proposal of small radii by using the
592 // // Metropolis-Hastings algorithm for asymmetric proposal densities
593 // double old_radius =
594 // std::sqrt(std::pow(particle_positions[0][0]-cyl_x,2) +
595 // std::pow(particle_positions[0][1]-cyl_y,2));
596 // double new_radius =
597 // std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2));
598 // auto const bf = std::min(1.0,
599 // exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius);
600
601 // Metropolis-Hastings algorithm for asymmetric proposal density
602 if (m_uniform_real_distribution(m_generator) < bf) {
603 // accept
606 return true;
607 }
608 // reject: restore original particle properties
610 return false;
611}
612
613/**
614 * Cleans the list of empty pids and searches for empty pid in the system
615 */
617 // Clean-up the list of empty pids
619
620 auto particle_ids = get_particle_ids_parallel();
621 std::sort(particle_ids.begin(), particle_ids.end());
622 auto pid1 = -1;
623 for (auto pid2 : particle_ids) {
624 for (int pid = pid1 + 1; pid < pid2; ++pid) {
626 }
627 pid1 = pid2;
628 }
629}
630
632 auto &system = System::get_system();
633 auto const obs = system.calculate_energy();
634 auto pot = obs->accumulate(-obs->kinetic[0]);
635 boost::mpi::broadcast(m_comm, pot, 0);
636 return pot;
637}
638
639Particle *ReactionAlgorithm::get_real_particle(int p_id) const {
640 assert(p_id >= 0);
641 auto const &system = System::get_system();
642 auto ptr = system.cell_structure->get_local_particle(p_id);
643 if (ptr != nullptr and ptr->is_ghost()) {
644 ptr = nullptr;
645 }
646 assert(boost::mpi::all_reduce(m_comm, static_cast<int>(ptr != nullptr),
647 std::plus<>()) == 1);
648 return ptr;
649}
650
651Particle *ReactionAlgorithm::get_local_particle(int p_id) const {
652 assert(p_id >= 0);
653 auto const &system = System::get_system();
654 auto ptr = system.cell_structure->get_local_particle(p_id);
655 assert(boost::mpi::all_reduce(
656 m_comm, static_cast<int>(ptr != nullptr and not ptr->is_ghost()),
657 std::plus<>()) == 1);
658 return ptr;
659}
660
661} // namespace ReactionMethods
Vector implementation and trait types for boost qvm interoperability.
std::optional< std::vector< int > > get_short_range_neighbors(System::System const &system, int const pid, double const distance)
Get ids of particles that are within a certain distance of another particle.
Definition cells.cpp:119
This file contains everything related to the global cell structure / cell system.
void add_reaction(std::shared_ptr< SingleReaction > const &new_reaction)
Adds a reaction to the reaction system.
double calculate_potential_energy() const
Compute the system potential energy.
void make_reaction_attempt(::ReactionMethods::SingleReaction const &reaction, ParticleChanges &bookkeeping)
Carry out a chemical reaction and save the old system state.
Utils::Vector3d get_random_position_in_box()
Writes a random position inside the central box into the provided array.
void set_slab_constraint(double slab_start_z, double slab_end_z)
int i_random(int maxint)
draws a random integer from the uniform distribution in the range [0,maxint-1]
void set_cyl_constraint(double center_x, double center_y, double radius)
void setup_bookkeeping_of_empty_pids()
Cleans the list of empty pids and searches for empty pid in the system.
void displacement_mc_move(int type, int n_particles)
Carry out displacement MC moves for particles of a given type.
void update_volume()
Automatically sets the volume which is used by the reaction ensemble to the volume of a cuboid box.
double make_reaction_mc_move_attempt(int reaction_id, double bf, double E_pot_old, double E_pot_new)
Accept or reject a reaction MC move made by create_new_trial_state based on a probability acceptance ...
void restore_old_system_state()
Restore last valid system state.
bool all_reactant_particles_exist(SingleReaction const &current_reaction) const
Checks whether all particles exist for the provided reaction.
auto & make_new_system_state()
Open new handle for system state tracking.
void clear_old_system_state()
Clear last valid system state.
void delete_particle(int p_id)
Deletes the particle with the given p_id and stores the id if the deletion created a hole in the part...
std::vector< std::shared_ptr< SingleReaction > > reactions
std::optional< double > create_new_trial_state(int reaction_id)
Carry out a reaction MC move and calculate the new potential energy.
bool make_displacement_mc_move_attempt(int type, int n_particles)
Attempt displacement MC moves for particles of a given type.
void on_particle_change()
Called every time a particle property changes.
void on_particle_local_change()
Called every time a particle local property changes.
std::shared_ptr< BoxGeometry > box_geo
This file contains the defaults for ESPResSo.
System & get_system()
bool contains(InputIt first, InputIt last, T const &value)
Check whether an iterator range contains a value.
Definition contains.hpp:36
void serialize(Archive &ar, std::tuple< T... > &pack, unsigned int const)
Serialize std::tuple.
auto constexpr new_part
auto constexpr any_type
std::vector< int > get_particle_ids_parallel()
void make_new_particle(int p_id, Utils::Vector3d const &pos)
Create a new particle and attach it to a cell.
int number_of_particles_with_type(int type)
void set_particle_pos(int p_id, Utils::Vector3d const &pos)
Move particle to a new position.
void init_type_map(int type)
int get_random_p_id(int type, int random_index_in_type_map)
Find a particle of given type and return its id.
void remove_particle(int p_id)
Remove particle with a given identity.
int get_maximal_particle_id()
Get maximal particle id.
void on_particle_type_change(int p_id, int old_type, int new_type)
Particles creation and deletion.
Struct holding all information for one particle.
Definition Particle.hpp:395
std::vector< std::tuple< int, Utils::Vector3d, Utils::Vector3d > > moved