19#ifndef REACTION_METHODS_REACTION_ALGORITHM_HPP
20#define REACTION_METHODS_REACTION_ALGORITHM_HPP
37#include <unordered_map>
50 boost::mpi::communicator
const &m_comm;
54 boost::mpi::communicator
const &comm,
int seed,
double kT,
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) {
61 throw std::domain_error(
"Invalid value for 'kT'");
64 throw std::domain_error(
"Invalid value for 'exclusion_range'");
72 std::vector<std::shared_ptr<SingleReaction>>
reactions;
97 if (new_volume <= 0.) {
98 throw std::domain_error(
"Invalid value for 'volume'");
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));
114 max_exclusion_range =
115 std::max(max_exclusion_range, 2. * exclusion_radius);
118 m_max_exclusion_range = max_exclusion_range;
125 if (m_reaction_constraint != ReactionConstraint::SLAB_Z) {
126 throw std::runtime_error(
"no slab constraint is currently active");
128 return {m_slab_start_z, m_slab_end_z};
133 void add_reaction(std::shared_ptr<SingleReaction>
const &new_reaction);
148 std::vector<std::tuple<int, int>>
hidden{};
149 std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>>
moved{};
157 throw std::runtime_error(
"No chemical reaction is currently under way");
195 double E_pot_old,
double E_pot_new);
229 std::uniform_int_distribution<int> uniform_int_dist(0, maxint - 1);
230 return uniform_int_dist(m_generator);
238 std::mt19937 m_generator;
239 std::normal_distribution<double> m_normal_distribution;
240 std::uniform_real_distribution<double> m_uniform_real_distribution;
242 std::unordered_map<int, int>
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);
251 auto get_random_velocity_vector() {
253 m_normal_distribution(m_generator),
254 m_normal_distribution(m_generator)};
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.;
266 Particle *get_real_particle(
int p_id)
const;
267 Particle *get_local_particle(
int p_id)
const;
Vector implementation and trait types for boost qvm interoperability.
Base class for reaction ensemble methods.
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.
int m_tried_configurational_MC_moves
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)
bool particle_inside_exclusion_range_touched
int i_random(int maxint)
draws a random integer from the uniform distribution in the range [0,maxint-1]
virtual ~ReactionAlgorithm()=default
std::unordered_map< int, double > charges_of_types
void set_cyl_constraint(double center_x, double center_y, double radius)
std::unordered_map< int, double > exclusion_radius_per_type
void setup_bookkeeping_of_empty_pids()
Cleans the list of empty pids and searches for empty pid in the system.
void delete_reaction(int reaction_id)
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.
auto get_exclusion_range() const
bool all_reactant_particles_exist(SingleReaction const ¤t_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.
std::shared_ptr< ParticleChanges > m_system_changes
void clear_old_system_state()
Clear last valid system state.
bool is_reaction_under_way() const
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...
int m_accepted_configurational_MC_moves
std::vector< std::shared_ptr< SingleReaction > > reactions
double get_acceptance_rate_configurational_moves() const
Utils::Vector2d get_slab_constraint_parameters() const
double exclusion_range
Hard sphere radius.
void set_volume(double new_volume)
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.
bool neighbor_search_order_n
auto const & get_old_system_state() const
std::vector< int > m_empty_p_ids_smaller_than_max_seen_particle
Communicator communicator
This file contains the defaults for ESPResSo.
Random number generation using Philox.
Struct holding all information for one particle.
std::vector< int > created
std::unordered_map< int, int > old_particle_numbers
std::vector< std::tuple< int, int > > hidden
std::vector< std::tuple< int, int > > changed
std::vector< std::tuple< int, Utils::Vector3d, Utils::Vector3d > > moved