From 801857de88af93e89e66e4597f9d8d353b7bb7fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 7 Feb 2023 19:36:14 +0100 Subject: [PATCH 1/5] Make reaction methods run in parallel MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rewrite reaction methods without MpiCallbacks. Add sanity checks. Rewrite some of the functionality as Python code and expose the `calculate_acceptance_probability()` function in Python to make prototyping of new Monte Carlo methods easier. Reduce overhead from the cell system when altering particle properties. Co-authored-by: Pablo Miguel Blanco Andrés --- src/core/particle_node.cpp | 94 ++-- src/core/particle_node.hpp | 8 +- src/core/reaction_methods/CMakeLists.txt | 8 +- .../reaction_methods/ConstantpHEnsemble.cpp | 50 -- .../reaction_methods/ConstantpHEnsemble.hpp | 12 +- .../reaction_methods/ReactionAlgorithm.cpp | 493 +++++++++--------- .../reaction_methods/ReactionAlgorithm.hpp | 182 ++++--- .../reaction_methods/ReactionEnsemble.cpp | 45 -- .../reaction_methods/ReactionEnsemble.hpp | 13 +- src/core/reaction_methods/SingleReaction.hpp | 4 + src/core/reaction_methods/WidomInsertion.cpp | 51 -- src/core/reaction_methods/WidomInsertion.hpp | 25 +- .../reaction_methods/tests/CMakeLists.txt | 4 - .../tests/ConstantpHEnsemble_test.cpp | 92 ---- .../tests/ReactionAlgorithm_test.cpp | 157 +++--- .../tests/ReactionEnsemble_test.cpp | 152 ------ .../tests/SingleReaction_test.cpp | 6 +- .../tests/particle_tracking_test.cpp | 2 - .../tests/reaction_methods_utils_test.cpp | 8 +- src/core/reaction_methods/utils.cpp | 68 ++- src/core/reaction_methods/utils.hpp | 18 +- src/python/espressomd/reaction_methods.py | 267 +++++++--- .../particle_data/ParticleHandle.cpp | 4 +- .../reaction_methods/CMakeLists.txt | 6 +- .../reaction_methods/ConstantpHEnsemble.hpp | 30 +- .../reaction_methods/ReactionAlgorithm.cpp | 248 +++++++++ .../reaction_methods/ReactionAlgorithm.hpp | 128 +---- .../reaction_methods/ReactionEnsemble.hpp | 19 +- .../reaction_methods/SingleReaction.hpp | 32 +- .../reaction_methods/WidomInsertion.hpp | 39 +- src/script_interface/system/System.cpp | 15 +- src/script_interface/tests/CMakeLists.txt | 4 + .../tests/ConstantpHEnsemble_test.cpp | 135 +++++ .../tests/ReactionEnsemble_test.cpp | 266 ++++++++++ testsuite/python/constant_pH.py | 5 +- testsuite/python/constant_pH_stats.py | 5 +- testsuite/python/reaction_ensemble.py | 13 +- .../python/reaction_methods_interface.py | 36 +- 38 files changed, 1594 insertions(+), 1150 deletions(-) delete mode 100644 src/core/reaction_methods/ConstantpHEnsemble.cpp delete mode 100644 src/core/reaction_methods/ReactionEnsemble.cpp delete mode 100644 src/core/reaction_methods/WidomInsertion.cpp delete mode 100644 src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp delete mode 100644 src/core/reaction_methods/tests/ReactionEnsemble_test.cpp create mode 100644 src/script_interface/reaction_methods/ReactionAlgorithm.cpp create mode 100644 src/script_interface/tests/ConstantpHEnsemble_test.cpp create mode 100644 src/script_interface/tests/ReactionEnsemble_test.cpp diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp index 00b99fb73ba..42d43d385ca 100644 --- a/src/core/particle_node.cpp +++ b/src/core/particle_node.cpp @@ -26,6 +26,7 @@ #include "communication.hpp" #include "event.hpp" #include "grid.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "partCfg_global.hpp" #include "particle_data.hpp" @@ -35,6 +36,7 @@ #include #include +#include #include #include #include @@ -86,15 +88,26 @@ static void mpi_synchronize_max_seen_pid_local() { REGISTER_CALLBACK(mpi_synchronize_max_seen_pid_local) void init_type_map(int type) { - type_list_enable = true; - if (type < 0) + if (type < 0) { throw std::runtime_error("Types may not be negative"); + } + ::type_list_enable = true; + make_particle_type_exist_local(type); - auto &map_for_type = particle_type_map[type]; - map_for_type.clear(); - for (auto const &p : partCfg()) { - if (p.type() == type) - map_for_type.insert(p.id()); + std::vector local_pids; + for (auto const &p : ::cell_structure.local_particles()) { + if (p.type() == type) { + local_pids.emplace_back(p.id()); + } + } + + std::vector> global_pids; + boost::mpi::all_gather(::comm_cart, local_pids, global_pids); + ::particle_type_map[type].clear(); + for (auto const &vec : global_pids) { + for (auto const &p_id : vec) { + ::particle_type_map[type].insert(p_id); + } } } @@ -110,11 +123,25 @@ static void add_id_to_type_map(int p_id, int type) { it->second.insert(p_id); } -void on_particle_type_change_head(int p_id, int old_type, int new_type) { - if (type_list_enable) { - assert(::this_node == 0); - if (old_type != new_type) { - remove_id_from_map(p_id, old_type); +void on_particle_type_change_parallel(int p_id, int old_type, int new_type) { + if (::type_list_enable) { + if (old_type == type_tracking::any_type) { + for (auto &kv : ::particle_type_map) { + auto it = kv.second.find(p_id); + if (it != kv.second.end()) { + kv.second.erase(it); +#ifndef NDEBUG + if (auto p = ::cell_structure.get_local_particle(p_id)) { + assert(p->type() == kv.first); + } +#endif + break; + } + } + } else if (old_type != type_tracking::new_part) { + if (old_type != new_type) { + remove_id_from_map(p_id, old_type); + } } add_id_to_type_map(p_id, new_type); } @@ -478,26 +505,21 @@ void remove_particle(int p_id) { } void remove_particle_parallel(int p_id) { - { - auto track_particle_types = ::type_list_enable; - boost::mpi::broadcast(::comm_cart, track_particle_types, 0); - if (track_particle_types) { - auto p = ::cell_structure.get_local_particle(p_id); - auto p_type = -1; - if (p != nullptr and not p->is_ghost()) { - if (this_node == 0) { - p_type = p->type(); - } else { - ::comm_cart.send(0, 42, p->type()); - } - } else if (this_node == 0) { - ::comm_cart.recv(boost::mpi::any_source, 42, p_type); - } + if (::type_list_enable) { + auto p = ::cell_structure.get_local_particle(p_id); + auto p_type = -1; + if (p != nullptr and not p->is_ghost()) { if (this_node == 0) { - assert(p_type != -1); - remove_id_from_map(p_id, p_type); + p_type = p->type(); + } else { + ::comm_cart.send(0, 42, p->type()); } + } else if (this_node == 0) { + ::comm_cart.recv(boost::mpi::any_source, 42, p_type); } + assert(not(this_node == 0 and p_type == -1)); + boost::mpi::broadcast(::comm_cart, p_type, 0); + remove_id_from_map(p_id, p_type); } if (this_node == 0) { @@ -575,7 +597,10 @@ int get_random_p_id(int type, int random_index_in_type_map) { if (random_index_in_type_map + 1 > it->second.size()) throw std::runtime_error("The provided index exceeds the number of " "particle types listed in the particle_type_map"); - return *std::next(it->second.begin(), random_index_in_type_map); + // there is no guarantee of order across MPI ranks + auto p_id = *std::next(it->second.begin(), random_index_in_type_map); + boost::mpi::broadcast(::comm_cart, p_id, 0); + return p_id; } int number_of_particles_with_type(int type) { @@ -605,6 +630,15 @@ std::vector get_particle_ids() { return ids; } +std::vector get_particle_ids_parallel() { + if (rebuild_needed()) { + build_particle_node_parallel(); + } + auto pids = Utils::keys(particle_node); + boost::mpi::broadcast(::comm_cart, pids, 0); + return pids; +} + int get_maximal_particle_id() { if (particle_node.empty()) build_particle_node(); diff --git a/src/core/particle_node.hpp b/src/core/particle_node.hpp index a519d4e5d78..fa60f258465 100644 --- a/src/core/particle_node.hpp +++ b/src/core/particle_node.hpp @@ -36,6 +36,11 @@ #include #include +namespace type_tracking { +auto constexpr any_type = -2; +auto constexpr new_part = -3; +} // namespace type_tracking + /** * @brief Get particle data. * @@ -102,7 +107,7 @@ void remove_all_particles(); void init_type_map(int type); void on_particle_type_change(int p_id, int type); -void on_particle_type_change_head(int p_id, int old_type, int new_type); +void on_particle_type_change_parallel(int p_id, int old_type, int new_type); /** Find a particle of given type and return its id */ int get_random_p_id(int type, int random_index_in_type_map); @@ -130,6 +135,7 @@ int get_particle_node(int p_id); * @return Sorted ids of all existing particles. */ std::vector get_particle_ids(); +std::vector get_particle_ids_parallel(); /** * @brief Get maximal particle id. diff --git a/src/core/reaction_methods/CMakeLists.txt b/src/core/reaction_methods/CMakeLists.txt index 40a571f51f4..88a046e9b79 100644 --- a/src/core/reaction_methods/CMakeLists.txt +++ b/src/core/reaction_methods/CMakeLists.txt @@ -18,12 +18,8 @@ # target_sources( - espresso_core - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/ReactionEnsemble.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/ConstantpHEnsemble.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/WidomInsertion.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp) + espresso_core PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp) if(ESPRESSO_BUILD_TESTS) add_subdirectory(tests) diff --git a/src/core/reaction_methods/ConstantpHEnsemble.cpp b/src/core/reaction_methods/ConstantpHEnsemble.cpp deleted file mode 100644 index 05fbe81287e..00000000000 --- a/src/core/reaction_methods/ConstantpHEnsemble.cpp +++ /dev/null @@ -1,50 +0,0 @@ -/* - * Copyright (C) 2010-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 . - */ - -#include "reaction_methods/ConstantpHEnsemble.hpp" - -#include "particle_data.hpp" - -#include "utils.hpp" - -#include -#include -#include - -namespace ReactionMethods { - -/** - * Calculates the expression in the acceptance probability of the constant pH - * method. - */ -double ConstantpHEnsemble::calculate_acceptance_probability( - SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, - std::map const &old_particle_numbers) const { - auto const beta = 1.0 / kT; - auto const pKa = -current_reaction.nu_bar * log10(current_reaction.gamma); - auto const ln_bf = (E_pot_new - E_pot_old) - current_reaction.nu_bar / beta * - log(10) * - (m_constant_pH - pKa); - const double factorial_expr = calculate_factorial_expression_cpH( - current_reaction, old_particle_numbers); - auto const bf = factorial_expr * exp(-beta * ln_bf); - return bf; -} - -} // namespace ReactionMethods diff --git a/src/core/reaction_methods/ConstantpHEnsemble.hpp b/src/core/reaction_methods/ConstantpHEnsemble.hpp index a63f9966216..ea526aa0613 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.hpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.hpp @@ -41,17 +41,13 @@ namespace ReactionMethods { class ConstantpHEnsemble : public ReactionAlgorithm { public: ConstantpHEnsemble( - int seed, double kT, double exclusion_range, double constant_pH, + boost::mpi::communicator const &comm, int seed, double kT, + double exclusion_range, double constant_pH, const std::unordered_map &exclusion_radius_per_type) - : ReactionAlgorithm(seed, kT, exclusion_range, exclusion_radius_per_type), + : ReactionAlgorithm(comm, seed, kT, exclusion_range, + exclusion_radius_per_type), m_constant_pH(constant_pH) {} double m_constant_pH; - -protected: - double calculate_acceptance_probability( - SingleReaction const ¤t_reaction, double E_pot_old, - double E_pot_new, - std::map const &old_particle_numbers) const override; }; } // namespace ReactionMethods diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index be949fd0e26..f259432cc8e 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -21,18 +21,20 @@ #include "reaction_methods/ReactionAlgorithm.hpp" -#include "analysis/statistics.hpp" #include "cells.hpp" #include "energy.hpp" +#include "event.hpp" #include "grid.hpp" -#include "partCfg_global.hpp" -#include "particle_data.hpp" #include "particle_node.hpp" #include #include #include +#include +#include +#include + #include #include #include @@ -40,29 +42,21 @@ #include #include #include +#include #include #include #include #include #include -namespace ReactionMethods { - -/** - * Performs a randomly selected reaction in the reaction ensemble - */ -void ReactionAlgorithm::do_reaction(int reaction_steps) { - // calculate potential energy; only consider potential energy since we - // assume that the kinetic part drops out in the process of calculating - // ensemble averages (kinetic part may be separated and crossed out) - auto current_E_pot = mpi_calculate_potential_energy(); - // Setup the list of empty pids for bookeeping - setup_bookkeeping_of_empty_pids(); - for (int i = 0; i < reaction_steps; i++) { - int reaction_id = i_random(static_cast(reactions.size())); - generic_oneway_reaction(*reactions[reaction_id], current_E_pot); - } +namespace boost::serialization { +template +void serialize(Archive &ar, std::tuple &tuple, const unsigned int) { + std::apply([&](auto &...item) { ((ar & item), ...); }, tuple); } +} // namespace boost::serialization + +namespace ReactionMethods { /** * Adds a reaction to the reaction system @@ -81,63 +75,46 @@ void ReactionAlgorithm::add_reaction( reactions.push_back(new_reaction); } -void ParticleChangeRecorder::restore_original_state() const { - auto const restore_state = [](StoredParticleProperty const &item) { -#ifdef ELECTROSTATICS - set_particle_q(item.p_id, item.charge); -#endif - set_particle_type(item.p_id, item.type); - }; - // 1) delete created product particles - std::for_each(m_created.begin(), m_created.end(), m_deleter); - // 2) restore previously hidden reactant particles - std::for_each(m_hidden.begin(), m_hidden.end(), restore_state); - // 3) restore previously changed reactant particles - std::for_each(m_changed.begin(), m_changed.end(), restore_state); -} - -void ParticleChangeRecorder::delete_hidden_particles() const { - for (auto const &item : m_hidden) { - // change back type otherwise the bookkeeping algorithm breaks - set_particle_type(item.p_id, item.type); - m_deleter(item.p_id); - } -} - /** - * Checks whether all necessary variables for the reaction ensemble have been - * set. + * @details This method tries to keep the cell system overhead to a minimum. + * Event callbacks are only called once after all particles are updated, + * except for particle deletion (the cell structure is still reinitialized + * after each deletion). */ -void ReactionAlgorithm::check_reaction_method() const { - if (reactions.empty()) { - throw std::runtime_error("Reaction system not initialized"); - } - +void ReactionAlgorithm::restore_old_system_state() { + auto const &old_state = get_old_system_state(); + // restore the properties of changed and hidden particles + for (auto const &state : {old_state.changed, old_state.hidden}) { + for (auto const &[p_id, p_type] : state) { + on_particle_type_change_parallel(p_id, type_tracking::any_type, p_type); + if (auto p = get_local_particle(p_id)) { + p->type() = p_type; #ifdef ELECTROSTATICS - // check for the existence of default charges for all types that take part in - // the reactions - - for (const auto ¤t_reaction : reactions) { - // check for reactants - for (int reactant_type : current_reaction->reactant_types) { - auto it = charges_of_types.find(reactant_type); - if (it == charges_of_types.end()) { - std::string message = std::string("Forgot to assign charge to type ") + - std::to_string(reactant_type); - throw std::runtime_error(message); + p->q() = charges_of_types.at(p_type); +#endif } } - // check for products - for (int product_type : current_reaction->product_types) { - auto it = charges_of_types.find(product_type); - if (it == charges_of_types.end()) { - std::string message = std::string("Forgot to assign charge to type ") + - std::to_string(product_type); - throw std::runtime_error(message); - } + } + // delete created particles + for (auto const p_id : old_state.created) { + delete_particle(p_id); + } + // restore original positions and velocities + for (auto const &[p_id, pos, vel] : old_state.moved) { + if (auto p = get_local_particle(p_id)) { + p->v() = vel; + auto folded_pos = pos; + auto image_box = Utils::Vector3i{}; + fold_position(folded_pos, image_box, box_geo); + p->pos() = folded_pos; + p->image_box() = image_box; } } -#endif + if (not old_state.moved.empty()) { + ::cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); + } + on_particle_change(); + clear_old_system_state(); } /** @@ -163,11 +140,8 @@ bool ReactionAlgorithm::all_reactant_particles_exist( return enough_particles; } -/** - *Performs a trial reaction move - */ -ParticleChangeRecorder -ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction) { +void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, + ParticleChanges &bookkeeping) { // create or hide particles of types with corresponding types in reaction auto const n_product_types = reaction.product_types.size(); auto const n_reactant_types = reaction.reactant_types.size(); @@ -175,18 +149,25 @@ ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction) { auto const random_index = i_random(number_of_particles_with_type(type)); return get_random_p_id(type, random_index); }; - ParticleChangeRecorder tracker{[this](int p_id) { delete_particle(p_id); }}; for (int i = 0; i < std::min(n_product_types, n_reactant_types); i++) { auto const n_product_coef = reaction.product_coefficients[i]; auto const n_reactant_coef = reaction.reactant_coefficients[i]; // change std::min(reactant_coefficients(i),product_coefficients(i)) many // particles of reactant_types(i) to product_types(i) - auto const type = reaction.reactant_types[i]; + auto const old_type = reaction.reactant_types[i]; + auto const new_type = reaction.product_types[i]; for (int j = 0; j < std::min(n_product_coef, n_reactant_coef); j++) { - auto const p_id = get_random_p_id_of_type(type); - tracker.save_changed_particle({p_id, type, charges_of_types[type]}); - replace_particle(p_id, reaction.product_types[i]); + auto const p_id = get_random_p_id_of_type(old_type); + on_particle_type_change_parallel(p_id, old_type, new_type); + if (auto p = get_local_particle(p_id)) { + p->type() = new_type; +#ifdef ELECTROSTATICS + p->q() = charges_of_types.at(new_type); +#endif + } + bookkeeping.changed.emplace_back(p_id, old_type); } + on_particle_change(); // create product_coefficients(i)-reactant_coefficients(i) many product // particles iff product_coefficients(i)-reactant_coefficients(i)>0, // iff product_coefficients(i)-reactant_coefficients(i)<0, hide this number @@ -196,107 +177,102 @@ ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction) { auto const type = reaction.product_types[i]; for (int j = 0; j < delta_n; j++) { auto const p_id = create_particle(type); - check_exclusion_range(p_id); - tracker.save_created_particle(p_id); + check_exclusion_range(p_id, type); + bookkeeping.created.emplace_back(p_id); } + on_particle_change(); } else if (delta_n < 0) { auto const type = reaction.reactant_types[i]; for (int j = 0; j < -delta_n; j++) { auto const p_id = get_random_p_id_of_type(type); - tracker.save_hidden_particle({p_id, type, charges_of_types[type]}); - check_exclusion_range(p_id); - hide_particle(p_id); + bookkeeping.hidden.emplace_back(p_id, type); + check_exclusion_range(p_id, type); + hide_particle(p_id, type); } + on_particle_change(); } } // create or hide particles of types with noncorresponding replacement types for (auto i = std::min(n_product_types, n_reactant_types); i < std::max(n_product_types, n_reactant_types); i++) { if (n_product_types < n_reactant_types) { - auto const type = reaction.reactant_types[i]; // hide superfluous reactant_types particles + auto const type = reaction.reactant_types[i]; for (int j = 0; j < reaction.reactant_coefficients[i]; j++) { auto const p_id = get_random_p_id_of_type(type); - tracker.save_hidden_particle({p_id, type, charges_of_types[type]}); - check_exclusion_range(p_id); - hide_particle(p_id); + bookkeeping.hidden.emplace_back(p_id, type); + check_exclusion_range(p_id, type); + hide_particle(p_id, type); } + on_particle_change(); } else { // create additional product_types particles + auto const type = reaction.product_types[i]; for (int j = 0; j < reaction.product_coefficients[i]; j++) { - auto const p_id = create_particle(reaction.product_types[i]); - check_exclusion_range(p_id); - tracker.save_created_particle(p_id); + auto const p_id = create_particle(type); + check_exclusion_range(p_id, type); + bookkeeping.created.emplace_back(p_id); } + on_particle_change(); } } - - return tracker; } -std::map ReactionAlgorithm::save_old_particle_numbers( - SingleReaction const ¤t_reaction) const { - std::map old_particle_numbers; +std::unordered_map +ReactionAlgorithm::get_particle_numbers(SingleReaction const &reaction) const { + std::unordered_map particle_numbers; // reactants - for (int type : current_reaction.reactant_types) { - old_particle_numbers[type] = number_of_particles_with_type(type); + for (int type : reaction.reactant_types) { + particle_numbers[type] = ::number_of_particles_with_type(type); } - // products - for (int type : current_reaction.product_types) { - old_particle_numbers[type] = number_of_particles_with_type(type); + for (int type : reaction.product_types) { + particle_numbers[type] = ::number_of_particles_with_type(type); } - return old_particle_numbers; + return particle_numbers; } -void ReactionAlgorithm::generic_oneway_reaction( - SingleReaction ¤t_reaction, double &E_pot_old) { - - current_reaction.tried_moves += 1; +std::optional +ReactionAlgorithm::generic_oneway_reaction_part_1(int reaction_id) { + auto &reaction = *reactions[reaction_id]; + reaction.tried_moves++; particle_inside_exclusion_range_touched = false; - if (!all_reactant_particles_exist(current_reaction)) { - // makes sure, no incomplete reaction is performed -> only need to consider - // rollback of complete reactions - return; + if (!all_reactant_particles_exist(reaction)) { + // make sure that no incomplete reaction is performed -> only need to + // consider rollback of complete reactions + return {}; } - - // find reacting molecules in reactants and save their properties for later - // re-creation if step is not accepted - auto const old_particle_numbers = save_old_particle_numbers(current_reaction); - auto const change_tracker = make_reaction_attempt(current_reaction); - - auto const E_pot_new = (particle_inside_exclusion_range_touched) - ? std::numeric_limits::max() - : mpi_calculate_potential_energy(); - - auto const bf = calculate_acceptance_probability( - current_reaction, E_pot_old, E_pot_new, old_particle_numbers); - - std::vector exponential = {exp(-1.0 / kT * (E_pot_new - E_pot_old))}; - current_reaction.accumulator_potential_energy_difference_exponential( - exponential); - - if (m_uniform_real_distribution(m_generator) < bf) { - // accept - // delete hidden reactant particles (remark: don't delete changed particles) - change_tracker.delete_hidden_particles(); - current_reaction.accepted_moves += 1; - E_pot_old = E_pot_new; // Update the system energy - } else { - // reject - change_tracker.restore_original_state(); + auto &bookkeeping = make_new_system_state(); + bookkeeping.reaction_id = reaction_id; + bookkeeping.old_particle_numbers = get_particle_numbers(reaction); + make_reaction_attempt(reaction, bookkeeping); + auto E_pot_new = std::numeric_limits::max(); + if (not particle_inside_exclusion_range_touched) { + E_pot_new = calculate_potential_energy(); } + return {E_pot_new}; } -/** - * Replaces a particle with the given particle id to be of a certain type. This - * especially means that the particle type and the particle charge are changed. - */ -void ReactionAlgorithm::replace_particle(int p_id, int desired_type) const { - set_particle_type(p_id, desired_type); -#ifdef ELECTROSTATICS - set_particle_q(p_id, charges_of_types.at(desired_type)); -#endif +double ReactionAlgorithm::generic_oneway_reaction_part_2(int reaction_id, + double bf, + double E_pot_old, + double E_pot_new) { + auto const exponential = std::exp(-(E_pot_new - E_pot_old) / kT); + auto &reaction = *reactions[reaction_id]; + reaction.accumulator_potential_energy_difference_exponential( + std::vector{exponential}); + if (get_random_uniform_number() >= bf) { + // reject trial move: restore previous state, energy is unchanged + restore_old_system_state(); + return E_pot_old; + } + // accept trial move: delete hidden particles and return new system energy + for (auto const &[p_id, p_type] : get_old_system_state().hidden) { + delete_particle(p_id); + } + reaction.accepted_moves++; + clear_old_system_state(); + return E_pot_new; } /** @@ -315,62 +291,81 @@ void ReactionAlgorithm::replace_particle(int p_id, int desired_type) const { * there would be a need for a rule for such "collision" reactions (a reaction * like the one above). */ -void ReactionAlgorithm::hide_particle(int p_id) const { - set_particle_type(p_id, non_interacting_type); +void ReactionAlgorithm::hide_particle(int p_id, int p_type) const { + on_particle_type_change_parallel(p_id, p_type, non_interacting_type); + if (auto p = get_local_particle(p_id)) { + p->type() = non_interacting_type; #ifdef ELECTROSTATICS - set_particle_q(p_id, 0.0); + p->q() = 0.; #endif + } } /** * Check if the inserted particle is too close to neighboring particles. */ -void ReactionAlgorithm::check_exclusion_range(int inserted_particle_id) { - - auto const &inserted_particle = get_particle_data(inserted_particle_id); +void ReactionAlgorithm::check_exclusion_range(int p_id, int p_type) { /* Check the exclusion radius of the inserted particle */ - if (exclusion_radius_per_type.count(inserted_particle.type()) != 0) { - if (exclusion_radius_per_type[inserted_particle.type()] == 0.) { + if (exclusion_radius_per_type.count(p_type) != 0) { + if (exclusion_radius_per_type[p_type] == 0.) { return; } } + auto p1_ptr = get_real_particle(p_id); + std::vector particle_ids; if (neighbor_search_order_n) { - particle_ids = get_particle_ids(); + auto all_ids = get_particle_ids_parallel(); /* remove the inserted particle id */ - particle_ids.erase(std::remove(particle_ids.begin(), particle_ids.end(), - inserted_particle_id), - particle_ids.end()); + all_ids.erase(std::remove(all_ids.begin(), all_ids.end(), p_id), + all_ids.end()); + particle_ids = all_ids; } else { - particle_ids = mpi_get_short_range_neighbors(inserted_particle.id(), - m_max_exclusion_range); + auto const local_ids = + mpi_get_short_range_neighbors_local(p_id, m_max_exclusion_range, true); + assert(p1_ptr == nullptr or !!local_ids); + if (local_ids) { + particle_ids = std::move(*local_ids); + } } - /* Check if the inserted particle within the exclusion radius of any other - * particle */ - for (auto const &particle_id : particle_ids) { - auto const &p = get_particle_data(particle_id); - double excluded_distance; - if (exclusion_radius_per_type.count(inserted_particle.type()) == 0 || - exclusion_radius_per_type.count(p.type()) == 0) { - excluded_distance = exclusion_range; - } else if (exclusion_radius_per_type[p.type()] == 0.) { - continue; - } else { - excluded_distance = exclusion_radius_per_type[inserted_particle.type()] + - exclusion_radius_per_type[p.type()]; + if (p1_ptr != nullptr) { + auto &p1 = *p1_ptr; + + /* Check if the inserted particle within the exclusion radius of any other + * particle */ + for (auto const p2_id : particle_ids) { + if (auto const p2_ptr = ::cell_structure.get_local_particle(p2_id)) { + auto const &p2 = *p2_ptr; + double excluded_distance; + if (exclusion_radius_per_type.count(p_type) == 0 || + exclusion_radius_per_type.count(p2.type()) == 0) { + excluded_distance = exclusion_range; + } else if (exclusion_radius_per_type[p2.type()] == 0.) { + continue; + } else { + excluded_distance = exclusion_radius_per_type[p_type] + + exclusion_radius_per_type[p2.type()]; + } + + auto const d_min = ::box_geo.get_mi_vector(p2.pos(), p1.pos()).norm(); + + if (d_min < excluded_distance) { + particle_inside_exclusion_range_touched = true; + break; + } + } } - - auto const d_min = - box_geo.get_mi_vector(p.pos(), inserted_particle.pos()).norm(); - - if (d_min < excluded_distance) { - particle_inside_exclusion_range_touched = true; - break; + if (m_comm.rank() != 0) { + m_comm.send(0, 1, particle_inside_exclusion_range_touched); } + } else if (m_comm.rank() == 0) { + m_comm.recv(boost::mpi::any_source, 1, + particle_inside_exclusion_range_touched); } + boost::mpi::broadcast(m_comm, particle_inside_exclusion_range_touched, 0); } /** @@ -380,10 +375,13 @@ void ReactionAlgorithm::check_exclusion_range(int inserted_particle_id) { * avoid the id range becoming excessively huge. */ void ReactionAlgorithm::delete_particle(int p_id) { - auto const old_max_seen_id = get_maximal_particle_id(); + if (p_id < 0) { + throw std::domain_error("Invalid particle id: " + std::to_string(p_id)); + } + auto const old_max_seen_id = ::get_maximal_particle_id_parallel(); if (p_id == old_max_seen_id) { // last particle, just delete - remove_particle(p_id); + remove_particle_parallel(p_id); // remove all saved empty p_ids which are greater than the max_seen_particle // this is needed in order to avoid the creation of holes for (auto p_id_iter = m_empty_p_ids_smaller_than_max_seen_particle.begin(); @@ -395,7 +393,7 @@ void ReactionAlgorithm::delete_particle(int p_id) { ++p_id_iter; } } else if (p_id <= old_max_seen_id) { - remove_particle(p_id); + remove_particle_parallel(p_id); m_empty_p_ids_smaller_than_max_seen_particle.push_back(p_id); } else { throw std::runtime_error( @@ -463,7 +461,7 @@ Utils::Vector3d ReactionAlgorithm::get_random_position_in_box() { /** * Creates a particle at the end of the observed particle id range. */ -int ReactionAlgorithm::create_particle(int desired_type) { +int ReactionAlgorithm::create_particle(int p_type) { int p_id; if (!m_empty_p_ids_smaller_than_max_seen_particle.empty()) { auto p_id_iter = std::min_element( @@ -472,37 +470,27 @@ int ReactionAlgorithm::create_particle(int desired_type) { p_id = *p_id_iter; m_empty_p_ids_smaller_than_max_seen_particle.erase(p_id_iter); } else { - p_id = get_maximal_particle_id() + 1; + p_id = ::get_maximal_particle_id_parallel() + 1; } - // we use mass=1 for all particles, think about adapting this - auto const new_pos = get_random_position_in_box(); - mpi_make_new_particle(p_id, new_pos); - move_particle(p_id, new_pos, std::sqrt(kT)); - set_particle_type(p_id, desired_type); + // create random velocity vector according to Maxwell-Boltzmann distribution + auto pos = get_random_position_in_box(); + auto vel = get_random_velocity_vector(); + + ::mpi_make_new_particle_local(p_id, pos); + if (auto p = get_local_particle(p_id)) { + p->v() = std::sqrt(kT / p->mass()) * vel; + p->type() = p_type; #ifdef ELECTROSTATICS - set_particle_q(p_id, charges_of_types[desired_type]); + p->q() = charges_of_types.at(p_type); #endif + } + on_particle_type_change_parallel(p_id, ::type_tracking::new_part, p_type); return p_id; } -void ReactionAlgorithm::move_particle(int p_id, Utils::Vector3d const &new_pos, - double velocity_prefactor) { - mpi_set_particle_pos(p_id, new_pos); - // create random velocity vector according to Maxwell-Boltzmann distribution - Utils::Vector3d vel; - vel[0] = velocity_prefactor * m_normal_distribution(m_generator); - vel[1] = velocity_prefactor * m_normal_distribution(m_generator); - vel[2] = velocity_prefactor * m_normal_distribution(m_generator); - set_particle_v(p_id, vel); -} - -std::vector> -ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) { - - std::vector> old_state; - old_state.reserve(n_particles); - +void ReactionAlgorithm::displacement_mc_move(int type, int n_particles) { + auto &bookkeeping = make_new_system_state(); // draw particle ids at random without replacement int p_id = -1; std::vector drawn_pids{p_id}; @@ -513,42 +501,42 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) { p_id = get_random_p_id(type, random_index); } drawn_pids.emplace_back(p_id); - // store original state - auto const &p = get_particle_data(p_id); - old_state.emplace_back(std::tuple{ - p_id, p.pos(), p.v()}); // write new position and new velocity - auto const prefactor = std::sqrt(kT / p.mass()); + typename decltype(ParticleChanges::moved)::value_type old_state; auto const new_pos = get_random_position_in_box(); - move_particle(p_id, new_pos, prefactor); - check_exclusion_range(p_id); + auto vel = get_random_velocity_vector(); + if (auto p = get_real_particle(p_id)) { + old_state = {p_id, p->pos(), p->v()}; + p->v() = std::sqrt(kT / p->mass()) * vel; + if (m_comm.rank() != 0) { + m_comm.send(0, 42, old_state); + } + } else if (m_comm.rank() == 0) { + m_comm.recv(boost::mpi::any_source, 42, old_state); + } + boost::mpi::broadcast(m_comm, old_state, 0); + bookkeeping.moved.emplace_back(old_state); + mpi_set_particle_pos_local(p_id, new_pos); + + check_exclusion_range(p_id, type); if (particle_inside_exclusion_range_touched) { break; } } - - return old_state; } -/** - * Performs a global MC move for particles of a given type. - * Particles are selected without replacement. - * @param type Type of particles to move. - * @param n_part Number of particles of that type to move. - * @returns true if all moves were accepted. - */ -bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type, - int n_part) { +bool ReactionAlgorithm::make_displacement_mc_move_attempt(int type, + int n_particles) { if (type < 0) { throw std::domain_error("Parameter 'type_mc' must be >= 0"); } - if (n_part < 0) { + if (n_particles < 0) { throw std::domain_error( "Parameter 'particle_number_to_be_changed' must be >= 0"); } - if (n_part == 0) { + if (n_particles == 0) { // reject return false; } @@ -556,24 +544,20 @@ bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type, m_tried_configurational_MC_moves += 1; particle_inside_exclusion_range_touched = false; - auto const n_particles_of_type = number_of_particles_with_type(type); - if (n_part > n_particles_of_type) { + auto const n_particles_of_type = ::number_of_particles_with_type(type); + if (n_particles > n_particles_of_type) { // reject return false; } - auto const E_pot_old = mpi_calculate_potential_energy(); - - auto const original_state = generate_new_particle_positions(type, n_part); - + auto const E_pot_old = calculate_potential_energy(); + displacement_mc_move(type, n_particles); auto const E_pot_new = (particle_inside_exclusion_range_touched) ? std::numeric_limits::max() - : mpi_calculate_potential_energy(); - - auto const beta = 1.0 / kT; + : calculate_potential_energy(); // Metropolis algorithm since proposal density is symmetric - auto const bf = std::min(1.0, exp(-beta * (E_pot_new - E_pot_old))); + auto const bf = std::min(1., std::exp(-(E_pot_new - E_pot_old) / kT)); // // correct for enhanced proposal of small radii by using the // // Metropolis-Hastings algorithm for asymmetric proposal densities @@ -589,13 +573,11 @@ bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type, if (m_uniform_real_distribution(m_generator) < bf) { // accept m_accepted_configurational_MC_moves += 1; + clear_old_system_state(); return true; } // reject: restore original particle properties - for (auto const &item : original_state) { - set_particle_v(std::get<0>(item), std::get<2>(item)); - mpi_set_particle_pos(std::get<0>(item), std::get<1>(item)); - } + restore_old_system_state(); return false; } @@ -606,7 +588,7 @@ void ReactionAlgorithm::setup_bookkeeping_of_empty_pids() { // Clean-up the list of empty pids m_empty_p_ids_smaller_than_max_seen_particle.clear(); - auto particle_ids = get_particle_ids(); + auto particle_ids = get_particle_ids_parallel(); std::sort(particle_ids.begin(), particle_ids.end()); auto pid1 = -1; for (auto pid2 : particle_ids) { @@ -617,4 +599,31 @@ void ReactionAlgorithm::setup_bookkeeping_of_empty_pids() { } } +double ReactionAlgorithm::calculate_potential_energy() const { + auto const obs = calculate_energy(); + auto pot = obs->accumulate(-obs->kinetic[0]); + boost::mpi::broadcast(m_comm, pot, 0); + return pot; +} + +Particle *ReactionAlgorithm::get_real_particle(int p_id) const { + assert(p_id >= 0); + auto ptr = ::cell_structure.get_local_particle(p_id); + if (ptr != nullptr and ptr->is_ghost()) { + ptr = nullptr; + } + assert(boost::mpi::all_reduce(m_comm, static_cast(ptr != nullptr), + std::plus<>()) == 1); + return ptr; +} + +Particle *ReactionAlgorithm::get_local_particle(int p_id) const { + assert(p_id >= 0); + auto ptr = ::cell_structure.get_local_particle(p_id); + assert(boost::mpi::all_reduce( + m_comm, static_cast(ptr != nullptr and not ptr->is_ghost()), + std::plus<>()) == 1); + return ptr; +} + } // namespace ReactionMethods diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp index 20c01cdb67b..9d61219a3aa 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.hpp +++ b/src/core/reaction_methods/ReactionAlgorithm.hpp @@ -23,13 +23,14 @@ #include "SingleReaction.hpp" +#include "Particle.hpp" #include "random.hpp" #include #include -#include #include +#include #include #include #include @@ -37,53 +38,23 @@ #include #include -namespace ReactionMethods { - -/** @brief Bookkeeping for particle changes. */ -struct ParticleChangeRecorder { - explicit ParticleChangeRecorder(std::function &&deleter) - : m_deleter{deleter} {} - - struct StoredParticleProperty { - int p_id; - int type; - double charge; - }; - - /** @brief Record particle creation. */ - void save_created_particle(int p_id) { m_created.push_back(p_id); } - - /** @brief Record particle state before it is hidden. */ - void save_hidden_particle(StoredParticleProperty &&item) { - m_hidden.push_back(item); - } - - /** @brief Record particle state before a property change. */ - void save_changed_particle(StoredParticleProperty &&item) { - m_changed.push_back(item); - } +namespace boost::mpi { +class communicator; +} // namespace boost::mpi - /** @brief Delete hidden particles from the system. */ - void delete_hidden_particles() const; - - /** @brief Restore original system state. */ - void restore_original_state() const; - -private: - std::function m_deleter; - std::vector m_created; - std::vector m_hidden; - std::vector m_changed; -}; +namespace ReactionMethods { /** Base class for reaction ensemble methods */ class ReactionAlgorithm { +private: + boost::mpi::communicator const &m_comm; public: ReactionAlgorithm( - int seed, double kT, double exclusion_range, + boost::mpi::communicator const &comm, int seed, double kT, + double exclusion_range, std::unordered_map const &exclusion_radius_per_type) - : kT{kT}, exclusion_range{exclusion_range}, + : m_comm{comm}, kT{kT}, exclusion_range{exclusion_range}, m_generator(Random::mt19937(std::seed_seq({seed, seed, seed}))), m_normal_distribution(0.0, 1.0), m_uniform_real_distribution(0.0, 1.0) { if (kT < 0.) { @@ -99,7 +70,7 @@ class ReactionAlgorithm { virtual ~ReactionAlgorithm() = default; std::vector> reactions; - std::map charges_of_types; + std::unordered_map charges_of_types; double kT; /** * Hard sphere radius. If particles are closer than this value, @@ -147,8 +118,6 @@ class ReactionAlgorithm { m_max_exclusion_range = max_exclusion_range; } - void do_reaction(int reaction_steps); - void check_reaction_method() const; void remove_constraint() { m_reaction_constraint = ReactionConstraint::NONE; } void set_cyl_constraint(double center_x, double center_y, double radius); void set_slab_constraint(double slab_start_z, double slab_end_z); @@ -166,37 +135,86 @@ class ReactionAlgorithm { reactions.erase(reactions.begin() + reaction_id); } - bool displacement_move_for_particles_of_type(int type, int n_part); - bool particle_inside_exclusion_range_touched = false; bool neighbor_search_order_n = true; protected: std::vector m_empty_p_ids_smaller_than_max_seen_particle; + +public: + struct ParticleChanges { + std::vector created{}; + std::vector> changed{}; + std::vector> hidden{}; + std::vector> moved{}; + std::unordered_map old_particle_numbers{}; + int reaction_id{-1}; + }; + + bool is_reaction_under_way() const { return m_system_changes != nullptr; } + auto const &get_old_system_state() const { + if (not is_reaction_under_way()) { + throw std::runtime_error("No chemical reaction is currently under way"); + } + return *m_system_changes; + } + +protected: + /** @brief Restore last valid system state. */ + void restore_old_system_state(); + /** @brief Clear last valid system state. */ + void clear_old_system_state() { + assert(is_reaction_under_way()); + m_system_changes = nullptr; + } + /** @brief Open new handle for system state tracking. */ + auto &make_new_system_state() { + assert(not is_reaction_under_way()); + m_system_changes = std::make_shared(); + return *m_system_changes; + } + std::shared_ptr m_system_changes; + +public: /** - * @brief Carry out a generic one-way chemical reaction. - * - * Generic one way reaction of the type - * A+B+...+G +... --> K+...X + Z +... - * You need to use 2A --> B instead of A+A --> B since - * in the latter you assume distinctness of the particles, however both - * ways to describe the reaction are equivalent in the thermodynamic limit - * (large particle numbers). Furthermore, it is crucial for the function - * in which order you provide the reactant and product types since particles - * will be replaced correspondingly! If there are less reactants than - * products, new product particles are created randomly in the box. - * Matching particles simply change the types. If there are more reactants - * than products, old reactant particles are deleted. - * - * @param[in,out] current_reaction The reaction to attempt. - * @param[in,out] E_pot_old The current potential energy. + * Attempt a reaction move and calculate the new potential energy. + * Particles are selected without replacement. + * @returns Potential energy of the system after the move. + */ + std::optional generic_oneway_reaction_part_1(int reaction_id); + /** + * Accept or reject moves made by @ref generic_oneway_reaction_part_1 based + * on a probability acceptance @c bf. + * @returns Potential energy of the system after the move was accepted or + * rejected. + */ + double generic_oneway_reaction_part_2(int reaction_id, double bf, + double E_pot_old, double E_pot_new); + /** + * Attempt a global MC move for particles of a given type. + * Particles are selected without replacement. + * @param type Type of particles to move. + * @param n_particles Number of particles of that type to move. + * @returns true if all moves were accepted. + */ + bool make_displacement_mc_move_attempt(int type, int n_particles); + /** + * Perform a global MC move for particles of a given type. + * Particles are selected without replacement. + * @param type Type of particles to move. + * @param n_particles Number of particles of that type to move. */ - void generic_oneway_reaction(SingleReaction ¤t_reaction, - double &E_pot_old); + void displacement_mc_move(int type, int n_particles); - ParticleChangeRecorder make_reaction_attempt(SingleReaction const &reaction); - std::vector> - generate_new_particle_positions(int type, int n_particles); + /** @brief Compute the system potential energy. */ + double calculate_potential_energy() const; + +protected: + /** @brief Carry out a chemical reaction and save the old system state. */ + void make_reaction_attempt(::ReactionMethods::SingleReaction const &reaction, + ParticleChanges &bookkeeping); + +public: /** * @brief draws a random integer from the uniform distribution in the range * [0,maxint-1] @@ -207,29 +225,30 @@ class ReactionAlgorithm { std::uniform_int_distribution uniform_int_dist(0, maxint - 1); return uniform_int_dist(m_generator); } + +protected: bool all_reactant_particles_exist(SingleReaction const ¤t_reaction) const; - virtual double - calculate_acceptance_probability(SingleReaction const &, double, double, - std::map const &) const { - return -10.; - } - private: std::mt19937 m_generator; std::normal_distribution m_normal_distribution; std::uniform_real_distribution m_uniform_real_distribution; - std::map - save_old_particle_numbers(SingleReaction const ¤t_reaction) const; + std::unordered_map + get_particle_numbers(SingleReaction const &reaction) const; - void replace_particle(int p_id, int desired_type) const; - int create_particle(int desired_type); - void hide_particle(int p_id) const; - void check_exclusion_range(int inserted_particle_id); - void move_particle(int p_id, Utils::Vector3d const &new_pos, - double velocity_prefactor); + int create_particle(int p_type); + void hide_particle(int p_id, int p_type) const; + void check_exclusion_range(int p_id, int p_type); + auto get_random_uniform_number() { + return m_uniform_real_distribution(m_generator); + } + auto get_random_velocity_vector() { + return Utils::Vector3d{m_normal_distribution(m_generator), + m_normal_distribution(m_generator), + m_normal_distribution(m_generator)}; + } enum class ReactionConstraint { NONE, CYL_Z, SLAB_Z }; ReactionConstraint m_reaction_constraint = ReactionConstraint::NONE; @@ -240,6 +259,9 @@ class ReactionAlgorithm { double m_slab_end_z = -10.0; double m_max_exclusion_range = 0.; + Particle *get_real_particle(int p_id) const; + Particle *get_local_particle(int p_id) const; + protected: Utils::Vector3d get_random_position_in_box(); }; diff --git a/src/core/reaction_methods/ReactionEnsemble.cpp b/src/core/reaction_methods/ReactionEnsemble.cpp deleted file mode 100644 index 06faefe36e2..00000000000 --- a/src/core/reaction_methods/ReactionEnsemble.cpp +++ /dev/null @@ -1,45 +0,0 @@ -/* - * Copyright (C) 2010-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 . - */ - -#include "reaction_methods/ReactionEnsemble.hpp" - -#include "utils.hpp" - -#include -#include - -namespace ReactionMethods { - -/** - * Calculates the expression in the acceptance probability in the reaction - * ensemble - */ -double ReactionEnsemble::calculate_acceptance_probability( - SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, - std::map const &old_particle_numbers) const { - const double factorial_expr = - calculate_factorial_expression(current_reaction, old_particle_numbers); - - const double beta = 1.0 / kT; - // calculate Boltzmann factor - return std::pow(volume, current_reaction.nu_bar) * current_reaction.gamma * - factorial_expr * exp(-beta * (E_pot_new - E_pot_old)); -} - -} // namespace ReactionMethods diff --git a/src/core/reaction_methods/ReactionEnsemble.hpp b/src/core/reaction_methods/ReactionEnsemble.hpp index d7242bd0dc5..f6ef81dd7c1 100644 --- a/src/core/reaction_methods/ReactionEnsemble.hpp +++ b/src/core/reaction_methods/ReactionEnsemble.hpp @@ -21,7 +21,7 @@ #include "reaction_methods/ReactionAlgorithm.hpp" -#include +#include namespace ReactionMethods { @@ -37,16 +37,11 @@ namespace ReactionMethods { class ReactionEnsemble : public ReactionAlgorithm { public: ReactionEnsemble( - int seed, double kT, double exclusion_radius, + boost::mpi::communicator const &comm, int seed, double kT, + double exclusion_radius, const std::unordered_map &exclusion_radius_per_type) - : ReactionAlgorithm(seed, kT, exclusion_radius, + : ReactionAlgorithm(comm, seed, kT, exclusion_radius, exclusion_radius_per_type) {} - -protected: - double calculate_acceptance_probability( - SingleReaction const ¤t_reaction, double E_pot_old, - double E_pot_new, - std::map const &old_particle_numbers) const override; }; } // namespace ReactionMethods diff --git a/src/core/reaction_methods/SingleReaction.hpp b/src/core/reaction_methods/SingleReaction.hpp index 395544ab1ec..1bad5e90be5 100644 --- a/src/core/reaction_methods/SingleReaction.hpp +++ b/src/core/reaction_methods/SingleReaction.hpp @@ -22,6 +22,7 @@ #include #include +#include #include namespace ReactionMethods { @@ -40,6 +41,9 @@ struct SingleReaction { throw std::invalid_argument( "products: number of types and coefficients have to match"); } + if (gamma <= 0.) { + throw std::domain_error("gamma needs to be a strictly positive value"); + } this->reactant_types = reactant_types; this->reactant_coefficients = reactant_coefficients; this->product_types = product_types; diff --git a/src/core/reaction_methods/WidomInsertion.cpp b/src/core/reaction_methods/WidomInsertion.cpp deleted file mode 100644 index 4b17f036688..00000000000 --- a/src/core/reaction_methods/WidomInsertion.cpp +++ /dev/null @@ -1,51 +0,0 @@ -/* - * Copyright (C) 2010-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 . - */ - -#include "reaction_methods/WidomInsertion.hpp" - -#include "energy.hpp" - -#include - -namespace ReactionMethods { - -double WidomInsertion::calculate_particle_insertion_potential_energy( - SingleReaction ¤t_reaction) { - - if (!all_reactant_particles_exist(current_reaction)) - throw std::runtime_error("Trying to remove some non-existing particles " - "from the system via the inverse Widom scheme."); - - auto const E_pot_old = mpi_calculate_potential_energy(); - - // Setup the list of empty pids for bookeeping - setup_bookkeeping_of_empty_pids(); - - // make reaction attempt and immediately reverse it - auto const change_tracker = make_reaction_attempt(current_reaction); - auto const E_pot_new = mpi_calculate_potential_energy(); - change_tracker.restore_original_state(); - - // calculate the particle insertion potential energy - auto const E_pot_insertion = E_pot_new - E_pot_old; - - return E_pot_insertion; -} - -} // namespace ReactionMethods diff --git a/src/core/reaction_methods/WidomInsertion.hpp b/src/core/reaction_methods/WidomInsertion.hpp index 772f54e9e07..d176b84a12b 100644 --- a/src/core/reaction_methods/WidomInsertion.hpp +++ b/src/core/reaction_methods/WidomInsertion.hpp @@ -21,6 +21,7 @@ #include "ReactionAlgorithm.hpp" +#include #include namespace ReactionMethods { @@ -29,12 +30,28 @@ namespace ReactionMethods { class WidomInsertion : public ReactionAlgorithm { public: WidomInsertion( - int seed, double kT, double exclusion_range, + boost::mpi::communicator const &comm, int seed, double kT, + double exclusion_range, const std::unordered_map &exclusion_radius_per_type) - : ReactionAlgorithm(seed, kT, exclusion_range, + : ReactionAlgorithm(comm, seed, kT, exclusion_range, exclusion_radius_per_type) {} - double calculate_particle_insertion_potential_energy( - SingleReaction ¤t_reaction); + + double calculate_particle_insertion_potential_energy(int reaction_id) { + auto &reaction = *reactions[reaction_id]; + if (not all_reactant_particles_exist(reaction)) { + throw std::runtime_error("Trying to remove some non-existing particles " + "from the system via the inverse Widom scheme."); + } + + // make reaction attempt and immediately reverse it + setup_bookkeeping_of_empty_pids(); + auto const E_pot_old = calculate_potential_energy(); + make_reaction_attempt(reaction, make_new_system_state()); + auto const E_pot_new = calculate_potential_energy(); + restore_old_system_state(); + + return E_pot_new - E_pot_old; + } }; } // namespace ReactionMethods diff --git a/src/core/reaction_methods/tests/CMakeLists.txt b/src/core/reaction_methods/tests/CMakeLists.txt index 4993e5b8267..73b6c318afa 100644 --- a/src/core/reaction_methods/tests/CMakeLists.txt +++ b/src/core/reaction_methods/tests/CMakeLists.txt @@ -19,12 +19,8 @@ include(unit_test) unit_test(NAME SingleReaction_test SRC SingleReaction_test.cpp DEPENDS espresso::core) -unit_test(NAME ConstantpHEnsemble_test SRC ConstantpHEnsemble_test.cpp DEPENDS - espresso::core Boost::mpi MPI::MPI_CXX) unit_test(NAME ReactionAlgorithm_test SRC ReactionAlgorithm_test.cpp DEPENDS espresso::core Boost::mpi MPI::MPI_CXX) -unit_test(NAME ReactionEnsemble_test SRC ReactionEnsemble_test.cpp DEPENDS - espresso::core Boost::mpi MPI::MPI_CXX) unit_test(NAME particle_tracking_test SRC particle_tracking_test.cpp DEPENDS espresso::core Boost::mpi MPI::MPI_CXX) unit_test(NAME reaction_methods_utils_test SRC reaction_methods_utils_test.cpp diff --git a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp deleted file mode 100644 index 126cf851f7e..00000000000 --- a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp +++ /dev/null @@ -1,92 +0,0 @@ -/* - * Copyright (C) 2021-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 . - */ - -/* Unit tests for the Constant pH Ensemble. */ - -#define BOOST_TEST_NO_MAIN -#define BOOST_TEST_MODULE Constant pH Ensemble test -#define BOOST_TEST_ALTERNATIVE_INIT_API -#define BOOST_TEST_DYN_LINK -#include - -#include "reaction_methods/ConstantpHEnsemble.hpp" -#include "reaction_methods/utils.hpp" - -#include "communication.hpp" - -#include - -#include -#include -#include -#include -#include - -// Check the Monte Carlo algorithm where moves depend on the system -// configuration, energy and pH. -BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { - using namespace ReactionMethods; - class ConstantpHEnsembleTest : public ConstantpHEnsemble { - public: - using ConstantpHEnsemble::calculate_acceptance_probability; - using ConstantpHEnsemble::ConstantpHEnsemble; - }; - constexpr double tol = 100 * std::numeric_limits::epsilon(); - - ConstantpHEnsembleTest r_algo(42, 20., 0., 1., {}); - - // exception if no reaction was added - BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error); - - // create a reaction A -> B + C - int const type_A = 0; - int const type_B = 1; - int const type_C = 2; - SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {1, 1}); - - // check acceptance probability - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - for (int k = 0; k < 3; ++k) { - // system contains i x A, j x B, and k x C - auto const p_numbers = - std::map{{type_A, i}, {type_B, j}, {type_C, k}}; - auto const energy = static_cast(i + 1); - auto const f_expr = - calculate_factorial_expression_cpH(reaction, p_numbers); - // acceptance = f_expr * exp(- E / T + nu_bar * log(10) * (pH - nu_bar * - // pKa)) - auto const acceptance_ref = - f_expr * std::exp(energy / r_algo.kT + - std::log(10.) * (r_algo.m_constant_pH + - std::log10(reaction.gamma))); - auto const acceptance = r_algo.calculate_acceptance_probability( - reaction, energy, 0., p_numbers); - BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); - } - } - } -} - -int main(int argc, char **argv) { - auto mpi_env = std::make_shared(argc, argv); - Communication::init(mpi_env); - - return boost::unit_test::unit_test_main(init_unit_test, argc, argv); -} diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp index 52d424a2f35..a7cd1a55e62 100644 --- a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp +++ b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp @@ -17,21 +17,23 @@ * along with this program. If not, see . */ -/* Unit tests for the Reaction Algorithm. */ - #define BOOST_TEST_NO_MAIN #define BOOST_TEST_MODULE ReactionAlgorithm test #define BOOST_TEST_ALTERNATIVE_INIT_API #define BOOST_TEST_DYN_LINK #include +#include "config/config.hpp" + #include "reaction_methods/ReactionAlgorithm.hpp" +#include "reaction_methods/SingleReaction.hpp" #include "EspressoSystemStandAlone.hpp" #include "Particle.hpp" +#include "cells.hpp" #include "communication.hpp" -#include "particle_data.hpp" #include "particle_node.hpp" +#include "unit_tests/ParticleFactory.hpp" #include @@ -41,6 +43,7 @@ #include #include #include +#include #include #include @@ -49,20 +52,27 @@ namespace espresso { static std::unique_ptr system; } // namespace espresso +namespace Testing { +class ReactionAlgorithm : public ReactionMethods::ReactionAlgorithm { +public: + using Base = ReactionMethods::ReactionAlgorithm; + using Base::clear_old_system_state; + using Base::displacement_mc_move; + using Base::get_old_system_state; + using Base::get_random_position_in_box; + using Base::make_displacement_mc_move_attempt; + using Base::ReactionAlgorithm; +}; +} // namespace Testing + // Check the base class for all Monte Carlo algorithms. -BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { - using namespace ReactionMethods; - class ReactionAlgorithmTest : public ReactionAlgorithm { - public: - using ReactionAlgorithm::calculate_acceptance_probability; - using ReactionAlgorithm::generate_new_particle_positions; - using ReactionAlgorithm::get_random_position_in_box; - using ReactionAlgorithm::ReactionAlgorithm; - }; - constexpr double tol = 100 * std::numeric_limits::epsilon(); +BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) { + using ReactionMethods::SingleReaction; + auto constexpr tol = 100. * std::numeric_limits::epsilon(); + auto const comm = boost::mpi::communicator(); // check acceptance rate - ReactionAlgorithmTest r_algo(42, 1., 0., {}); + auto r_algo = Testing::ReactionAlgorithm(comm, 42, 1., 0., {}); for (int tried_moves = 1; tried_moves < 5; ++tried_moves) { for (int accepted_moves = 0; accepted_moves < 5; ++accepted_moves) { r_algo.m_tried_configurational_MC_moves = tried_moves; @@ -74,19 +84,12 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { } } - // exception if no reaction was added - BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error); - // create a reaction A -> 3 B + 4 C int const type_A = 0; int const type_B = 1; int const type_C = 2; SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); - - // track particles - init_type_map(type_A); - init_type_map(type_B); - init_type_map(type_C); + r_algo.non_interacting_type = 5; // check reaction addition { @@ -106,26 +109,6 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { BOOST_CHECK_EQUAL(value.gamma, reaction.gamma); } - // check acceptance probability - { - double probability = - r_algo.calculate_acceptance_probability(reaction, -1., -1., {{1, 2}}); - BOOST_CHECK_EQUAL(probability, -10.); - } - -#ifdef ELECTROSTATICS - // exception if reactant types have no charge information - BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error); - r_algo.charges_of_types[0] = 1; - // exception if product types have no charge information - BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error); - r_algo.charges_of_types[1] = 1; - r_algo.charges_of_types[2] = 0; -#endif - - // sanity checks should now pass - r_algo.check_reaction_method(); - // check reaction removal { auto const new_gamma = 5.; @@ -162,25 +145,27 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { BOOST_CHECK(!r_algo.particle_inside_exclusion_range_touched); r_algo.particle_inside_exclusion_range_touched = false; r_algo.exclusion_range = box_l; - auto const bookkeeping = r_algo.generate_new_particle_positions(0, 2); + r_algo.displacement_mc_move(0, 2); + auto const &bookkeeping = r_algo.get_old_system_state(); BOOST_CHECK(r_algo.particle_inside_exclusion_range_touched); // check moves and bookkeeping - for (auto const &item : bookkeeping) { - auto const pid = std::get<0>(item); + for (auto const &[pid, old_pos, old_vel] : bookkeeping.moved) { BOOST_REQUIRE(pid == 0 or pid == 1); auto const ref_old_pos = ref_positions[pid].first; auto const ref_old_vel = ref_positions[pid].second; - auto const &p = get_particle_data(pid); - auto const &new_pos = p.pos(); - auto const &new_vel = p.v(); - BOOST_CHECK_EQUAL(std::get<1>(item), ref_old_pos); - BOOST_CHECK_EQUAL(std::get<2>(item), ref_old_vel); - BOOST_CHECK_GE(new_pos, Utils::Vector3d::broadcast(0.)); - BOOST_CHECK_LE(new_pos, Utils::Vector3d::broadcast(box_l)); - BOOST_CHECK_GE((new_pos - ref_old_pos).norm(), 0.1); - BOOST_CHECK_GE((new_vel - ref_old_vel).norm(), 10.); + if (auto const p = ::cell_structure.get_local_particle(pid)) { + auto const &new_pos = p->pos(); + auto const &new_vel = p->v(); + BOOST_CHECK_EQUAL(old_pos, ref_old_pos); + BOOST_CHECK_EQUAL(old_vel, ref_old_vel); + BOOST_CHECK_GE(new_pos, Utils::Vector3d::broadcast(0.)); + BOOST_CHECK_LE(new_pos, Utils::Vector3d::broadcast(box_l)); + BOOST_CHECK_GE((new_pos - ref_old_pos).norm(), 0.1); + BOOST_CHECK_GE((new_vel - ref_old_vel).norm(), 10.); + } } // cleanup + r_algo.clear_old_system_state(); remove_particle(0); remove_particle(1); } @@ -190,15 +175,15 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { // set up particles auto const box_l = 1.; std::vector ref_positions{{0.1, 0.2, 0.3}, - {0.4, 0.5, 0.6}}; + {0.6, 0.7, 0.8}}; mpi_make_new_particle(0, ref_positions[0]); mpi_make_new_particle(1, ref_positions[1]); set_particle_type(0, type_A); set_particle_type(1, type_A); // check runtime errors when a MC move cannot be physically performed - auto displacement_move = - std::bind(&ReactionAlgorithm::displacement_move_for_particles_of_type, - &r_algo, std::placeholders::_1, std::placeholders::_2); + auto displacement_move = std::bind( + &Testing::ReactionAlgorithm::make_displacement_mc_move_attempt, &r_algo, + std::placeholders::_1, std::placeholders::_2); BOOST_REQUIRE(!displacement_move(type_C, 1)); BOOST_REQUIRE(!displacement_move(type_B, 2)); BOOST_REQUIRE(!displacement_move(type_A, 0)); @@ -207,24 +192,24 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { // force all MC moves to be rejected by picking particles inside // their exclusion radius r_algo.exclusion_range = box_l; - r_algo.particle_inside_exclusion_range_touched = false; - BOOST_REQUIRE(!r_algo.displacement_move_for_particles_of_type(type_A, 2)); + BOOST_REQUIRE(!displacement_move(type_A, 2)); // check none of the particles moved for (auto const pid : {0, 1}) { auto const ref_old_pos = ref_positions[pid]; - auto const &p = get_particle_data(pid); - auto const &new_pos = p.pos(); - BOOST_CHECK_LE((new_pos - ref_old_pos).norm(), tol); + if (auto const p = ::cell_structure.get_local_particle(pid)) { + auto const &new_pos = p->pos(); + BOOST_CHECK_LE((new_pos - ref_old_pos).norm(), tol); + } } // force a MC move to be accepted by using a constant Hamiltonian r_algo.exclusion_range = 0.; - r_algo.particle_inside_exclusion_range_touched = false; - BOOST_REQUIRE(r_algo.displacement_move_for_particles_of_type(type_A, 1)); + BOOST_REQUIRE(displacement_move(type_A, 1)); std::vector distances(2); // check that only one particle moved for (auto const pid : {0, 1}) { - auto const &p = get_particle_data(pid); - distances[pid] = (ref_positions[pid] - p.pos()).norm(); + if (auto const p = ::cell_structure.get_local_particle(pid)) { + distances[pid] = (ref_positions[pid] - p->pos()).norm(); + } } BOOST_CHECK_LE(std::min(distances[0], distances[1]), tol); BOOST_CHECK_GE(std::max(distances[0], distances[1]), 0.1); @@ -245,6 +230,13 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { auto const pos = r_algo.get_random_position_in_box(); BOOST_TEST(pos <= box_l, boost::test_tools::per_element()); BOOST_TEST(pos >= origin, boost::test_tools::per_element()); + // nodes must agree on random numbers + auto pos_head_node = Utils::Vector3d::broadcast(0.); + if (comm.rank() == 0) { + pos_head_node = pos; + } + boost::mpi::broadcast(comm, pos_head_node, 0); + BOOST_TEST(pos == pos_head_node, boost::test_tools::per_element()); } // slab case auto const start_z{0.2}, end_z{0.6}; @@ -286,10 +278,8 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { } { - // domain error if negative exclusion_range is provided - - BOOST_CHECK_THROW(ReactionAlgorithmTest r_algo(40, 1., -1, {}), + BOOST_CHECK_THROW(Testing::ReactionAlgorithm(comm, 40, 1., -1, {}), std::domain_error); // domain error if a negative value is provided in exclusion_radius_per_type @@ -298,11 +288,10 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { exclusion_radius_per_type[type_B] = -1; BOOST_CHECK_THROW( - ReactionAlgorithmTest r_algo(40, 1., 1, exclusion_radius_per_type), + Testing::ReactionAlgorithm(comm, 40, 1., 1, exclusion_radius_per_type), std::domain_error); - auto const box_l = Utils::Vector3d{1, 1, 1}; - espresso::system->set_box_l(box_l); + espresso::system->set_box_l({1., 1., 1.}); // set up particle mpi_make_new_particle(0, {0.5, 0.5, 0.5}); @@ -311,39 +300,35 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { set_particle_type(1, type_B); exclusion_radius_per_type[type_A] = 0.1; exclusion_radius_per_type[type_B] = 1; - ReactionAlgorithmTest r_algo(40, 1., 0, exclusion_radius_per_type); + auto r_algo = + Testing::ReactionAlgorithm(comm, 40, 1., 0, exclusion_radius_per_type); // the new position will always be in the excluded range since the sum of // the radii of both particle types is larger than box length. The exclusion // range value should be ignored - - r_algo.generate_new_particle_positions(type_B, 1); - + r_algo.displacement_mc_move(type_B, 1); BOOST_REQUIRE(r_algo.particle_inside_exclusion_range_touched); + r_algo.clear_old_system_state(); // the new position will never be in the excluded range because the // exclusion_radius of the particle is 0 - r_algo.exclusion_radius_per_type[type_B] = 0; r_algo.particle_inside_exclusion_range_touched = false; - r_algo.generate_new_particle_positions(type_B, 1); - + r_algo.displacement_mc_move(type_B, 1); BOOST_REQUIRE(!r_algo.particle_inside_exclusion_range_touched); + r_algo.clear_old_system_state(); + // the new position will never accepted since the value in exclusion_range // will be used if the particle does not have a defined excluded radius - r_algo.exclusion_range = 1; r_algo.exclusion_radius_per_type = {{type_A, 0}}; - r_algo.generate_new_particle_positions(type_B, 1); - + r_algo.displacement_mc_move(type_B, 1); BOOST_REQUIRE(r_algo.particle_inside_exclusion_range_touched); - - // + r_algo.clear_old_system_state(); } } int main(int argc, char **argv) { espresso::system = std::make_unique(argc, argv); - return boost::unit_test::unit_test_main(init_unit_test, argc, argv); } diff --git a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp deleted file mode 100644 index acfa7bf0866..00000000000 --- a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp +++ /dev/null @@ -1,152 +0,0 @@ -/* - * Copyright (C) 2021-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 . - */ - -/* Unit tests for the Reaction Ensemble. */ - -#define BOOST_TEST_NO_MAIN -#define BOOST_TEST_MODULE Reaction Ensemble test -#define BOOST_TEST_ALTERNATIVE_INIT_API -#define BOOST_TEST_DYN_LINK -#include - -#include "reaction_methods/ReactionEnsemble.hpp" -#include "reaction_methods/utils.hpp" - -#include "EspressoSystemStandAlone.hpp" -#include "communication.hpp" -#include "particle_data.hpp" -#include "particle_node.hpp" - -#include - -#include - -#include -#include -#include -#include -#include - -namespace espresso { -// ESPResSo system instance -static std::unique_ptr system; -} // namespace espresso - -// Check the Monte Carlo algorithm where moves depend on the system -// configuration and energy. -BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { - using namespace ReactionMethods; - class ReactionEnsembleTest : public ReactionEnsemble { - public: - using ReactionEnsemble::calculate_acceptance_probability; - using ReactionEnsemble::generic_oneway_reaction; - using ReactionEnsemble::ReactionEnsemble; - }; - auto constexpr tol = 100 * std::numeric_limits::epsilon(); - - // check basic interface - { - ReactionEnsembleTest r_algo(42, 20., 0., {}); - r_algo.set_volume(10.); - - // exception if no reaction was added - BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error); - - // create a reaction A -> 3 B + 4 C - int const type_A = 0; - int const type_B = 1; - int const type_C = 2; - SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); - - // check acceptance probability - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - for (int k = 0; k < 3; ++k) { - // system contains i x A, j x B, and k x C - auto const p_numbers = - std::map{{type_A, i}, {type_B, j}, {type_C, k}}; - auto const energy = static_cast(i + 1); - auto const f_expr = - calculate_factorial_expression(reaction, p_numbers); - // acceptance = V^{nu_bar} * gamma * f_expr * exp(- E / T) - auto const acceptance_ref = std::pow(r_algo.volume, reaction.nu_bar) * - reaction.gamma * f_expr * - std::exp(energy / r_algo.kT); - auto const acceptance = r_algo.calculate_acceptance_probability( - reaction, energy, 0., p_numbers); - BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); - } - } - } - } - - // check that the system energy is updated after a successful reaction - { - ReactionEnsembleTest test_reaction(42, 1., 0., {}); - test_reaction.set_volume(1.); - - // create a generic identity exchange reaction D <-> E - int const type_D = 0; - int const type_E = 1; - - test_reaction.charges_of_types[type_D] = 0; - test_reaction.charges_of_types[type_E] = 0; - - // track particles - init_type_map(type_D); - init_type_map(type_E); - - auto const gamma = 1e100; // reaction completely shifted to product creation - ReactionMethods::SingleReaction reaction(gamma, {type_D}, {1}, {type_E}, - {1}); - - // resize system box - auto const box_l = Utils::Vector3d{0.5, 0.4, 0.7}; - espresso::system->set_box_l(box_l); - - // create a D particle in the system - auto const pid = 0; - auto const ref_position = Utils::Vector3d{0.1, 0.2, 0.3}; - mpi_make_new_particle(pid, ref_position); - set_particle_type(pid, type_D); - - // sentinel value to check energies get updated - double energy = 1.0; - - // for an ideal system with gamma ~ inf, the reaction is always accepted - test_reaction.generic_oneway_reaction(reaction, energy); - - // the potential energy of the new state has to be 0 since it is an ideal - // system - double const energy_ref = 0.0; - BOOST_CHECK_CLOSE(energy, energy_ref, tol); - - // the reaction was updated - BOOST_CHECK_EQUAL(reaction.tried_moves, 1); - BOOST_CHECK_EQUAL(reaction.accepted_moves, 1); - } -} - -int main(int argc, char **argv) { - auto mpi_env = std::make_shared(argc, argv); - espresso::system = std::make_unique(argc, argv); - Communication::init(mpi_env); - - return boost::unit_test::unit_test_main(init_unit_test, argc, argv); -} diff --git a/src/core/reaction_methods/tests/SingleReaction_test.cpp b/src/core/reaction_methods/tests/SingleReaction_test.cpp index 6d9c6a43e5b..aeae253cec1 100644 --- a/src/core/reaction_methods/tests/SingleReaction_test.cpp +++ b/src/core/reaction_methods/tests/SingleReaction_test.cpp @@ -17,8 +17,6 @@ * along with this program. If not, see . */ -/* Unit tests for the SingleReaction. */ - #define BOOST_TEST_MODULE SingleReaction test #define BOOST_TEST_DYN_LINK #include @@ -27,8 +25,8 @@ #include "reaction_methods/utils.hpp" #include -#include #include +#include // Check a simple chemical reaction, the Monte Carlo acceptance rate // and the configurational move probability for a given system state. @@ -63,7 +61,7 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) { for (int k = 0; k < 3; ++k) { // system contains i x A, j x B, and k x C auto const p_numbers = - std::map{{type_A, i}, {type_B, j}, {type_C, k}}; + std::unordered_map{{type_A, i}, {type_B, j}, {type_C, k}}; auto const val = calculate_factorial_expression(reaction, p_numbers); auto const ref = g(i, -1) * g(j, 3) * g(k, 4); BOOST_CHECK_CLOSE(val, ref, 5 * tol); diff --git a/src/core/reaction_methods/tests/particle_tracking_test.cpp b/src/core/reaction_methods/tests/particle_tracking_test.cpp index 79aecaaacee..c3a79e2396b 100644 --- a/src/core/reaction_methods/tests/particle_tracking_test.cpp +++ b/src/core/reaction_methods/tests/particle_tracking_test.cpp @@ -17,8 +17,6 @@ * along with this program. If not, see . */ -/* Unit tests for the particle tracking mechanism. */ - #define BOOST_TEST_NO_MAIN #define BOOST_TEST_MODULE Particle tracking test #define BOOST_TEST_ALTERNATIVE_INIT_API diff --git a/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp b/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp index 579fd206d34..e2d57e08f1f 100644 --- a/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp +++ b/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp @@ -17,8 +17,6 @@ * along with this program. If not, see . */ -/* Unit tests for the ReactionMethods utility functions. */ - #define BOOST_TEST_MODULE ReactionMethods utility functions test #define BOOST_TEST_DYN_LINK #include @@ -27,12 +25,10 @@ #include #include -#include -#include BOOST_AUTO_TEST_CASE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i_test) { using namespace ReactionMethods; - constexpr double tol = 100 * std::numeric_limits::epsilon(); + constexpr double tol = 100. * std::numeric_limits::epsilon(); auto const reaction_ensemble_combinations = [](int N, int nu) { return (N + nu < 0) ? 0. : std::tgamma(N + 1) / std::tgamma(N + nu + 1); @@ -42,7 +38,7 @@ BOOST_AUTO_TEST_CASE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i_test) { for (int nu = -4; nu <= 4; ++nu) { auto const val = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N0, nu); auto const ref = reaction_ensemble_combinations(N0, nu); - BOOST_CHECK_CLOSE(val, ref, 10 * tol); + BOOST_CHECK_CLOSE(val, ref, 10. * tol); } } } diff --git a/src/core/reaction_methods/utils.cpp b/src/core/reaction_methods/utils.cpp index 47136c82ed7..6e29ee17813 100644 --- a/src/core/reaction_methods/utils.cpp +++ b/src/core/reaction_methods/utils.cpp @@ -17,74 +17,68 @@ * along with this program. If not, see . */ +#include "reaction_methods/SingleReaction.hpp" + #include "reaction_methods/utils.hpp" #include -#include +#include namespace ReactionMethods { double factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(int Ni0, int nu_i) { - double value = 1.0; + auto value = 1.; if (nu_i) { if (nu_i > 0) { for (int i = 1; i <= nu_i; i++) { - value /= Ni0 + i; + value *= static_cast(Ni0 + i); } + value = 1. / value; } else { - auto const abs_nu_i = std::abs(nu_i); - for (int i = 0; i < abs_nu_i; i++) { - value *= Ni0 - i; + for (int i = 0; i < -nu_i; i++) { + value *= static_cast(Ni0 - i); } } } return value; } -double -calculate_factorial_expression(SingleReaction const ¤t_reaction, - std::map const &old_particle_numbers) { - double factorial_expr = 1.0; +double calculate_factorial_expression( + SingleReaction const &reaction, + std::unordered_map const &particle_numbers) { + auto value = 1.; // factorial contribution of reactants - for (int i = 0; i < current_reaction.reactant_types.size(); i++) { - int nu_i = -1 * current_reaction.reactant_coefficients[i]; - int N_i0 = old_particle_numbers.at(current_reaction.reactant_types[i]); - factorial_expr *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i( - N_i0, nu_i); // zeta = 1 (see @cite smith94c) - // since we only perform one reaction - // at one call of the function + for (int i = 0; i < reaction.reactant_types.size(); i++) { + auto const nu_i = -1 * reaction.reactant_coefficients[i]; + auto const N_i0 = particle_numbers.at(reaction.reactant_types[i]); + value *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); } // factorial contribution of products - for (int i = 0; i < current_reaction.product_types.size(); i++) { - int nu_i = current_reaction.product_coefficients[i]; - int N_i0 = old_particle_numbers.at(current_reaction.product_types[i]); - factorial_expr *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i( - N_i0, nu_i); // zeta = 1 (see @cite smith94c) - // since we only perform one reaction - // at one call of the function + for (int i = 0; i < reaction.product_types.size(); i++) { + auto const nu_i = reaction.product_coefficients[i]; + auto const N_i0 = particle_numbers.at(reaction.product_types[i]); + value *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); } - return factorial_expr; + return value; } double calculate_factorial_expression_cpH( - SingleReaction const ¤t_reaction, - std::map const &old_particle_numbers) { - double factorial_expr = 1.0; + SingleReaction const &reaction, + std::unordered_map const &particle_numbers) { + auto value = 1.; // factorial contribution of reactants { - int nu_i = -1 * current_reaction.reactant_coefficients[0]; - int N_i0 = old_particle_numbers.at(current_reaction.reactant_types[0]); - factorial_expr *= - factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); + auto const nu_i = -1 * reaction.reactant_coefficients[0]; + auto const N_i0 = particle_numbers.at(reaction.reactant_types[0]); + value *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); } // factorial contribution of products { - int nu_i = current_reaction.product_coefficients[0]; - int N_i0 = old_particle_numbers.at(current_reaction.product_types[0]); - factorial_expr *= - factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); + auto const nu_i = reaction.product_coefficients[0]; + auto const N_i0 = particle_numbers.at(reaction.product_types[0]); + value *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); } - return factorial_expr; + return value; } } // namespace ReactionMethods diff --git a/src/core/reaction_methods/utils.hpp b/src/core/reaction_methods/utils.hpp index bf033f1bf7c..f49019ed257 100644 --- a/src/core/reaction_methods/utils.hpp +++ b/src/core/reaction_methods/utils.hpp @@ -19,9 +19,9 @@ #ifndef REACTION_METHODS_UTILS_HPP #define REACTION_METHODS_UTILS_HPP -#include "reaction_methods/ReactionAlgorithm.hpp" +#include "reaction_methods/SingleReaction.hpp" -#include +#include namespace ReactionMethods { @@ -31,19 +31,21 @@ namespace ReactionMethods { * * See @cite smith94c. */ -double -calculate_factorial_expression(SingleReaction const ¤t_reaction, - std::map const &old_particle_numbers); +double calculate_factorial_expression( + SingleReaction const &reaction, + std::unordered_map const &particle_numbers); /** * Calculates the factorial expression which occurs in the constant pH method * with symmetric proposal probability. * - * See @cite landsgesell17b + * See @cite landsgesell17b for details. + * zeta = 1 (see @cite smith94c) since we only perform one reaction + * at one call of the function. */ double calculate_factorial_expression_cpH( - SingleReaction const ¤t_reaction, - std::map const &old_particle_numbers); + SingleReaction const &reaction, + std::unordered_map const &particle_numbers); /** * Calculates the factorial expression which occurs in the reaction ensemble diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py index 65bcb7db666..c02445736e7 100644 --- a/src/python/espressomd/reaction_methods.py +++ b/src/python/espressomd/reaction_methods.py @@ -17,27 +17,32 @@ import numpy as np import warnings +import math +import sys from .script_interface import ScriptInterfaceHelper, script_interface_register +from .code_features import has_features from . import utils @script_interface_register class SingleReaction(ScriptInterfaceHelper): _so_name = "ReactionMethods::SingleReaction" - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" + _so_bind_methods = ("get_acceptance_rate",) def __init__(self, **kwargs): super().__init__(**kwargs) if not 'sip' in kwargs: utils.check_valid_keys(self.valid_keys(), kwargs.keys()) + self.reaction_types = tuple( + sorted(set(tuple(self.reactant_types) + tuple(self.product_types)))) def valid_keys(self): - return ("reactant_types", "reactant_coefficients", - "product_types", "product_coefficients", "gamma") + return self.required_keys() def required_keys(self): - return ("reactant_types", "reactant_coefficients", - "product_types", "product_coefficients", "gamma") + return {"reactant_types", "reactant_coefficients", "gamma", + "product_types", "product_coefficients"} def make_backward_reaction(self): return SingleReaction( @@ -47,7 +52,6 @@ def make_backward_reaction(self): product_coefficients=self.reactant_coefficients) -@script_interface_register class ReactionAlgorithm(ScriptInterfaceHelper): """ @@ -67,7 +71,7 @@ class ReactionAlgorithm(ScriptInterfaceHelper): Thermal energy of the system in simulation units exclusion_range : :obj:`float` Minimal distance from any particle, within which new particles will not - be inserted. + be inserted. seed : :obj:`int` Initial counter value (or seed) of the Mersenne Twister RNG. exclusion_radius_per_type : :obj:`dict`, optional @@ -201,14 +205,6 @@ class ReactionAlgorithm(ScriptInterfaceHelper): get_non_interacting_type() Returns the type which is used for hiding particle - reaction() - Performs randomly selected reactions. - - Parameters - ---------- - reaction_steps : :obj:`int`, optional - The number of reactions to be performed at once, defaults to 1. - displacement_mc_move_for_particles_of_type() Performs displacement Monte Carlo moves for particles of a given type. New positions of the displaced particles are chosen from the whole box @@ -254,29 +250,15 @@ class ReactionAlgorithm(ScriptInterfaceHelper): Parameters ---------- reaction_id : :obj:`int` - Reaction id + Identifier of the reaction to modify. + Will be multiplied by 2 internally! gamma : :obj:`float` - New reaction constant - - delete_reaction() - Delete a reaction from the set of used reactions - (the forward and backward reaction). - The ``reaction_id`` which is assigned to a reaction - depends on the order in which :meth:`add_reaction` was called. - The 0th reaction has ``reaction_id=0``, the next added - reaction needs to be addressed with ``reaction_id=1``, etc. - After the deletion of a reaction subsequent reactions - take the ``reaction_id`` of the deleted reaction. - - Parameters - ---------- - reaction_id : :obj:`int` - Reaction id + New reaction constant for the forward reaction. """ _so_name = "ReactionMethods::ReactionAlgorithm" - _so_creation_policy = "LOCAL" + _so_creation_policy = "GLOBAL" _so_bind_methods = ("remove_constraint", "get_wall_constraints_in_z_direction", "set_wall_constraints_in_z_direction", @@ -286,21 +268,22 @@ class ReactionAlgorithm(ScriptInterfaceHelper): "get_acceptance_rate_reaction", "set_non_interacting_type", "get_non_interacting_type", - "reaction", "displacement_mc_move_for_particles_of_type", - "check_reaction_method", "change_reaction_constant", - "delete_reaction", "delete_particle", ) def __init__(self, **kwargs): + if self._so_name == ReactionAlgorithm._so_name: + raise RuntimeError( + f"Base class '{self.__class__.__name__}' cannot be instantiated") if 'exclusion_radius' in kwargs: raise KeyError( 'the keyword `exclusion_radius` is obsolete. Currently, the equivalent keyword is `exclusion_range`') super().__init__(**kwargs) if not 'sip' in kwargs: utils.check_valid_keys(self.valid_keys(), kwargs.keys()) + self._rebuild_reaction_cache() def valid_keys(self): return {"kT", "exclusion_range", "seed", @@ -338,6 +321,8 @@ def add_reaction(self, **kwargs): """ default_charges = kwargs.pop("default_charges") neutrality_check = kwargs.pop("check_for_electroneutrality", True) + if not isinstance(default_charges, dict): + raise TypeError("Argument 'default_charges' needs to be a dict") forward_reaction = SingleReaction(**kwargs) backward_reaction = forward_reaction.make_backward_reaction() if neutrality_check: @@ -345,17 +330,57 @@ def add_reaction(self, **kwargs): type2charge=default_charges, reaction=forward_reaction) + old_default_charges = self.default_charges + for ptype, charge in default_charges.items(): + if ptype in old_default_charges: + if abs(old_default_charges[ptype] - charge) > 1e-10: + raise ValueError( + f"Cannot override charge on particle type {ptype} (from {old_default_charges[ptype]} to {charge})") + for ptype, charge in default_charges.items(): + if ptype not in old_default_charges: + self.call_method( + "set_charge_of_type", + type=ptype, + charge=charge) + self.call_method("add_reaction", reaction=forward_reaction) self.call_method("add_reaction", reaction=backward_reaction) + self.check_reaction_method() + self._rebuild_reaction_cache() - for ptype, charge in default_charges.items(): - self.call_method("set_charge_of_type", type=ptype, charge=charge) - self.call_method("check_reaction_method") + def delete_reaction(self, **kwargs): + """ + Delete a reaction from the set of used reactions + (the forward and backward reaction). + The ``reaction_id`` which is assigned to a reaction + depends on the order in which :meth:`add_reaction` was called. + The 0th reaction has ``reaction_id=0``, the next added + reaction needs to be addressed with ``reaction_id=1``, etc. + After the deletion of a reaction subsequent reactions + take the ``reaction_id`` of the deleted reaction. + + Parameters + ---------- + reaction_id : :obj:`int` + Reaction id + + """ + self.call_method("delete_reaction", **kwargs) + self._rebuild_reaction_cache() + + def check_reaction_method(self): + if len(self.reactions) == 0: + raise RuntimeError("Reaction system not initialized") + + # charges of all reactive types need to be known + if has_features("ELECTROSTATICS"): + for reaction in self.reactions: + for p_type in reaction.reaction_types: + if p_type not in self.default_charges: + raise RuntimeError( + f"Forgot to assign charge to type {p_type}") def _check_charge_neutrality(self, type2charge, reaction): - if not isinstance(type2charge, dict): - raise ValueError( - "No dictionary for relation between types and default charges provided.") charges = np.array(list(type2charge.values())) if np.count_nonzero(charges) == 0: # all particles have zero charge @@ -374,6 +399,102 @@ def _check_charge_neutrality(self, type2charge, reaction): if abs(net_charge_change) / min_abs_nonzero_charge > 1e-10: raise ValueError("Reaction system is not charge neutral") + def _rebuild_reaction_cache(self): + """ + Maintain a local cache of the Python representation of the + script interface reactions. This helps reducing the overhead + when calculating the acceptance probability. + """ + self._reactions_cache = [x for x in self.reactions] + + def reaction(self, reaction_steps): + """ + Performs randomly selected reactions. + + Parameters + ---------- + reaction_steps : :obj:`int`, optional + The number of reactions to be performed at once, defaults to 1. + + """ + self.call_method("setup_bookkeeping_of_empty_pids") + E_pot = self.call_method("potential_energy") + for _ in range(reaction_steps): + reaction_id = self.call_method("get_random_reaction_index") + E_pot = self.generic_oneway_reaction(reaction_id, E_pot) + + def calculate_acceptance_probability(self, reaction_id, E_pot_diff): + """ + Calculate the acceptance probability of a Monte Carlo move. + + Parameters + ---------- + reaction_id : :obj:`int` + Identifier of the reaction that was carried out in the move. + E_pot_diff : :obj:`float` + The potential energy difference for the move. + + Returns + ------- + :obj:`float` + The acceptance probability. + + """ + factorial_expr = self.call_method("calculate_factorial_expression") + reaction = self._reactions_cache[reaction_id] + return self.get_volume()**reaction.nu_bar * reaction.gamma * \ + factorial_expr * math.exp(-E_pot_diff / self.kT) + + def generic_oneway_reaction(self, reaction_id, E_pot_old): + """ + Carry out a generic one-way chemical reaction of the type + `A+B+...+G +... --> K+...X + Z +...` + + You need to use `2A --> B` instead of `A+A --> B` since + in the latter you assume distinctness of the particles, however both + ways to describe the reaction are equivalent in the thermodynamic limit + (large particle numbers). Furthermore, the order of the reactant and + product types matters since particles will be replaced in that order! + If there are less reactants than products, new product particles are + created randomly in the box. Reactants get their type and charge + changed to the corresponding type and charge of the products. + If there are more reactants than products, excess reactant particles + are deleted. + + Parameters + ---------- + reaction_id : :obj:`int` + Identifier of the reaction to attempt. + E_pot_old : :obj:`float` + The current potential energy. + + Returns + ------- + E_pot_new : :obj:`float` + The potential energy after the move. + + """ + try: + E_pot_new = self.call_method( + "generic_oneway_reaction_part_1", reaction_id=reaction_id) + if E_pot_new is None: + return E_pot_old + E_pot_diff = E_pot_new - E_pot_old + bf = self.calculate_acceptance_probability(reaction_id, E_pot_diff) + return self.call_method("generic_oneway_reaction_part_2", + reaction_id=reaction_id, bf=bf, + E_pot_new=E_pot_new, E_pot_old=E_pot_old) + except BaseException as err: + tb = sys.exc_info()[2] + raise RuntimeError( + "An exception was raised by a chemical reaction; the particle " + "state tracking is no longer guaranteed to be correct! -- " + f"{err}").with_traceback(tb) + + def _check_reaction_index(self, reaction_index): + if reaction_index < 0 or reaction_index >= len(self.reactions): + raise IndexError(f"No reaction with id {reaction_index}") + def get_status(self): """ Returns the status of the reaction ensemble in a dictionary containing @@ -381,19 +502,19 @@ def get_status(self): """ - self.call_method("check_reaction_method") + self.check_reaction_method() reactions_list = [] - for core_reaction in self.reactions: - reaction = {"reactant_coefficients": core_reaction.reactant_coefficients, - "reactant_types": core_reaction.reactant_types, - "product_types": core_reaction.product_types, - "product_coefficients": core_reaction.product_coefficients, - "reactant_types": core_reaction.reactant_types, - "gamma": core_reaction.gamma} + for reaction in self.reactions: + reaction = {"reactant_coefficients": reaction.reactant_coefficients, + "reactant_types": reaction.reactant_types, + "product_types": reaction.product_types, + "product_coefficients": reaction.product_coefficients, + "reactant_types": reaction.reactant_types, + "gamma": reaction.gamma} reactions_list.append(reaction) - return {"reactions": reactions_list, "kT": self.kT, + return {"reactions": reactions_list, "kT": self.kT, "exclusion_range": self.exclusion_range} @@ -404,7 +525,6 @@ class ReactionEnsemble(ReactionAlgorithm): """ _so_name = "ReactionMethods::ReactionEnsemble" - _so_creation_policy = "LOCAL" @script_interface_register @@ -423,7 +543,6 @@ class ConstantpHEnsemble(ReactionAlgorithm): """ _so_name = "ReactionMethods::ConstantpHEnsemble" - _so_creation_policy = "LOCAL" def valid_keys(self): return {"kT", "exclusion_range", "seed", @@ -432,6 +551,29 @@ def valid_keys(self): def required_keys(self): return {"kT", "exclusion_range", "seed", "constant_pH"} + def calculate_acceptance_probability(self, reaction_id, E_pot_diff): + """ + Calculate the acceptance probability of a Monte Carlo move. + + Parameters + ---------- + reaction_id : :obj:`int` + Identifier of the reaction that was carried out in the move. + E_pot_diff : :obj:`float` + The potential energy difference for the move. + + Returns + ------- + :obj:`float` + The acceptance probability. + + """ + factorial_expr = self.call_method("calculate_factorial_expression") + reaction = self._reactions_cache[reaction_id] + ln_bf = E_pot_diff - reaction.nu_bar * self.kT * math.log(10.) * ( + self.constant_pH + reaction.nu_bar * math.log10(reaction.gamma)) + return factorial_expr * math.exp(-ln_bf / self.kT) + def add_reaction(self, *args, **kwargs): warn_msg = ( "arguments 'reactant_coefficients' and 'product_coefficients' " @@ -479,7 +621,6 @@ class WidomInsertion(ReactionAlgorithm): """ _so_name = "ReactionMethods::WidomInsertion" - _so_creation_policy = "LOCAL" def required_keys(self): return {"kT", "seed"} @@ -495,7 +636,7 @@ def calculate_particle_insertion_potential_energy(self, **kwargs): """ Measures the potential energy when particles are inserted in the system following the reaction provided in ``reaction_id``. Please - define the insertion moves first by calling the method + define the insertion moves by calling the method :meth:`~ReactionAlgorithm.add_reaction` (with only product types specified). @@ -507,16 +648,19 @@ def calculate_particle_insertion_potential_energy(self, **kwargs): Parameters ---------- reaction_id : :obj:`int` - Reaction identifier. - """ - # make inverse widom scheme (deletion of particles) inaccessible. - # The deletion reactions are the odd reaction_ids + Reaction identifier. Will be multiplied by 2 internally to + skip reverse reactions, i.e. deletion reactions! + + Returns + ------- + :obj:`float` + The particle insertion potential energy. + """ return self.call_method( "calculate_particle_insertion_potential_energy", **kwargs) - def calculate_excess_chemical_potential( - self, **kwargs): + def calculate_excess_chemical_potential(self, **kwargs): """ Given a set of samples of the particle insertion potential energy, calculates the excess chemical potential and its statistical error. @@ -534,6 +678,7 @@ def calculate_excess_chemical_potential( Mean excess chemical potential. error : :obj:`float` Standard error of the mean. + """ def do_block_analysis(samples, N_blocks): diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp index f9248d13a60..d392084113a 100644 --- a/src/script_interface/particle_data/ParticleHandle.cpp +++ b/src/script_interface/particle_data/ParticleHandle.cpp @@ -167,9 +167,7 @@ ParticleHandle::ParticleHandle() { error_msg("type", "must be an integer >= 0")); } make_particle_type_exist_local(new_type); - if (context()->is_head_node()) { - on_particle_type_change_head(m_pid, old_type, new_type); - } + on_particle_type_change_parallel(m_pid, old_type, new_type); set_particle_property(&Particle::type, value); }, [this]() { return get_particle_data(m_pid).type(); }}, diff --git a/src/script_interface/reaction_methods/CMakeLists.txt b/src/script_interface/reaction_methods/CMakeLists.txt index 74036344736..7d17e33e001 100644 --- a/src/script_interface/reaction_methods/CMakeLists.txt +++ b/src/script_interface/reaction_methods/CMakeLists.txt @@ -17,5 +17,7 @@ # along with this program. If not, see . # -target_sources(espresso_script_interface - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/initialize.cpp) +target_sources( + espresso_script_interface + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/initialize.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp) diff --git a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp index 3d6e943417a..4156a98ef89 100644 --- a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp +++ b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp @@ -26,6 +26,7 @@ #include "core/reaction_methods/ConstantpHEnsemble.hpp" #include "core/reaction_methods/ReactionAlgorithm.hpp" +#include "core/reaction_methods/utils.hpp" #include #include @@ -38,6 +39,10 @@ class ConstantpHEnsemble : public ReactionAlgorithm { std::shared_ptr<::ReactionMethods::ReactionAlgorithm> RE() override { return m_re; } + std::shared_ptr<::ReactionMethods::ReactionAlgorithm> const + RE() const override { + return m_re; + } ConstantpHEnsemble() { add_parameters({ @@ -50,21 +55,32 @@ class ConstantpHEnsemble : public ReactionAlgorithm { } void do_construct(VariantMap const ¶ms) override { - m_re = std::make_shared<::ReactionMethods::ConstantpHEnsemble>( - get_value(params, "seed"), get_value(params, "kT"), - get_value(params, "exclusion_range"), - get_value(params, "constant_pH"), - get_value_or>( - params, "exclusion_radius_per_type", {})); + context()->parallel_try_catch([&]() { + m_re = std::make_shared<::ReactionMethods::ConstantpHEnsemble>( + context()->get_comm(), get_value(params, "seed"), + get_value(params, "kT"), + get_value(params, "exclusion_range"), + get_value(params, "constant_pH"), + get_value_or>( + params, "exclusion_radius_per_type", {})); + }); do_set_parameter("search_algorithm", Variant{get_value_or( params, "search_algorithm", "order_n")}); } +protected: + double calculate_factorial_expression( + ::ReactionMethods::SingleReaction const &reaction, + std::unordered_map const &particle_numbers) const override { + return ::ReactionMethods::calculate_factorial_expression_cpH( + reaction, particle_numbers); + } + private: std::shared_ptr<::ReactionMethods::ConstantpHEnsemble> m_re; }; } /* namespace ReactionMethods */ } /* namespace ScriptInterface */ -#endif \ No newline at end of file +#endif diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp new file mode 100644 index 00000000000..fbe51a7a707 --- /dev/null +++ b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp @@ -0,0 +1,248 @@ +/* + * Copyright (C) 2021-2023 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 . + */ + +#include "ReactionAlgorithm.hpp" + +#include "script_interface/ScriptInterface.hpp" + +#include "SingleReaction.hpp" + +#include "config/config.hpp" + +#include "core/cells.hpp" +#include "core/communication.hpp" +#include "core/energy.hpp" +#include "core/event.hpp" +#include "core/particle_node.hpp" +#include "core/reaction_methods/ReactionAlgorithm.hpp" +#include "core/reaction_methods/utils.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace ScriptInterface { +namespace ReactionMethods { + +ReactionAlgorithm::ReactionAlgorithm() { + add_parameters( + {{"reactions", AutoParameter::read_only, + [this]() { + std::vector out; + for (auto const &e : m_reactions) { + out.emplace_back(e); + } + return out; + }}, + {"kT", AutoParameter::read_only, [this]() { return RE()->get_kT(); }}, + {"search_algorithm", + [this](Variant const &v) { + context()->parallel_try_catch([&]() { + auto const key = get_value(v); + if (key == "order_n") { + RE()->neighbor_search_order_n = true; + } else if (key == "parallel") { + RE()->neighbor_search_order_n = false; + } else { + throw std::invalid_argument("Unknown search algorithm '" + key + + "'"); + } + }); + }, + [this]() { + if (RE()->neighbor_search_order_n) { + return std::string("order_n"); + } + return std::string("parallel"); + }}, + {"particle_inside_exclusion_range_touched", + [this](Variant const &v) { + RE()->particle_inside_exclusion_range_touched = get_value(v); + }, + [this]() { return RE()->particle_inside_exclusion_range_touched; }}, + {"default_charges", AutoParameter::read_only, + [this]() { + return make_unordered_map_of_variants(RE()->charges_of_types); + }}, + {"exclusion_range", AutoParameter::read_only, + [this]() { return RE()->get_exclusion_range(); }}, + {"exclusion_radius_per_type", + [this](Variant const &v) { + context()->parallel_try_catch([&]() { + RE()->set_exclusion_radius_per_type( + get_value>(v)); + }); + }, + [this]() { + return make_unordered_map_of_variants( + RE()->exclusion_radius_per_type); + }}}); +} + +static auto get_real_particle(boost::mpi::communicator const &comm, int p_id) { + assert(p_id >= 0); + auto ptr = ::cell_structure.get_local_particle(p_id); + if (ptr != nullptr and ptr->is_ghost()) { + ptr = nullptr; + } + assert(boost::mpi::all_reduce(comm, static_cast(ptr != nullptr), + std::plus<>()) == 1); + return ptr; +} + +Variant ReactionAlgorithm::do_call_method(std::string const &name, + VariantMap const ¶ms) { + if (name == "calculate_factorial_expression") { + if (context()->is_head_node()) { + auto &bookkeeping = RE()->get_old_system_state(); + auto &old_particle_numbers = bookkeeping.old_particle_numbers; + auto &reaction = *m_reactions[bookkeeping.reaction_id]->get_reaction(); + return calculate_factorial_expression(reaction, old_particle_numbers); + } + return {}; + } + if (name == "get_random_reaction_index") { + return RE()->i_random(static_cast(RE()->reactions.size())); + } + if (name == "potential_energy") { + return RE()->calculate_potential_energy(); + } + if (name == "generic_oneway_reaction_part_1") { + auto const reaction_id = get_value(params, "reaction_id"); + Variant result{}; + context()->parallel_try_catch([&]() { + auto const optional = RE()->generic_oneway_reaction_part_1(reaction_id); + if (optional) { + result = *optional; + } + }); + return result; + } + if (name == "generic_oneway_reaction_part_2") { + auto const bf = get_value(params, "bf"); + auto const E_pot_old = get_value(params, "E_pot_old"); + auto const E_pot_new = get_value(params, "E_pot_new"); + auto const reaction_id = get_value(params, "reaction_id"); + Variant result; + context()->parallel_try_catch([&]() { + result = RE()->generic_oneway_reaction_part_2(reaction_id, bf, E_pot_old, + E_pot_new); + }); + return result; + } + if (name == "setup_bookkeeping_of_empty_pids") { + RE()->setup_bookkeeping_of_empty_pids(); + } else if (name == "remove_constraint") { + RE()->remove_constraint(); + } else if (name == "set_cylindrical_constraint_in_z_direction") { + context()->parallel_try_catch([&]() { + RE()->set_cyl_constraint(get_value(params, "center_x"), + get_value(params, "center_y"), + get_value(params, "radius")); + }); + } else if (name == "set_wall_constraints_in_z_direction") { + context()->parallel_try_catch([&]() { + RE()->set_slab_constraint(get_value(params, "slab_start_z"), + get_value(params, "slab_end_z")); + }); + } else if (name == "get_wall_constraints_in_z_direction") { + if (context()->is_head_node()) { + return RE()->get_slab_constraint_parameters(); + } + return {}; + } else if (name == "set_volume") { + context()->parallel_try_catch( + [&]() { RE()->set_volume(get_value(params, "volume")); }); + } else if (name == "get_volume") { + return RE()->get_volume(); + } else if (name == "get_acceptance_rate_reaction") { + auto const index = get_value(params, "reaction_id"); + context()->parallel_try_catch([&]() { + if (index < 0 or index >= static_cast(m_reactions.size())) { + throw std::out_of_range("No reaction with id " + std::to_string(index)); + } + }); + return m_reactions[index]->get_reaction()->get_acceptance_rate(); + } else if (name == "set_non_interacting_type") { + auto const type = get_value(params, "type"); + context()->parallel_try_catch([&]() { + if (type < 0) { + throw std::domain_error("Invalid type: " + std::to_string(type)); + } + }); + RE()->non_interacting_type = type; + ::init_type_map(type); + } else if (name == "get_non_interacting_type") { + return RE()->non_interacting_type; + } else if (name == "displacement_mc_move_for_particles_of_type") { + auto const type = get_value(params, "type_mc"); + auto const n_particles = + get_value_or(params, "particle_number_to_be_changed", 1); + auto result = false; + context()->parallel_try_catch([&]() { + result = RE()->make_displacement_mc_move_attempt(type, n_particles); + }); + return result; + } else if (name == "delete_particle") { + context()->parallel_try_catch( + [&]() { RE()->delete_particle(get_value(params, "p_id")); }); + } else if (name == "delete_reaction") { + auto const reaction_id = get_value(params, "reaction_id"); + auto const index = get_reaction_index(reaction_id); + // delete forward and backward reactions + delete_reaction(index + 1); + delete_reaction(index + 0); + } else if (name == "add_reaction") { + auto const reaction = + get_value>(params, "reaction"); + m_reactions.push_back(reaction); + RE()->add_reaction(reaction->get_reaction()); + } else if (name == "change_reaction_constant") { + auto const gamma = get_value(params, "gamma"); + auto const reaction_id = get_value(params, "reaction_id"); + context()->parallel_try_catch([&]() { + if (reaction_id % 2 == 1) { + throw std::invalid_argument("Only forward reactions can be selected"); + } + if (gamma <= 0.) { + throw std::domain_error("gamma needs to be a strictly positive value"); + } + }); + auto const index = get_reaction_index(reaction_id); + m_reactions[index + 0]->get_reaction()->gamma = gamma; + m_reactions[index + 1]->get_reaction()->gamma = 1. / gamma; + } else if (name == "set_charge_of_type") { + auto const type = get_value(params, "type"); + auto const charge = get_value(params, "charge"); + RE()->charges_of_types[type] = charge; + } else if (context()->is_head_node()) { + throw std::runtime_error("unknown method '" + name + "()'"); + } + return {}; +} + +} /* namespace ReactionMethods */ +} /* namespace ScriptInterface */ diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp index 8029552df9c..ae3fa822fa3 100644 --- a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp +++ b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp @@ -23,8 +23,10 @@ #include "SingleReaction.hpp" #include "script_interface/ScriptInterface.hpp" +#include "script_interface/particle_data/ParticleHandle.hpp" #include "core/reaction_methods/ReactionAlgorithm.hpp" +#include "core/reaction_methods/utils.hpp" #include #include @@ -47,119 +49,31 @@ class ReactionAlgorithm : public AutoParameters { */ int get_reaction_index(int reaction_id) const { auto const index = 2 * reaction_id; - if (index < 0 or index >= static_cast(m_reactions.size())) { - throw std::out_of_range("This reaction is not present"); - } + context()->parallel_try_catch([&]() { + if (index < 0 or index >= static_cast(m_reactions.size())) { + throw std::out_of_range("No reaction with id " + + std::to_string(reaction_id)); + } + }); return index; } public: virtual std::shared_ptr<::ReactionMethods::ReactionAlgorithm> RE() = 0; + virtual std::shared_ptr<::ReactionMethods::ReactionAlgorithm> const + RE() const = 0; - ReactionAlgorithm() { - add_parameters( - {{"reactions", AutoParameter::read_only, - [this]() { - std::vector out; - for (auto const &e : m_reactions) { - out.emplace_back(e); - } - return out; - }}, - {"kT", AutoParameter::read_only, [this]() { return RE()->get_kT(); }}, - {"search_algorithm", - [this](Variant const &v) { - auto const key = get_value(v); - if (key == "order_n") { - RE()->neighbor_search_order_n = true; - } else if (key == "parallel") { - RE()->neighbor_search_order_n = false; - } else { - throw std::invalid_argument("Unknown search algorithm '" + key + - "'"); - } - }, - [this]() { - if (RE()->neighbor_search_order_n) { - return std::string("order_n"); - } - return std::string("parallel"); - }}, - {"exclusion_range", AutoParameter::read_only, - [this]() { return RE()->get_exclusion_range(); }}, - {"exclusion_radius_per_type", - [this](Variant const &v) { - RE()->set_exclusion_radius_per_type( - get_value>(v)); - }, - [this]() { - return make_unordered_map_of_variants( - RE()->exclusion_radius_per_type); - }}}); - } + ReactionAlgorithm(); Variant do_call_method(std::string const &name, - VariantMap const ¶meters) override { - if (name == "remove_constraint") { - RE()->remove_constraint(); - } else if (name == "set_cylindrical_constraint_in_z_direction") { - RE()->set_cyl_constraint(get_value(parameters, "center_x"), - get_value(parameters, "center_y"), - get_value(parameters, "radius")); - } else if (name == "set_wall_constraints_in_z_direction") { - RE()->set_slab_constraint(get_value(parameters, "slab_start_z"), - get_value(parameters, "slab_end_z")); - } else if (name == "get_wall_constraints_in_z_direction") { - return RE()->get_slab_constraint_parameters(); - } else if (name == "set_volume") { - RE()->set_volume(get_value(parameters, "volume")); - } else if (name == "get_volume") { - return RE()->get_volume(); - } else if (name == "get_acceptance_rate_reaction") { - auto const index = get_value(parameters, "reaction_id"); - if (index < 0 or index >= static_cast(m_reactions.size())) { - throw std::out_of_range("This reaction is not present"); - } - return m_reactions[index]->get_reaction()->get_acceptance_rate(); - } else if (name == "set_non_interacting_type") { - RE()->non_interacting_type = get_value(parameters, "type"); - } else if (name == "get_non_interacting_type") { - return RE()->non_interacting_type; - } else if (name == "reaction") { - RE()->do_reaction(get_value_or(parameters, "reaction_steps", 1)); - } else if (name == "displacement_mc_move_for_particles_of_type") { - return RE()->displacement_move_for_particles_of_type( - get_value(parameters, "type_mc"), - get_value_or(parameters, "particle_number_to_be_changed", 1)); - } else if (name == "check_reaction_method") { - RE()->check_reaction_method(); - } else if (name == "delete_particle") { - RE()->delete_particle(get_value(parameters, "p_id")); - } else if (name == "delete_reaction") { - auto const reaction_id = get_value(parameters, "reaction_id"); - auto const index = get_reaction_index(reaction_id); - // delete forward and backward reactions - delete_reaction(index + 1); - delete_reaction(index + 0); - } else if (name == "add_reaction") { - auto const reaction = - get_value>(parameters, "reaction"); - m_reactions.push_back(reaction); - RE()->add_reaction(reaction->get_reaction()); - } else if (name == "change_reaction_constant") { - auto const gamma = get_value(parameters, "gamma"); - auto const reaction_id = get_value(parameters, "reaction_id"); - auto const index = get_reaction_index(reaction_id); - m_reactions[index + 0]->get_reaction()->gamma = gamma; - m_reactions[index + 1]->get_reaction()->gamma = 1. / gamma; - } else if (name == "set_charge_of_type") { - auto const type = get_value(parameters, "type"); - auto const charge = get_value(parameters, "charge"); - RE()->charges_of_types[type] = charge; - } else { - throw std::runtime_error(("unknown method '" + name + "()'").c_str()); - } - return {}; + VariantMap const ¶ms) override; + +protected: + virtual double calculate_factorial_expression( + ::ReactionMethods::SingleReaction const &reaction, + std::unordered_map const &particle_numbers) const { + return ::ReactionMethods::calculate_factorial_expression(reaction, + particle_numbers); } private: @@ -169,10 +83,12 @@ class ReactionAlgorithm : public AutoParameters { } std::string get_internal_state() const override { + assert(not RE()->is_reaction_under_way()); throw std::runtime_error("Reaction methods do not support checkpointing"); } }; + } /* namespace ReactionMethods */ } /* namespace ScriptInterface */ -#endif \ No newline at end of file +#endif diff --git a/src/script_interface/reaction_methods/ReactionEnsemble.hpp b/src/script_interface/reaction_methods/ReactionEnsemble.hpp index a520c95f6af..6d2ea320dd2 100644 --- a/src/script_interface/reaction_methods/ReactionEnsemble.hpp +++ b/src/script_interface/reaction_methods/ReactionEnsemble.hpp @@ -38,13 +38,20 @@ class ReactionEnsemble : public ReactionAlgorithm { std::shared_ptr<::ReactionMethods::ReactionAlgorithm> RE() override { return m_re; } + std::shared_ptr<::ReactionMethods::ReactionAlgorithm> const + RE() const override { + return m_re; + } void do_construct(VariantMap const ¶ms) override { - m_re = std::make_shared<::ReactionMethods::ReactionEnsemble>( - get_value(params, "seed"), get_value(params, "kT"), - get_value(params, "exclusion_range"), - get_value_or>( - params, "exclusion_radius_per_type", {})); + context()->parallel_try_catch([&]() { + m_re = std::make_shared<::ReactionMethods::ReactionEnsemble>( + context()->get_comm(), get_value(params, "seed"), + get_value(params, "kT"), + get_value(params, "exclusion_range"), + get_value_or>( + params, "exclusion_radius_per_type", {})); + }); do_set_parameter("search_algorithm", Variant{get_value_or( params, "search_algorithm", "order_n")}); @@ -56,4 +63,4 @@ class ReactionEnsemble : public ReactionAlgorithm { } /* namespace ReactionMethods */ } /* namespace ScriptInterface */ -#endif \ No newline at end of file +#endif diff --git a/src/script_interface/reaction_methods/SingleReaction.hpp b/src/script_interface/reaction_methods/SingleReaction.hpp index c821ef87734..264b9007f92 100644 --- a/src/script_interface/reaction_methods/SingleReaction.hpp +++ b/src/script_interface/reaction_methods/SingleReaction.hpp @@ -20,9 +20,12 @@ #ifndef SCRIPT_INTERFACE_REACTION_METHODS_SINGLE_REACTION_HPP #define SCRIPT_INTERFACE_REACTION_METHODS_SINGLE_REACTION_HPP -#include "core/reaction_methods/SingleReaction.hpp" #include "script_interface/ScriptInterface.hpp" -#include + +#include "core/reaction_methods/SingleReaction.hpp" + +#include +#include #include namespace ScriptInterface { @@ -32,6 +35,7 @@ class SingleReaction : public AutoParameters { public: SingleReaction() { add_parameters({ + {"nu_bar", AutoParameter::read_only, [this]() { return m_sr->nu_bar; }}, {"gamma", AutoParameter::read_only, [this]() { return m_sr->gamma; }}, {"reactant_types", AutoParameter::read_only, [this]() { return m_sr->reactant_types; }}, @@ -45,22 +49,32 @@ class SingleReaction : public AutoParameters { } void do_construct(VariantMap const ¶ms) override { - m_sr = std::make_shared<::ReactionMethods::SingleReaction>( - get_value(params, "gamma"), - get_value>(params, "reactant_types"), - get_value>(params, "reactant_coefficients"), - get_value>(params, "product_types"), - get_value>(params, "product_coefficients")); + context()->parallel_try_catch([&]() { + m_sr = std::make_shared<::ReactionMethods::SingleReaction>( + get_value(params, "gamma"), + get_value>(params, "reactant_types"), + get_value>(params, "reactant_coefficients"), + get_value>(params, "product_types"), + get_value>(params, "product_coefficients")); + }); } std::shared_ptr<::ReactionMethods::SingleReaction> get_reaction() { return m_sr; } + Variant do_call_method(std::string const &name, VariantMap const &) override { + if (name == "get_acceptance_rate") { + return m_sr->get_acceptance_rate(); + } + return {}; + } + private: std::shared_ptr<::ReactionMethods::SingleReaction> m_sr; }; + } /* namespace ReactionMethods */ } /* namespace ScriptInterface */ -#endif \ No newline at end of file +#endif diff --git a/src/script_interface/reaction_methods/WidomInsertion.hpp b/src/script_interface/reaction_methods/WidomInsertion.hpp index f90c0aef834..bd660f64187 100644 --- a/src/script_interface/reaction_methods/WidomInsertion.hpp +++ b/src/script_interface/reaction_methods/WidomInsertion.hpp @@ -39,37 +39,50 @@ class WidomInsertion : public ReactionAlgorithm { std::shared_ptr<::ReactionMethods::ReactionAlgorithm> RE() override { return m_re; } + std::shared_ptr<::ReactionMethods::ReactionAlgorithm> const + RE() const override { + return m_re; + } WidomInsertion() { add_parameters({{"search_algorithm", - [](Variant const &) { - throw std::runtime_error( - "No search algorithm for WidomInsertion"); + [this](Variant const &) { + if (context()->is_head_node()) { + throw std::runtime_error( + "No search algorithm for WidomInsertion"); + } }, []() { return none; }}}); } void do_construct(VariantMap const ¶ms) override { - m_re = std::make_shared<::ReactionMethods::WidomInsertion>( - get_value(params, "seed"), get_value(params, "kT"), 0., - std::unordered_map{}); + context()->parallel_try_catch([&]() { + m_re = std::make_shared<::ReactionMethods::WidomInsertion>( + context()->get_comm(), get_value(params, "seed"), + get_value(params, "kT"), 0., + std::unordered_map{}); + }); } Variant do_call_method(std::string const &name, - VariantMap const ¶meters) override { + VariantMap const ¶ms) override { if (name == "calculate_particle_insertion_potential_energy") { - auto const reaction_id = get_value(parameters, "reaction_id"); - auto const index = get_reaction_index(reaction_id); - auto &reaction = *m_reactions[index]->get_reaction(); - return m_re->calculate_particle_insertion_potential_energy(reaction); + Variant result; + context()->parallel_try_catch([&]() { + auto const reaction_id = get_value(params, "reaction_id"); + auto const index = get_reaction_index(reaction_id); + result = m_re->calculate_particle_insertion_potential_energy(index); + }); + return result; } - return ReactionAlgorithm::do_call_method(name, parameters); + return ReactionAlgorithm::do_call_method(name, params); } private: std::shared_ptr<::ReactionMethods::WidomInsertion> m_re; }; + } /* namespace ReactionMethods */ } /* namespace ScriptInterface */ -#endif \ No newline at end of file +#endif diff --git a/src/script_interface/system/System.cpp b/src/script_interface/system/System.cpp index 5b00ffb1b49..20e490baea1 100644 --- a/src/script_interface/system/System.cpp +++ b/src/script_interface/system/System.cpp @@ -96,20 +96,15 @@ Variant System::do_call_method(std::string const &name, return {}; } if (name == "setup_type_map") { - if (context()->is_head_node()) { - auto const types = get_value>(parameters, "type_list"); - for (auto const type : types) { - init_type_map(type); - } + auto const types = get_value>(parameters, "type_list"); + for (auto const type : types) { + ::init_type_map(type); } return {}; } if (name == "number_of_particles") { - if (context()->is_head_node()) { - auto const type = get_value(parameters, "type"); - return number_of_particles_with_type(type); - } - return {}; + auto const type = get_value(parameters, "type"); + return ::number_of_particles_with_type(type); } if (name == "velocity_difference") { auto const pos1 = get_value(parameters, "pos1"); diff --git a/src/script_interface/tests/CMakeLists.txt b/src/script_interface/tests/CMakeLists.txt index 29827d64772..55a699da218 100644 --- a/src/script_interface/tests/CMakeLists.txt +++ b/src/script_interface/tests/CMakeLists.txt @@ -56,3 +56,7 @@ unit_test(NAME Constraints_test SRC Constraints_test.cpp DEPENDS espresso::script_interface espresso::core) unit_test(NAME Actors_test SRC Actors_test.cpp DEPENDS espresso::script_interface espresso::core) +unit_test(NAME ConstantpHEnsemble_test SRC ConstantpHEnsemble_test.cpp DEPENDS + espresso::core espresso::script_interface Boost::mpi MPI::MPI_CXX) +unit_test(NAME ReactionEnsemble_test SRC ReactionEnsemble_test.cpp DEPENDS + espresso::core espresso::script_interface Boost::mpi MPI::MPI_CXX) diff --git a/src/script_interface/tests/ConstantpHEnsemble_test.cpp b/src/script_interface/tests/ConstantpHEnsemble_test.cpp new file mode 100644 index 00000000000..ec3d485f7a9 --- /dev/null +++ b/src/script_interface/tests/ConstantpHEnsemble_test.cpp @@ -0,0 +1,135 @@ +/* + * Copyright (C) 2021-2023 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 . + */ + +#define BOOST_TEST_NO_MAIN +#define BOOST_TEST_MODULE ConstantpHEnsemble test +#define BOOST_TEST_ALTERNATIVE_INIT_API +#define BOOST_TEST_DYN_LINK +#include + +#include "core/reaction_methods/ConstantpHEnsemble.hpp" +#include "core/reaction_methods/SingleReaction.hpp" +#include "script_interface/LocalContext.hpp" +#include "script_interface/Variant.hpp" +#include "script_interface/reaction_methods/ConstantpHEnsemble.hpp" + +#include "core/EspressoSystemStandAlone.hpp" +#include "core/Particle.hpp" +#include "core/particle_node.hpp" +#include "core/unit_tests/ParticleFactory.hpp" + +#include + +#include +#include +#include +#include + +using ::ReactionMethods::SingleReaction; + +namespace espresso { +// ESPResSo system instance +static std::unique_ptr system; +} // namespace espresso + +namespace ScriptInterface::Testing { +class ConstantpHEnsemble + : public ScriptInterface::ReactionMethods::ConstantpHEnsemble { +public: + auto calculate_acceptance_probability( + SingleReaction const &reaction, double E_pot_diff, + std::unordered_map const &old_particle_numbers) const { + auto const factorial_expr = + ::ReactionMethods::calculate_factorial_expression_cpH( + reaction, old_particle_numbers); + auto const pH = + std::dynamic_pointer_cast<::ReactionMethods::ConstantpHEnsemble>(RE()) + ->m_constant_pH; + auto const ln_bf = + E_pot_diff - reaction.nu_bar * RE()->kT * std::log(10.) * + (pH + reaction.nu_bar * std::log10(reaction.gamma)); + return factorial_expr * std::exp(-ln_bf / RE()->kT); + } +}; +} // namespace ScriptInterface::Testing + +// Check the Monte Carlo algorithm where moves depend on the system +// configuration, energy and pH. +BOOST_FIXTURE_TEST_CASE(ConstantpHEnsemble_test, ParticleFactory) { + using namespace ScriptInterface; + auto constexpr tol = 100. * std::numeric_limits::epsilon(); + + Utils::Factory factory; + factory.register_new( + "Testing::ConstantpHEnsemble"); + + auto const comm = boost::mpi::communicator(); + auto const make_algo = [&factory, + &comm](int seed, double kT, double pH, + double exclusion_range, + std::unordered_map const &radii) { + auto ctx = std::make_shared(factory, comm); + VariantMap params{}; + params["seed"] = seed; + params["kT"] = kT; + params["constant_pH"] = pH; + params["exclusion_range"] = exclusion_range; + params["exclusion_radius_per_type"] = make_unordered_map_of_variants(radii); + auto &&sp = ctx->make_shared_local("Testing::ConstantpHEnsemble", params); + return std::dynamic_pointer_cast(sp); + }; + + auto const constant_pH = 1.; + auto &&r_algo_si = make_algo(42, 20., constant_pH, 0., {}); + auto &r_algo = *r_algo_si->RE(); + r_algo.non_interacting_type = 5; + + // create a reaction A -> B + C + int const type_A = 0; + int const type_B = 1; + int const type_C = 2; + SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {1, 1}); + + // check acceptance probability + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // system contains i x A, j x B, and k x C + auto const p_numbers = + std::unordered_map{{type_A, i}, {type_B, j}, {type_C, k}}; + auto const energy = -static_cast(i + 1); + auto const f_expr = + calculate_factorial_expression_cpH(reaction, p_numbers); + // bf = f_expr * exp(- E / kT + nu_bar * log(10) * (pH - nu_bar * pKa)) + auto const acceptance_ref = + f_expr * std::exp(-energy / r_algo.kT + + std::log(10.) * + (constant_pH + std::log10(reaction.gamma))); + auto const acceptance = r_algo_si->calculate_acceptance_probability( + reaction, energy, p_numbers); + BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5. * tol); + } + } + } +} + +int main(int argc, char **argv) { + espresso::system = std::make_unique(argc, argv); + return boost::unit_test::unit_test_main(init_unit_test, argc, argv); +} diff --git a/src/script_interface/tests/ReactionEnsemble_test.cpp b/src/script_interface/tests/ReactionEnsemble_test.cpp new file mode 100644 index 00000000000..664d819fcbf --- /dev/null +++ b/src/script_interface/tests/ReactionEnsemble_test.cpp @@ -0,0 +1,266 @@ +/* + * Copyright (C) 2021-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 . + */ + +#define BOOST_TEST_NO_MAIN +#define BOOST_TEST_MODULE ReactionEnsemble test +#define BOOST_TEST_ALTERNATIVE_INIT_API +#define BOOST_TEST_DYN_LINK +#include + +#include "core/reaction_methods/ReactionEnsemble.hpp" +#include "core/reaction_methods/SingleReaction.hpp" +#include "script_interface/LocalContext.hpp" +#include "script_interface/None.hpp" +#include "script_interface/Variant.hpp" +#include "script_interface/reaction_methods/ReactionEnsemble.hpp" + +#include "core/EspressoSystemStandAlone.hpp" +#include "core/Particle.hpp" +#include "core/particle_node.hpp" +#include "core/unit_tests/ParticleFactory.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using ::ReactionMethods::SingleReaction; + +namespace espresso { +// ESPResSo system instance +static std::unique_ptr system; +} // namespace espresso + +namespace ScriptInterface::Testing { +class ReactionEnsemble + : public ScriptInterface::ReactionMethods::ReactionEnsemble { +public: + auto calculate_acceptance_probability( + SingleReaction const &reaction, double E_pot_diff, + std::unordered_map const &old_particle_numbers) const { + auto const factorial_expr = + ::ReactionMethods::calculate_factorial_expression(reaction, + old_particle_numbers); + return std::pow(RE()->get_volume(), reaction.nu_bar) * reaction.gamma * + factorial_expr * std::exp(-E_pot_diff / RE()->kT); + } +}; +} // namespace ScriptInterface::Testing + +// Check the Monte Carlo algorithm where moves depend on the system +// configuration and energy. +BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { + using namespace ScriptInterface; + auto constexpr tol = 100. * std::numeric_limits::epsilon(); + + Utils::Factory factory; + factory.register_new("Testing::ReactionEnsemble"); + factory.register_new( + "SingleReaction"); + + auto const comm = boost::mpi::communicator(); + auto const make_algo = [&factory, + &comm](int seed, double kT, double exclusion_range, + std::unordered_map const &radii) { + using namespace ScriptInterface; + auto ctx = std::make_shared(factory, comm); + VariantMap params{}; + params["seed"] = seed; + params["kT"] = kT; + params["exclusion_range"] = exclusion_range; + params["exclusion_radius_per_type"] = make_unordered_map_of_variants(radii); + auto &&sp = ctx->make_shared_local("Testing::ReactionEnsemble", params); + return std::dynamic_pointer_cast(sp); + }; + auto const make_reaction = + [&factory, &comm](double gamma, std::vector const &reactant_types, + std::vector const &reactant_coefficients, + std::vector const &product_types, + std::vector const &product_coefficients) { + using namespace ScriptInterface; + auto ctx = + std::make_shared(factory, comm); + VariantMap params{}; + params["gamma"] = gamma; + params["reactant_types"] = make_vector_of_variants(reactant_types); + params["reactant_coefficients"] = + make_vector_of_variants(reactant_coefficients); + params["product_types"] = make_vector_of_variants(product_types); + params["product_coefficients"] = + make_vector_of_variants(product_coefficients); + auto &&si_obj = ctx->make_shared_local("SingleReaction", params); + return std::dynamic_pointer_cast< + ScriptInterface::ReactionMethods::SingleReaction>(si_obj); + }; + + // check basic interface + { + auto &&r_algo_si = make_algo(42, 20., 0., {}); + auto &r_algo = *r_algo_si->RE(); + r_algo.set_volume(10.); + r_algo.non_interacting_type = 5; + + // create a reaction A -> 3 B + 4 C + int const type_A = 0; + int const type_B = 1; + int const type_C = 2; + SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); + + // check acceptance probability + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // system contains i x A, j x B, and k x C + auto const p_numbers = std::unordered_map{ + {type_A, i}, {type_B, j}, {type_C, k}}; + auto const energy = static_cast(i + 1); + auto const f_expr = + calculate_factorial_expression(reaction, p_numbers); + // acceptance = V^{nu_bar} * gamma * f_expr * exp(- E / T) + auto const acceptance_ref = std::pow(r_algo.volume, reaction.nu_bar) * + reaction.gamma * f_expr * + std::exp(-energy / r_algo.kT); + auto const acceptance = r_algo_si->calculate_acceptance_probability( + reaction, energy, p_numbers); + BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); + } + } + } + } + + // check that the system energy is updated after a successful reaction + { + auto &&r_algo_si = make_algo(42, 1., 0., {}); + auto &r_algo = *r_algo_si->RE(); + r_algo.set_volume(1.); + r_algo.non_interacting_type = 5; + + // create a generic identity exchange reaction D <-> E + int const type_D = 0; + int const type_E = 1; + + r_algo.charges_of_types[type_D] = 0; + r_algo.charges_of_types[type_E] = 0; + + // track particles + init_type_map(type_D); + init_type_map(type_E); + + // set up a reaction completely shifted to product creation + auto const gamma = 1e100; + auto &&reaction_si = make_reaction(gamma, {type_D}, {1}, {type_E}, {1}); + auto &reaction = *reaction_si->get_reaction(); + r_algo_si->do_call_method("add_reaction", {{"reaction", reaction_si}}); + auto const reaction_id = 0; + + // resize system box + auto const box_l = Utils::Vector3d{0.5, 0.4, 0.7}; + espresso::system->set_box_l(box_l); + + // without reactants, no reaction will take place + auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + BOOST_REQUIRE(not result.has_value()); + + // the reaction was updated + BOOST_CHECK_EQUAL(reaction.tried_moves, 1); + BOOST_CHECK_EQUAL(reaction.accepted_moves, 0); + + // add reactants in the system + auto const ref_position = Utils::Vector3d{0.1, 0.2, 0.3}; + create_particle(ref_position, 1, type_D); + create_particle(ref_position, 2, type_D); + + { + // for an ideal system with gamma ~ inf, the reaction is always accepted; + // the potential energy of the new state has to be 0 since it is an ideal + // system + auto const energy_ref = 0.; + + auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + BOOST_REQUIRE(result.has_value()); + auto const energy_move = *result; + + // verify bookkeeping + auto const &bookkeeping = r_algo.get_old_system_state(); + BOOST_REQUIRE_EQUAL(bookkeeping.changed.size(), 1ul); + BOOST_REQUIRE_EQUAL(bookkeeping.created.size(), 0ul); + BOOST_REQUIRE_EQUAL(bookkeeping.hidden.size(), 0ul); + BOOST_REQUIRE_EQUAL(bookkeeping.moved.size(), 0ul); + BOOST_REQUIRE_EQUAL(bookkeeping.reaction_id, reaction_id); + BOOST_REQUIRE_EQUAL(bookkeeping.old_particle_numbers.at(type_D), 2); + BOOST_REQUIRE_EQUAL(bookkeeping.old_particle_numbers.at(type_E), 0); + BOOST_REQUIRE_EQUAL(std::get<1>(bookkeeping.changed[0]), type_D); + + auto const bf = r_algo_si->calculate_acceptance_probability( + reaction, energy_move, {{type_D, 1}, {type_E, 0}}); + + auto const energy_end = r_algo.generic_oneway_reaction_part_2( + reaction_id, bf, 0., energy_move); + BOOST_CHECK_CLOSE(energy_end, energy_ref, tol); + + // verify bookkeeping was cleared + BOOST_CHECK_THROW(r_algo.get_old_system_state(), std::runtime_error); + + // the reaction was updated + BOOST_CHECK_EQUAL(reaction.tried_moves, 2); + BOOST_CHECK_EQUAL(reaction.accepted_moves, 1); + } + { + // attempt a second reaction + auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + BOOST_REQUIRE(result.has_value()); + + // verify bookkeeping + auto const &bookkeeping = r_algo.get_old_system_state(); + BOOST_REQUIRE_EQUAL(bookkeeping.changed.size(), 1ul); + BOOST_REQUIRE_EQUAL(bookkeeping.created.size(), 0ul); + BOOST_REQUIRE_EQUAL(bookkeeping.hidden.size(), 0ul); + BOOST_REQUIRE_EQUAL(bookkeeping.moved.size(), 0ul); + BOOST_REQUIRE_EQUAL(bookkeeping.reaction_id, reaction_id); + BOOST_REQUIRE_EQUAL(bookkeeping.old_particle_numbers.at(type_D), 1); + BOOST_REQUIRE_EQUAL(bookkeeping.old_particle_numbers.at(type_E), 1); + BOOST_REQUIRE_EQUAL(std::get<1>(bookkeeping.changed[0]), type_D); + + // force move to be rejected + auto const energy_reject = + r_algo.generic_oneway_reaction_part_2(reaction_id, 0., 0.2, 0.1); + BOOST_CHECK_CLOSE(energy_reject, 0.2, tol); + + // the reaction was updated + BOOST_CHECK_EQUAL(reaction.tried_moves, 3); + BOOST_CHECK_EQUAL(reaction.accepted_moves, 1); + + // verify bookkeeping was cleared + BOOST_CHECK_THROW(r_algo.get_old_system_state(), std::runtime_error); + } + } +} + +int main(int argc, char **argv) { + espresso::system = std::make_unique(argc, argv); + return boost::unit_test::unit_test_main(init_unit_test, argc, argv); +} diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py index d4d5fba590b..4b9cb5ab644 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH.py @@ -56,16 +56,13 @@ def test_ideal_alpha(self): exclusion_range=1, seed=44, constant_pH=pH) + RE.set_non_interacting_type(type=max(types.values()) + 1) RE.add_reaction( gamma=10**(-pKa), reactant_types=[types["HA"]], product_types=[types["A-"], types["H+"]], default_charges=charges_dict) - # Set the hidden particle type to the lowest possible number to speed - # up the simulation - RE.set_non_interacting_type(type=max(types.values()) + 1) - # equilibration RE.reaction(reaction_steps=800) diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py index 18e1f33b900..ad0a6ae4afd 100644 --- a/testsuite/python/constant_pH_stats.py +++ b/testsuite/python/constant_pH_stats.py @@ -54,6 +54,7 @@ class Test(ut.TestCase): RE = espressomd.reaction_methods.ConstantpHEnsemble( kT=1., exclusion_range=1., seed=44, constant_pH=pH, search_algorithm="parallel") + RE.set_non_interacting_type(type=max(types.values()) + 1) @classmethod def setUpClass(cls): @@ -77,10 +78,6 @@ def test_ideal_titration_curve(self): types = self.types system = self.system - # Set the hidden particle type to the lowest possible number to speed - # up the simulation - RE.set_non_interacting_type(type=max(types.values()) + 1) - # chemical warmup - get close to chemical equilibrium before we start # sampling RE.reaction(reaction_steps=40 * N0) diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py index 40b5b1b5846..7f4936b3a98 100644 --- a/testsuite/python/reaction_ensemble.py +++ b/testsuite/python/reaction_ensemble.py @@ -70,6 +70,7 @@ class ReactionEnsembleTest(ut.TestCase): RE = espressomd.reaction_methods.ReactionEnsemble( seed=12, kT=temperature, exclusion_range=exclusion_range, search_algorithm="parallel") + RE.set_non_interacting_type(type=max(types.values()) + 1) @classmethod def setUpClass(cls): @@ -94,10 +95,6 @@ def test_ideal_titration_curve(self): RE = ReactionEnsembleTest.RE target_alpha = ReactionEnsembleTest.target_alpha - # Set the hidden particle type to the lowest possible number to speed - # up the simulation - RE.set_non_interacting_type(type=max(types.values()) + 1) - # chemical warmup - get close to chemical equilibrium before we start # sampling RE.reaction(reaction_steps=5 * N0) @@ -137,6 +134,14 @@ def test_ideal_titration_curve(self): rate1 = RE.get_acceptance_rate_reaction(reaction_id=1) self.assertAlmostEqual(rate0, 0.85, delta=0.05) self.assertAlmostEqual(rate1, 0.85, delta=0.05) + self.assertAlmostEqual( + rate0, + RE.reactions[0].get_acceptance_rate(), + delta=1e-10) + self.assertAlmostEqual( + rate1, + RE.reactions[1].get_acceptance_rate(), + delta=1e-10) if __name__ == "__main__": diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py index 6e1f07be643..2f677d08f2e 100644 --- a/testsuite/python/reaction_methods_interface.py +++ b/testsuite/python/reaction_methods_interface.py @@ -108,6 +108,10 @@ def count_by_type(types): type_mc=0, particle_number_to_be_changed=0)) self.assertFalse(method.displacement_mc_move_for_particles_of_type( type_mc=0, particle_number_to_be_changed=100000)) + method.particle_inside_exclusion_range_touched = True + self.assertTrue(method.particle_inside_exclusion_range_touched) + method.particle_inside_exclusion_range_touched = False + self.assertFalse(method.particle_inside_exclusion_range_touched) # check constraints method.set_wall_constraints_in_z_direction( @@ -135,13 +139,21 @@ def count_by_type(types): reactions = method.reactions self.assertEqual(len(reactions), 2) check_reaction_parameters(method.reactions, reaction_parameters) + self.assertIsNone(reactions[0].call_method("unknown")) - # check reactions after parameter change + # check reactions after unsuccessful parameter change + with self.assertRaises(ValueError): + method.change_reaction_constant(reaction_id=0, gamma=0.) + check_reaction_parameters(method.reactions, reaction_parameters) + check_reaction_parameters(method._reactions_cache, reaction_parameters) + + # check reactions after successful parameter change new_gamma = 634. reaction_forward['gamma'] = new_gamma reaction_backward['gamma'] = 1. / new_gamma method.change_reaction_constant(reaction_id=0, gamma=new_gamma) check_reaction_parameters(method.reactions, reaction_parameters) + check_reaction_parameters(method._reactions_cache, reaction_parameters) status = method.get_status() self.assertAlmostEqual( status['reactions'][0]['gamma'], @@ -179,7 +191,6 @@ def test_interface(self): params = {'exclusion_range': 0.8, 'exclusion_radius_per_type': {1: 0.1}} - # reaction ensemble with self.subTest(msg="reaction ensemble"): method = espressomd.reaction_methods.ReactionEnsemble( kT=1.4, seed=12, search_algorithm="order_n", **params) @@ -224,6 +235,10 @@ def test_exceptions(self): method.add_reaction(**{**reaction_params, 'reactant_types': []}) with self.assertRaisesRegex(ValueError, f'products: {err_msg}'): method.add_reaction(**{**reaction_params, 'product_types': []}) + with self.assertRaisesRegex(ValueError, 'gamma'): + method.add_reaction(**{**reaction_params, 'gamma': 0.}) + with self.assertRaisesRegex(ValueError, 'gamma'): + method.add_reaction(**{**reaction_params, 'gamma': -2.}) # check charge conservation err_msg = 'Reaction system is not charge neutral' @@ -233,15 +248,20 @@ def test_exceptions(self): with self.assertRaisesRegex(ValueError, err_msg): method.add_reaction(default_charges={2: 1, 3: 0, 4: 1 + 1e-10}, **single_reaction_params) + with self.assertRaisesRegex(TypeError, "needs to be a dict"): + method.add_reaction(default_charges=(1, 2), + **single_reaction_params) # check invalid reaction id exceptions # (note: reactions id = 2 * reactions index) self.assertEqual(len(method.reactions), 2) for i in [-2, -1, 1, 2, 3]: - with self.assertRaisesRegex(IndexError, 'This reaction is not present'): + with self.assertRaisesRegex(IndexError, f"No reaction with id {i}"): method.delete_reaction(reaction_id=i) - with self.assertRaisesRegex(IndexError, 'This reaction is not present'): + with self.assertRaisesRegex(IndexError, f"No reaction with id {2*i}"): method.get_acceptance_rate_reaction(reaction_id=2 * i) + with self.assertRaisesRegex(ValueError, "Only forward reactions can be selected"): + method.change_reaction_constant(reaction_id=1, gamma=1.) # check constraint exceptions set_cyl_constraint = method.set_cylindrical_constraint_in_z_direction @@ -280,6 +300,8 @@ def test_exceptions(self): # check exceptions for missing particles with self.assertRaisesRegex(RuntimeError, "Particle id is greater than the max seen particle id"): method.delete_particle(p_id=0) + with self.assertRaisesRegex(ValueError, "Invalid particle id: -2"): + method.delete_particle(p_id=-2) with self.assertRaisesRegex(RuntimeError, "Trying to remove some non-existing particles from the system via the inverse Widom scheme"): widom.calculate_particle_insertion_potential_energy(reaction_id=0) with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"): @@ -312,6 +334,12 @@ def test_exceptions(self): with self.assertRaisesRegex(ValueError, "Parameter 'type_mc' must be >= 0"): method.displacement_mc_move_for_particles_of_type( type_mc=-1, particle_number_to_be_changed=1) + with self.assertRaisesRegex(RuntimeError, "No chemical reaction is currently under way"): + method.call_method("calculate_factorial_expression") + with self.assertRaisesRegex(RuntimeError, "cannot be instantiated"): + espressomd.reaction_methods.ReactionAlgorithm() + with self.assertRaisesRegex(ValueError, "Invalid type: -1"): + method.set_non_interacting_type(type=-1) # check invalid exclusion ranges and radii with self.assertRaisesRegex(ValueError, "Invalid value for 'exclusion_range'"): From 82cda68d475f28c81e555f942f3719cad0facb2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 9 Feb 2023 22:53:34 +0100 Subject: [PATCH 2/5] Remove particle_data.cpp Rewrite the EspressoSystemStandAlone class and all particle property setters, getters and events without MpiCallbacks. Run unit tests in parallel. --- src/core/CMakeLists.txt | 1 - src/core/EspressoSystemStandAlone.cpp | 36 +- src/core/bond_breakage/bond_breakage.cpp | 1 - src/core/collision.cpp | 10 +- src/core/constraints/ShapeBasedConstraint.hpp | 2 +- src/core/electrostatics/scafacos_impl.cpp | 1 - src/core/event.cpp | 2 +- src/core/grid.cpp | 1 - src/core/magnetostatics/dlc.cpp | 1 - src/core/magnetostatics/scafacos_impl.cpp | 1 - .../nonbonded_interaction_data.cpp | 19 +- .../nonbonded_interaction_data.hpp | 5 - src/core/particle_data.cpp | 249 -------------- src/core/particle_data.hpp | 53 --- src/core/particle_node.cpp | 111 ++---- src/core/particle_node.hpp | 13 +- .../reaction_methods/ReactionAlgorithm.cpp | 20 +- .../reaction_methods/tests/CMakeLists.txt | 2 +- .../tests/ReactionAlgorithm_test.cpp | 12 +- .../tests/particle_tracking_test.cpp | 2 +- src/core/thermostat.hpp | 4 + src/core/unit_tests/CMakeLists.txt | 10 +- .../EspressoSystemInterface_test.cpp | 58 ++-- .../EspressoSystemStandAlone_test.cpp | 322 +++++++++--------- src/core/unit_tests/ParticleFactory.hpp | 39 ++- src/core/unit_tests/VerletCriterion_test.cpp | 1 - src/core/unit_tests/Verlet_list_test.cpp | 165 +++++---- .../interactions/NonBondedInteraction.hpp | 3 +- .../interactions/NonBondedInteractions.hpp | 4 +- .../particle_data/ParticleHandle.cpp | 12 +- .../particle_data/ParticleList.cpp | 14 +- src/script_interface/tests/CMakeLists.txt | 6 +- 32 files changed, 401 insertions(+), 779 deletions(-) delete mode 100644 src/core/particle_data.cpp delete mode 100644 src/core/particle_data.hpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 2d12bcb2c01..04465777b7d 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -41,7 +41,6 @@ add_library( integrate.cpp npt.cpp partCfg_global.cpp - particle_data.cpp particle_node.cpp polymer.cpp pressure.cpp diff --git a/src/core/EspressoSystemStandAlone.cpp b/src/core/EspressoSystemStandAlone.cpp index cfa2aa0c5e6..3dee8e90b87 100644 --- a/src/core/EspressoSystemStandAlone.cpp +++ b/src/core/EspressoSystemStandAlone.cpp @@ -20,7 +20,6 @@ #include "config/config.hpp" #include "EspressoSystemStandAlone.hpp" -#include "MpiCallbacks.hpp" #include "communication.hpp" #include "event.hpp" #include "grid.hpp" @@ -47,48 +46,21 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) { #ifdef VIRTUAL_SITES set_virtual_sites(std::make_shared()); #endif - - // initialize the MpiCallbacks loop (blocking on worker nodes) - mpi_loop(); -} - -static void mpi_set_box_length_local(Utils::Vector3d const &value) { - set_box_length(value); -} - -static void mpi_set_node_grid_local(Utils::Vector3i const &value) { - set_node_grid(value); -} - -static void mpi_set_time_step_local(double const &value) { - set_time_step(value); } -REGISTER_CALLBACK(mpi_set_box_length_local) -REGISTER_CALLBACK(mpi_set_node_grid_local) -REGISTER_CALLBACK(mpi_set_time_step_local) - void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const { - if (!head_node) - return; - mpi_call_all(mpi_set_box_length_local, box_l); + set_box_length(box_l); } void EspressoSystemStandAlone::set_node_grid( Utils::Vector3i const &node_grid) const { - if (!head_node) - return; - mpi_call_all(mpi_set_node_grid_local, node_grid); + ::set_node_grid(node_grid); } void EspressoSystemStandAlone::set_time_step(double time_step) const { - if (!head_node) - return; - mpi_call_all(mpi_set_time_step_local, time_step); + ::set_time_step(time_step); } void EspressoSystemStandAlone::set_skin(double new_skin) const { - if (!head_node) - return; - mpi_call_all(mpi_set_skin_local, new_skin); + mpi_set_skin_local(new_skin); } diff --git a/src/core/bond_breakage/bond_breakage.cpp b/src/core/bond_breakage/bond_breakage.cpp index caee1750393..b0248047ed1 100644 --- a/src/core/bond_breakage/bond_breakage.cpp +++ b/src/core/bond_breakage/bond_breakage.cpp @@ -23,7 +23,6 @@ #include "communication.hpp" #include "errorhandling.hpp" #include "event.hpp" -#include "particle_data.hpp" #include diff --git a/src/core/collision.cpp b/src/core/collision.cpp index 0c199d7d940..41f47f6d738 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -193,7 +193,7 @@ void Collision_parameters::initialize() { throw std::domain_error("Collision detection particle type for virtual " "sites needs to be >=0"); } - make_particle_type_exist_local(collision_params.vs_particle_type); + make_particle_type_exist(collision_params.vs_particle_type); } if (collision_params.mode == CollisionModeType::GLUE_TO_SURF) { @@ -201,25 +201,25 @@ void Collision_parameters::initialize() { throw std::domain_error("Collision detection particle type for virtual " "sites needs to be >=0"); } - make_particle_type_exist_local(collision_params.vs_particle_type); + make_particle_type_exist(collision_params.vs_particle_type); if (collision_params.part_type_to_be_glued < 0) { throw std::domain_error("Collision detection particle type to be glued " "needs to be >=0"); } - make_particle_type_exist_local(collision_params.part_type_to_be_glued); + make_particle_type_exist(collision_params.part_type_to_be_glued); if (collision_params.part_type_to_attach_vs_to < 0) { throw std::domain_error("Collision detection particle type to attach " "the virtual site to needs to be >=0"); } - make_particle_type_exist_local(collision_params.part_type_to_attach_vs_to); + make_particle_type_exist(collision_params.part_type_to_attach_vs_to); if (collision_params.part_type_after_glueing < 0) { throw std::domain_error("Collision detection particle type after gluing " "needs to be >=0"); } - make_particle_type_exist_local(collision_params.part_type_after_glueing); + make_particle_type_exist(collision_params.part_type_after_glueing); } on_short_range_ia_change(); diff --git a/src/core/constraints/ShapeBasedConstraint.hpp b/src/core/constraints/ShapeBasedConstraint.hpp index 351536bc19d..a2c9a9cb920 100644 --- a/src/core/constraints/ShapeBasedConstraint.hpp +++ b/src/core/constraints/ShapeBasedConstraint.hpp @@ -77,7 +77,7 @@ class ShapeBasedConstraint : public Constraint { void set_type(const int &type) { part_rep.type() = type; - make_particle_type_exist_local(type); + make_particle_type_exist(type); } Utils::Vector3d total_force() const; diff --git a/src/core/electrostatics/scafacos_impl.cpp b/src/core/electrostatics/scafacos_impl.cpp index 3a7ce9b1021..365eaeeb2a0 100644 --- a/src/core/electrostatics/scafacos_impl.cpp +++ b/src/core/electrostatics/scafacos_impl.cpp @@ -31,7 +31,6 @@ #include "event.hpp" #include "grid.hpp" #include "integrate.hpp" -#include "particle_data.hpp" #include "tuning.hpp" #include diff --git a/src/core/event.cpp b/src/core/event.cpp index 8cce1321c82..bf1a3588d16 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -84,7 +84,7 @@ void on_program_start() { cells_re_init(CellStructureType::CELL_STRUCTURE_REGULAR); /* make sure interaction 0<->0 always exists */ - make_particle_type_exist_local(0); + make_particle_type_exist(0); } void on_integration_start(double time_step) { diff --git a/src/core/grid.cpp b/src/core/grid.cpp index d06ec670a40..db541fb5a34 100644 --- a/src/core/grid.cpp +++ b/src/core/grid.cpp @@ -28,7 +28,6 @@ #include "communication.hpp" #include "event.hpp" -#include "particle_data.hpp" #include #include diff --git a/src/core/magnetostatics/dlc.cpp b/src/core/magnetostatics/dlc.cpp index 96d8231585d..fb6811b936a 100644 --- a/src/core/magnetostatics/dlc.cpp +++ b/src/core/magnetostatics/dlc.cpp @@ -35,7 +35,6 @@ #include "communication.hpp" #include "errorhandling.hpp" #include "grid.hpp" -#include "particle_data.hpp" #include #include diff --git a/src/core/magnetostatics/scafacos_impl.cpp b/src/core/magnetostatics/scafacos_impl.cpp index c71f69af5f1..21bdf3ba540 100644 --- a/src/core/magnetostatics/scafacos_impl.cpp +++ b/src/core/magnetostatics/scafacos_impl.cpp @@ -29,7 +29,6 @@ #include "cells.hpp" #include "communication.hpp" #include "grid.hpp" -#include "particle_data.hpp" #include #include diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index 8b1925ba306..0644de67fb8 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -23,7 +23,6 @@ */ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "communication.hpp" #include "electrostatics/coulomb.hpp" #include "event.hpp" @@ -50,7 +49,7 @@ static double min_global_cut = INACTIVE_CUTOFF; * general low-level functions *****************************************/ -void mpi_realloc_ia_params_local(int new_size) { +static void realloc_ia_params(int new_size) { auto const old_size = ::max_seen_particle_type; if (new_size <= old_size) return; @@ -76,8 +75,6 @@ void mpi_realloc_ia_params_local(int new_size) { ::nonbonded_ia_params = std::move(new_params); } -REGISTER_CALLBACK(mpi_realloc_ia_params_local) - static double recalc_maximal_cutoff(const IA_parameters &data) { auto max_cut_current = INACTIVE_CUTOFF; @@ -165,19 +162,7 @@ double maximal_cutoff_nonbonded() { return max_cut_nonbonded; } -bool is_new_particle_type(int type) { - return (type + 1) > max_seen_particle_type; -} - -void make_particle_type_exist(int type) { - if (is_new_particle_type(type)) - mpi_call_all(mpi_realloc_ia_params_local, type + 1); -} - -void make_particle_type_exist_local(int type) { - if (is_new_particle_type(type)) - mpi_realloc_ia_params_local(type + 1); -} +void make_particle_type_exist(int type) { realloc_ia_params(type + 1); } void set_min_global_cut(double min_global_cut) { ::min_global_cut = min_global_cut; diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index a6370f456d5..2fd5ad27f67 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp @@ -387,17 +387,12 @@ inline IA_parameters &get_ia_param(int i, int j) { return *::nonbonded_ia_params[get_ia_param_key(i, j)]; } -void mpi_realloc_ia_params_local(int new_size); - -bool is_new_particle_type(int type); /** Make sure that ia_params is large enough to cover interactions * for this particle type. The interactions are initialized with values * such that no physical interaction occurs. */ void make_particle_type_exist(int type); -void make_particle_type_exist_local(int type); - /** Check if a non-bonded interaction is defined */ inline bool checkIfInteraction(IA_parameters const &data) { return data.max_cut != INACTIVE_CUTOFF; diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp deleted file mode 100644 index f61db4ba881..00000000000 --- a/src/core/particle_data.cpp +++ /dev/null @@ -1,249 +0,0 @@ -/* - * Copyright (C) 2010-2022 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * 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 . - */ - -#include "particle_data.hpp" - -#include "Particle.hpp" -#include "cells.hpp" -#include "communication.hpp" -#include "event.hpp" -#include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "particle_node.hpp" - -#include "config/config.hpp" - -#include - -#include -#include - -#include - -constexpr auto some_tag = 42; - -namespace { -/** - * @brief A generic particle update. - * - * Here the sub-struct structure of Particle is - * used: the specification of the data member to update - * consists of two parts, the pointer to the substruct @p s - * and a pointer to a member of that substruct @p m. - * - * @tparam S Substruct type of Particle - * @tparam s Pointer to a member of Particle - * @tparam T Type of the member to update, must be serializable - * @tparam m Pointer to the member. - */ -template -struct UpdateParticle { - T value; - - void operator()(Particle &p) const { (p.*s).*m = value; } - - template void serialize(Archive &ar, long int) { ar &value; } -}; - -template -using UpdateProperty = UpdateParticle; - -template -using UpdateLocalProperty = UpdateParticle; -template -using UpdatePosition = UpdateParticle; -template -using UpdateMomentum = UpdateParticle; -template -using UpdateForce = UpdateParticle; - -using Prop = ParticleProperties; - -// clang-format off -using UpdatePropertyMessage = boost::variant - < UpdateProperty -#ifdef ELECTROSTATICS - , UpdateProperty -#endif -#ifdef EXTERNAL_FORCES - , UpdateProperty -#endif - >; - -using UpdatePositionMessage = boost::variant - < UpdatePosition >; - -using UpdateMomentumMessage = boost::variant - < UpdateMomentum >; - -using UpdateForceMessage = boost::variant - < UpdateForce >; - -/** - * @brief Top-level message. - * - * A message is either updates a property, - * or a position, or ... - * New messages can be added here, if they - * fulfill the type requirements, namely: - * They either have an integer member id indicating the particle - * that should be updated with an operator()(const Particle&) - * that is called with the particle, or a tree of - * variants with leaves that have such an operator() member. - */ -using UpdateMessage = boost::variant - < UpdatePropertyMessage - , UpdatePositionMessage - , UpdateMomentumMessage - , UpdateForceMessage - >; -// clang-format on - -/** - * @brief Meta-function to detect the message type from - * the particle substruct. - */ -template struct message_type; - -template <> struct message_type { - using type = UpdatePropertyMessage; -}; - -template <> struct message_type { - using type = UpdatePositionMessage; -}; - -template <> struct message_type { - using type = UpdateMomentumMessage; -}; - -template <> struct message_type { - using type = UpdateForceMessage; -}; - -template -using message_type_t = typename message_type::type; - -/** - * @brief Visitor for message evaluation. - * - * This visitor either recurses into the active type - * if it is a variant type, or otherwise calls - * operator()(Particle&) on the active type with - * the particle that ought to be updated. This construction - * allows to nest message variants to allow for sub- - * categories. Those are mostly used here to differentiate - * the updates for the substructs of Particle. - */ -struct UpdateVisitor : public boost::static_visitor { - explicit UpdateVisitor(int id) : id(id) {} - - const int id; - - /* Recurse into sub-variants */ - template - void operator()(const boost::variant &msg) const { - boost::apply_visitor(*this, msg); - } - /* Plain messages are just called. */ - template void operator()(const Message &msg) const { - assert(cell_structure.get_local_particle(id)); - msg(*cell_structure.get_local_particle(id)); - } -}; -} // namespace - -static void mpi_send_update_message_local(int node, int id) { - if (node == comm_cart.rank()) { - UpdateMessage msg{}; - comm_cart.recv(0, some_tag, msg); - boost::apply_visitor(UpdateVisitor(id), msg); - } - - on_particle_change(); -} - -REGISTER_CALLBACK(mpi_send_update_message_local) - -/** - * @brief Send a particle update message. - * - * This sends the message to the node that is responsible for the particle, - * where - * @p msg is called with the particle as argument. The message then performs the - * change to the particle that is encoded in it. The mechanism to call a functor - * based on the active type of a variant is called visitation. Here we can use - * @c UpdateVisitor with a single templated call operator on the visitor because - * all message (eventually) provide the same interface. Overall this is - * logically equivalent to nested switch statements over the message types, - * where the case statements play the role of the call of the messages in our - * case. A general introduction can be found in the documentation of - * boost::variant. - * - * @param id Id of the particle to update - * @param msg The message - */ -static void mpi_send_update_message(int id, const UpdateMessage &msg) { - auto const pnode = get_particle_node(id); - - mpi_call(mpi_send_update_message_local, pnode, id); - - /* If the particle is remote, send the - * message to the target, otherwise we - * can just apply the update directly. */ - if (pnode != comm_cart.rank()) { - comm_cart.send(pnode, some_tag, msg); - } else { - boost::apply_visitor(UpdateVisitor(id), msg); - } - - on_particle_change(); -} - -template -void mpi_update_particle(int id, const T &value) { - using MessageType = message_type_t; - MessageType msg = UpdateParticle{value}; - mpi_send_update_message(id, msg); -} - -template -void mpi_update_particle_property(int id, const T &value) { - mpi_update_particle(id, value); -} - -void set_particle_v(int part, Utils::Vector3d const &v) { - mpi_update_particle(part, v); -} - -#ifdef ELECTROSTATICS -void set_particle_q(int part, double q) { - mpi_update_particle_property(part, q); -} -#else -const constexpr double ParticleProperties::q; -#endif - -void set_particle_type(int p_id, int type) { - make_particle_type_exist(type); - on_particle_type_change(p_id, type); - mpi_update_particle_property(p_id, type); -} diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp deleted file mode 100644 index 0799203b009..00000000000 --- a/src/core/particle_data.hpp +++ /dev/null @@ -1,53 +0,0 @@ -/* - * Copyright (C) 2010-2022 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * 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 . - */ -#ifndef CORE_PARTICLE_DATA_HPP -#define CORE_PARTICLE_DATA_HPP -/** \file - * Particles property update. - * - * This file contains everything related to updating particle properties. - */ - -#include "config/config.hpp" - -#include - -/** Call only on the head node: set particle velocity. - * @param part the particle. - * @param v its new velocity. - */ -void set_particle_v(int part, Utils::Vector3d const &v); - -#ifdef ELECTROSTATICS -/** Call only on the head node: set particle charge. - * @param part the particle. - * @param q its new charge. - */ -void set_particle_q(int part, double q); -#endif - -/** Call only on the head node: set particle type. - * @param p_id the particle. - * @param type its new type. - */ -void set_particle_type(int p_id, int type); - -#endif diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp index 42d43d385ca..a5d186fb148 100644 --- a/src/core/particle_node.cpp +++ b/src/core/particle_node.cpp @@ -28,7 +28,6 @@ #include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "partCfg_global.hpp" -#include "particle_data.hpp" #include #include @@ -92,7 +91,7 @@ void init_type_map(int type) { throw std::runtime_error("Types may not be negative"); } ::type_list_enable = true; - make_particle_type_exist_local(type); + make_particle_type_exist(type); std::vector local_pids; for (auto const &p : ::cell_structure.local_particles()) { @@ -123,7 +122,7 @@ static void add_id_to_type_map(int p_id, int type) { it->second.insert(p_id); } -void on_particle_type_change_parallel(int p_id, int old_type, int new_type) { +void on_particle_type_change(int p_id, int old_type, int new_type) { if (::type_list_enable) { if (old_type == type_tracking::any_type) { for (auto &kv : ::particle_type_map) { @@ -147,21 +146,6 @@ void on_particle_type_change_parallel(int p_id, int old_type, int new_type) { } } -void on_particle_type_change(int p_id, int type) { - if (type_list_enable) { - assert(::this_node == 0); - // check if the particle exists already and the type is changed, then remove - // it from the list which contains it - auto const &cur_par = get_particle_data(p_id); - int prev_type = cur_par.type(); - if (prev_type != type) { - // particle existed before so delete it from the list - remove_id_from_map(p_id, prev_type); - } - add_id_to_type_map(p_id, type); - } -} - namespace { /* Limit cache to 100 MiB */ std::size_t const max_cache_size = (100ul * 1048576ul) / sizeof(Particle); @@ -386,6 +370,29 @@ int get_particle_node(int p_id) { return needle->second; } +int get_particle_node_parallel(int p_id) { + if (p_id < 0) { + throw std::domain_error("Invalid particle id: " + std::to_string(p_id)); + } + + if (rebuild_needed()) { + build_particle_node_parallel(); + } + + if (this_node != 0) { + return -1; + } + + auto const needle = particle_node.find(p_id); + + // Check if particle has a node, if not, we assume it does not exist. + if (needle == particle_node.end()) { + throw std::runtime_error("Particle node for id " + std::to_string(p_id) + + " not found!"); + } + return needle->second; +} + void clear_particle_node() { ::max_seen_pid = -1; particle_node.clear(); @@ -462,49 +469,14 @@ void particle_checks(int p_id, Utils::Vector3d const &pos) { #endif // __FAST_MATH__ } -static void mpi_remove_particle_local(int p_id) { - cell_structure.remove_particle(p_id); - on_particle_change(); - mpi_synchronize_max_seen_pid_local(); -} - -REGISTER_CALLBACK(mpi_remove_particle_local) - -static void mpi_remove_all_particles_local() { - cell_structure.remove_all_particles(); +void remove_all_particles() { + ::cell_structure.remove_all_particles(); on_particle_change(); clear_particle_node(); clear_particle_type_map(); } -REGISTER_CALLBACK(mpi_remove_all_particles_local) - -void remove_all_particles() { mpi_call_all(mpi_remove_all_particles_local); } - void remove_particle(int p_id) { - auto const &cur_par = get_particle_data(p_id); - if (type_list_enable) { - // remove particle from its current type_list - int type = cur_par.type(); - remove_id_from_map(p_id, type); - } - - particle_node[p_id] = -1; - mpi_call_all(mpi_remove_particle_local, p_id); - particle_node.erase(p_id); - if (p_id == max_seen_pid) { - --max_seen_pid; - // if there is a gap (i.e. there is no particle with id max_seen_pid - 1, - // then the cached value is invalidated and has to be recomputed (slow) - if (particle_node.count(max_seen_pid) == 0 or - particle_node[max_seen_pid] == -1) { - max_seen_pid = calculate_max_seen_id(); - } - } - mpi_call_all(mpi_synchronize_max_seen_pid_local); -} - -void remove_particle_parallel(int p_id) { if (::type_list_enable) { auto p = ::cell_structure.get_local_particle(p_id); auto p_type = -1; @@ -525,7 +497,9 @@ void remove_particle_parallel(int p_id) { if (this_node == 0) { particle_node[p_id] = -1; } - mpi_remove_particle_local(p_id); + ::cell_structure.remove_particle(p_id); + on_particle_change(); + mpi_synchronize_max_seen_pid_local(); if (this_node == 0) { particle_node.erase(p_id); if (p_id == ::max_seen_pid) { @@ -541,7 +515,7 @@ void remove_particle_parallel(int p_id) { mpi_synchronize_max_seen_pid_local(); } -void mpi_make_new_particle_local(int p_id, Utils::Vector3d const &pos) { +void make_new_particle(int p_id, Utils::Vector3d const &pos) { if (rebuild_needed()) { build_particle_node_parallel(); } @@ -559,14 +533,7 @@ void mpi_make_new_particle_local(int p_id, Utils::Vector3d const &pos) { mpi_synchronize_max_seen_pid_local(); } -REGISTER_CALLBACK(mpi_make_new_particle_local) - -void mpi_make_new_particle(int p_id, Utils::Vector3d const &pos) { - particle_checks(p_id, pos); - mpi_call_all(mpi_make_new_particle_local, p_id, pos); -} - -void mpi_set_particle_pos_local(int p_id, Utils::Vector3d const &pos) { +void set_particle_pos(int p_id, Utils::Vector3d const &pos) { auto const has_moved = maybe_move_particle(p_id, pos); ::cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); on_particle_change(); @@ -579,13 +546,6 @@ void mpi_set_particle_pos_local(int p_id, Utils::Vector3d const &pos) { } } -REGISTER_CALLBACK(mpi_set_particle_pos_local) - -void mpi_set_particle_pos(int p_id, Utils::Vector3d const &pos) { - particle_checks(p_id, pos); - mpi_call_all(mpi_set_particle_pos_local, p_id, pos); -} - int get_random_p_id(int type, int random_index_in_type_map) { auto it = particle_type_map.find(type); if (it == particle_type_map.end()) { @@ -640,13 +600,6 @@ std::vector get_particle_ids_parallel() { } int get_maximal_particle_id() { - if (particle_node.empty()) - build_particle_node(); - - return max_seen_pid; -} - -int get_maximal_particle_id_parallel() { if (rebuild_needed()) { build_particle_node_parallel(); } diff --git a/src/core/particle_node.hpp b/src/core/particle_node.hpp index fa60f258465..d6577cbf2e8 100644 --- a/src/core/particle_node.hpp +++ b/src/core/particle_node.hpp @@ -75,16 +75,13 @@ std::size_t fetch_cache_max_size(); */ void clear_particle_node(); -void mpi_set_particle_pos_local(int p_id, Utils::Vector3d const &pos); - /** * @brief Create a new particle and attach it to a cell. * Also call @ref on_particle_change. * @param p_id The identity of the particle to create. * @param pos The particle position. */ -void mpi_make_new_particle(int p_id, Utils::Vector3d const &pos); -void mpi_make_new_particle_local(int p_id, Utils::Vector3d const &pos); +void make_new_particle(int p_id, Utils::Vector3d const &pos); void particle_checks(int p_id, Utils::Vector3d const &pos); /** @@ -93,21 +90,19 @@ void particle_checks(int p_id, Utils::Vector3d const &pos); * @param p_id The identity of the particle to move. * @param pos The new particle position. */ -void mpi_set_particle_pos(int p_id, Utils::Vector3d const &pos); +void set_particle_pos(int p_id, Utils::Vector3d const &pos); /** Remove particle with a given identity. Also removes all bonds to the * particle. * @param p_id identity of the particle to remove */ void remove_particle(int p_id); -void remove_particle_parallel(int p_id); /** Remove all particles. */ void remove_all_particles(); void init_type_map(int type); -void on_particle_type_change(int p_id, int type); -void on_particle_type_change_parallel(int p_id, int old_type, int new_type); +void on_particle_type_change(int p_id, int old_type, int new_type); /** Find a particle of given type and return its id */ int get_random_p_id(int type, int random_index_in_type_map); @@ -128,6 +123,7 @@ bool particle_exists(int p_id); * @return The MPI rank the particle is on. */ int get_particle_node(int p_id); +int get_particle_node_parallel(int p_id); /** * @brief Get all particle ids. @@ -141,7 +137,6 @@ std::vector get_particle_ids_parallel(); * @brief Get maximal particle id. */ int get_maximal_particle_id(); -int get_maximal_particle_id_parallel(); /** * @brief Get number of particles. diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index f259432cc8e..78cedbab7ba 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -86,7 +86,7 @@ void ReactionAlgorithm::restore_old_system_state() { // restore the properties of changed and hidden particles for (auto const &state : {old_state.changed, old_state.hidden}) { for (auto const &[p_id, p_type] : state) { - on_particle_type_change_parallel(p_id, type_tracking::any_type, p_type); + on_particle_type_change(p_id, type_tracking::any_type, p_type); if (auto p = get_local_particle(p_id)) { p->type() = p_type; #ifdef ELECTROSTATICS @@ -158,7 +158,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, auto const new_type = reaction.product_types[i]; for (int j = 0; j < std::min(n_product_coef, n_reactant_coef); j++) { auto const p_id = get_random_p_id_of_type(old_type); - on_particle_type_change_parallel(p_id, old_type, new_type); + on_particle_type_change(p_id, old_type, new_type); if (auto p = get_local_particle(p_id)) { p->type() = new_type; #ifdef ELECTROSTATICS @@ -292,7 +292,7 @@ double ReactionAlgorithm::generic_oneway_reaction_part_2(int reaction_id, * like the one above). */ void ReactionAlgorithm::hide_particle(int p_id, int p_type) const { - on_particle_type_change_parallel(p_id, p_type, non_interacting_type); + on_particle_type_change(p_id, p_type, non_interacting_type); if (auto p = get_local_particle(p_id)) { p->type() = non_interacting_type; #ifdef ELECTROSTATICS @@ -378,10 +378,10 @@ void ReactionAlgorithm::delete_particle(int p_id) { if (p_id < 0) { throw std::domain_error("Invalid particle id: " + std::to_string(p_id)); } - auto const old_max_seen_id = ::get_maximal_particle_id_parallel(); + auto const old_max_seen_id = ::get_maximal_particle_id(); if (p_id == old_max_seen_id) { // last particle, just delete - remove_particle_parallel(p_id); + remove_particle(p_id); // remove all saved empty p_ids which are greater than the max_seen_particle // this is needed in order to avoid the creation of holes for (auto p_id_iter = m_empty_p_ids_smaller_than_max_seen_particle.begin(); @@ -393,7 +393,7 @@ void ReactionAlgorithm::delete_particle(int p_id) { ++p_id_iter; } } else if (p_id <= old_max_seen_id) { - remove_particle_parallel(p_id); + remove_particle(p_id); m_empty_p_ids_smaller_than_max_seen_particle.push_back(p_id); } else { throw std::runtime_error( @@ -470,14 +470,14 @@ int ReactionAlgorithm::create_particle(int p_type) { p_id = *p_id_iter; m_empty_p_ids_smaller_than_max_seen_particle.erase(p_id_iter); } else { - p_id = ::get_maximal_particle_id_parallel() + 1; + p_id = ::get_maximal_particle_id() + 1; } // create random velocity vector according to Maxwell-Boltzmann distribution auto pos = get_random_position_in_box(); auto vel = get_random_velocity_vector(); - ::mpi_make_new_particle_local(p_id, pos); + ::make_new_particle(p_id, pos); if (auto p = get_local_particle(p_id)) { p->v() = std::sqrt(kT / p->mass()) * vel; p->type() = p_type; @@ -485,7 +485,7 @@ int ReactionAlgorithm::create_particle(int p_type) { p->q() = charges_of_types.at(p_type); #endif } - on_particle_type_change_parallel(p_id, ::type_tracking::new_part, p_type); + on_particle_type_change(p_id, ::type_tracking::new_part, p_type); return p_id; } @@ -516,7 +516,7 @@ void ReactionAlgorithm::displacement_mc_move(int type, int n_particles) { } boost::mpi::broadcast(m_comm, old_state, 0); bookkeeping.moved.emplace_back(old_state); - mpi_set_particle_pos_local(p_id, new_pos); + ::set_particle_pos(p_id, new_pos); check_exclusion_range(p_id, type); if (particle_inside_exclusion_range_touched) { diff --git a/src/core/reaction_methods/tests/CMakeLists.txt b/src/core/reaction_methods/tests/CMakeLists.txt index 73b6c318afa..7dd67a9e634 100644 --- a/src/core/reaction_methods/tests/CMakeLists.txt +++ b/src/core/reaction_methods/tests/CMakeLists.txt @@ -20,7 +20,7 @@ include(unit_test) unit_test(NAME SingleReaction_test SRC SingleReaction_test.cpp DEPENDS espresso::core) unit_test(NAME ReactionAlgorithm_test SRC ReactionAlgorithm_test.cpp DEPENDS - espresso::core Boost::mpi MPI::MPI_CXX) + espresso::core Boost::mpi MPI::MPI_CXX NUM_PROC 2) unit_test(NAME particle_tracking_test SRC particle_tracking_test.cpp DEPENDS espresso::core Boost::mpi MPI::MPI_CXX) unit_test(NAME reaction_methods_utils_test SRC reaction_methods_utils_test.cpp diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp index a7cd1a55e62..c7ed0c5292d 100644 --- a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp +++ b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp @@ -135,8 +135,8 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) { std::vector> ref_positions{ {{0.1, 0.2, 0.3}, {+10., +20., +30.}}, {{0.4, 0.5, 0.6}, {-10., -20., -30.}}}; - mpi_make_new_particle(0, ref_positions[0].first); - mpi_make_new_particle(1, ref_positions[1].first); + ::make_new_particle(0, ref_positions[0].first); + ::make_new_particle(1, ref_positions[1].first); set_particle_v(0, ref_positions[0].second); set_particle_v(1, ref_positions[1].second); set_particle_type(0, type_A); @@ -176,8 +176,8 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) { auto const box_l = 1.; std::vector ref_positions{{0.1, 0.2, 0.3}, {0.6, 0.7, 0.8}}; - mpi_make_new_particle(0, ref_positions[0]); - mpi_make_new_particle(1, ref_positions[1]); + ::make_new_particle(0, ref_positions[0]); + ::make_new_particle(1, ref_positions[1]); set_particle_type(0, type_A); set_particle_type(1, type_A); // check runtime errors when a MC move cannot be physically performed @@ -294,8 +294,8 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) { espresso::system->set_box_l({1., 1., 1.}); // set up particle - mpi_make_new_particle(0, {0.5, 0.5, 0.5}); - mpi_make_new_particle(1, {0.7, 0.7, 0.7}); + ::make_new_particle(0, {0.5, 0.5, 0.5}); + ::make_new_particle(1, {0.7, 0.7, 0.7}); set_particle_type(0, type_A); set_particle_type(1, type_B); exclusion_radius_per_type[type_A] = 0.1; diff --git a/src/core/reaction_methods/tests/particle_tracking_test.cpp b/src/core/reaction_methods/tests/particle_tracking_test.cpp index c3a79e2396b..7d19a401e3b 100644 --- a/src/core/reaction_methods/tests/particle_tracking_test.cpp +++ b/src/core/reaction_methods/tests/particle_tracking_test.cpp @@ -26,7 +26,7 @@ #include "unit_tests/ParticleFactory.hpp" #include "communication.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index 14e294146bc..45e0726b313 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -331,6 +331,7 @@ struct StokesianThermostat : public BaseThermostat { * @param thermostat The thermostat object name */ #define NEW_THERMOSTAT(thermostat) \ + void mpi_##thermostat##_set_rng_seed(uint32_t seed); \ void thermostat##_set_rng_seed(uint32_t seed); \ void thermostat##_set_rng_counter(uint64_t seed); @@ -376,11 +377,14 @@ void mpi_set_langevin_gamma_rot(Thermostat::GammaType const &gamma); void mpi_set_thermo_virtual(bool thermo_virtual); void mpi_set_temperature(double temperature); +void mpi_set_temperature_local(double temperature); void mpi_set_thermo_switch(int thermo_switch); +void mpi_set_thermo_switch_local(int thermo_switch); #ifdef NPT void mpi_set_nptiso_gammas(double gamma0, double gammav); +void mpi_set_nptiso_gammas_local(double gamma0, double gammav); #endif #endif diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index 9be4e36a851..2ece0659ee7 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -26,11 +26,15 @@ unit_test(NAME RuntimeError_test SRC RuntimeError_test.cpp DEPENDS Boost::serialization) unit_test(NAME RuntimeErrorCollector_test SRC RuntimeErrorCollector_test.cpp DEPENDS espresso::core Boost::mpi MPI::MPI_CXX NUM_PROC 2) -unit_test(NAME EspressoSystemStandAlone_test SRC +unit_test(NAME EspressoSystemStandAlone_parallel_test SRC EspressoSystemStandAlone_test.cpp DEPENDS espresso::core Boost::mpi - MPI::MPI_CXX NUM_PROC 2) + MPI::MPI_CXX NUM_PROC 4) +unit_test(NAME EspressoSystemStandAlone_serial_test SRC + EspressoSystemStandAlone_test.cpp DEPENDS espresso::core Boost::mpi + MPI::MPI_CXX NUM_PROC 1) unit_test(NAME EspressoSystemInterface_test SRC - EspressoSystemInterface_test.cpp DEPENDS espresso::core Boost::mpi) + EspressoSystemInterface_test.cpp DEPENDS espresso::core Boost::mpi + NUM_PROC 2) unit_test(NAME MpiCallbacks_test SRC MpiCallbacks_test.cpp DEPENDS espresso::utils Boost::mpi MPI::MPI_CXX NUM_PROC 2) unit_test(NAME ParticleIterator_test SRC ParticleIterator_test.cpp DEPENDS diff --git a/src/core/unit_tests/EspressoSystemInterface_test.cpp b/src/core/unit_tests/EspressoSystemInterface_test.cpp index ca7e269e29a..f2d97586dac 100644 --- a/src/core/unit_tests/EspressoSystemInterface_test.cpp +++ b/src/core/unit_tests/EspressoSystemInterface_test.cpp @@ -23,12 +23,13 @@ #define BOOST_TEST_DYN_LINK #include +#include "ParticleFactory.hpp" + #include "EspressoSystemInterface.hpp" +#include "Particle.hpp" +#include "cells.hpp" #include "communication.hpp" -#include "particle_data.hpp" #include "particle_node.hpp" -#include "virtual_sites.hpp" -#include "virtual_sites/VirtualSitesOff.hpp" #include "config/config.hpp" @@ -66,39 +67,53 @@ boost::test_tools::assertion_result has_gpu(boost::unit_test::test_unit_id) { return has_compatible_gpu; } -BOOST_AUTO_TEST_CASE(check_with_gpu, *boost::unit_test::precondition(has_gpu)) { +BOOST_FIXTURE_TEST_CASE(check_with_gpu, ParticleFactory, + *boost::unit_test::precondition(has_gpu)) { + auto const rank = boost::mpi::communicator().rank(); + check_uninitialized_device_pointers(); BOOST_CHECK_EQUAL(espresso::system.npart_gpu(), 0); - auto const pid = 1; - mpi_make_new_particle(pid, {0., 0., 0.}); - set_particle_type(pid, 0); + auto const p_id = 1; + create_particle({0., 0., 0.}, p_id, 0); BOOST_CHECK_EQUAL(espresso::system.npart_gpu(), 0); espresso::system.update(); BOOST_CHECK_EQUAL(espresso::system.npart_gpu(), 0); espresso::system.requestParticleStructGpu(); espresso::system.update(); - BOOST_CHECK_EQUAL(espresso::system.npart_gpu(), 1); + BOOST_CHECK_EQUAL(espresso::system.npart_gpu(), (rank == 0) ? 1 : 0); // check position split BOOST_TEST(espresso::system.hasRGpu()); espresso::system.requestRGpu(); espresso::system.update(); - BOOST_TEST(espresso::system.rGpuBegin() != nullptr); + if (rank == 0) { + BOOST_TEST(espresso::system.rGpuBegin() != nullptr); + } else { + BOOST_TEST(espresso::system.rGpuBegin() == nullptr); + } // check force split BOOST_TEST(espresso::system.hasFGpu()); espresso::system.requestFGpu(); espresso::system.update(); - BOOST_TEST(espresso::system.fGpuBegin() != nullptr); + if (rank == 0) { + BOOST_TEST(espresso::system.fGpuBegin() != nullptr); + } else { + BOOST_TEST(espresso::system.fGpuBegin() == nullptr); + } // check torque split #ifdef ROTATION BOOST_CHECK(espresso::system.hasTorqueGpu()); espresso::system.requestTorqueGpu(); espresso::system.update(); - BOOST_TEST(espresso::system.torqueGpuBegin() != nullptr); + if (rank == 0) { + BOOST_TEST(espresso::system.torqueGpuBegin() != nullptr); + } else { + BOOST_TEST(espresso::system.torqueGpuBegin() == nullptr); + } #else BOOST_CHECK(!espresso::system.hasTorqueGpu()); BOOST_CHECK_THROW(espresso::system.requestTorqueGpu(), std::runtime_error); @@ -109,7 +124,11 @@ BOOST_AUTO_TEST_CASE(check_with_gpu, *boost::unit_test::precondition(has_gpu)) { BOOST_CHECK(espresso::system.hasQGpu()); espresso::system.requestQGpu(); espresso::system.update(); - BOOST_TEST(espresso::system.qGpuBegin() != nullptr); + if (rank == 0) { + BOOST_TEST(espresso::system.qGpuBegin() != nullptr); + } else { + BOOST_TEST(espresso::system.qGpuBegin() == nullptr); + } #else BOOST_CHECK(!espresso::system.hasQGpu()); BOOST_CHECK_THROW(espresso::system.requestQGpu(), std::runtime_error); @@ -120,14 +139,18 @@ BOOST_AUTO_TEST_CASE(check_with_gpu, *boost::unit_test::precondition(has_gpu)) { BOOST_CHECK(espresso::system.hasDipGpu()); espresso::system.requestDipGpu(); espresso::system.update(); - BOOST_TEST(espresso::system.dipGpuBegin() != nullptr); + if (rank == 0) { + BOOST_TEST(espresso::system.dipGpuBegin() != nullptr); + } else { + BOOST_TEST(espresso::system.dipGpuBegin() == nullptr); + } #else BOOST_CHECK(!espresso::system.hasDipGpu()); BOOST_CHECK_THROW(espresso::system.requestDipGpu(), std::runtime_error); #endif // clear device memory - remove_particle(pid); + remove_particle(p_id); espresso::system.update(); BOOST_CHECK_EQUAL(espresso::system.npart_gpu(), 0); } @@ -154,15 +177,8 @@ BOOST_AUTO_TEST_CASE(check_without_cuda) { int main(int argc, char **argv) { auto mpi_env = mpi_init(argc, argv); - // this unit test only runs on 1 core - assert(boost::mpi::communicator().size() == 1); - // initialize the MpiCallbacks framework Communication::init(mpi_env); -#ifdef VIRTUAL_SITES - set_virtual_sites(std::make_shared()); -#endif - mpi_loop(); espresso::system.init(); diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 06e303e624f..b653ef1b69d 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -27,7 +27,6 @@ namespace utf = boost::unit_test; #include "ParticleFactory.hpp" #include "EspressoSystemStandAlone.hpp" -#include "MpiCallbacks.hpp" #include "Particle.hpp" #include "accumulators/TimeSeries.hpp" #include "bonded_interactions/bonded_interaction_utils.hpp" @@ -43,7 +42,6 @@ namespace utf = boost::unit_test; #include "integrate.hpp" #include "nonbonded_interactions/lj.hpp" #include "observables/ParticleVelocities.hpp" -#include "particle_data.hpp" #include "particle_node.hpp" #include @@ -52,14 +50,15 @@ namespace utf = boost::unit_test; #include #include +#include #include -#include #include #include #include #include #include +#include #include #include #include @@ -69,103 +68,42 @@ namespace espresso { static std::unique_ptr system; } // namespace espresso -/** Decorator to run a unit test only on the head node. */ -struct if_head_node { - boost::test_tools::assertion_result operator()(utf::test_unit_id) { - return world.rank() == 0; - } - -private: - boost::mpi::communicator world; -}; - -static void mpi_create_bonds_local(int harm_bond_id, int fene_bond_id) { - // set up a harmonic bond - auto const harm_bond = HarmonicBond(200.0, 0.3, 1.0); - auto const harm_bond_ia = std::make_shared(harm_bond); - bonded_ia_params.insert(harm_bond_id, harm_bond_ia); - // set up a FENE bond - auto const fene_bond = FeneBond(300.0, 1.0, 0.3); - auto const fene_bond_ia = std::make_shared(fene_bond); - bonded_ia_params.insert(fene_bond_id, fene_bond_ia); -} - -REGISTER_CALLBACK(mpi_create_bonds_local) - -static void mpi_create_bonds(int harm_bond_id, int fene_bond_id) { - mpi_call_all(mpi_create_bonds_local, harm_bond_id, fene_bond_id); +static void remove_translational_motion() { + Galilei{}.kill_particle_motion(false); } -static void mpi_add_bond_local(int p_id, int bond_id, - std::vector const &partner_ids) { +static auto copy_particle_to_head_node(boost::mpi::communicator const &comm, + int p_id) { + boost::optional result{}; auto p = ::cell_structure.get_local_particle(p_id); - if (p != nullptr and not p->is_ghost()) { - p->bonds().insert(BondView(bond_id, partner_ids)); + if (p and not p->is_ghost()) { + if (comm.rank() == 0) { + result = *p; + } else { + comm.send(0, p_id, *p); + } } - on_particle_change(); -} - -REGISTER_CALLBACK(mpi_add_bond_local) - -static void mpi_add_bond(int p_id, int bond_id, - std::vector const &partner_ids) { - mpi_call_all(mpi_add_bond_local, p_id, bond_id, partner_ids); -} - -static void mpi_remove_translational_motion_local() { - Galilei{}.kill_particle_motion(false); -} - -REGISTER_CALLBACK(mpi_remove_translational_motion_local) - -static void mpi_remove_translational_motion() { - mpi_call_all(mpi_remove_translational_motion_local); -} - -static auto mpi_calculate_energy() { - return mpi_call(Communication::Result::main_rank, calculate_energy); -} - -#ifdef P3M -static void mpi_set_tuned_p3m_local(double prefactor) { - auto p3m = P3MParameters{false, - 0.0, - 3.5, - Utils::Vector3i::broadcast(8), - Utils::Vector3d::broadcast(0.5), - 5, - 0.654, - 1e-3}; - auto solver = - std::make_shared(std::move(p3m), prefactor, 1, false, true); - ::Coulomb::add_actor(solver); -} - -REGISTER_CALLBACK(mpi_set_tuned_p3m_local) - -static void mpi_set_tuned_p3m(double prefactor) { - mpi_call_all(mpi_set_tuned_p3m_local, prefactor); -} -#endif // P3M - -#ifdef LENNARD_JONES -void mpi_set_lj_local(int key, double eps, double sig, double cut, - double offset, double min, double shift) { - LJ_Parameters lj{eps, sig, cut, offset, min, shift}; - ::nonbonded_ia_params[key]->lj = lj; - on_non_bonded_ia_change(); + if (comm.rank() == 0 and not result) { + Particle p{}; + comm.recv(boost::mpi::any_source, p_id, p); + result = p; + } + return result; } -REGISTER_CALLBACK(mpi_set_lj_local) -#endif // LENNARD_JONES - -BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, - *utf::precondition(if_head_node())) { +BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { constexpr auto tol = 100. * std::numeric_limits::epsilon(); + auto const comm = boost::mpi::communicator(); + auto const rank = comm.rank(); + auto const n_nodes = comm.size(); - auto const box_l = 8.; + auto const box_l = 12.; auto const box_center = box_l / 2.; + auto const time_step = 0.001; + auto const skin = 0.4; espresso::system->set_box_l(Utils::Vector3d::broadcast(box_l)); + espresso::system->set_time_step(time_step); + espresso::system->set_skin(skin); // particle properties auto const pid1 = 9; @@ -183,25 +121,20 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, create_particle(start_positions.at(pid1), pid1, type_a); create_particle(start_positions.at(pid2), pid2, type_b); create_particle(start_positions.at(pid3), pid3, type_b); - - { - boost::mpi::communicator world; - auto const n_nodes = world.size(); - if (n_nodes % 2 == 0) { - BOOST_REQUIRE_EQUAL(get_particle_node(pid1), 0); - BOOST_REQUIRE_GE(get_particle_node(pid2), 1); - BOOST_REQUIRE_GE(get_particle_node(pid3), 1); - } + if (n_nodes % 2 == 0) { + BOOST_REQUIRE_EQUAL(get_particle_node_parallel(pid1), rank ? -1 : 0); + BOOST_REQUIRE_GE(get_particle_node_parallel(pid2), rank ? -1 : 1); + BOOST_REQUIRE_GE(get_particle_node_parallel(pid3), rank ? -1 : 1); } auto const reset_particle_positions = [&start_positions]() { for (auto const &kv : start_positions) { - mpi_set_particle_pos(kv.first, kv.second); + set_particle_pos(kv.first, kv.second); } }; // check accumulators - { + if (n_nodes == 1) { auto const pids = std::vector{pid2}; auto obs = std::make_shared(pids); auto acc = Accumulators::TimeSeries(obs, 1); @@ -211,7 +144,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, BOOST_REQUIRE_EQUAL_COLLECTIONS(obs_shape.begin(), obs_shape.end(), ref_shape.begin(), ref_shape.end()); - mpi_remove_translational_motion(); + remove_translational_motion(); for (int i = 0; i < 5; ++i) { set_particle_v(pid2, {static_cast(i), 0., 0.}); @@ -229,14 +162,16 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, // check kinetic energy { - mpi_remove_translational_motion(); + remove_translational_motion(); for (int i = 0; i < 5; ++i) { set_particle_v(pid2, {static_cast(i), 0., 0.}); - auto const &p = get_particle_data(pid2); - auto const kinetic_energy = 0.5 * p.mass() * p.v().norm2(); - auto const obs_energy = mpi_calculate_energy(); - BOOST_CHECK_CLOSE(obs_energy->kinetic[0], kinetic_energy, tol); - BOOST_CHECK_CLOSE(mpi_observable_compute_energy(), kinetic_energy, tol); + auto const obs_energy = calculate_energy(); + auto const p_opt = copy_particle_to_head_node(comm, pid2); + if (rank == 0) { + auto const &p = *p_opt; + auto const kinetic_energy = 0.5 * p.mass() * p.v().norm2(); + BOOST_CHECK_CLOSE(obs_energy->kinetic[0], kinetic_energy, tol); + } } } @@ -256,8 +191,10 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, make_particle_type_exist(std::max(type_a, type_b)); auto const key_ab = get_ia_param_key(type_a, type_b); auto const key_bb = get_ia_param_key(type_b, type_b); - mpi_call_all(mpi_set_lj_local, key_ab, eps, sig, cut, offset, min, shift); - mpi_call_all(mpi_set_lj_local, key_bb, eps, sig, cut, offset, min, shift); + LJ_Parameters lj{eps, sig, cut, offset, min, shift}; + ::nonbonded_ia_params[key_ab]->lj = lj; + ::nonbonded_ia_params[key_bb]->lj = lj; + on_non_bonded_ia_change(); // matrix indices and reference energy value auto const max_type = type_b + 1; @@ -268,12 +205,14 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift); // measure energies - auto const obs_energy = mpi_calculate_energy(); - for (int i = 0; i < n_pairs; ++i) { - auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.; - auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.; - BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 500. * tol); - BOOST_CHECK_CLOSE(obs_energy->non_bonded_intra[i], ref_intra, 500. * tol); + auto const obs_energy = calculate_energy(); + if (rank == 0) { + for (int i = 0; i < n_pairs; ++i) { + auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.; + auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.; + BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 1e-10); + BOOST_CHECK_CLOSE(obs_energy->non_bonded_intra[i], ref_intra, 1e-10); + } } } #endif // LENNARD_JONES @@ -286,53 +225,79 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, auto const harm_bond_id = 0; auto const none_bond_id = 1; auto const fene_bond_id = 2; - mpi_create_bonds(harm_bond_id, fene_bond_id); + { + // set up a harmonic bond + auto const bond = HarmonicBond(200.0, 0.3, 1.0); + auto const bond_ia = std::make_shared(bond); + bonded_ia_params.insert(harm_bond_id, bond_ia); + } + { + // set up a FENE bond + auto const bond = FeneBond(300.0, 1.0, 0.3); + auto const bond_ia = std::make_shared(bond); + bonded_ia_params.insert(fene_bond_id, bond_ia); + } auto const &harm_bond = *boost::get(bonded_ia_params.at(harm_bond_id).get()); auto const &fene_bond = *boost::get(bonded_ia_params.at(fene_bond_id).get()); - mpi_add_bond(pid2, harm_bond_id, {pid1}); - mpi_add_bond(pid2, fene_bond_id, {pid3}); + insert_particle_bond(pid2, harm_bond_id, {pid1}); + insert_particle_bond(pid2, fene_bond_id, {pid3}); // measure energies - auto const obs_energy = mpi_calculate_energy(); - auto const none_energy = 0.0; - auto const harm_energy = 0.5 * harm_bond.k * Utils::sqr(harm_bond.r - dist); - auto const fene_energy = - -0.5 * fene_bond.k * Utils::sqr(fene_bond.drmax) * - std::log(1.0 - Utils::sqr((dist - fene_bond.r0) / fene_bond.drmax)); - BOOST_CHECK_CLOSE(obs_energy->bonded[none_bond_id], none_energy, 0.0); - BOOST_CHECK_CLOSE(obs_energy->bonded[harm_bond_id], harm_energy, 40. * tol); - BOOST_CHECK_CLOSE(obs_energy->bonded[fene_bond_id], fene_energy, 40. * tol); + auto const obs_energy = calculate_energy(); + if (rank == 0) { + auto const none_energy = 0.0; + auto const harm_energy = + 0.5 * harm_bond.k * Utils::sqr(harm_bond.r - dist); + auto const fene_energy = + -0.5 * fene_bond.k * Utils::sqr(fene_bond.drmax) * + std::log(1.0 - Utils::sqr((dist - fene_bond.r0) / fene_bond.drmax)); + BOOST_CHECK_CLOSE(obs_energy->bonded[none_bond_id], none_energy, 0.0); + BOOST_CHECK_CLOSE(obs_energy->bonded[harm_bond_id], harm_energy, 1e-10); + BOOST_CHECK_CLOSE(obs_energy->bonded[fene_bond_id], fene_energy, 1e-10); + } } // check electrostatics #ifdef P3M { // add charges - set_particle_q(pid1, 1.); - set_particle_q(pid2, -1.); + set_particle_property(pid1, &Particle::q, +1.); + set_particle_property(pid2, &Particle::q, -1.); // set up P3M auto const prefactor = 2.; - mpi_set_tuned_p3m(prefactor); + auto p3m = P3MParameters{false, + 0.0, + 3.5, + Utils::Vector3i::broadcast(12), + Utils::Vector3d::broadcast(0.5), + 5, + 0.615, + 1e-3}; + auto solver = + std::make_shared(std::move(p3m), prefactor, 1, false, true); + ::Coulomb::add_actor(solver); // measure energies auto const step = 0.02; - auto const pos1 = get_particle_data(pid1).pos(); + auto const pos1 = start_positions.at(pid1); Utils::Vector3d pos2{box_center, box_center - 0.1, 1.0}; for (int i = 0; i < 10; ++i) { // move particle pos2[0] += step; - mpi_set_particle_pos(pid2, pos2); + set_particle_pos(pid2, pos2); auto const r = (pos2 - pos1).norm(); // check P3M energy - auto const obs_energy = mpi_calculate_energy(); - // at very short distances, the real-space contribution to - // the energy is much larger than the k-space contribution - auto const energy_ref = -prefactor / r; - auto const energy_p3m = obs_energy->coulomb[0] + obs_energy->coulomb[1]; - BOOST_CHECK_CLOSE(energy_p3m, energy_ref, 0.01); + auto const obs_energy = calculate_energy(); + if (rank == 0) { + // at very short distances, the real-space contribution to + // the energy is much larger than the k-space contribution + auto const energy_ref = -prefactor / r; + auto const energy_p3m = obs_energy->coulomb[0] + obs_energy->coulomb[1]; + BOOST_CHECK_CLOSE(energy_p3m, energy_ref, 0.002); + } } } #endif // P3M @@ -340,37 +305,38 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, // check integration { // set up velocity-Verlet integrator - auto const time_step = 0.001; - auto const skin = 0.4; - espresso::system->set_time_step(time_step); - espresso::system->set_skin(skin); set_integ_switch(INTEG_METHOD_NVT); // reset system - mpi_remove_translational_motion(); + remove_translational_motion(); reset_particle_positions(); // recalculate forces without propagating the system - mpi_integrate(0, 0); + integrate(0, 0); // particles are arranged in a right triangle - auto const &p1 = get_particle_data(pid1); - auto const &p2 = get_particle_data(pid2); - auto const &p3 = get_particle_data(pid3); - // forces are symmetric - BOOST_CHECK_CLOSE(p1.force()[0], -p2.force()[0], tol); - BOOST_CHECK_CLOSE(p3.force()[1], -p2.force()[1], tol); - // periodic image contributions to the electrostatic force are negligible - BOOST_CHECK_LE(std::abs(p1.force()[1]), tol); - BOOST_CHECK_LE(std::abs(p1.force()[2]), tol); - BOOST_CHECK_LE(std::abs(p2.force()[2]), tol); - // zero long-range contribution for uncharged particles - BOOST_CHECK_EQUAL(p3.force()[0], 0.); - BOOST_CHECK_EQUAL(p3.force()[2], 0.); - // velocities are not propagated - BOOST_CHECK_EQUAL(p1.v().norm(), 0.); - BOOST_CHECK_EQUAL(p2.v().norm(), 0.); - BOOST_CHECK_EQUAL(p3.v().norm(), 0.); + auto const p1_opt = copy_particle_to_head_node(comm, pid1); + auto const p2_opt = copy_particle_to_head_node(comm, pid2); + auto const p3_opt = copy_particle_to_head_node(comm, pid3); + if (rank == 0) { + auto const &p1 = *p1_opt; + auto const &p2 = *p2_opt; + auto const &p3 = *p3_opt; + // forces are symmetric + BOOST_CHECK_CLOSE(p1.force()[0], -p2.force()[0], tol); + BOOST_CHECK_CLOSE(p3.force()[1], -p2.force()[1], tol); + // periodic image contributions to the electrostatic force are negligible + BOOST_CHECK_LE(std::abs(p1.force()[1]), tol); + BOOST_CHECK_LE(std::abs(p1.force()[2]), tol); + BOOST_CHECK_LE(std::abs(p2.force()[2]), tol); + // zero long-range contribution for uncharged particles + BOOST_CHECK_EQUAL(p3.force()[0], 0.); + BOOST_CHECK_EQUAL(p3.force()[2], 0.); + // velocities are not propagated + BOOST_CHECK_EQUAL(p1.v().norm(), 0.); + BOOST_CHECK_EQUAL(p2.v().norm(), 0.); + BOOST_CHECK_EQUAL(p3.v().norm(), 0.); + } // check integrated trajectory; the time step is chosen // small enough so that particles don't travel too far @@ -381,19 +347,35 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, for (int i = 0; i < 10; ++i) { std::unordered_map expected; for (auto pid : pids) { - auto p = get_particle_data(pid); - p.v() += 0.5 * time_step * p.force() / p.mass(); - p.pos() += time_step * p.v(); - expected[pid] = p.pos(); + auto p_opt = copy_particle_to_head_node(comm, pid); + if (rank == 0) { + auto &p = *p_opt; + p.v() += 0.5 * time_step * p.force() / p.mass(); + p.pos() += time_step * p.v(); + expected[pid] = p.pos(); + } } - mpi_integrate(1, 0); + integrate(1, 0); for (auto pid : pids) { - auto const &p = get_particle_data(pid); - BOOST_CHECK_LE((p.pos() - expected[pid]).norm(), tol); - assert((p.pos() - pos_com).norm() < 0.5); + auto const p_opt = copy_particle_to_head_node(comm, pid); + if (rank == 0) { + auto &p = *p_opt; + BOOST_CHECK_LE((p.pos() - expected[pid]).norm(), tol); + assert((p.pos() - pos_com).norm() < 0.5); + } } } } + + // check exceptions + { + BOOST_CHECK_THROW(get_particle_node_parallel(-1), std::domain_error); + if (rank == 0) { + BOOST_CHECK_THROW(get_particle_node_parallel(12345), std::runtime_error); + } else { + get_particle_node_parallel(12345); + } + } } int main(int argc, char **argv) { diff --git a/src/core/unit_tests/ParticleFactory.hpp b/src/core/unit_tests/ParticleFactory.hpp index 696435549e7..697d8d3d7cf 100644 --- a/src/core/unit_tests/ParticleFactory.hpp +++ b/src/core/unit_tests/ParticleFactory.hpp @@ -19,7 +19,9 @@ #ifndef REACTION_ENSEMBLE_TESTS_PARTICLE_FACTORY_HPP #define REACTION_ENSEMBLE_TESTS_PARTICLE_FACTORY_HPP -#include "particle_data.hpp" +#include "BondList.hpp" +#include "cells.hpp" +#include "event.hpp" #include "particle_node.hpp" #include @@ -36,10 +38,37 @@ struct ParticleFactory { } } - void create_particle(Utils::Vector3d const &pos, int pid, int type) { - mpi_make_new_particle(pid, pos); - set_particle_type(pid, type); - particle_cache.emplace_back(pid); + void create_particle(Utils::Vector3d const &pos, int p_id, int type) { + ::make_new_particle(p_id, pos); + set_particle_property(p_id, &Particle::type, type); + on_particle_type_change(p_id, type_tracking::new_part, type); + particle_cache.emplace_back(p_id); + } + + void set_particle_type(int p_id, int type) const { + set_particle_property(p_id, &Particle::type, type); + on_particle_type_change(p_id, type_tracking::any_type, type); + } + + void set_particle_v(int p_id, Utils::Vector3d const &vel) const { + set_particle_property(p_id, &Particle::v, vel); + } + + void insert_particle_bond(int p_id, int bond_id, + std::vector const &partner_ids) const { + auto p = ::cell_structure.get_local_particle(p_id); + if (p != nullptr and not p->is_ghost()) { + p->bonds().insert(BondView(bond_id, partner_ids)); + } + on_particle_change(); + } + + template + void set_particle_property(int p_id, T &(Particle::*setter)(), + T const &value) const { + if (auto p = ::cell_structure.get_local_particle(p_id)) { + (p->*setter)() = value; + } } private: diff --git a/src/core/unit_tests/VerletCriterion_test.cpp b/src/core/unit_tests/VerletCriterion_test.cpp index e84e9b595a3..d46103a11a1 100644 --- a/src/core/unit_tests/VerletCriterion_test.cpp +++ b/src/core/unit_tests/VerletCriterion_test.cpp @@ -25,7 +25,6 @@ #include "cell_system/CellStructure.hpp" #include "config/config.hpp" #include "nonbonded_interactions/VerletCriterion.hpp" -#include "particle_data.hpp" BOOST_AUTO_TEST_CASE(VerletCriterion_test) { auto constexpr skin = 0.4; diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp index 660463ba75f..78e4c420f1a 100644 --- a/src/core/unit_tests/Verlet_list_test.cpp +++ b/src/core/unit_tests/Verlet_list_test.cpp @@ -32,10 +32,10 @@ namespace utf = boost::unit_test; namespace bdata = boost::unit_test::data; +#include "ParticleFactory.hpp" + #include "EspressoSystemStandAlone.hpp" -#include "MpiCallbacks.hpp" #include "Particle.hpp" -#include "ParticleFactory.hpp" #include "cells.hpp" #include "communication.hpp" #include "event.hpp" @@ -61,34 +61,31 @@ namespace espresso { static std::unique_ptr system; } // namespace espresso -/** Decorator to run a unit test only on the head node. */ -struct if_head_node { - boost::test_tools::assertion_result operator()(utf::test_unit_id) { - return world.rank() == 0; +static auto copy_particle_to_head_node(boost::mpi::communicator const &comm, + int p_id) { + boost::optional result{}; + auto p = ::cell_structure.get_local_particle(p_id); + if (p and not p->is_ghost()) { + if (comm.rank() == 0) { + result = *p; + } else { + comm.send(0, p_id, *p); + } } - -private: - boost::mpi::communicator world; -}; - -#ifdef EXTERNAL_FORCES -static void mpi_set_particle_ext_force_local(int p_id, - Utils::Vector3d const &ext_force) { - if (auto p = ::cell_structure.get_local_particle(p_id)) { - p->ext_force() = ext_force; + if (comm.rank() == 0 and not result) { + Particle p{}; + comm.recv(boost::mpi::any_source, p_id, p); + result = p; } - on_particle_change(); + return result; } -REGISTER_CALLBACK(mpi_set_particle_ext_force_local) -#endif // EXTERNAL_FORCES - namespace Testing { /** * Helper class to setup an integrator and particle properties such that the * particle is on a collision course with another particle on the x-axis. */ -struct IntegratorHelper { +struct IntegratorHelper : public ParticleFactory { IntegratorHelper() = default; virtual ~IntegratorHelper() = default; /** Set integrator parameters. */ @@ -102,35 +99,25 @@ struct IntegratorHelper { } }; -void mpi_set_integrator_sd_local() { - register_integrator(SteepestDescentParameters(0., 0.01, 100.)); - set_integ_switch(INTEG_METHOD_STEEPEST_DESCENT); -} - -REGISTER_CALLBACK(mpi_set_integrator_sd_local) - #ifdef EXTERNAL_FORCES struct : public IntegratorHelper { void set_integrator() const override { - mpi_set_thermo_switch(THERMO_OFF); - mpi_call_all(mpi_set_integrator_sd_local); + mpi_set_thermo_switch_local(THERMO_OFF); + register_integrator(SteepestDescentParameters(0., 0.01, 100.)); + set_integ_switch(INTEG_METHOD_STEEPEST_DESCENT); } void set_particle_properties(int pid) const override { - mpi_call_all(mpi_set_particle_ext_force_local, pid, - Utils::Vector3d{20., 0., 0.}); + set_particle_property(pid, &Particle::ext_force, + Utils::Vector3d{20., 0., 0.}); } char const *name() const override { return "SteepestDescent"; } } steepest_descent; #endif // EXTERNAL_FORCES -void mpi_set_integrator_vv_local() { set_integ_switch(INTEG_METHOD_NVT); } - -REGISTER_CALLBACK(mpi_set_integrator_vv_local) - struct : public IntegratorHelper { void set_integrator() const override { - mpi_set_thermo_switch(THERMO_OFF); - mpi_call_all(mpi_set_integrator_vv_local); + mpi_set_thermo_switch_local(THERMO_OFF); + set_integ_switch(INTEG_METHOD_NVT); } void set_particle_properties(int pid) const override { set_particle_v(pid, {20., 0., 0.}); @@ -139,20 +126,14 @@ struct : public IntegratorHelper { } velocity_verlet; #ifdef NPT -void mpi_set_integrator_npt_local() { - ::nptiso = NptIsoParameters(1., 1e9, {true, true, true}, true); - set_integ_switch(INTEG_METHOD_NPT_ISO); -} - -REGISTER_CALLBACK(mpi_set_integrator_npt_local) - struct : public IntegratorHelper { void set_integrator() const override { - mpi_call_all(mpi_set_integrator_npt_local); - mpi_set_temperature(1.); - npt_iso_set_rng_seed(0); - mpi_set_thermo_switch(thermo_switch | THERMO_NPT_ISO); - mpi_set_nptiso_gammas(0., 0.); // disable box volume change + ::nptiso = NptIsoParameters(1., 1e9, {true, true, true}, true); + set_integ_switch(INTEG_METHOD_NPT_ISO); + mpi_set_temperature_local(1.); + mpi_npt_iso_set_rng_seed(0); + mpi_set_thermo_switch_local(thermo_switch | THERMO_NPT_ISO); + mpi_set_nptiso_gammas_local(0., 0.); // disable box volume change } void set_particle_properties(int pid) const override { set_particle_v(pid, {20., 0., 0.}); @@ -163,15 +144,6 @@ struct : public IntegratorHelper { } // namespace Testing -void mpi_set_lj_local(int key, double eps, double sig, double cut, - double offset, double min, double shift) { - LJ_Parameters lj{eps, sig, cut, offset, min, shift}; - ::nonbonded_ia_params[key]->lj = lj; - on_non_bonded_ia_change(); -} - -REGISTER_CALLBACK(mpi_set_lj_local) - inline double get_dist_from_last_verlet_update(Particle const &p) { return (p.pos() - p.pos_at_last_verlet_update()).norm(); } @@ -192,11 +164,12 @@ auto const propagators = #endif // EXTERNAL_FORCES }; -BOOST_TEST_DECORATOR(*utf::precondition(if_head_node())) BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, bdata::make(node_grids) * bdata::make(propagators), node_grid, integration_helper) { constexpr auto tol = 100. * std::numeric_limits::epsilon(); + auto const comm = boost::mpi::communicator(); + auto const rank = comm.rank(); auto const box_l = 8.; espresso::system->set_box_l(Utils::Vector3d::broadcast(box_l)); @@ -219,7 +192,9 @@ BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, auto const cut = r_off + 1e-3; make_particle_type_exist(1); auto const key = get_ia_param_key(0, 1); - mpi_call_all(mpi_set_lj_local, key, eps, sig, cut, offset, min, shift); + LJ_Parameters lj{eps, sig, cut, offset, min, shift}; + ::nonbonded_ia_params[key]->lj = lj; + on_non_bonded_ia_change(); // set up velocity-Verlet integrator auto const time_step = 0.01; @@ -237,9 +212,10 @@ BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, create_particle({2. - 0.10, 1., 1.}, pid1, 0); create_particle({4. + 0.15, 1., 1.}, pid2, 1); create_particle({2. - 0.10, 5., 1.}, pid3, 0); - BOOST_REQUIRE_EQUAL(get_particle_node(pid1), 0); - BOOST_REQUIRE_EQUAL(get_particle_node(pid2), 2); - BOOST_REQUIRE_EQUAL(get_particle_node(pid3), (node_grid[0] == 4) ? 0 : 1); + BOOST_REQUIRE_EQUAL(get_particle_node_parallel(pid1), rank ? -1 : 0); + BOOST_REQUIRE_EQUAL(get_particle_node_parallel(pid2), rank ? -1 : 2); + BOOST_REQUIRE_EQUAL(get_particle_node_parallel(pid3), + rank ? -1 : ((node_grid[0] == 4) ? 0 : 1)); // check that particles in different cells will eventually interact during // normal integration (the integration helper sets a collision course) @@ -248,44 +224,61 @@ BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, // integrate until both particles are closer than cutoff { - mpi_integrate(11, 0); - auto const &p1 = get_particle_data(pid1); - auto const &p2 = get_particle_data(pid2); - BOOST_REQUIRE_LT(get_dist_from_pair(p1, p2), cut); + integrate(11, 0); + auto const p1_opt = copy_particle_to_head_node(comm, pid1); + auto const p2_opt = copy_particle_to_head_node(comm, pid2); + if (rank == 0) { + auto const &p1 = *p1_opt; + auto const &p2 = *p2_opt; + BOOST_REQUIRE_LT(get_dist_from_pair(p1, p2), cut); + } } // check forces and Verlet update { - mpi_integrate(1, 0); - auto const &p1 = get_particle_data(pid1); + integrate(1, 0); + auto const p1_opt = copy_particle_to_head_node(comm, pid1); +#ifdef EXTERNAL_FORCES + auto const p2_opt = copy_particle_to_head_node(comm, pid2); +#endif // EXTERNAL_FORCES + if (rank == 0) { + auto const &p1 = *p1_opt; + auto const &p2 = *p2_opt; #ifdef EXTERNAL_FORCES - auto const &p2 = get_particle_data(pid2); - BOOST_CHECK_CLOSE(p1.force()[0] - p1.ext_force()[0], 480., 1e-9); + BOOST_CHECK_CLOSE(p1.force()[0] - p1.ext_force()[0], 480., 1e-9); #endif // EXTERNAL_FORCES - BOOST_CHECK_CLOSE(p1.force()[1], 0., tol); - BOOST_CHECK_CLOSE(p1.force()[2], 0., tol); + BOOST_CHECK_CLOSE(p1.force()[1], 0., tol); + BOOST_CHECK_CLOSE(p1.force()[2], 0., tol); #ifdef EXTERNAL_FORCES - BOOST_TEST(p1.force() - p1.ext_force() == -p2.force(), - boost::test_tools::per_element()); + BOOST_TEST(p1.force() - p1.ext_force() == -p2.force(), + boost::test_tools::per_element()); #endif // EXTERNAL_FORCES - BOOST_CHECK_LT(get_dist_from_last_verlet_update(p1), skin / 2.); + BOOST_CHECK_LT(get_dist_from_last_verlet_update(p1), skin / 2.); + } } } // check that particles in different cells will interact when manually - // placed next to each other (@c mpi_set_particle_pos() resorts particles) + // placed next to each other (@c set_particle_pos() resorts particles) { - mpi_set_particle_pos(pid3, {4. + 0.10, 1., 1.0}); + ::set_particle_pos(pid3, {4. + 0.10, 1., 1.0}); { - auto const &p2 = get_particle_data(pid2); - auto const &p3 = get_particle_data(pid3); - BOOST_REQUIRE_LT(get_dist_from_pair(p2, p3), cut); - BOOST_CHECK_GT(get_dist_from_last_verlet_update(p3), skin / 2.); + auto const p2_opt = copy_particle_to_head_node(comm, pid2); + auto const p3_opt = copy_particle_to_head_node(comm, pid3); + if (rank == 0) { + auto const &p2 = *p2_opt; + auto const &p3 = *p3_opt; + BOOST_REQUIRE_LT(get_dist_from_pair(p2, p3), cut); + BOOST_CHECK_GT(get_dist_from_last_verlet_update(p3), skin / 2.); + } } { - mpi_integrate(0, 0); - auto const &p3 = get_particle_data(pid3); - BOOST_CHECK_LT(get_dist_from_last_verlet_update(p3), skin / 2.); + integrate(0, 0); + auto const p3_opt = copy_particle_to_head_node(comm, pid3); + if (rank == 0) { + auto const &p3 = *p3_opt; + BOOST_CHECK_LT(get_dist_from_last_verlet_update(p3), skin / 2.); + } } } } diff --git a/src/script_interface/interactions/NonBondedInteraction.hpp b/src/script_interface/interactions/NonBondedInteraction.hpp index 6ed08ce4906..14cd3953902 100644 --- a/src/script_interface/interactions/NonBondedInteraction.hpp +++ b/src/script_interface/interactions/NonBondedInteraction.hpp @@ -899,8 +899,7 @@ class NonBondedInteractionHandle auto const types = get_value>(params.at("_types")); m_types[0] = std::min(types[0], types[1]); m_types[1] = std::max(types[0], types[1]); - // make type exist - mpi_realloc_ia_params_local(m_types[1] + 1); + make_particle_type_exist(m_types[1]); // create interface objects auto const key = get_ia_param_key(m_types[0], m_types[1]); m_interaction = ::nonbonded_ia_params[key]; diff --git a/src/script_interface/interactions/NonBondedInteractions.hpp b/src/script_interface/interactions/NonBondedInteractions.hpp index 945ac1be2f1..99182b573e2 100644 --- a/src/script_interface/interactions/NonBondedInteractions.hpp +++ b/src/script_interface/interactions/NonBondedInteractions.hpp @@ -69,7 +69,7 @@ class NonBondedInteractions : public ObjectHandle { void do_construct(VariantMap const ¶ms) override { auto const size = ::max_seen_particle_type; - make_particle_type_exist_local(size); + make_particle_type_exist(size); for (int i = 0; i < size; i++) { for (int j = i; j < size; j++) { auto const key = Utils::upper_triangular(i, j, size); @@ -89,7 +89,7 @@ class NonBondedInteractions : public ObjectHandle { } if (name == "insert") { auto const types = get_value>(params.at("key")); - make_particle_type_exist_local(std::max(types[0], types[1])); + make_particle_type_exist(std::max(types[0], types[1])); auto const key = get_ia_param_key(std::min(types[0], types[1]), std::max(types[0], types[1])); auto obj_ptr = get_value>( diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp index d392084113a..4de2ff356a0 100644 --- a/src/script_interface/particle_data/ParticleHandle.cpp +++ b/src/script_interface/particle_data/ParticleHandle.cpp @@ -166,14 +166,14 @@ ParticleHandle::ParticleHandle() { throw std::domain_error( error_msg("type", "must be an integer >= 0")); } - make_particle_type_exist_local(new_type); - on_particle_type_change_parallel(m_pid, old_type, new_type); + make_particle_type_exist(new_type); + on_particle_type_change(m_pid, old_type, new_type); set_particle_property(&Particle::type, value); }, [this]() { return get_particle_data(m_pid).type(); }}, {"pos", [this](Variant const &value) { - mpi_set_particle_pos_local(m_pid, get_value(value)); + set_particle_pos(m_pid, get_value(value)); }, [this]() { auto const p = get_particle_data(m_pid); @@ -581,7 +581,7 @@ Variant ParticleHandle::do_call_method(std::string const &name, if (name == "remove_particle") { context()->parallel_try_catch([&]() { std::ignore = get_real_particle(context()->get_comm(), m_pid); - remove_particle_parallel(m_pid); + remove_particle(m_pid); }); #ifdef VIRTUAL_SITES_RELATIVE } else if (name == "vs_relate_to") { @@ -711,7 +711,7 @@ void ParticleHandle::do_construct(VariantMap const ¶ms) { return params.count(key) == 1; }; m_pid = (has_param("id")) ? get_value(params, "id") - : get_maximal_particle_id_parallel() + 1; + : get_maximal_particle_id() + 1; #ifndef NDEBUG if (!has_param("id")) { @@ -754,7 +754,7 @@ void ParticleHandle::do_construct(VariantMap const ¶ms) { #endif // ROTATION // create a default-constructed particle - mpi_make_new_particle_local(m_pid, pos); + make_new_particle(m_pid, pos); context()->parallel_try_catch([&]() { // set particle properties (filter out read-only and deferred properties) diff --git a/src/script_interface/particle_data/ParticleList.cpp b/src/script_interface/particle_data/ParticleList.cpp index 9552ec16a9f..900c3bbd8c6 100644 --- a/src/script_interface/particle_data/ParticleList.cpp +++ b/src/script_interface/particle_data/ParticleList.cpp @@ -220,6 +220,13 @@ Variant ParticleList::do_call_method(std::string const &name, return {}; } #endif // EXCLUSIONS + if (name == "get_highest_particle_id") { + return get_maximal_particle_id(); + } + if (name == "clear") { + remove_all_particles(); + return {}; + } if (not context()->is_head_node()) { return {}; } @@ -232,9 +239,7 @@ Variant ParticleList::do_call_method(std::string const &name, if (name == "particle_exists") { return particle_exists(get_value(params, "p_id")); } - if (name == "clear") { - remove_all_particles(); - } else if (name == "add_particle") { + if (name == "add_particle") { assert(params.count("bonds") == 0); auto obj = context()->make_shared("Particles::ParticleHandle", params); auto &p_handle = dynamic_cast(*obj); @@ -245,9 +250,6 @@ Variant ParticleList::do_call_method(std::string const &name, #endif // EXCLUSIONS return p_handle.get_parameter("id"); } - if (name == "get_highest_particle_id") { - return get_maximal_particle_id(); - } return {}; } diff --git a/src/script_interface/tests/CMakeLists.txt b/src/script_interface/tests/CMakeLists.txt index 55a699da218..35b14fed086 100644 --- a/src/script_interface/tests/CMakeLists.txt +++ b/src/script_interface/tests/CMakeLists.txt @@ -57,6 +57,8 @@ unit_test(NAME Constraints_test SRC Constraints_test.cpp DEPENDS unit_test(NAME Actors_test SRC Actors_test.cpp DEPENDS espresso::script_interface espresso::core) unit_test(NAME ConstantpHEnsemble_test SRC ConstantpHEnsemble_test.cpp DEPENDS - espresso::core espresso::script_interface Boost::mpi MPI::MPI_CXX) + espresso::core espresso::script_interface Boost::mpi MPI::MPI_CXX + NUM_PROC 2) unit_test(NAME ReactionEnsemble_test SRC ReactionEnsemble_test.cpp DEPENDS - espresso::core espresso::script_interface Boost::mpi MPI::MPI_CXX) + espresso::core espresso::script_interface Boost::mpi MPI::MPI_CXX + NUM_PROC 2) From 4cfc4cec6f21bbbe8e9c36a334221f11e971c983 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 13 Feb 2023 18:15:31 +0100 Subject: [PATCH 3/5] core: Give MC functions more descriptive names MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Pablo Miguel Blanco Andrés --- .../reaction_methods/ReactionAlgorithm.cpp | 10 +++---- .../reaction_methods/ReactionAlgorithm.hpp | 28 +++++++++++-------- src/python/espressomd/reaction_methods.py | 7 +++-- .../reaction_methods/ReactionAlgorithm.cpp | 10 +++---- .../tests/ReactionEnsemble_test.cpp | 10 +++---- .../python/reaction_methods_interface.py | 4 +++ 6 files changed, 39 insertions(+), 30 deletions(-) diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index 78cedbab7ba..f723af7a90d 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -233,7 +233,7 @@ ReactionAlgorithm::get_particle_numbers(SingleReaction const &reaction) const { } std::optional -ReactionAlgorithm::generic_oneway_reaction_part_1(int reaction_id) { +ReactionAlgorithm::create_new_trial_state(int reaction_id) { auto &reaction = *reactions[reaction_id]; reaction.tried_moves++; particle_inside_exclusion_range_touched = false; @@ -253,10 +253,10 @@ ReactionAlgorithm::generic_oneway_reaction_part_1(int reaction_id) { return {E_pot_new}; } -double ReactionAlgorithm::generic_oneway_reaction_part_2(int reaction_id, - double bf, - double E_pot_old, - double E_pot_new) { +double ReactionAlgorithm::make_reaction_mc_move_attempt(int reaction_id, + double bf, + double E_pot_old, + double E_pot_new) { auto const exponential = std::exp(-(E_pot_new - E_pot_old) / kT); auto &reaction = *reactions[reaction_id]; reaction.accumulator_potential_energy_difference_exponential( diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp index 9d61219a3aa..e617dbdebca 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.hpp +++ b/src/core/reaction_methods/ReactionAlgorithm.hpp @@ -177,39 +177,43 @@ class ReactionAlgorithm { public: /** - * Attempt a reaction move and calculate the new potential energy. + * Carry out a reaction MC move and calculate the new potential energy. * Particles are selected without replacement. + * The previous state of the system is cached. * @returns Potential energy of the system after the move. */ - std::optional generic_oneway_reaction_part_1(int reaction_id); + std::optional create_new_trial_state(int reaction_id); /** - * Accept or reject moves made by @ref generic_oneway_reaction_part_1 based - * on a probability acceptance @c bf. + * Accept or reject a reaction MC move made by @ref create_new_trial_state + * based on a probability acceptance @c bf. + * The previous state of the system is either restored from the cache if + * the move is rejected, or cleared from the cache if the move is accepted. * @returns Potential energy of the system after the move was accepted or * rejected. */ - double generic_oneway_reaction_part_2(int reaction_id, double bf, - double E_pot_old, double E_pot_new); + double make_reaction_mc_move_attempt(int reaction_id, double bf, + double E_pot_old, double E_pot_new); /** - * Attempt a global MC move for particles of a given type. + * Attempt displacement MC moves for particles of a given type. * Particles are selected without replacement. * @param type Type of particles to move. * @param n_particles Number of particles of that type to move. * @returns true if all moves were accepted. */ bool make_displacement_mc_move_attempt(int type, int n_particles); + + /** @brief Compute the system potential energy. */ + double calculate_potential_energy() const; + +protected: /** - * Perform a global MC move for particles of a given type. + * Carry out displacement MC moves for particles of a given type. * Particles are selected without replacement. * @param type Type of particles to move. * @param n_particles Number of particles of that type to move. */ void displacement_mc_move(int type, int n_particles); - /** @brief Compute the system potential energy. */ - double calculate_potential_energy() const; - -protected: /** @brief Carry out a chemical reaction and save the old system state. */ void make_reaction_attempt(::ReactionMethods::SingleReaction const &reaction, ParticleChanges &bookkeeping); diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py index c02445736e7..7d6a45fe510 100644 --- a/src/python/espressomd/reaction_methods.py +++ b/src/python/espressomd/reaction_methods.py @@ -476,12 +476,12 @@ def generic_oneway_reaction(self, reaction_id, E_pot_old): """ try: E_pot_new = self.call_method( - "generic_oneway_reaction_part_1", reaction_id=reaction_id) + "create_new_trial_state", reaction_id=reaction_id) if E_pot_new is None: return E_pot_old E_pot_diff = E_pot_new - E_pot_old bf = self.calculate_acceptance_probability(reaction_id, E_pot_diff) - return self.call_method("generic_oneway_reaction_part_2", + return self.call_method("make_reaction_mc_move_attempt", reaction_id=reaction_id, bf=bf, E_pot_new=E_pot_new, E_pot_old=E_pot_old) except BaseException as err: @@ -515,7 +515,8 @@ def get_status(self): reactions_list.append(reaction) return {"reactions": reactions_list, "kT": self.kT, - "exclusion_range": self.exclusion_range} + "exclusion_range": self.exclusion_range, + "exclusion_radius_per_type": self.exclusion_radius_per_type} @script_interface_register diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp index fbe51a7a707..d11c875a4bf 100644 --- a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp +++ b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp @@ -130,26 +130,26 @@ Variant ReactionAlgorithm::do_call_method(std::string const &name, if (name == "potential_energy") { return RE()->calculate_potential_energy(); } - if (name == "generic_oneway_reaction_part_1") { + if (name == "create_new_trial_state") { auto const reaction_id = get_value(params, "reaction_id"); Variant result{}; context()->parallel_try_catch([&]() { - auto const optional = RE()->generic_oneway_reaction_part_1(reaction_id); + auto const optional = RE()->create_new_trial_state(reaction_id); if (optional) { result = *optional; } }); return result; } - if (name == "generic_oneway_reaction_part_2") { + if (name == "make_reaction_mc_move_attempt") { auto const bf = get_value(params, "bf"); auto const E_pot_old = get_value(params, "E_pot_old"); auto const E_pot_new = get_value(params, "E_pot_new"); auto const reaction_id = get_value(params, "reaction_id"); Variant result; context()->parallel_try_catch([&]() { - result = RE()->generic_oneway_reaction_part_2(reaction_id, bf, E_pot_old, - E_pot_new); + result = RE()->make_reaction_mc_move_attempt(reaction_id, bf, E_pot_old, + E_pot_new); }); return result; } diff --git a/src/script_interface/tests/ReactionEnsemble_test.cpp b/src/script_interface/tests/ReactionEnsemble_test.cpp index 664d819fcbf..dc29e554ef5 100644 --- a/src/script_interface/tests/ReactionEnsemble_test.cpp +++ b/src/script_interface/tests/ReactionEnsemble_test.cpp @@ -182,7 +182,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { espresso::system->set_box_l(box_l); // without reactants, no reaction will take place - auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + auto const result = r_algo.create_new_trial_state(reaction_id); BOOST_REQUIRE(not result.has_value()); // the reaction was updated @@ -200,7 +200,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { // system auto const energy_ref = 0.; - auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + auto const result = r_algo.create_new_trial_state(reaction_id); BOOST_REQUIRE(result.has_value()); auto const energy_move = *result; @@ -218,7 +218,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { auto const bf = r_algo_si->calculate_acceptance_probability( reaction, energy_move, {{type_D, 1}, {type_E, 0}}); - auto const energy_end = r_algo.generic_oneway_reaction_part_2( + auto const energy_end = r_algo.make_reaction_mc_move_attempt( reaction_id, bf, 0., energy_move); BOOST_CHECK_CLOSE(energy_end, energy_ref, tol); @@ -231,7 +231,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { } { // attempt a second reaction - auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + auto const result = r_algo.create_new_trial_state(reaction_id); BOOST_REQUIRE(result.has_value()); // verify bookkeeping @@ -247,7 +247,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { // force move to be rejected auto const energy_reject = - r_algo.generic_oneway_reaction_part_2(reaction_id, 0., 0.2, 0.1); + r_algo.make_reaction_mc_move_attempt(reaction_id, 0., 0.2, 0.1); BOOST_CHECK_CLOSE(energy_reject, 0.2, tol); // the reaction was updated diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py index 2f677d08f2e..8f2bfcf4276 100644 --- a/testsuite/python/reaction_methods_interface.py +++ b/testsuite/python/reaction_methods_interface.py @@ -93,6 +93,7 @@ def count_by_type(types): list(method.exclusion_radius_per_type.keys()), [2]) self.assertAlmostEqual( method.exclusion_radius_per_type[2], 0.2, delta=1e-10) + exclusion_radius_per_type = {2: 0.2} self.assertAlmostEqual( method.get_volume(), self.system.volume(), delta=1e-10) method.set_volume(volume=1.) @@ -125,6 +126,9 @@ def count_by_type(types): status = method.get_status() self.assertEqual(status['kT'], kT) self.assertEqual(status['exclusion_range'], exclusion_range) + self.assertEqual( + status['exclusion_radius_per_type'], + exclusion_radius_per_type) self.assertEqual(len(status['reactions']), 2) for reaction_flat, params in zip( status['reactions'], reaction_parameters): From 881d214633957faf044803583f4a78b6ea2da853 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 14 Feb 2023 00:17:34 +0100 Subject: [PATCH 4/5] python: Rename MC method reaction argument MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit For consistency with the integrator.run() method. Co-authored-by: Pablo Miguel Blanco Andrés --- doc/tutorials/constant_pH/constant_pH.ipynb | 4 ++-- maintainer/benchmarks/mc_acid_base_reservoir.py | 6 +++--- samples/grand_canonical.py | 6 +++--- samples/reaction_ensemble_complex_reaction.py | 4 ++-- samples/reaction_methods.py | 2 +- src/python/espressomd/reaction_methods.py | 6 +++--- testsuite/python/constant_pH.py | 4 ++-- testsuite/python/constant_pH_stats.py | 4 ++-- testsuite/python/reaction_bookkeeping.py | 2 +- testsuite/python/reaction_complex.py | 4 ++-- testsuite/python/reaction_ensemble.py | 4 ++-- 11 files changed, 23 insertions(+), 23 deletions(-) diff --git a/doc/tutorials/constant_pH/constant_pH.ipynb b/doc/tutorials/constant_pH/constant_pH.ipynb index 9f78b64d30f..bcfd5bdf993 100644 --- a/doc/tutorials/constant_pH/constant_pH.ipynb +++ b/doc/tutorials/constant_pH/constant_pH.ipynb @@ -686,7 +686,7 @@ "source": [ "```python\n", "def equilibrate_reaction(reaction_steps=1):\n", - " RE.reaction(reaction_steps=reaction_steps)\n", + " RE.reaction(steps=reaction_steps)\n", "```" ] }, @@ -743,7 +743,7 @@ " if USE_WCA and np.random.random() < prob_integration:\n", " system.integrator.run(integration_steps)\n", " # we should do at least one reaction attempt per reactive particle\n", - " RE.reaction(reaction_steps=reaction_steps) \n", + " RE.reaction(steps=reaction_steps)\n", " num_As[i] = system.number_of_particles(type=type_A)\n", "```" ] diff --git a/maintainer/benchmarks/mc_acid_base_reservoir.py b/maintainer/benchmarks/mc_acid_base_reservoir.py index 0557e53b3ee..97e465e8c70 100644 --- a/maintainer/benchmarks/mc_acid_base_reservoir.py +++ b/maintainer/benchmarks/mc_acid_base_reservoir.py @@ -259,7 +259,7 @@ def calc_donnan_coefficient(c_acid, I_res, charge=-1): def equilibrate_reaction(reaction_steps=1): - RE.reaction(reaction_steps=reaction_steps) + RE.reaction(steps=reaction_steps) def report_progress(system, i, next_i): @@ -295,7 +295,7 @@ def report_progress(system, i, next_i): if MC_STEPS_PER_SAMPLE > 0: tick_MC = time.time() - RE.reaction(reaction_steps=MC_STEPS_PER_SAMPLE) + RE.reaction(steps=MC_STEPS_PER_SAMPLE) tock_MC = time.time() t_MC = (tock_MC - tick_MC) / MC_STEPS_PER_SAMPLE @@ -332,7 +332,7 @@ def report_progress(system, i, next_i): for i in range(NUM_SAMPLES): if RUN_INTEGRATION: system.integrator.run(INTEGRATION_STEPS_PER_SAMPLE) - RE.reaction(reaction_steps=MC_STEPS_PER_SAMPLE) + RE.reaction(steps=MC_STEPS_PER_SAMPLE) n_A = system.number_of_particles(type=TYPES['A']) n_As.append(n_A) n_All = len(system.part) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 150050e849f..0d581d863a3 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -101,7 +101,7 @@ # up the simulation RE.set_non_interacting_type(type=max(types) + 1) -RE.reaction(reaction_steps=10000) +RE.reaction(steps=10000) p3m = espressomd.electrostatics.P3M(prefactor=2.0, accuracy=1e-3) system.actors.add(p3m) @@ -134,14 +134,14 @@ system.thermostat.set_langevin(kT=temperature, gamma=.5, seed=42) # MC warmup -RE.reaction(reaction_steps=1000) +RE.reaction(steps=1000) n_int_cycles = 10000 n_int_steps = 600 num_As = [] deviation = None for i in range(n_int_cycles): - RE.reaction(reaction_steps=10) + RE.reaction(steps=10) system.integrator.run(steps=n_int_steps) num_As.append(system.number_of_particles(type=1)) if i > 2 and i % 50 == 0: diff --git a/samples/reaction_ensemble_complex_reaction.py b/samples/reaction_ensemble_complex_reaction.py index cd727a58773..93c2c2805dd 100644 --- a/samples/reaction_ensemble_complex_reaction.py +++ b/samples/reaction_ensemble_complex_reaction.py @@ -100,10 +100,10 @@ RE.set_non_interacting_type(type=max(types) + 1) # warmup -RE.reaction(reaction_steps=200) +RE.reaction(steps=200) for i in range(200): - RE.reaction(reaction_steps=10) + RE.reaction(steps=10) for _type in types: numbers[_type].append(system.number_of_particles(type=_type)) diff --git a/samples/reaction_methods.py b/samples/reaction_methods.py index 220e2675b00..9e1416a3e0d 100644 --- a/samples/reaction_methods.py +++ b/samples/reaction_methods.py @@ -107,7 +107,7 @@ for i in range(10000): - RE.reaction(reaction_steps=1) + RE.reaction(steps=1) if i % 100 == 0: print(f"HA {system.number_of_particles(type=types['HA'])}", f"A- {system.number_of_particles(type=types['A-'])}", diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py index 7d6a45fe510..f539d776187 100644 --- a/src/python/espressomd/reaction_methods.py +++ b/src/python/espressomd/reaction_methods.py @@ -407,19 +407,19 @@ def _rebuild_reaction_cache(self): """ self._reactions_cache = [x for x in self.reactions] - def reaction(self, reaction_steps): + def reaction(self, steps): """ Performs randomly selected reactions. Parameters ---------- - reaction_steps : :obj:`int`, optional + steps : :obj:`int`, optional The number of reactions to be performed at once, defaults to 1. """ self.call_method("setup_bookkeeping_of_empty_pids") E_pot = self.call_method("potential_energy") - for _ in range(reaction_steps): + for _ in range(steps): reaction_id = self.call_method("get_random_reaction_index") E_pot = self.generic_oneway_reaction(reaction_id, E_pot) diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py index 4b9cb5ab644..afbbf692d6a 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH.py @@ -64,12 +64,12 @@ def test_ideal_alpha(self): default_charges=charges_dict) # equilibration - RE.reaction(reaction_steps=800) + RE.reaction(steps=800) # sampling alphas = [] for _ in range(80): - RE.reaction(reaction_steps=15) + RE.reaction(steps=15) num_H = system.number_of_particles(type=types["H+"]) num_HA = system.number_of_particles(type=types["HA"]) num_A = system.number_of_particles(type=types["A-"]) diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py index ad0a6ae4afd..69057f61479 100644 --- a/testsuite/python/constant_pH_stats.py +++ b/testsuite/python/constant_pH_stats.py @@ -80,14 +80,14 @@ def test_ideal_titration_curve(self): # chemical warmup - get close to chemical equilibrium before we start # sampling - RE.reaction(reaction_steps=40 * N0) + RE.reaction(steps=40 * N0) average_NH = 0.0 average_NHA = 0.0 average_NA = 0.0 num_samples = 1000 for _ in range(num_samples): - RE.reaction(reaction_steps=10) + RE.reaction(steps=10) average_NH += system.number_of_particles(type=types["H+"]) average_NHA += system.number_of_particles(type=types["HA"]) average_NA += system.number_of_particles(type=types["A-"]) diff --git a/testsuite/python/reaction_bookkeeping.py b/testsuite/python/reaction_bookkeeping.py index 7ded2e9b8e3..5e02f220a07 100644 --- a/testsuite/python/reaction_bookkeeping.py +++ b/testsuite/python/reaction_bookkeeping.py @@ -88,7 +88,7 @@ def setUpClass(cls): def test_reaction_bookeeping(self): self.widom.calculate_particle_insertion_potential_energy(reaction_id=0) - self.cph.reaction(reaction_steps=100) + self.cph.reaction(steps=100) # Measure the chemical potential for _ in range(50): diff --git a/testsuite/python/reaction_complex.py b/testsuite/python/reaction_complex.py index f4dbf0993ab..9f955270614 100644 --- a/testsuite/python/reaction_complex.py +++ b/testsuite/python/reaction_complex.py @@ -84,11 +84,11 @@ def test_equilibrium_concentrations(self): RE.set_non_interacting_type(type=max(types) + 1) # warmup - RE.reaction(reaction_steps=50) + RE.reaction(steps=50) numbers = {type_A: [], type_B: [], type_C: [], type_D: [], type_E: []} for _ in range(40): - RE.reaction(reaction_steps=6) + RE.reaction(steps=6) for key in types: numbers[key].append(self.system.number_of_particles(type=key)) diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py index 7f4936b3a98..0e3427de867 100644 --- a/testsuite/python/reaction_ensemble.py +++ b/testsuite/python/reaction_ensemble.py @@ -97,14 +97,14 @@ def test_ideal_titration_curve(self): # chemical warmup - get close to chemical equilibrium before we start # sampling - RE.reaction(reaction_steps=5 * N0) + RE.reaction(steps=5 * N0) average_NH = 0.0 average_NHA = 0.0 average_NA = 0.0 num_samples = 300 for _ in range(num_samples): - RE.reaction(reaction_steps=10) + RE.reaction(steps=10) average_NH += system.number_of_particles(type=types["H+"]) average_NHA += system.number_of_particles(type=types["HA"]) average_NA += system.number_of_particles(type=types["A-"]) From a716f36123ca94336fff935d98e42316608d670f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 23 Feb 2023 15:37:00 +0100 Subject: [PATCH 5/5] python: Add pylint check W0109 (duplicate-key) Duplicated dict keys will silently override one another. --- .pylintrc | 1 + src/python/espressomd/reaction_methods.py | 14 ++++---------- 2 files changed, 5 insertions(+), 10 deletions(-) diff --git a/.pylintrc b/.pylintrc index 9932b4f2f97..c6667ebc3aa 100644 --- a/.pylintrc +++ b/.pylintrc @@ -67,6 +67,7 @@ disable=all # multiple time (only on the command line, not in the configuration file where # it should appear only once). See also the "--disable" option for examples. enable=dangerous-default-value, # W0102 + duplicate-key, # W0109 wildcard-import, # W0401 assert-on-tuple, # W0199 unused-import, # W0611 diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py index f539d776187..7cf5eb9ca17 100644 --- a/src/python/espressomd/reaction_methods.py +++ b/src/python/espressomd/reaction_methods.py @@ -503,16 +503,10 @@ def get_status(self): """ self.check_reaction_method() - reactions_list = [] - - for reaction in self.reactions: - reaction = {"reactant_coefficients": reaction.reactant_coefficients, - "reactant_types": reaction.reactant_types, - "product_types": reaction.product_types, - "product_coefficients": reaction.product_coefficients, - "reactant_types": reaction.reactant_types, - "gamma": reaction.gamma} - reactions_list.append(reaction) + property_keys = {"reactant_coefficients", "reactant_types", + "product_coefficients", "product_types", "gamma"} + reactions_list = [{key: getattr(reaction, key) for key in property_keys} + for reaction in self.reactions] return {"reactions": reactions_list, "kT": self.kT, "exclusion_range": self.exclusion_range,