ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
core/reaction_methods/ReactionAlgorithm.hpp
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#ifndef REACTION_METHODS_REACTION_ALGORITHM_HPP
20#define REACTION_METHODS_REACTION_ALGORITHM_HPP
21
22#include "config/config.hpp"
23
24#include "SingleReaction.hpp"
25
26#include "Particle.hpp"
27#include "random.hpp"
28
29#include <utils/Vector.hpp>
30
31#include <functional>
32#include <memory>
33#include <optional>
34#include <random>
35#include <stdexcept>
36#include <tuple>
37#include <unordered_map>
38#include <utility>
39#include <vector>
40
41namespace boost::mpi {
42class communicator;
43} // namespace boost::mpi
44
45namespace ReactionMethods {
46
47/** Base class for reaction ensemble methods */
49private:
50 boost::mpi::communicator const &m_comm;
51
52public:
54 boost::mpi::communicator const &comm, int seed, double kT,
55 double exclusion_range,
56 std::unordered_map<int, double> const &exclusion_radius_per_type)
57 : m_comm{comm}, kT{kT}, exclusion_range{exclusion_range},
58 m_generator(Random::mt19937(std::seed_seq({seed, seed, seed}))),
59 m_normal_distribution(0.0, 1.0), m_uniform_real_distribution(0.0, 1.0) {
60 if (kT < 0.) {
61 throw std::domain_error("Invalid value for 'kT'");
62 }
63 if (exclusion_range < 0.) {
64 throw std::domain_error("Invalid value for 'exclusion_range'");
65 }
68 }
69
70 virtual ~ReactionAlgorithm() = default;
71
72 std::vector<std::shared_ptr<SingleReaction>> reactions;
73 std::unordered_map<int, double> charges_of_types;
74 double kT;
75 /**
76 * Hard sphere radius. If particles are closer than this value,
77 * it is assumed that their interaction energy gets approximately
78 * infinite, therefore these configurations do not contribute
79 * to the partition function and ensemble averages.
80 */
82 std::unordered_map<int, double> exclusion_radius_per_type;
83 double volume;
85
89 return static_cast<double>(m_accepted_configurational_MC_moves) /
90 static_cast<double>(m_tried_configurational_MC_moves);
91 }
92
93 auto get_kT() const { return kT; }
94 auto get_exclusion_range() const { return exclusion_range; }
95 auto get_volume() const { return volume; }
96 void set_volume(double new_volume) {
97 if (new_volume <= 0.) {
98 throw std::domain_error("Invalid value for 'volume'");
99 }
100 volume = new_volume;
101 }
102 void update_volume();
103 void
104 set_exclusion_radius_per_type(std::unordered_map<int, double> const &map) {
105 auto max_exclusion_range = exclusion_range;
106 for (auto const &item : map) {
107 auto const type = item.first;
108 auto const exclusion_radius = item.second;
109 if (exclusion_radius < 0.) {
110 throw std::domain_error("Invalid excluded_radius value for type " +
111 std::to_string(type) + ": radius " +
112 std::to_string(exclusion_radius));
113 }
114 max_exclusion_range =
115 std::max(max_exclusion_range, 2. * exclusion_radius);
116 }
118 m_max_exclusion_range = max_exclusion_range;
119 }
120
121 void remove_constraint() { m_reaction_constraint = ReactionConstraint::NONE; }
122 void set_cyl_constraint(double center_x, double center_y, double radius);
123 void set_slab_constraint(double slab_start_z, double slab_end_z);
125 if (m_reaction_constraint != ReactionConstraint::SLAB_Z) {
126 throw std::runtime_error("no slab constraint is currently active");
127 }
128 return {m_slab_start_z, m_slab_end_z};
129 }
130
132 void delete_particle(int p_id);
133 void add_reaction(std::shared_ptr<SingleReaction> const &new_reaction);
134 void delete_reaction(int reaction_id) {
135 reactions.erase(reactions.begin() + reaction_id);
136 }
137
140
141protected:
143
144public:
146 std::vector<int> created{};
147 std::vector<std::tuple<int, int>> changed{};
148 std::vector<std::tuple<int, int>> hidden{};
149 std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>> moved{};
150 std::unordered_map<int, int> old_particle_numbers{};
151 int reaction_id{-1};
152 };
153
154 bool is_reaction_under_way() const { return m_system_changes != nullptr; }
155 auto const &get_old_system_state() const {
156 if (not is_reaction_under_way()) {
157 throw std::runtime_error("No chemical reaction is currently under way");
158 }
159 return *m_system_changes;
160 }
161
162protected:
163 /** @brief Restore last valid system state. */
165 /** @brief Clear last valid system state. */
167 assert(is_reaction_under_way());
168 m_system_changes = nullptr;
169 }
170 /** @brief Open new handle for system state tracking. */
172 assert(not is_reaction_under_way());
173 m_system_changes = std::make_shared<ParticleChanges>();
174 return *m_system_changes;
175 }
176 std::shared_ptr<ParticleChanges> m_system_changes;
177
178public:
179 /**
180 * Carry out a reaction MC move and calculate the new potential energy.
181 * Particles are selected without replacement.
182 * The previous state of the system is cached.
183 * @returns Potential energy of the system after the move.
184 */
185 std::optional<double> create_new_trial_state(int reaction_id);
186 /**
187 * Accept or reject a reaction MC move made by @ref create_new_trial_state
188 * based on a probability acceptance @c bf.
189 * The previous state of the system is either restored from the cache if
190 * the move is rejected, or cleared from the cache if the move is accepted.
191 * @returns Potential energy of the system after the move was accepted or
192 * rejected.
193 */
194 double make_reaction_mc_move_attempt(int reaction_id, double bf,
195 double E_pot_old, double E_pot_new);
196 /**
197 * Attempt displacement MC moves for particles of a given type.
198 * Particles are selected without replacement.
199 * @param type Type of particles to move.
200 * @param n_particles Number of particles of that type to move.
201 * @returns true if all moves were accepted.
202 */
203 bool make_displacement_mc_move_attempt(int type, int n_particles);
204
205 /** @brief Compute the system potential energy. */
206 double calculate_potential_energy() const;
207
208protected:
209 /**
210 * Carry out displacement MC moves for particles of a given type.
211 * Particles are selected without replacement.
212 * @param type Type of particles to move.
213 * @param n_particles Number of particles of that type to move.
214 */
215 void displacement_mc_move(int type, int n_particles);
216
217 /** @brief Carry out a chemical reaction and save the old system state. */
219 ParticleChanges &bookkeeping);
220
221public:
222 /**
223 * @brief draws a random integer from the uniform distribution in the range
224 * [0,maxint-1]
225 *
226 * @param maxint range.
227 */
228 int i_random(int maxint) {
229 std::uniform_int_distribution<int> uniform_int_dist(0, maxint - 1);
230 return uniform_int_dist(m_generator);
231 }
232
233protected:
234 bool
235 all_reactant_particles_exist(SingleReaction const &current_reaction) const;
236
237private:
238 std::mt19937 m_generator;
239 std::normal_distribution<double> m_normal_distribution;
240 std::uniform_real_distribution<double> m_uniform_real_distribution;
241
242 std::unordered_map<int, int>
243 get_particle_numbers(SingleReaction const &reaction) const;
244
245 int create_particle(int p_type);
246 void hide_particle(int p_id, int p_type) const;
247 void check_exclusion_range(int p_id, int p_type);
248 auto get_random_uniform_number() {
249 return m_uniform_real_distribution(m_generator);
250 }
251 auto get_random_velocity_vector() {
252 return Utils::Vector3d{m_normal_distribution(m_generator),
253 m_normal_distribution(m_generator),
254 m_normal_distribution(m_generator)};
255 }
256
257 enum class ReactionConstraint { NONE, CYL_Z, SLAB_Z };
258 ReactionConstraint m_reaction_constraint = ReactionConstraint::NONE;
259 double m_cyl_radius = -10.0;
260 double m_cyl_x = -10.0;
261 double m_cyl_y = -10.0;
262 double m_slab_start_z = -10.0;
263 double m_slab_end_z = -10.0;
264 double m_max_exclusion_range = 0.;
265
266 Particle *get_real_particle(int p_id) const;
267 Particle *get_local_particle(int p_id) const;
268
269protected:
271};
272
273} // namespace ReactionMethods
274#endif
Vector implementation and trait types for boost qvm interoperability.
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]
virtual ~ReactionAlgorithm()=default
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.
void set_exclusion_radius_per_type(std::unordered_map< int, double > const &map)
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.
ReactionAlgorithm(boost::mpi::communicator const &comm, int seed, double kT, double exclusion_range, std::unordered_map< int, double > const &exclusion_radius_per_type)
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.
Communicator communicator
This file contains the defaults for ESPResSo.
Random number generation using Philox.
Struct holding all information for one particle.
Definition Particle.hpp:395
std::vector< std::tuple< int, Utils::Vector3d, Utils::Vector3d > > moved