From 7add1c62323aef028c028a3c2722d696be1cfb7d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 11 Mar 2021 15:24:35 +0100 Subject: [PATCH 01/23] core: Move SingleReaction constructor --- src/core/reaction_ensemble.cpp | 29 ++--------------------------- src/core/reaction_ensemble.hpp | 28 +++++++++++++++++++++++----- 2 files changed, 25 insertions(+), 32 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 3e783fed4b4..6ebcc80c55b 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -102,16 +102,8 @@ void ReactionAlgorithm::add_reaction( const std::vector &_reactant_coefficients, const std::vector &_product_types, const std::vector &_product_coefficients) { - SingleReaction new_reaction; - - new_reaction.gamma = gamma; - new_reaction.reactant_types = _reactant_types; - new_reaction.reactant_coefficients = _reactant_coefficients; - new_reaction.product_types = _product_types; - new_reaction.product_coefficients = _product_coefficients; - - new_reaction.nu_bar = calculate_nu_bar(new_reaction.reactant_coefficients, - new_reaction.product_coefficients); + SingleReaction new_reaction(gamma, _reactant_types, _reactant_coefficients, + _product_types, _product_coefficients); // make ESPResSo count the particle numbers which take part in the reactions for (int reactant_type : new_reaction.reactant_types) @@ -493,23 +485,6 @@ void ReactionAlgorithm::generic_oneway_reaction(int reaction_id) { on_end_reaction(accepted_state); } -/** - * Calculates the change in particle numbers for the given reaction - */ -int ReactionAlgorithm::calculate_nu_bar( - std::vector &reactant_coefficients, - std::vector &product_coefficients) { - // should only be used at when defining a new reaction - int nu_bar = 0; - for (int reactant_coefficient : reactant_coefficients) { - nu_bar -= reactant_coefficient; - } - for (int product_coefficient : product_coefficients) { - nu_bar += product_coefficient; - } - return nu_bar; -} - /** * 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. diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 35c1ad5818c..041593feb9f 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -27,8 +27,10 @@ #include #include +#include #include #include +#include #include #include #include @@ -38,6 +40,26 @@ namespace ReactionEnsemble { struct SingleReaction { + SingleReaction() = default; + SingleReaction(double gamma, std::vector const &reactant_types, + std::vector const &reactant_coefficients, + std::vector const &product_types, + std::vector const &product_coefficients) { + std::copy(reactant_types.begin(), reactant_types.end(), + std::back_inserter(this->reactant_types)); + std::copy(reactant_coefficients.begin(), reactant_coefficients.end(), + std::back_inserter(this->reactant_coefficients)); + std::copy(product_types.begin(), product_types.end(), + std::back_inserter(this->product_types)); + std::copy(product_coefficients.begin(), product_coefficients.end(), + std::back_inserter(this->product_coefficients)); + this->gamma = gamma; + nu_bar = std::accumulate(product_coefficients.begin(), + product_coefficients.end(), 0) - + std::accumulate(reactant_coefficients.begin(), + reactant_coefficients.end(), 0); + } + // strict input to the algorithm std::vector reactant_types; std::vector reactant_coefficients; @@ -45,7 +67,7 @@ struct SingleReaction { std::vector product_coefficients; double gamma = {}; // calculated values that are stored for performance reasons - int nu_bar = {}; + int nu_bar = {}; ///< change in particle numbers for the reaction Utils::Accumulator accumulator_exponentials = Utils::Accumulator(1); int tried_moves = 0; int accepted_moves = 0; @@ -213,10 +235,6 @@ class ReactionAlgorithm { std::map save_old_particle_numbers(int reaction_id); - int calculate_nu_bar( - std::vector &reactant_coefficients, - std::vector &product_coefficients); // should only be used when - // defining a new reaction void replace_particle(int p_id, int desired_type); int create_particle(int desired_type); void hide_particle(int p_id, int previous_type); From b62a88c725c8894193137dd86f175c52741c435e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 11 Mar 2021 15:34:10 +0100 Subject: [PATCH 02/23] core: ReactionEnsemble: make everything const --- src/core/reaction_ensemble.cpp | 52 ++++++++++----------- src/core/reaction_ensemble.hpp | 84 ++++++++++++++++++---------------- 2 files changed, 70 insertions(+), 66 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 6ebcc80c55b..59d1d0906bc 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -120,7 +120,7 @@ void ReactionAlgorithm::add_reaction( * Checks whether all necessary variables for the reaction ensemble have been * set. */ -void ReactionAlgorithm::check_reaction_ensemble() { +void ReactionAlgorithm::check_reaction_ensemble() const { if (reactions.empty()) { throw std::runtime_error("Reaction system not initialized"); } @@ -192,7 +192,7 @@ double factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(int Ni0, int nu_i) { /** * Checks whether all particles exist for the provided reaction. */ -bool ReactionAlgorithm::all_reactant_particles_exist(int reaction_id) { +bool ReactionAlgorithm::all_reactant_particles_exist(int reaction_id) const { bool enough_particles = true; for (int i = 0; i < reactions[reaction_id].reactant_types.size(); i++) { int current_number = @@ -222,7 +222,7 @@ void ReactionAlgorithm::append_particle_property_of_random_particle( *Performs a trial reaction move */ void ReactionAlgorithm::make_reaction_attempt( - SingleReaction ¤t_reaction, + SingleReaction const ¤t_reaction, std::vector &changed_particles_properties, std::vector &p_ids_created_particles, std::vector &hidden_particles_properties) { @@ -296,7 +296,7 @@ void ReactionAlgorithm::make_reaction_attempt( * when a reaction attempt is rejected. */ void ReactionAlgorithm::restore_properties( - std::vector &property_list, + std::vector const &property_list, const int number_of_saved_properties) { // this function restores all properties of all particles provided in the // property list, the format of the property list is (p_id,charge,type) @@ -318,10 +318,10 @@ void ReactionAlgorithm::restore_properties( * ensemble */ double ReactionEnsemble::calculate_acceptance_probability( - SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, - std::map &old_particle_numbers, int dummy_old_state_index, + SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, + std::map const &old_particle_numbers, int dummy_old_state_index, int dummy_new_state_index, - bool dummy_only_make_configuration_changing_move) { + bool dummy_only_make_configuration_changing_move) const { const double factorial_expr = calculate_factorial_expression(current_reaction, old_particle_numbers); @@ -818,10 +818,10 @@ void WangLandauReactionEnsemble::add_new_CV_potential_energy( * Returns the flattened index of the multidimensional Wang-Landau histogram */ int WangLandauReactionEnsemble::get_flattened_index_wang_landau( - std::vector ¤t_state, - std::vector &collective_variables_minimum_values, - std::vector &collective_variables_maximum_values, - std::vector &delta_collective_variables_values, + std::vector const ¤t_state, + std::vector const &collective_variables_minimum_values, + std::vector const &collective_variables_maximum_values, + std::vector const &delta_collective_variables_values, int nr_collective_variables) { int index = -10; // negative number is not allowed as index and therefore // indicates error @@ -917,7 +917,7 @@ int WangLandauReactionEnsemble:: * grid which starts at 0 */ double WangLandauReactionEnsemble::get_minimum_CV_value_on_delta_CV_spaced_grid( - double min_CV_value, double delta_CV) { + double min_CV_value, double delta_CV) const { // assume grid has it s origin at 0 double minimum_CV_value_on_delta_CV_spaced_grid = floor(min_CV_value / delta_CV) * delta_CV; @@ -950,7 +950,7 @@ double WangLandauReactionEnsemble::calculate_delta_degree_of_association( /** * Initializes the Wang-Landau histogram. */ -int WangLandauReactionEnsemble::get_num_needed_bins() { +int WangLandauReactionEnsemble::get_num_needed_bins() const { int needed_bins = 1; for (const auto ¤t_collective_variable : collective_variables) { needed_bins = needed_bins * (int((current_collective_variable->CV_maximum - @@ -1044,9 +1044,9 @@ int WangLandauReactionEnsemble::initialize_wang_landau() { * Modify Boltzmann factor according to Wang-Landau algorithm in @cite yan02b. */ double WangLandauReactionEnsemble::calculate_acceptance_probability( - SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, - std::map &old_particle_numbers, int old_state_index, - int new_state_index, bool only_make_configuration_changing_move) { + SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, + std::map const &old_particle_numbers, int old_state_index, + int new_state_index, bool only_make_configuration_changing_move) const { double beta = 1.0 / temperature; double bf; if (do_not_sample_reaction_partition_function || @@ -1156,7 +1156,7 @@ void WangLandauReactionEnsemble::update_wang_landau_potential_and_histogram( /** *Determines whether we can reduce the Wang-Landau parameter */ -bool WangLandauReactionEnsemble::can_refine_wang_landau_one_over_t() { +bool WangLandauReactionEnsemble::can_refine_wang_landau_one_over_t() const { double minimum_required_value = 0.80 * average_list_of_allowed_entries( histogram); // This is an additional constraint to sample @@ -1214,7 +1214,7 @@ void WangLandauReactionEnsemble::refine_wang_landau_parameter_one_over_t() { *Determine whether the desired number of refinements was achieved. */ bool WangLandauReactionEnsemble:: - achieved_desired_number_of_refinements_one_over_t() { + achieved_desired_number_of_refinements_one_over_t() const { if (wang_landau_parameter < final_wang_landau_parameter) { printf("Achieved desired number of refinements\n"); return true; @@ -1539,10 +1539,10 @@ int ConstantpHEnsemble::do_reaction(int reaction_steps) { * method. */ double ConstantpHEnsemble::calculate_acceptance_probability( - SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, - std::map &dummy_old_particle_numbers, int dummy_old_state_index, - int dummy_new_state_index, - bool dummy_only_make_configuration_changing_move) { + SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, + std::map const &dummy_old_particle_numbers, + int dummy_old_state_index, int dummy_new_state_index, + bool dummy_only_make_configuration_changing_move) const { auto const beta = 1.0 / temperature; 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 * @@ -1603,13 +1603,13 @@ WidomInsertion::measure_excess_chemical_potential(int reaction_id) { * See @cite smith94c. */ double -calculate_factorial_expression(SingleReaction ¤t_reaction, - std::map &old_particle_numbers) { +calculate_factorial_expression(SingleReaction const ¤t_reaction, + std::map const &old_particle_numbers) { double factorial_expr = 1.0; // 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[current_reaction.reactant_types[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 @@ -1618,7 +1618,7 @@ calculate_factorial_expression(SingleReaction ¤t_reaction, // 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[current_reaction.product_types[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 diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 041593feb9f..56daa2f8da4 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -71,10 +71,10 @@ struct SingleReaction { Utils::Accumulator accumulator_exponentials = Utils::Accumulator(1); int tried_moves = 0; int accepted_moves = 0; - double get_acceptance_rate() { + double get_acceptance_rate() const { return static_cast(accepted_moves) / static_cast(tried_moves); - }; + } }; struct StoredParticleProperty { @@ -87,9 +87,9 @@ struct CollectiveVariable { double CV_minimum = {}; double CV_maximum = {}; double delta_CV = {}; - virtual double determine_current_state() = 0; // use pure virtual, otherwise - // this will be used in vector - // of collective variables + // use pure virtual, otherwise this will be used in vector of collective + // variables + virtual double determine_current_state() const = 0; virtual ~CollectiveVariable() = default; }; @@ -97,7 +97,7 @@ class WangLandauReactionEnsemble; struct EnergyCollectiveVariable : public CollectiveVariable { std::string energy_boundaries_filename; - double determine_current_state() override { + double determine_current_state() const override { return calculate_current_potential_energy_of_system(); } void @@ -107,7 +107,7 @@ struct EnergyCollectiveVariable : public CollectiveVariable { struct DegreeOfAssociationCollectiveVariable : public CollectiveVariable { std::vector corresponding_acid_types; int associated_type; - double determine_current_state() override { + double determine_current_state() const override { return calculate_degree_of_association(); } @@ -117,7 +117,7 @@ struct DegreeOfAssociationCollectiveVariable : public CollectiveVariable { * is needed since you may use multiple degrees of association as collective * variable for the Wang-Landau algorithm. */ - double calculate_degree_of_association() { + double calculate_degree_of_association() const { int total_number_of_corresponding_acid = 0; for (int corresponding_acid_type : corresponding_acid_types) { int num_of_current_type = @@ -170,14 +170,14 @@ class ReactionAlgorithm { int m_accepted_configurational_MC_moves = 0; int m_tried_configurational_MC_moves = 0; - double get_acceptance_rate_configurational_moves() { + double get_acceptance_rate_configurational_moves() const { return static_cast(m_accepted_configurational_MC_moves) / static_cast(m_tried_configurational_MC_moves); } void set_cuboid_reaction_ensemble_volume(); virtual int do_reaction(int reaction_steps); - void check_reaction_ensemble(); + void check_reaction_ensemble() const; int delete_particle(int p_id); void add_reaction(double gamma, const std::vector &_reactant_types, @@ -206,15 +206,16 @@ class ReactionAlgorithm { virtual void on_mc_rejection_directly_after_entry(int &old_state_index){}; virtual void on_mc_accept(int &new_state_index){}; virtual void on_mc_reject(int &old_state_index){}; - virtual int on_mc_use_WL_get_new_state() { return -10; }; + virtual int on_mc_use_WL_get_new_state() { return -10; } void make_reaction_attempt( - SingleReaction ¤t_reaction, + SingleReaction const ¤t_reaction, std::vector &changed_particles_properties, std::vector &p_ids_created_particles, std::vector &hidden_particles_properties); - void restore_properties(std::vector &property_list, - int number_of_saved_properties); + void + restore_properties(std::vector const &property_list, + int number_of_saved_properties); /** * @brief draws a random integer from the uniform distribution in the range @@ -226,7 +227,7 @@ class ReactionAlgorithm { std::uniform_int_distribution uniform_int_dist(0, maxint - 1); return uniform_int_dist(m_generator); } - bool all_reactant_particles_exist(int reaction_id); + bool all_reactant_particles_exist(int reaction_id) const; private: std::mt19937 m_generator; @@ -243,11 +244,12 @@ class ReactionAlgorithm { int type, std::vector &list_of_particles); virtual double calculate_acceptance_probability( - SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, - std::map &old_particle_numbers, int old_state_index, - int new_state_index, bool only_make_configuration_changing_move) { + SingleReaction const ¤t_reaction, double E_pot_old, + double E_pot_new, std::map const &old_particle_numbers, + int old_state_index, int new_state_index, + bool only_make_configuration_changing_move) const { return -10; - }; + } void add_types_to_index(std::vector &type_list); Utils::Vector3d get_random_position_in_box(); @@ -270,10 +272,10 @@ class ReactionEnsemble : public ReactionAlgorithm { private: double calculate_acceptance_probability( - SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, - std::map &old_particle_numbers, int dummy_old_state_index, - int dummy_new_state_index, - bool dummy_only_make_configuration_changing_move) override; + SingleReaction const ¤t_reaction, double E_pot_old, + double E_pot_new, std::map const &old_particle_numbers, + int dummy_old_state_index, int dummy_new_state_index, + bool dummy_only_make_configuration_changing_move) const override; }; /** Wang-Landau reaction ensemble method */ @@ -323,9 +325,10 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { void on_attempted_reaction(int &new_state_index) override; void on_end_reaction(int &accepted_state) override; double calculate_acceptance_probability( - SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, - std::map &old_particle_numbers, int old_state_index, - int new_state_index, bool only_make_configuration_changing_move) override; + SingleReaction const ¤t_reaction, double E_pot_old, + double E_pot_new, std::map const &old_particle_numbers, + int old_state_index, int new_state_index, + bool only_make_configuration_changing_move) const override; void on_mc_rejection_directly_after_entry(int &old_state_index) override; void on_mc_accept(int &new_state_index) override; void on_mc_reject(int &old_state_index) override; @@ -351,30 +354,30 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { int collective_variable_index_energy_observable); // needed for energy int get_flattened_index_wang_landau( - std::vector ¤t_state, - std::vector &collective_variables_minimum_values, - std::vector &collective_variables_maximum_values, - std::vector &delta_collective_variables_values, + std::vector const ¤t_state, + std::vector const &collective_variables_minimum_values, + std::vector const &collective_variables_maximum_values, + std::vector const &delta_collective_variables_values, int nr_collective_variables); // collective variable int get_flattened_index_wang_landau_of_current_state(); void update_wang_landau_potential_and_histogram( int index_of_state_after_acceptance_or_rejection); int m_WL_tries = 0; - bool can_refine_wang_landau_one_over_t(); + bool can_refine_wang_landau_one_over_t() const; bool m_system_is_in_1_over_t_regime = false; - bool achieved_desired_number_of_refinements_one_over_t(); + bool achieved_desired_number_of_refinements_one_over_t() const; void refine_wang_landau_parameter_one_over_t(); int initialize_wang_landau(); // has to be called (at least) after the last // collective variable is added double calculate_delta_degree_of_association( DegreeOfAssociationCollectiveVariable ¤t_collective_variable); - int get_num_needed_bins(); + int get_num_needed_bins() const; void invalidate_bins(); void reset_histogram(); double get_minimum_CV_value_on_delta_CV_spaced_grid(double min_CV_value, - double delta_CV); + double delta_CV) const; }; /** @@ -398,10 +401,10 @@ class ConstantpHEnsemble : public ReactionAlgorithm { private: double calculate_acceptance_probability( - SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, - std::map &dummy_old_particle_numbers, int dummy_old_state_index, - int dummy_new_state_index, - bool dummy_only_make_configuration_changing_move) override; + SingleReaction const ¤t_reaction, double E_pot_old, + double E_pot_new, std::map const &dummy_old_particle_numbers, + int dummy_old_state_index, int dummy_new_state_index, + bool dummy_only_make_configuration_changing_move) const override; int get_random_valid_p_id(); }; @@ -416,8 +419,9 @@ class WidomInsertion : public ReactionAlgorithm { // utility functions // /////////////////////// -double calculate_factorial_expression(SingleReaction ¤t_reaction, - std::map &old_particle_numbers); +double +calculate_factorial_expression(SingleReaction const ¤t_reaction, + std::map const &old_particle_numbers); double factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(int Ni0, int nu_i); From 974e8e9f21e2d8db49a76d9ecfbdbed1400878b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 11 Mar 2021 15:57:01 +0100 Subject: [PATCH 03/23] tests: Unit test ReactionEnsemble classes --- src/core/unit_tests/CMakeLists.txt | 3 + .../reaction_ensemble_classes_test.cpp | 235 ++++++++++++++++++ .../reaction_ensemble_utils_test.cpp | 4 +- 3 files changed, 241 insertions(+), 1 deletion(-) create mode 100644 src/core/unit_tests/reaction_ensemble_classes_test.cpp diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index e5a8bb4e837..754601ac1da 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -48,3 +48,6 @@ unit_test(NAME random_test SRC random_test.cpp DEPENDS EspressoUtils Random123) unit_test(NAME BondList_test SRC BondList_test.cpp DEPENDS EspressoCore) unit_test(NAME reaction_ensemble_utils_test SRC reaction_ensemble_utils_test.cpp DEPENDS EspressoCore) +unit_test(NAME reaction_ensemble_classes_test SRC + reaction_ensemble_classes_test.cpp DEPENDS EspressoCore Boost::mpi + MPI::MPI_CXX) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp new file mode 100644 index 00000000000..c6853ac8de5 --- /dev/null +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -0,0 +1,235 @@ +/* + * Copyright (C) 2021 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 ReactionEnsemble classes. */ + +#define BOOST_TEST_NO_MAIN +#define BOOST_TEST_MODULE ReactionEnsemble classes test +#define BOOST_TEST_ALTERNATIVE_INIT_API +#define BOOST_TEST_DYN_LINK +#include + +#include "reaction_ensemble.hpp" + +#include "communication.hpp" +#include "particle_data.hpp" + +#include + +#include +#include +#include +#include + +/** Fixture to create particles during a test and remove them at the end. */ +struct ParticleFactory { + ParticleFactory() = default; + ~ParticleFactory() { + for (auto pid : particle_cache) { + remove_particle(pid); + } + } + void create_particle(int pid, int type) { + double pos[3] = {0., 0., 0.}; + place_particle(pid, pos); + set_particle_type(pid, type); + particle_cache.emplace_back(pid); + } + +private: + std::vector particle_cache; +}; + +BOOST_FIXTURE_TEST_CASE(particle_type_map_test, ParticleFactory) { + constexpr double tol = 100 * std::numeric_limits::epsilon(); + + // particle properties + int const type = 10; + int const pid = 1; + + // exception for untracked particle ids + BOOST_CHECK_THROW(number_of_particles_with_type(type), std::runtime_error); + + // exception for negative particle ids + BOOST_CHECK_THROW(init_type_map(-10), std::runtime_error); + + // check particle counting + init_type_map(type); + BOOST_CHECK_EQUAL(number_of_particles_with_type(type), 0); + create_particle(pid, type); + BOOST_CHECK_EQUAL(number_of_particles_with_type(type), 1); + + // exception for random index that exceeds the number of particles + BOOST_CHECK_THROW(get_random_p_id(type, 10), std::runtime_error); + + // check particle selection + BOOST_CHECK_EQUAL(get_random_p_id(type, 0), pid); +} + +BOOST_FIXTURE_TEST_CASE(DegreeOfAssociationCollectiveVariable_test, + ParticleFactory) { + using namespace ReactionEnsemble; + constexpr double tol = 100 * std::numeric_limits::epsilon(); + + // particle types + int const type_A = 0; + int const type_AH = 1; + int const type_AH2 = 2; + init_type_map(type_A); + init_type_map(type_AH); + init_type_map(type_AH2); + + // collective variable + DegreeOfAssociationCollectiveVariable doa_cv{}; + doa_cv.corresponding_acid_types = {type_A, type_AH, type_AH2}; + doa_cv.associated_type = type_A; + + // exception if no base is present + BOOST_CHECK_THROW(doa_cv.determine_current_state(), std::runtime_error); + + // add base + create_particle(0, type_A); + BOOST_CHECK_CLOSE(doa_cv.determine_current_state(), 1.0, tol); + + // add acid + create_particle(1, type_AH); + BOOST_CHECK_CLOSE(doa_cv.determine_current_state(), 0.5, tol); + + // add acid + create_particle(2, type_AH); + create_particle(3, type_AH2); + BOOST_CHECK_CLOSE(doa_cv.determine_current_state(), 0.25, tol); +} + +BOOST_AUTO_TEST_CASE(SingleReaction_test) { + using namespace ReactionEnsemble; + constexpr double tol = 100 * std::numeric_limits::epsilon(); + + // check derived parameter + SingleReaction reaction(2., {0}, {1}, {1, 2}, {3, 4}); + BOOST_CHECK_EQUAL(reaction.nu_bar, 6); + + // check acceptance rate + for (int tried_moves = 1; tried_moves < 5; ++tried_moves) { + for (int accepted_moves = 0; accepted_moves < 5; ++accepted_moves) { + reaction.tried_moves = tried_moves; + reaction.accepted_moves = accepted_moves; + auto const ref_rate = static_cast(accepted_moves) / + static_cast(tried_moves); + BOOST_CHECK_CLOSE(reaction.get_acceptance_rate(), ref_rate, tol); + } + } +} + +BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { + using namespace ReactionEnsemble; + constexpr double tol = 100 * std::numeric_limits::epsilon(); + + // check acceptance rate + ReactionAlgorithm r_algo(42); + 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; + r_algo.m_accepted_configurational_MC_moves = accepted_moves; + auto const ref_rate = static_cast(accepted_moves) / + static_cast(tried_moves); + BOOST_CHECK_CLOSE(r_algo.get_acceptance_rate_configurational_moves(), + ref_rate, tol); + } + } + + // exception if no reaction was added + BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), std::runtime_error); + + // check reaction addition + { + SingleReaction const ref(2., {0}, {1}, {1, 2}, {3, 4}); + r_algo.add_reaction(ref.gamma, ref.reactant_types, + ref.reactant_coefficients, ref.product_types, + ref.product_coefficients); + BOOST_CHECK_EQUAL(r_algo.reactions.size(), 1ul); + auto const &reaction = r_algo.reactions[0]; + BOOST_TEST(reaction.reactant_types == ref.reactant_types, + boost::test_tools::per_element()); + BOOST_TEST(reaction.reactant_coefficients == ref.reactant_coefficients, + boost::test_tools::per_element()); + BOOST_TEST(reaction.product_types == ref.product_types, + boost::test_tools::per_element()); + BOOST_TEST(reaction.product_coefficients == ref.product_coefficients, + boost::test_tools::per_element()); + BOOST_CHECK_EQUAL(reaction.gamma, ref.gamma); + } + + // exception if temperature is negative + BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), std::runtime_error); + + // set temperature + r_algo.temperature = 10; + +#ifdef ELECTROSTATICS + // exception if reactant types have no charge information + BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), 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_ensemble(), 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_ensemble(); + + // check reaction removal + { + r_algo.add_reaction(5, {1}, {1}, {2}, {1}); + BOOST_CHECK_EQUAL(r_algo.reactions.size(), 2ul); + BOOST_CHECK_EQUAL(r_algo.reactions[1].gamma, 5); + r_algo.delete_reaction(1); + BOOST_CHECK_EQUAL(r_algo.reactions.size(), 1ul); + BOOST_CHECK_EQUAL(r_algo.reactions[0].gamma, 2); + r_algo.delete_reaction(0); + BOOST_CHECK_EQUAL(r_algo.reactions.size(), 0ul); + } + + // exception if deleting a non-existent particle + BOOST_CHECK_THROW(r_algo.delete_particle(5), std::runtime_error); +} + +BOOST_AUTO_TEST_CASE(WidomInsertion_test) { + using namespace ReactionEnsemble; + constexpr double tol = 100 * std::numeric_limits::epsilon(); + + // check acceptance rate + WidomInsertion r_algo(42); + + SingleReaction const ref(2., {0}, {1}, {1, 2}, {3, 4}); + r_algo.add_reaction(ref.gamma, ref.reactant_types, ref.reactant_coefficients, + ref.product_types, ref.product_coefficients); + + // exception if not enough particles + BOOST_CHECK_THROW(r_algo.measure_excess_chemical_potential(0), + std::runtime_error); +} + +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/unit_tests/reaction_ensemble_utils_test.cpp b/src/core/unit_tests/reaction_ensemble_utils_test.cpp index 3f78f7691ef..b8e47dfd219 100644 --- a/src/core/unit_tests/reaction_ensemble_utils_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_utils_test.cpp @@ -17,7 +17,9 @@ * along with this program. If not, see . */ -#define BOOST_TEST_MODULE ReactionEnsemble test +/* Unit tests for the ReactionEnsemble utility functions. */ + +#define BOOST_TEST_MODULE ReactionEnsemble utility functions test #define BOOST_TEST_DYN_LINK #include From 3191b974929835dcc8bc24846cfd18fb23118d5f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 12 Mar 2021 12:53:38 +0100 Subject: [PATCH 04/23] tests: Check factorial expression --- .../reaction_ensemble_classes_test.cpp | 14 ++++++++++++ .../reaction_ensemble_utils_test.cpp | 22 ++++++++++--------- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index c6853ac8de5..8cc9b3519e3 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -135,6 +135,20 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) { BOOST_CHECK_CLOSE(reaction.get_acceptance_rate(), ref_rate, tol); } } + + // check factorial expression + constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // i adduct #0, j product #1, k product #2 + auto const p_numbers = std::map{{0, i}, {1, j}, {2, 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); + } + } + } } BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { diff --git a/src/core/unit_tests/reaction_ensemble_utils_test.cpp b/src/core/unit_tests/reaction_ensemble_utils_test.cpp index b8e47dfd219..a234bfab15c 100644 --- a/src/core/unit_tests/reaction_ensemble_utils_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_utils_test.cpp @@ -25,6 +25,7 @@ #include "reaction_ensemble.hpp" +#include #include #include #include @@ -59,14 +60,15 @@ BOOST_AUTO_TEST_CASE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i_test) { using namespace ReactionEnsemble; constexpr double tol = 100 * std::numeric_limits::epsilon(); - BOOST_CHECK_CLOSE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(2, 0), 1.0, - tol); - BOOST_CHECK_CLOSE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(2, 0), 1.0, - tol); - BOOST_CHECK_CLOSE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(1, 2), - 1.0 / 6.0, tol); - BOOST_CHECK_CLOSE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(1, -2), - 0.0, tol); - BOOST_CHECK_CLOSE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(2, -2), - 2.0, tol); + auto const reaction_ensemble_combinations = [](int N, int nu) { + return (N + nu < 0) ? 0. : std::tgamma(N + 1) / std::tgamma(N + nu + 1); + }; + + for (int N0 = 0; N0 < 6; ++N0) { + 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, 5 * tol); + } + } } From acf5786d2b5b1dfde6858523d3954af9df47e260 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 12 Mar 2021 13:08:40 +0100 Subject: [PATCH 05/23] core: Make acceptance probability methods protected This way the methods can be unit tested. --- src/core/reaction_ensemble.hpp | 35 +++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 56daa2f8da4..345f9b3ffea 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -229,6 +229,15 @@ class ReactionAlgorithm { } bool all_reactant_particles_exist(int reaction_id) const; +protected: + virtual double calculate_acceptance_probability( + SingleReaction const ¤t_reaction, double E_pot_old, + double E_pot_new, std::map const &old_particle_numbers, + int old_state_index, int new_state_index, + bool only_make_configuration_changing_move) const { + return -10; + } + private: std::mt19937 m_generator; std::normal_distribution m_normal_distribution; @@ -243,14 +252,6 @@ class ReactionAlgorithm { void append_particle_property_of_random_particle( int type, std::vector &list_of_particles); - virtual double calculate_acceptance_probability( - SingleReaction const ¤t_reaction, double E_pot_old, - double E_pot_new, std::map const &old_particle_numbers, - int old_state_index, int new_state_index, - bool only_make_configuration_changing_move) const { - return -10; - } - void add_types_to_index(std::vector &type_list); Utils::Vector3d get_random_position_in_box(); }; @@ -270,7 +271,7 @@ class ReactionEnsemble : public ReactionAlgorithm { public: ReactionEnsemble(int seed) : ReactionAlgorithm(seed) {} -private: +protected: double calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, std::map const &old_particle_numbers, @@ -318,17 +319,19 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { void write_wang_landau_results_to_file( const std::string &full_path_to_output_filename); +protected: + double calculate_acceptance_probability( + SingleReaction const ¤t_reaction, double E_pot_old, + double E_pot_new, std::map const &old_particle_numbers, + int old_state_index, int new_state_index, + bool only_make_configuration_changing_move) const override; + private: void on_reaction_entry(int &old_state_index) override; void on_reaction_rejection_directly_after_entry(int &old_state_index) override; void on_attempted_reaction(int &new_state_index) override; void on_end_reaction(int &accepted_state) override; - double calculate_acceptance_probability( - SingleReaction const ¤t_reaction, double E_pot_old, - double E_pot_new, std::map const &old_particle_numbers, - int old_state_index, int new_state_index, - bool only_make_configuration_changing_move) const override; void on_mc_rejection_directly_after_entry(int &old_state_index) override; void on_mc_accept(int &new_state_index) override; void on_mc_reject(int &old_state_index) override; @@ -399,12 +402,14 @@ class ConstantpHEnsemble : public ReactionAlgorithm { double m_constant_pH = -10; int do_reaction(int reaction_steps) override; -private: +protected: double calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, std::map const &dummy_old_particle_numbers, int dummy_old_state_index, int dummy_new_state_index, bool dummy_only_make_configuration_changing_move) const override; + +private: int get_random_valid_p_id(); }; From 3d22a1d8e872b94ad1b643dc15200d2635d1d24f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 12 Mar 2021 13:49:50 +0100 Subject: [PATCH 06/23] tests: Check acceptance probability methods --- .../reaction_ensemble_classes_test.cpp | 142 +++++++++++++++++- 1 file changed, 141 insertions(+), 1 deletion(-) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index 8cc9b3519e3..3160122edb4 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -30,8 +30,11 @@ #include "communication.hpp" #include "particle_data.hpp" +#include + #include +#include #include #include #include @@ -153,10 +156,15 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) { BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { using namespace ReactionEnsemble; + class ReactionAlgorithmTest : public ReactionAlgorithm { + public: + using ReactionAlgorithm::calculate_acceptance_probability; + using ReactionAlgorithm::ReactionAlgorithm; + }; constexpr double tol = 100 * std::numeric_limits::epsilon(); // check acceptance rate - ReactionAlgorithm r_algo(42); + ReactionAlgorithmTest r_algo(42); 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; @@ -190,6 +198,14 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { BOOST_CHECK_EQUAL(reaction.gamma, ref.gamma); } + // check acceptance probability + { + SingleReaction const reaction(2., {0}, {1}, {1, 2}, {3, 4}); + double probability = r_algo.calculate_acceptance_probability( + reaction, -1., -1., {{1, 2}}, -1, -1, false); + BOOST_CHECK_EQUAL(probability, -10.); + } + // exception if temperature is negative BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), std::runtime_error); @@ -225,6 +241,130 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { BOOST_CHECK_THROW(r_algo.delete_particle(5), std::runtime_error); } +BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { + class ReactionEnsembleTest : public ReactionEnsemble::ReactionEnsemble { + public: + using ReactionEnsemble::ReactionEnsemble::calculate_acceptance_probability; + using ReactionEnsemble::ReactionEnsemble::ReactionEnsemble; + }; + using namespace ReactionEnsemble; + constexpr double tol = 100 * std::numeric_limits::epsilon(); + + ReactionEnsembleTest r_algo(42); + r_algo.volume = 10.; + r_algo.temperature = 20.; + + // exception if no reaction was added + BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), std::runtime_error); + + // check acceptance probability + constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i; + SingleReaction const reaction(2., {0}, {1}, {1, 2}, {3, 4}); + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // i adduct #0, j product #1, k product #2 + auto const p_numbers = std::map{{0, i}, {1, j}, {2, k}}; + auto const energy = static_cast(i + 1); + auto const f_expr = calculate_factorial_expression(reaction, p_numbers); + auto const ref = 2e6 * f_expr * std::exp(energy / 20.); + auto const val = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, -1, false); + BOOST_CHECK_CLOSE(val, ref, 5 * tol); + } + } + } +} + +BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { + using namespace ReactionEnsemble; + class WangLandauReactionEnsembleTest : public WangLandauReactionEnsemble { + public: + using WangLandauReactionEnsemble::calculate_acceptance_probability; + using WangLandauReactionEnsemble::WangLandauReactionEnsemble; + }; + constexpr double tol = 100 * std::numeric_limits::epsilon(); + + WangLandauReactionEnsembleTest r_algo(42); + r_algo.volume = 10.; + r_algo.temperature = 20.; + + // exception if no reaction was added + BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), std::runtime_error); + + // check acceptance probability + constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i; + SingleReaction const reaction(2., {0}, {1}, {1, 2}, {3, 4}); + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // i adduct #0, j product #1, k product #2 + auto const p_numbers = std::map{{0, i}, {1, j}, {2, k}}; + auto const energy = static_cast(i + 1); + auto const f_expr = calculate_factorial_expression(reaction, p_numbers); + auto const prefactor = 2e6 * f_expr; + auto const boltzmann = std::exp(energy / 20.); + double val; + r_algo.do_energy_reweighting = false; + val = r_algo.calculate_acceptance_probability(reaction, energy, 0., + p_numbers, -1, -1, false); + BOOST_CHECK_CLOSE(val, prefactor * boltzmann, 5 * tol); + val = r_algo.calculate_acceptance_probability(reaction, energy, 0., + p_numbers, -1, -1, true); + BOOST_CHECK_CLOSE(val, boltzmann, 5 * tol); + r_algo.do_energy_reweighting = true; + val = r_algo.calculate_acceptance_probability(reaction, energy, 0., + p_numbers, -1, -1, false); + BOOST_CHECK_CLOSE(val, prefactor, 5 * tol); + val = r_algo.calculate_acceptance_probability(reaction, energy, 0., + p_numbers, -1, -1, true); + BOOST_CHECK_CLOSE(val, 1., 5 * tol); + val = r_algo.calculate_acceptance_probability(reaction, energy, 0., + p_numbers, -1, 0, true); + BOOST_CHECK_CLOSE(val, 10., 5 * tol); + val = r_algo.calculate_acceptance_probability(reaction, energy, 0., + p_numbers, 0, -1, true); + BOOST_CHECK_CLOSE(val, -10., 5 * tol); + } + } + } +} + +BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { + using namespace ReactionEnsemble; + 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); + r_algo.temperature = 20.; + r_algo.m_constant_pH = 1.; + + // exception if no reaction was added + BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), std::runtime_error); + + // check acceptance probability + constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i; + SingleReaction const reaction(2., {0}, {1}, {1, 2}, {1, 1}); + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // i adduct #0, j product #1, k product #2 + auto const p_numbers = std::map{{0, i}, {1, j}, {2, k}}; + auto const energy = static_cast(i + 1); + auto const ref = + std::exp(energy / 20. + std::log(10.) * (1. + std::log10(2.))); + auto const val = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, -1, false); + BOOST_CHECK_CLOSE(val, ref, 5 * tol); + } + } + } +} + BOOST_AUTO_TEST_CASE(WidomInsertion_test) { using namespace ReactionEnsemble; constexpr double tol = 100 * std::numeric_limits::epsilon(); From b032021816286e6ae11729f731809f39d33eacd4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 17 Mar 2021 11:12:38 +0100 Subject: [PATCH 07/23] core: Remove prototype without implementation --- src/core/reaction_ensemble.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 345f9b3ffea..aadb5070756 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -252,7 +252,6 @@ class ReactionAlgorithm { void append_particle_property_of_random_particle( int type, std::vector &list_of_particles); - void add_types_to_index(std::vector &type_list); Utils::Vector3d get_random_position_in_box(); }; From 62cdcb5cd54ea8838ac233699de2b126b84e7579 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 17 Mar 2021 12:06:29 +0100 Subject: [PATCH 08/23] core: Add missing includes --- src/core/reaction_ensemble.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 59d1d0906bc..553dcad8e9b 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -23,6 +23,7 @@ #include "energy.hpp" #include "grid.hpp" #include "partCfg_global.hpp" +#include "particle_data.hpp" #include "statistics.hpp" #include @@ -30,12 +31,18 @@ #include #include +#include #include #include #include #include #include #include +#include +#include +#include +#include +#include #include namespace ReactionEnsemble { From 959b42683c53a138b83b50d54f6cf1d2a9881a09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 17 Mar 2021 12:59:29 +0100 Subject: [PATCH 09/23] core: Fix -Wconversion --- src/core/reaction_ensemble.cpp | 159 +++++++++++++++++---------------- src/core/reaction_ensemble.hpp | 6 +- 2 files changed, 83 insertions(+), 82 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 553dcad8e9b..95fd85826d1 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -95,7 +95,7 @@ void EnergyCollectiveVariable::load_CV_boundaries( */ int ReactionAlgorithm::do_reaction(int reaction_steps) { for (int i = 0; i < reaction_steps; i++) { - int reaction_id = i_random(reactions.size()); + int reaction_id = i_random(static_cast(reactions.size())); generic_oneway_reaction(reaction_id); } return 0; @@ -274,8 +274,8 @@ void ReactionAlgorithm::make_reaction_attempt( } } // create or hide particles of types with noncorresponding replacement types - for (int i = std::min(current_reaction.product_types.size(), - current_reaction.reactant_types.size()); + for (auto i = std::min(current_reaction.product_types.size(), + current_reaction.reactant_types.size()); i < std::max(current_reaction.product_types.size(), current_reaction.reactant_types.size()); i++) { @@ -462,13 +462,11 @@ void ReactionAlgorithm::generic_oneway_reaction(int reaction_id) { std::vector to_be_deleted_hidden_types( len_hidden_particles_properties); for (int i = 0; i < len_hidden_particles_properties; i++) { - auto p_id = static_cast(hidden_particles_properties[i].p_id); + auto p_id = hidden_particles_properties[i].p_id; to_be_deleted_hidden_ids[i] = p_id; to_be_deleted_hidden_types[i] = hidden_particles_properties[i].type; - set_particle_type(p_id, - hidden_particles_properties[i] - .type); // change back type otherwise the - // bookkeeping algorithm is not working + // change back type otherwise the bookkeeping algorithm is not working + set_particle_type(p_id, hidden_particles_properties[i].type); } for (int i = 0; i < len_hidden_particles_properties; i++) { @@ -830,15 +828,10 @@ int WangLandauReactionEnsemble::get_flattened_index_wang_landau( std::vector const &collective_variables_maximum_values, std::vector const &delta_collective_variables_values, int nr_collective_variables) { - int index = -10; // negative number is not allowed as index and therefore - // indicates error - std::vector individual_indices(nr_collective_variables); // pre result - // individual_indices.resize(nr_collective_variables,-1); //initialize - // individual_indices to -1 - - // check for the current state to be an allowed state in the [range - // collective_variables_minimum_values:collective_variables_maximum_values], - // else return a negative index + + // check for the current state to be an allowed state in the range + // [collective_variables_minimum_values:collective_variables_maximum_values], + // else return a negative index (to indicate error) for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { if (current_state[CV_i] > collective_variables_maximum_values[CV_i] + @@ -850,26 +843,28 @@ int WangLandauReactionEnsemble::get_flattened_index_wang_landau( } } + std::vector individual_indices(nr_collective_variables); for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - if (CV_i == collective_variables.size() - 1 && - do_energy_reweighting) // for energy collective variable (simple - // truncating conversion desired) - individual_indices[CV_i] = static_cast( - (current_state[CV_i] - collective_variables_minimum_values[CV_i]) / - delta_collective_variables_values[CV_i]); - else // for degree of association collective variables (rounding conversion - // desired) - individual_indices[CV_i] = std::lround( - (current_state[CV_i] - collective_variables_minimum_values[CV_i]) / - delta_collective_variables_values[CV_i]); - if (individual_indices[CV_i] < 0 or - individual_indices[CV_i] >= - nr_subindices_of_collective_variable[CV_i]) { // sanity check + auto const value = + (current_state[CV_i] - collective_variables_minimum_values[CV_i]) / + delta_collective_variables_values[CV_i]; + int rounded_value; + if (CV_i == collective_variables.size() - 1 && do_energy_reweighting) { + // for energy collective variable (simple truncating conversion desired) + rounded_value = static_cast(value); + } else { + // for degree of association collective variables (rounding conversion) + rounded_value = static_cast(std::floor(value)); + } + if (rounded_value < 0 or + rounded_value >= nr_subindices_of_collective_variable[CV_i]) { return -10; } + individual_indices[CV_i] = rounded_value; } // get flattened index from individual_indices - index = 0; // this is already part of the algorithm to find the correct index + int index = 0; + // this is already part of the algorithm to find the correct index for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { int factor = 1; for (int j = CV_i + 1; j < nr_collective_variables; j++) { @@ -886,7 +881,8 @@ int WangLandauReactionEnsemble::get_flattened_index_wang_landau( */ int WangLandauReactionEnsemble:: get_flattened_index_wang_landau_of_current_state() { - int nr_collective_variables = collective_variables.size(); + auto const nr_collective_variables = + static_cast(collective_variables.size()); // get current state std::vector current_state(nr_collective_variables); for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) @@ -976,28 +972,28 @@ void WangLandauReactionEnsemble::invalidate_bins() { // prohibit them int empty_bins_in_memory = 0; - for (int flattened_index = 0; flattened_index < wang_landau_potential.size(); - flattened_index++) { - std::vector unraveled_index( + for (std::size_t flattened_index = 0; + flattened_index < wang_landau_potential.size(); flattened_index++) { + std::vector unraveled_index( nr_subindices_of_collective_variable.size()); Utils::unravel_index(nr_subindices_of_collective_variable.begin(), nr_subindices_of_collective_variable.end(), unraveled_index.begin(), flattened_index); // use unraveled index int EnergyCollectiveVariable_index = 0; - if (collective_variables.size() > 1) + if (collective_variables.size() > 1) { + // assume the energy collective variable to be the last added + // collective variable EnergyCollectiveVariable_index = - static_cast(collective_variables.size()) - - 1; // assume the energy collective - // variable to be the last added - // collective variable + static_cast(collective_variables.size()) - 1; + } double current_energy = - unraveled_index[EnergyCollectiveVariable_index] * + static_cast(unraveled_index[EnergyCollectiveVariable_index]) * collective_variables[EnergyCollectiveVariable_index]->delta_CV + collective_variables[EnergyCollectiveVariable_index]->CV_minimum; int flat_index_without_energy_CV = get_flattened_index_wang_landau_without_energy_collective_variable( - flattened_index, EnergyCollectiveVariable_index); + static_cast(flattened_index), EnergyCollectiveVariable_index); std::shared_ptr energy_CV = collective_variables[EnergyCollectiveVariable_index]; if (current_energy > @@ -1118,7 +1114,7 @@ double WangLandauReactionEnsemble::calculate_acceptance_probability( int WangLandauReactionEnsemble::do_reaction(int reaction_steps) { m_WL_tries += reaction_steps; for (int step = 0; step < reaction_steps; step++) { - int reaction_id = i_random(reactions.size()); + int reaction_id = i_random(static_cast(reactions.size())); generic_oneway_reaction(reaction_id); if (can_refine_wang_landau_one_over_t() && m_WL_tries % 10000 == 0) { // check for convergence @@ -1202,8 +1198,8 @@ void WangLandauReactionEnsemble::reset_histogram() { *Refine the Wang-Landau parameter using the 1/t rule. */ void WangLandauReactionEnsemble::refine_wang_landau_parameter_one_over_t() { - double monte_carlo_time = - static_cast(monte_carlo_trial_moves) / used_bins; + double monte_carlo_time = static_cast(monte_carlo_trial_moves) / + static_cast(used_bins); if (wang_landau_parameter / 2.0 <= 1.0 / monte_carlo_time || m_system_is_in_1_over_t_regime) { wang_landau_parameter = 1.0 / monte_carlo_time; @@ -1240,27 +1236,28 @@ void WangLandauReactionEnsemble::write_wang_landau_results_to_file( if (pFile == nullptr) { throw std::runtime_error("ERROR: Wang-Landau file could not be written\n"); } - for (int flattened_index = 0; flattened_index < wang_landau_potential.size(); - flattened_index++) { + for (std::size_t flattened_index = 0; + flattened_index < wang_landau_potential.size(); flattened_index++) { // unravel index if (std::abs(wang_landau_potential[flattened_index] - double_fill_value) > - 1) { // only output data if they are not equal to - // double_fill_value. This if ensures - // that for the energy observable not allowed energies (energies - // in the interval [global_E_min, global_E_max]) in the - // multidimensional Wang-Landau potential are printed out, since - // the range [E_min(nbar), E_max(nbar)] for each nbar may be a - // different one - std::vector unraveled_index( + 1) { + // only output data if they are not equal to double_fill_value. + // This ensures that for the energy observable not allowed energies + // (energies in the interval [global_E_min, global_E_max]) in the + // multidimensional Wang-Landau potential are printed out, since + // the range [E_min(nbar), E_max(nbar)] for each nbar may be a + // different one + std::vector unraveled_index( nr_subindices_of_collective_variable.size()); Utils::unravel_index(nr_subindices_of_collective_variable.begin(), nr_subindices_of_collective_variable.end(), unraveled_index.begin(), flattened_index); // use unraveled index - for (int i = 0; i < collective_variables.size(); i++) { - fprintf(pFile, "%f ", - unraveled_index[i] * collective_variables[i]->delta_CV + - collective_variables[i]->CV_minimum); + for (std::size_t i = 0; i < collective_variables.size(); i++) { + auto const value = static_cast(unraveled_index[i]) * + collective_variables[i]->delta_CV + + collective_variables[i]->CV_minimum; + fprintf(pFile, "%f ", value); } fprintf(pFile, "%f \n", wang_landau_potential[flattened_index]); } @@ -1314,19 +1311,20 @@ void WangLandauReactionEnsemble::write_out_preliminary_energy_run_results( } fprintf(pFile, "#nbar E_min E_max\n"); - for (int flattened_index = 0; flattened_index < wang_landau_potential.size(); - flattened_index++) { + for (std::size_t flattened_index = 0; + flattened_index < wang_landau_potential.size(); flattened_index++) { // unravel index - std::vector unraveled_index( + std::vector unraveled_index( nr_subindices_of_collective_variable.size()); Utils::unravel_index(nr_subindices_of_collective_variable.begin(), nr_subindices_of_collective_variable.end(), unraveled_index.begin(), flattened_index); // use unraveled index - for (int i = 0; i < collective_variables.size(); i++) { - fprintf(pFile, "%f ", - unraveled_index[i] * collective_variables[i]->delta_CV + - collective_variables[i]->CV_minimum); + for (std::size_t i = 0; i < collective_variables.size(); i++) { + auto const value = static_cast(unraveled_index[i]) * + collective_variables[i]->delta_CV + + collective_variables[i]->CV_minimum; + fprintf(pFile, "%f ", value); } fprintf(pFile, "%f %f \n", minimum_energies_at_flat_index[flattened_index], maximum_energies_at_flat_index[flattened_index]); @@ -1343,18 +1341,20 @@ int WangLandauReactionEnsemble:: get_flattened_index_wang_landau_without_energy_collective_variable( int flattened_index_with_EnergyCollectiveVariable, int CV_index_energy_observable) { - std::vector unraveled_index(nr_subindices_of_collective_variable.size()); - Utils::unravel_index(nr_subindices_of_collective_variable.begin(), - nr_subindices_of_collective_variable.end(), - unraveled_index.begin(), - flattened_index_with_EnergyCollectiveVariable); - // use unraveled index - const int nr_collective_variables = - static_cast(collective_variables.size()) - - 1; // forget the last collective variable (the energy collective variable) + std::vector unraveled_index( + nr_subindices_of_collective_variable.size()); + Utils::unravel_index( + nr_subindices_of_collective_variable.begin(), + nr_subindices_of_collective_variable.end(), unraveled_index.begin(), + static_cast(flattened_index_with_EnergyCollectiveVariable)); + // use unraveled index, but skip last collective variable + // (the energy collective variable) + auto const nr_collective_variables = + static_cast(collective_variables.size()) - 1; std::vector current_state(nr_collective_variables); for (int i = 0; i < nr_collective_variables; i++) { - current_state[i] = unraveled_index[i] * collective_variables[i]->delta_CV + + current_state[i] = static_cast(unraveled_index[i]) * + collective_variables[i]->delta_CV + collective_variables[i]->CV_minimum; } @@ -1534,8 +1534,9 @@ int ConstantpHEnsemble::do_reaction(int reaction_steps) { } } // randomly select a reaction to be performed - int reaction_id = list_of_reaction_ids_with_given_reactant_type[i_random( - list_of_reaction_ids_with_given_reactant_type.size())]; + int reaction_id = + list_of_reaction_ids_with_given_reactant_type[i_random(static_cast( + list_of_reaction_ids_with_given_reactant_type.size()))]; generic_oneway_reaction(reaction_id); } return 0; diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index aadb5070756..9aa13276757 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -132,8 +132,7 @@ struct DegreeOfAssociationCollectiveVariable : public CollectiveVariable { int num_of_associated_acid = number_of_particles_with_type(associated_type); double degree_of_association = static_cast(num_of_associated_acid) / - total_number_of_corresponding_acid; // cast to double because otherwise - // any fractional part is lost + static_cast(total_number_of_corresponding_acid); return degree_of_association; } }; @@ -451,7 +450,8 @@ double average_list_of_allowed_entries(const std::vector &rng) { } } if (counter_allowed_entries) { - return static_cast(result) / counter_allowed_entries; + return static_cast(result) / + static_cast(counter_allowed_entries); } return 0.0; } From b32e086c446779dd2e80f928efa91679f590c98f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 17 Mar 2021 16:48:47 +0100 Subject: [PATCH 10/23] tests: Check Wang-Landau CV and histogram --- .../reaction_ensemble_classes_test.cpp | 71 +++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index 3160122edb4..a115bbd3fe4 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -33,11 +33,19 @@ #include #include +#include +#include +#include +#include +#include #include +#include +#include #include #include #include +#include #include /** Fixture to create particles during a test and remove them at the end. */ @@ -328,6 +336,69 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { } } } + + // check collective variables + { + using namespace boost::adaptors; + // setup input + auto const delta_CV = 0.5; + auto const filename_input = std::string("wl_input.dat"); + auto const filename_output = std::string("wl_output.dat"); + std::ofstream wl_file; + wl_file.open(filename_input); + wl_file << "# header 1.0 0.5\n"; + wl_file << "0 1.0 4.0\n"; + wl_file << "1\t2.0\t5.0\n"; + wl_file.close(); + r_algo.do_energy_reweighting = false; + // add collective variable + r_algo.add_new_CV_potential_energy(filename_input, delta_CV); + BOOST_REQUIRE_EQUAL(r_algo.collective_variables.size(), 1ul); + BOOST_REQUIRE_EQUAL(r_algo.min_boundaries_energies.size(), 2ul); + BOOST_REQUIRE_EQUAL(r_algo.max_boundaries_energies.size(), 2ul); + BOOST_CHECK_EQUAL(r_algo.do_energy_reweighting, true); + BOOST_CHECK_EQUAL(r_algo.collective_variables[0]->delta_CV, delta_CV); + BOOST_CHECK_EQUAL(r_algo.collective_variables[0]->CV_minimum, 1.); + BOOST_CHECK_EQUAL(r_algo.collective_variables[0]->CV_maximum, 5.); + BOOST_CHECK_CLOSE(r_algo.min_boundaries_energies[0], 1., tol); + BOOST_CHECK_CLOSE(r_algo.min_boundaries_energies[1], 2., tol); + BOOST_CHECK_CLOSE(r_algo.max_boundaries_energies[0], 4., tol); + BOOST_CHECK_CLOSE(r_algo.max_boundaries_energies[1], 5., tol); + // check Wang-Landau histogram + { + r_algo.write_wang_landau_results_to_file(filename_output); + std::ifstream outfile; + outfile.open(filename_output); + BOOST_TEST_REQUIRE(!outfile.fail(), "output file not generated"); + std::istream_iterator start(outfile), end; + std::vector flat_array(start, end); + outfile.close(); + BOOST_REQUIRE_EQUAL(flat_array.size(), 14ul); + std::vector bin_edges; + std::vector histogram; + boost::range::push_back(bin_edges, flat_array | strided(2)); + boost::range::push_back( + histogram, flat_array | sliced(1, flat_array.size()) | strided(2)); + std::vector bin_edges_ref(7u, 0.); + std::vector histogram_ref(7u, 0.); + std::generate(bin_edges_ref.begin(), bin_edges_ref.end(), + [n = .5]() mutable { return n += .5; }); + BOOST_TEST(bin_edges == bin_edges_ref, boost::test_tools::per_element()); + BOOST_TEST(histogram == histogram_ref, boost::test_tools::per_element()); + } + // add collective variable (with different energies in WL object) + r_algo.min_boundaries_energies.emplace_back(-1.); + r_algo.max_boundaries_energies.emplace_back(10.); + r_algo.add_new_CV_potential_energy(filename_input, delta_CV); + BOOST_REQUIRE_EQUAL(r_algo.collective_variables.size(), 2ul); + BOOST_REQUIRE_EQUAL(r_algo.min_boundaries_energies.size(), 5ul); + BOOST_REQUIRE_EQUAL(r_algo.max_boundaries_energies.size(), 5ul); + BOOST_CHECK_EQUAL(r_algo.collective_variables[1]->CV_minimum, -1.); + BOOST_CHECK_EQUAL(r_algo.collective_variables[1]->CV_maximum, 10.); + // exception if file doesn't exist + BOOST_CHECK_THROW(r_algo.add_new_CV_potential_energy("unknown", 0.), + std::runtime_error); + } } BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { From 989395b839319beb2293816c3b00e8c1a76cf8b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 17 Mar 2021 17:58:26 +0100 Subject: [PATCH 11/23] core: Improve documentation of RE classes --- doc/sphinx/advanced_methods.rst | 4 +-- src/core/reaction_ensemble.cpp | 57 +++++++++++++-------------------- src/core/reaction_ensemble.hpp | 40 ++++++++++++++--------- 3 files changed, 48 insertions(+), 53 deletions(-) diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index 95a7f5ec087..1bbf4b3d56c 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -1693,9 +1693,7 @@ In practice, this constant is often used with the dimension of :math:`(c^{\ominu A simulation in the reaction ensemble consists of two types of moves: the *reaction move* and the *configuration move*. The configuration move changes the configuration -of the system. It is not performed by the Reaction Ensemble module, and can be -performed by a suitable molecular dynamics or a Monte Carlo scheme. The -``reactant_ensemble`` command takes care only of the reaction moves. +of the system. In the *forward* reaction, the appropriate number of reactants (given by :math:`\nu_i`) is removed from the system, and the concomitant number of products is inserted into the system. In the *backward* reaction, diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 95fd85826d1..686569a2b26 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -56,17 +56,17 @@ void EnergyCollectiveVariable::load_CV_boundaries( m_current_wang_landau_system.do_energy_reweighting = true; // load energy boundaries from file std::ifstream infile; - - infile.open( - energy_boundaries_filename); // file containing numbers in 3 columns + infile.open(energy_boundaries_filename); if (infile.fail()) throw std::runtime_error("ERROR: energy boundaries file for the specific " "system could not be read.\n"); // Note that you cannot change the other collective variables in the // pre-production run and the production run - // Note: the number of collective variables is unknown, therefore we cannot - // read in the file in an easier way + + // The data is formatted as space-separated floaing point values + // (the first line is a header). The min and max energies are stored in + // the last two columns. The first N columns are the collective variables. std::string line = ""; std::getline(infile, line); // dummy read to throw away header while (std::getline(infile, line)) { @@ -363,10 +363,9 @@ void WangLandauReactionEnsemble::on_reaction_entry(int &old_state_index) { void WangLandauReactionEnsemble::on_reaction_rejection_directly_after_entry( int &old_state_index) { - update_wang_landau_potential_and_histogram( - old_state_index); // increase the Wang-Landau potential and histogram at - // the current nbar (this case covers the cases nbar=0 - // or nbar=1) + // increase the Wang-Landau potential and histogram at the current nbar + // (this case covers the cases nbar=0 or nbar=1) + update_wang_landau_potential_and_histogram(old_state_index); } void WangLandauReactionEnsemble::on_attempted_reaction(int &new_state_index) { @@ -609,8 +608,7 @@ int ReactionAlgorithm::create_particle(int desired_type) { } Utils::Vector3d pos_vec; - // create random velocity vector according to Maxwell Boltzmann distribution - // for components + // create random velocity vector according to Maxwell-Boltzmann distribution double vel[3]; // we use mass=1 for all particles, think about adapting this vel[0] = std::sqrt(temperature) * m_normal_distribution(m_generator); @@ -664,7 +662,7 @@ int WangLandauReactionEnsemble::on_mc_use_WL_get_new_state() { } /** - * Performs a global mc move for a particle of the provided type. + * Performs a global MC move for a particle of the provided type. */ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type( int type, int particle_number_of_type_to_be_changed, bool use_wang_landau) { @@ -950,18 +948,15 @@ double WangLandauReactionEnsemble::calculate_delta_degree_of_association( return delta; } -/** - * Initializes the Wang-Landau histogram. - */ int WangLandauReactionEnsemble::get_num_needed_bins() const { int needed_bins = 1; + // add 1 for degrees of association-related part of histogram (think of only + // one acid particle) for (const auto ¤t_collective_variable : collective_variables) { needed_bins = needed_bins * (int((current_collective_variable->CV_maximum - current_collective_variable->CV_minimum) / current_collective_variable->delta_CV) + - 1); // plus 1 needed for degrees of association - // related part of histogram (think of only - // one acid particle) + 1); } return needed_bins; } @@ -1010,29 +1005,23 @@ void WangLandauReactionEnsemble::invalidate_bins() { static_cast(wang_landau_potential.size()) - empty_bins_in_memory; } -/** - * Initializes the current Wang-Landau system. - */ int WangLandauReactionEnsemble::initialize_wang_landau() { nr_subindices_of_collective_variable.resize(collective_variables.size(), 0); auto const new_CV_i = collective_variables.size() - 1; + // add 1 for collective variables which are of type degree of association nr_subindices_of_collective_variable[new_CV_i] = int((collective_variables[new_CV_i]->CV_maximum - collective_variables[new_CV_i]->CV_minimum) / collective_variables[new_CV_i]->delta_CV) + - 1; //+1 for collective variables which are of type degree of association + 1; - // construct (possibly higher dimensional) histogram over gamma (the room - // which should be equally sampled when the Wang-Landau algorithm has - // converged) + // construct (possibly higher dimensional) histogram and potential over + // gamma (the space which should be equally sampled when the Wang-Landau + // algorithm has converged) int needed_bins = get_num_needed_bins(); - histogram.resize(needed_bins, 0); // initialize new values with 0 - - // construct (possibly higher dimensional) wang_landau potential over gamma - // (the room which should be equally sampled when the Wang-Landau algorithm - // has converged) - wang_landau_potential.resize(needed_bins, 0); // initialize new values with 0 + histogram.resize(needed_bins, 0); + wang_landau_potential.resize(needed_bins, 0); used_bins = needed_bins; // initialize for 1/t wang_landau algorithm @@ -1593,14 +1582,12 @@ WidomInsertion::measure_excess_chemical_potential(int reaction_id) { exp(-1.0 / temperature * (E_pot_new - E_pot_old))}; current_reaction.accumulator_exponentials(exponential); + // calculate mean excess chemical potential and standard error of the mean std::pair result = std::make_pair( -temperature * log(current_reaction.accumulator_exponentials.mean()[0]), std::abs(-temperature / current_reaction.accumulator_exponentials.mean()[0] * - current_reaction.accumulator_exponentials - .std_error()[0])); // excess chemical potential; error - // excess chemical potential, - // determined via error propagation + current_reaction.accumulator_exponentials.std_error()[0])); return result; } diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 9aa13276757..8bec73358b8 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -104,9 +104,17 @@ struct EnergyCollectiveVariable : public CollectiveVariable { load_CV_boundaries(WangLandauReactionEnsemble &m_current_wang_landau_system); }; +/** Measure the degree of association of a chemical species. + * As an example, consider a polybasic acid A which can be protonated + * into species AH and AH2. The degree of association of species A is + * equal to n(A) / (n(A) + n(AH) + n(AH2)). + */ struct DegreeOfAssociationCollectiveVariable : public CollectiveVariable { + /** List of all conjugated species */ std::vector corresponding_acid_types; + /** Reference species from which the degree of association is measured */ int associated_type; + /** Calculate the degree of association of the reference species */ double determine_current_state() const override { return calculate_degree_of_association(); } @@ -297,16 +305,13 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { std::vector min_boundaries_energies; std::vector max_boundaries_energies; - std::vector - minimum_energies_at_flat_index; // only present in energy preparation run - std::vector - maximum_energies_at_flat_index; // only present in energy preparation run + // only present in energy preparation run + std::vector minimum_energies_at_flat_index; + // only present in energy preparation run + std::vector maximum_energies_at_flat_index; - int update_maximum_and_minimum_energies_at_current_state(); // use for - // preliminary - // energy - // reweighting - // runs + // use for preliminary energy reweighting runs + int update_maximum_and_minimum_energies_at_current_state(); int do_reaction(int reaction_steps) override; void write_out_preliminary_energy_run_results(const std::string &filename); @@ -336,13 +341,15 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { int on_mc_use_WL_get_new_state() override; std::vector histogram; - std::vector wang_landau_potential; // equals the logarithm to basis e - // of the degeneracy of the states + + // logarithm to basis e of the degeneracy of the states + std::vector wang_landau_potential; std::vector nr_subindices_of_collective_variable; - double wang_landau_parameter = 1.0; // equals the logarithm to basis e of the - // modification factor of the degeneracy of + + // logarithm to basis e of the modification factor of the degeneracy of // states when the state is visited + double wang_landau_parameter = 1.0; int int_fill_value = -10; double double_fill_value = -10.0; @@ -370,8 +377,11 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { bool achieved_desired_number_of_refinements_one_over_t() const; void refine_wang_landau_parameter_one_over_t(); - int initialize_wang_landau(); // has to be called (at least) after the last - // collective variable is added + /** + * Initialize the current Wang-Landau system. + * Has to be called (at least) after the last collective variable is added. + */ + int initialize_wang_landau(); double calculate_delta_degree_of_association( DegreeOfAssociationCollectiveVariable ¤t_collective_variable); int get_num_needed_bins() const; From a8d0bd55b640142ec7188e8e6a794c39e2ed0c77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 17 Mar 2021 18:11:03 +0100 Subject: [PATCH 12/23] test: Improve clarity of RE unit tests --- .../reaction_ensemble_classes_test.cpp | 146 ++++++++++++------ 1 file changed, 97 insertions(+), 49 deletions(-) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index a115bbd3fe4..711dc9b7191 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -67,6 +67,8 @@ struct ParticleFactory { std::vector particle_cache; }; +// Check the mechanism that tracks particles of a certain type and the +// function that selects a random particle in the pool of tracked particles. BOOST_FIXTURE_TEST_CASE(particle_type_map_test, ParticleFactory) { constexpr double tol = 100 * std::numeric_limits::epsilon(); @@ -93,6 +95,10 @@ BOOST_FIXTURE_TEST_CASE(particle_type_map_test, ParticleFactory) { BOOST_CHECK_EQUAL(get_random_p_id(type, 0), pid); } +// Check the collective variable that tracks the degree of dissociation +// of a polybasic acid. The state measures the concentration of the base +// (i.e. fully deprotonated species) divided by the concentration of the +// base and all of its the conjugated acids. BOOST_FIXTURE_TEST_CASE(DegreeOfAssociationCollectiveVariable_test, ParticleFactory) { using namespace ReactionEnsemble; @@ -128,12 +134,19 @@ BOOST_FIXTURE_TEST_CASE(DegreeOfAssociationCollectiveVariable_test, BOOST_CHECK_CLOSE(doa_cv.determine_current_state(), 0.25, tol); } +// Check a simple chemical reaction, the Monte Carlo acceptance rate +// and the configurational move probability for a given system state. BOOST_AUTO_TEST_CASE(SingleReaction_test) { using namespace ReactionEnsemble; constexpr double tol = 100 * std::numeric_limits::epsilon(); - // check derived parameter - SingleReaction reaction(2., {0}, {1}, {1, 2}, {3, 4}); + // create a reaction A -> 3 B + 4 C + int const type_A = 0; + int const type_B = 1; + int const type_C = 2; + SingleReaction reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); + + // check derived parameter nu_bar BOOST_CHECK_EQUAL(reaction.nu_bar, 6); // check acceptance rate @@ -152,8 +165,9 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) { for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { for (int k = 0; k < 3; ++k) { - // i adduct #0, j product #1, k product #2 - auto const p_numbers = std::map{{0, i}, {1, j}, {2, 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 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); @@ -162,6 +176,7 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) { } } +// Check the base class for all Monte Carlo algorithms. BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { using namespace ReactionEnsemble; class ReactionAlgorithmTest : public ReactionAlgorithm { @@ -189,11 +204,15 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { // check reaction addition { - SingleReaction const ref(2., {0}, {1}, {1, 2}, {3, 4}); + // 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 ref(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); r_algo.add_reaction(ref.gamma, ref.reactant_types, ref.reactant_coefficients, ref.product_types, ref.product_coefficients); - BOOST_CHECK_EQUAL(r_algo.reactions.size(), 1ul); + BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 1ul); auto const &reaction = r_algo.reactions[0]; BOOST_TEST(reaction.reactant_types == ref.reactant_types, boost::test_tools::per_element()); @@ -235,20 +254,22 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { // check reaction removal { - r_algo.add_reaction(5, {1}, {1}, {2}, {1}); - BOOST_CHECK_EQUAL(r_algo.reactions.size(), 2ul); - BOOST_CHECK_EQUAL(r_algo.reactions[1].gamma, 5); + r_algo.add_reaction(5., {1}, {1}, {2}, {1}); + BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 2ul); + BOOST_CHECK_EQUAL(r_algo.reactions[1].gamma, 5.); r_algo.delete_reaction(1); - BOOST_CHECK_EQUAL(r_algo.reactions.size(), 1ul); - BOOST_CHECK_EQUAL(r_algo.reactions[0].gamma, 2); + BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 1ul); + BOOST_CHECK_EQUAL(r_algo.reactions[0].gamma, 2.); r_algo.delete_reaction(0); - BOOST_CHECK_EQUAL(r_algo.reactions.size(), 0ul); + BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 0ul); } // exception if deleting a non-existent particle BOOST_CHECK_THROW(r_algo.delete_particle(5), std::runtime_error); } +// Check the Monte Carlo algorithm where moves depend on the system +// configuration and energy. BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { class ReactionEnsembleTest : public ReactionEnsemble::ReactionEnsemble { public: @@ -265,25 +286,33 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { // exception if no reaction was added BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), 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 constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i; - SingleReaction const reaction(2., {0}, {1}, {1, 2}, {3, 4}); for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { for (int k = 0; k < 3; ++k) { - // i adduct #0, j product #1, k product #2 - auto const p_numbers = std::map{{0, i}, {1, j}, {2, 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); - auto const ref = 2e6 * f_expr * std::exp(energy / 20.); - auto const val = r_algo.calculate_acceptance_probability( + auto const acceptance_ref = 2e6 * f_expr * std::exp(energy / 20.); + auto const acceptance = r_algo.calculate_acceptance_probability( reaction, energy, 0., p_numbers, -1, -1, false); - BOOST_CHECK_CLOSE(val, ref, 5 * tol); + BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); } } } } +// Check the Monte Carlo algorithm where moves depend on the system +// configuration and/or energy. BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { using namespace ReactionEnsemble; class WangLandauReactionEnsembleTest : public WangLandauReactionEnsemble { @@ -300,39 +329,45 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { // exception if no reaction was added BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), 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 constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i; - SingleReaction const reaction(2., {0}, {1}, {1, 2}, {3, 4}); for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { for (int k = 0; k < 3; ++k) { - // i adduct #0, j product #1, k product #2 - auto const p_numbers = std::map{{0, i}, {1, j}, {2, 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); auto const prefactor = 2e6 * f_expr; auto const boltzmann = std::exp(energy / 20.); - double val; + double acceptance; r_algo.do_energy_reweighting = false; - val = r_algo.calculate_acceptance_probability(reaction, energy, 0., - p_numbers, -1, -1, false); - BOOST_CHECK_CLOSE(val, prefactor * boltzmann, 5 * tol); - val = r_algo.calculate_acceptance_probability(reaction, energy, 0., - p_numbers, -1, -1, true); - BOOST_CHECK_CLOSE(val, boltzmann, 5 * tol); + acceptance = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, -1, false); + BOOST_CHECK_CLOSE(acceptance, prefactor * boltzmann, 5 * tol); + acceptance = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, -1, true); + BOOST_CHECK_CLOSE(acceptance, boltzmann, 5 * tol); r_algo.do_energy_reweighting = true; - val = r_algo.calculate_acceptance_probability(reaction, energy, 0., - p_numbers, -1, -1, false); - BOOST_CHECK_CLOSE(val, prefactor, 5 * tol); - val = r_algo.calculate_acceptance_probability(reaction, energy, 0., - p_numbers, -1, -1, true); - BOOST_CHECK_CLOSE(val, 1., 5 * tol); - val = r_algo.calculate_acceptance_probability(reaction, energy, 0., - p_numbers, -1, 0, true); - BOOST_CHECK_CLOSE(val, 10., 5 * tol); - val = r_algo.calculate_acceptance_probability(reaction, energy, 0., - p_numbers, 0, -1, true); - BOOST_CHECK_CLOSE(val, -10., 5 * tol); + acceptance = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, -1, false); + BOOST_CHECK_CLOSE(acceptance, prefactor, 5 * tol); + acceptance = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, -1, true); + BOOST_CHECK_CLOSE(acceptance, 1., 5 * tol); + acceptance = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, 0, true); + BOOST_CHECK_CLOSE(acceptance, 10., 5 * tol); + acceptance = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, 0, -1, true); + BOOST_CHECK_CLOSE(acceptance, -10., 5 * tol); } } } @@ -401,6 +436,8 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { } } +// Check the Monte Carlo algorithm where moves depend on the system +// configuration, energy and pH. BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { using namespace ReactionEnsemble; class ConstantpHEnsembleTest : public ConstantpHEnsemble { @@ -417,20 +454,26 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { // exception if no reaction was added BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), 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 constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i; - SingleReaction const reaction(2., {0}, {1}, {1, 2}, {1, 1}); for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { for (int k = 0; k < 3; ++k) { - // i adduct #0, j product #1, k product #2 - auto const p_numbers = std::map{{0, i}, {1, j}, {2, 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 ref = + auto const acceptance_ref = std::exp(energy / 20. + std::log(10.) * (1. + std::log10(2.))); - auto const val = r_algo.calculate_acceptance_probability( + auto const acceptance = r_algo.calculate_acceptance_probability( reaction, energy, 0., p_numbers, -1, -1, false); - BOOST_CHECK_CLOSE(val, ref, 5 * tol); + BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); } } } @@ -443,9 +486,14 @@ BOOST_AUTO_TEST_CASE(WidomInsertion_test) { // check acceptance rate WidomInsertion r_algo(42); - SingleReaction const ref(2., {0}, {1}, {1, 2}, {3, 4}); - r_algo.add_reaction(ref.gamma, ref.reactant_types, ref.reactant_coefficients, - ref.product_types, ref.product_coefficients); + // 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}); + r_algo.add_reaction(reaction.gamma, reaction.reactant_types, + reaction.reactant_coefficients, reaction.product_types, + reaction.product_coefficients); // exception if not enough particles BOOST_CHECK_THROW(r_algo.measure_excess_chemical_potential(0), From 3d282f8a3ef0d932bbbe8e4f5a7e8059dd71ff66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 18 Mar 2021 09:16:33 +0100 Subject: [PATCH 13/23] core: Simplify RE code --- src/core/reaction_ensemble.cpp | 84 +++++++++++++++------------------- src/core/reaction_ensemble.hpp | 27 +++++------ 2 files changed, 48 insertions(+), 63 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 686569a2b26..d7eed526a4c 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -31,6 +31,8 @@ #include #include +#include + #include #include #include @@ -51,9 +53,9 @@ namespace ReactionEnsemble { * variables under min_boundaries_energies, max_boundaries_energies */ void EnergyCollectiveVariable::load_CV_boundaries( - WangLandauReactionEnsemble &m_current_wang_landau_system) { + WangLandauReactionEnsemble &wl_system) { - m_current_wang_landau_system.do_energy_reweighting = true; + wl_system.do_energy_reweighting = true; // load energy boundaries from file std::ifstream infile; infile.open(energy_boundaries_filename); @@ -76,18 +78,14 @@ void EnergyCollectiveVariable::load_CV_boundaries( while (iss >> value) { values.push_back(value); } - m_current_wang_landau_system.min_boundaries_energies.push_back( - values[values.size() - 2]); - m_current_wang_landau_system.max_boundaries_energies.push_back( - values[values.size() - 1]); - } - - CV_minimum = *(std::min_element( - m_current_wang_landau_system.min_boundaries_energies.begin(), - m_current_wang_landau_system.min_boundaries_energies.end())); - CV_maximum = *(std::max_element( - m_current_wang_landau_system.max_boundaries_energies.begin(), - m_current_wang_landau_system.max_boundaries_energies.end())); + assert(values.size() >= 2); + auto it = values.end(); + wl_system.max_boundaries_energies.push_back(*(--it)); + wl_system.min_boundaries_energies.push_back(*(--it)); + } + + CV_minimum = *boost::range::min_element(wl_system.min_boundaries_energies); + CV_maximum = *boost::range::max_element(wl_system.max_boundaries_energies); } /** @@ -105,12 +103,12 @@ int ReactionAlgorithm::do_reaction(int reaction_steps) { * Adds a reaction to the reaction system */ void ReactionAlgorithm::add_reaction( - double gamma, const std::vector &_reactant_types, - const std::vector &_reactant_coefficients, - const std::vector &_product_types, - const std::vector &_product_coefficients) { - SingleReaction new_reaction(gamma, _reactant_types, _reactant_coefficients, - _product_types, _product_coefficients); + double gamma, const std::vector &reactant_types, + const std::vector &reactant_coefficients, + const std::vector &product_types, + const std::vector &product_coefficients) { + SingleReaction new_reaction(gamma, reactant_types, reactant_coefficients, + product_types, product_coefficients); // make ESPResSo count the particle numbers which take part in the reactions for (int reactant_type : new_reaction.reactant_types) @@ -810,10 +808,8 @@ void WangLandauReactionEnsemble::add_new_CV_potential_energy( std::make_shared(); new_collective_variable->energy_boundaries_filename = filename; new_collective_variable->delta_CV = delta_CV; - collective_variables.push_back(new_collective_variable); new_collective_variable->load_CV_boundaries(*this); - collective_variables[collective_variables.size() - 1] = - new_collective_variable; + collective_variables.emplace_back(new_collective_variable); initialize_wang_landau(); } @@ -948,15 +944,13 @@ double WangLandauReactionEnsemble::calculate_delta_degree_of_association( return delta; } -int WangLandauReactionEnsemble::get_num_needed_bins() const { - int needed_bins = 1; +long WangLandauReactionEnsemble::get_num_needed_bins() const { + long needed_bins = 1; // add 1 for degrees of association-related part of histogram (think of only // one acid particle) - for (const auto ¤t_collective_variable : collective_variables) { - needed_bins = needed_bins * (int((current_collective_variable->CV_maximum - - current_collective_variable->CV_minimum) / - current_collective_variable->delta_CV) + - 1); + for (const auto &cv : collective_variables) { + needed_bins *= + static_cast((cv->CV_maximum - cv->CV_minimum) / cv->delta_CV) + 1; } return needed_bins; } @@ -966,7 +960,7 @@ void WangLandauReactionEnsemble::invalidate_bins() { // allowed at the given degree of association, because the energy boundaries // prohibit them - int empty_bins_in_memory = 0; + long empty_bins_in_memory = 0; for (std::size_t flattened_index = 0; flattened_index < wang_landau_potential.size(); flattened_index++) { std::vector unraveled_index( @@ -997,29 +991,26 @@ void WangLandauReactionEnsemble::invalidate_bins() { energy_CV->delta_CV) { histogram[flattened_index] = int_fill_value; wang_landau_potential[flattened_index] = double_fill_value; - empty_bins_in_memory += 1; + empty_bins_in_memory++; } } used_bins = - static_cast(wang_landau_potential.size()) - empty_bins_in_memory; + static_cast(wang_landau_potential.size()) - empty_bins_in_memory; } int WangLandauReactionEnsemble::initialize_wang_landau() { nr_subindices_of_collective_variable.resize(collective_variables.size(), 0); - auto const new_CV_i = collective_variables.size() - 1; + auto const &cv = collective_variables.back(); // add 1 for collective variables which are of type degree of association - nr_subindices_of_collective_variable[new_CV_i] = - int((collective_variables[new_CV_i]->CV_maximum - - collective_variables[new_CV_i]->CV_minimum) / - collective_variables[new_CV_i]->delta_CV) + - 1; + nr_subindices_of_collective_variable.back() = + static_cast((cv->CV_maximum - cv->CV_minimum) / cv->delta_CV) + 1; // construct (possibly higher dimensional) histogram and potential over // gamma (the space which should be equally sampled when the Wang-Landau // algorithm has converged) - int needed_bins = get_num_needed_bins(); + auto const needed_bins = get_num_needed_bins(); histogram.resize(needed_bins, 0); wang_landau_potential.resize(needed_bins, 0); @@ -1039,12 +1030,12 @@ double WangLandauReactionEnsemble::calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, std::map const &old_particle_numbers, int old_state_index, int new_state_index, bool only_make_configuration_changing_move) const { + double beta = 1.0 / temperature; - double bf; - if (do_not_sample_reaction_partition_function || - only_make_configuration_changing_move) { - bf = 1.0; - } else { + double bf = 1.0; + + if (!(do_not_sample_reaction_partition_function || + only_make_configuration_changing_move)) { auto const factorial_expr = calculate_factorial_expression(current_reaction, old_particle_numbers); bf = std::pow(volume, current_reaction.nu_bar) * current_reaction.gamma * @@ -1053,9 +1044,8 @@ double WangLandauReactionEnsemble::calculate_acceptance_probability( if (!do_energy_reweighting) { bf *= exp(-beta * (E_pot_new - E_pot_old)); - } else { - // pass } + // Check whether the proposed state lies in the reaction coordinate space // gamma and add the Wang-Landau modification factor, this is a bit nasty // due to the energy collective variable case (memory layout of storage diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 8bec73358b8..ba28635bcb8 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -45,14 +45,10 @@ struct SingleReaction { std::vector const &reactant_coefficients, std::vector const &product_types, std::vector const &product_coefficients) { - std::copy(reactant_types.begin(), reactant_types.end(), - std::back_inserter(this->reactant_types)); - std::copy(reactant_coefficients.begin(), reactant_coefficients.end(), - std::back_inserter(this->reactant_coefficients)); - std::copy(product_types.begin(), product_types.end(), - std::back_inserter(this->product_types)); - std::copy(product_coefficients.begin(), product_coefficients.end(), - std::back_inserter(this->product_coefficients)); + this->reactant_types = reactant_types; + this->reactant_coefficients = reactant_coefficients; + this->product_types = product_types; + this->product_coefficients = product_coefficients; this->gamma = gamma; nu_bar = std::accumulate(product_coefficients.begin(), product_coefficients.end(), 0) - @@ -100,8 +96,7 @@ struct EnergyCollectiveVariable : public CollectiveVariable { double determine_current_state() const override { return calculate_current_potential_energy_of_system(); } - void - load_CV_boundaries(WangLandauReactionEnsemble &m_current_wang_landau_system); + void load_CV_boundaries(WangLandauReactionEnsemble &wl_system); }; /** Measure the degree of association of a chemical species. @@ -187,10 +182,10 @@ class ReactionAlgorithm { void check_reaction_ensemble() const; int delete_particle(int p_id); - void add_reaction(double gamma, const std::vector &_reactant_types, - const std::vector &_reactant_coefficients, - const std::vector &_product_types, - const std::vector &_product_coefficients); + void add_reaction(double gamma, const std::vector &reactant_types, + const std::vector &reactant_coefficients, + const std::vector &product_types, + const std::vector &product_coefficients); void delete_reaction(int reaction_id) { reactions.erase(reactions.begin() + reaction_id); } @@ -354,7 +349,7 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { int int_fill_value = -10; double double_fill_value = -10.0; - int used_bins = -10; // for 1/t algorithm + long used_bins = -10; // for 1/t algorithm int monte_carlo_trial_moves = 0; // for 1/t algorithm int get_flattened_index_wang_landau_without_energy_collective_variable( @@ -384,7 +379,7 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { int initialize_wang_landau(); double calculate_delta_degree_of_association( DegreeOfAssociationCollectiveVariable ¤t_collective_variable); - int get_num_needed_bins() const; + long get_num_needed_bins() const; void invalidate_bins(); void reset_histogram(); double get_minimum_CV_value_on_delta_CV_spaced_grid(double min_CV_value, From ccfd29069658eca72919595d78a0625b0bb10d6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 18 Mar 2021 09:25:24 +0100 Subject: [PATCH 14/23] tests: Rename Reaction Ensemble statistical tests Also do some cosmetic changes. --- testsuite/python/CMakeLists.txt | 4 +-- .../{constant_pH.py => constant_pH_stats.py} | 9 ++---- ...ction_ensemble.py => wang_landau_stats.py} | 30 +++++++++---------- 3 files changed, 19 insertions(+), 24 deletions(-) rename testsuite/python/{constant_pH.py => constant_pH_stats.py} (93%) rename testsuite/python/{wang_landau_reaction_ensemble.py => wang_landau_stats.py} (86%) diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index f13ae95562d..d71099be2ce 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -119,7 +119,7 @@ python_test(FILE rotational_dynamics.py MAX_NUM_PROC 1) python_test(FILE script_interface_object_params.py MAX_NUM_PROC 4) python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4) python_test(FILE widom_insertion.py MAX_NUM_PROC 1) -python_test(FILE constant_pH.py MAX_NUM_PROC 4 LABELS long) +python_test(FILE constant_pH_stats.py MAX_NUM_PROC 4 LABELS long) python_test(FILE writevtf.py MAX_NUM_PROC 4) python_test(FILE lb_stokes_sphere.py MAX_NUM_PROC 4 LABELS gpu long) python_test(FILE lb_pressure_tensor.py MAX_NUM_PROC 1 LABELS gpu long) @@ -171,7 +171,7 @@ python_test(FILE analyze_distance.py MAX_NUM_PROC 1) python_test(FILE analyze_acf.py MAX_NUM_PROC 1) python_test(FILE comfixed.py MAX_NUM_PROC 2) python_test(FILE rescale.py MAX_NUM_PROC 2) -python_test(FILE wang_landau_reaction_ensemble.py MAX_NUM_PROC 1 LABELS long) +python_test(FILE wang_landau_stats.py MAX_NUM_PROC 1 LABELS long) python_test(FILE array_properties.py MAX_NUM_PROC 4) python_test(FILE analyze_distribution.py MAX_NUM_PROC 1) python_test(FILE observable_profile.py MAX_NUM_PROC 4) diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH_stats.py similarity index 93% rename from testsuite/python/constant_pH.py rename to testsuite/python/constant_pH_stats.py index 878f5f83141..5ed6c79bf26 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH_stats.py @@ -17,17 +17,15 @@ # along with this program. If not, see . # -"""Testmodule for the Reaction Ensemble. -""" import unittest as ut import numpy as np import espressomd -from espressomd import reaction_ensemble +import espressomd.reaction_ensemble class ReactionEnsembleTest(ut.TestCase): - """Test the core implementation of the reaction ensemble.""" + """Test the core implementation of the constant pH reaction ensemble.""" N0 = 40 c0 = 0.00028 @@ -35,7 +33,6 @@ class ReactionEnsembleTest(ut.TestCase): type_A = 1 type_H = 5 temperature = 1.0 - # avoid extreme regions in the titration curve e.g. via the choice # choose target alpha not too far from 0.5 to get good statistics in a # small number of steps pKa_minus_pH = -0.2 @@ -47,7 +44,7 @@ class ReactionEnsembleTest(ut.TestCase): np.random.seed(69) # make reaction code fully deterministic system.cell_system.skin = 0.4 system.time_step = 0.01 - RE = reaction_ensemble.ConstantpHEnsemble( + RE = espressomd.reaction_ensemble.ConstantpHEnsemble( temperature=1.0, exclusion_radius=1, seed=44) @classmethod diff --git a/testsuite/python/wang_landau_reaction_ensemble.py b/testsuite/python/wang_landau_stats.py similarity index 86% rename from testsuite/python/wang_landau_reaction_ensemble.py rename to testsuite/python/wang_landau_stats.py index f534517230e..be4b7835aa1 100644 --- a/testsuite/python/wang_landau_reaction_ensemble.py +++ b/testsuite/python/wang_landau_stats.py @@ -22,18 +22,17 @@ import unittest as ut import espressomd -from espressomd.interactions import HarmonicBond -from espressomd import reaction_ensemble +import espressomd.interactions +import espressomd.reaction_ensemble -class ReactionEnsembleTest(ut.TestCase): +class WangLandauReactionEnsembleTest(ut.TestCase): - """Test the core implementation of the wang_landau reaction ensemble. + """Test the core implementation of the Wang-Landau reaction ensemble. - Create a harmonic bond between the two reacting particles. Therefore the - potential energy is quadratic in the elongation of the bond and - therefore the density of states is known as the one of the harmonic - oscillator + Create a harmonic bond between the two reacting particles. Therefore + the potential energy is quadratic in the elongation of the bond and + the density of states is known as the one of the harmonic oscillator. """ # System parameters @@ -53,18 +52,17 @@ class ReactionEnsembleTest(ut.TestCase): # Setup System # - N0 = 1 # number of titratable units K_diss = 0.0088 - system.part.add(id=0, pos=[0, 0, 0] * system.box_l, type=3) - system.part.add(id=1, pos=[1.0, 1.0, 1.0] * system.box_l / 2.0, type=1) - system.part.add(id=2, pos=np.random.random() * system.box_l, type=2) - system.part.add(id=3, pos=np.random.random() * system.box_l, type=2) + system.part.add(id=0, type=3, pos=[0, 0, 0] * system.box_l) + system.part.add(id=1, type=1, pos=[1.0, 1.0, 1.0] * system.box_l / 2.0) + system.part.add(id=2, type=2, pos=np.random.random() * system.box_l) + system.part.add(id=3, type=2, pos=np.random.random() * system.box_l) - h = HarmonicBond(r_0=0, k=1) + h = espressomd.interactions.HarmonicBond(r_0=0, k=1) system.bonded_inter[0] = h system.part[0].add_bond((h, 1)) - WLRE = reaction_ensemble.WangLandauReactionEnsemble( + WLRE = espressomd.reaction_ensemble.WangLandauReactionEnsemble( temperature=temperature, exclusion_radius=0, seed=86) WLRE.add_reaction( gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], @@ -106,7 +104,7 @@ def test_wang_landau_output(self): self.WLRE.reaction() for _ in range(2): self.WLRE.displacement_mc_move_for_particles_of_type(3) - except reaction_ensemble.WangLandauHasConverged: + except espressomd.reaction_ensemble.WangLandauHasConverged: break nbars, Epots, WL_potentials = np.loadtxt(self.file_output, unpack=True) From f40d61e3bb7d0ad2447be5ee25466859908ab651 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 18 Mar 2021 09:28:03 +0100 Subject: [PATCH 15/23] tests: Split Wang-Landau test Move the Wang-Landau output files checks to an interface test. The statistical test runs the Monte Carlo algorithm until convergence. --- testsuite/python/CMakeLists.txt | 1 + testsuite/python/wang_landau.py | 130 ++++++++++++++++++++++++++ testsuite/python/wang_landau_stats.py | 39 +------- 3 files changed, 134 insertions(+), 36 deletions(-) create mode 100644 testsuite/python/wang_landau.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index d71099be2ce..172525f1fd6 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -171,6 +171,7 @@ python_test(FILE analyze_distance.py MAX_NUM_PROC 1) python_test(FILE analyze_acf.py MAX_NUM_PROC 1) python_test(FILE comfixed.py MAX_NUM_PROC 2) python_test(FILE rescale.py MAX_NUM_PROC 2) +python_test(FILE wang_landau.py MAX_NUM_PROC 1) python_test(FILE wang_landau_stats.py MAX_NUM_PROC 1 LABELS long) python_test(FILE array_properties.py MAX_NUM_PROC 4) python_test(FILE analyze_distribution.py MAX_NUM_PROC 1) diff --git a/testsuite/python/wang_landau.py b/testsuite/python/wang_landau.py new file mode 100644 index 00000000000..a3e85b335ac --- /dev/null +++ b/testsuite/python/wang_landau.py @@ -0,0 +1,130 @@ +# +# Copyright (C) 2013-2019 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 . +# +"""Testmodule for the Wang-Landau Reaction Ensemble. +""" +import numpy as np +import unittest as ut + +import espressomd +import espressomd.interactions +import espressomd.reaction_ensemble + + +class WangLandauReactionEnsembleTest(ut.TestCase): + + """Test the interface of the Wang-Landau reaction ensemble.""" + + # System parameters + # + box_l = 6 * np.sqrt(2) + temperature = 1.0 + + # Integration parameters + # + system = espressomd.System(box_l=[box_l, box_l, box_l]) + np.random.seed(seed=42) + system.time_step = 0.01 + system.cell_system.skin = 0 + system.cell_system.set_n_square(use_verlet_lists=False) + + # + # Setup System + # + + K_diss = 0.0088 + + system.part.add(id=0, type=3, pos=[0, 0, 0] * system.box_l) + system.part.add(id=1, type=1, pos=[1.0, 1.0, 1.0] * system.box_l / 2.0) + system.part.add(id=2, type=2, pos=np.random.random() * system.box_l) + system.part.add(id=3, type=2, pos=np.random.random() * system.box_l) + + h = espressomd.interactions.HarmonicBond(r_0=0, k=1) + system.bonded_inter[0] = h + system.part[0].add_bond((h, 1)) + WLRE = espressomd.reaction_ensemble.WangLandauReactionEnsemble( + temperature=temperature, exclusion_radius=0, seed=86) + WLRE.add_reaction( + gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], + product_types=[1, 2], product_coefficients=[1, 1], + default_charges={0: 0, 1: -1, 2: +1}) + system.setup_type_map([0, 1, 2, 3]) + # initialize wang_landau + file_input = "energy_boundaries.dat" + file_output = "WL_potential_out.dat" + # generate preliminary_energy_run_results here, this should be done in a + # separate simulation without energy reweighting using the update energy + # functions + np.savetxt(file_input, np.transpose([[0, 1], [0, 0], [9, 9]]), + delimiter='\t', header="nbar E_potmin E_potmax") + + WLRE.add_collective_variable_degree_of_association( + associated_type=0, min=0, max=1, corresponding_acid_types=[0, 1]) + WLRE.set_wang_landau_parameters( + final_wang_landau_parameter=0.8 * 1e-2, + do_not_sample_reaction_partition_function=True, + full_path_to_output_filename=file_output) + + def test_01_energy_recording(self): + self.WLRE.update_maximum_and_minimum_energies_at_current_state() + self.WLRE.write_out_preliminary_energy_run_results() + nbars, E_mins, E_maxs = np.loadtxt( + "preliminary_energy_run_results", unpack=True) + np.testing.assert_almost_equal(nbars, [0, 1]) + np.testing.assert_almost_equal(E_mins, [27.0, -10]) + np.testing.assert_almost_equal(E_maxs, [27.0, -10]) + + def check_checkpoint(self, filename): + # write first checkpoint + self.WLRE.write_wang_landau_checkpoint() + old_checkpoint = np.loadtxt(filename) + + # modify old_checkpoint in memory and in file (this destroys the + # information contained in the checkpoint, but allows for testing of + # the functions) + modified_checkpoint = old_checkpoint + modified_checkpoint[0] = 1 + np.savetxt(filename, modified_checkpoint) + + # check whether changes are carried out correctly + self.WLRE.load_wang_landau_checkpoint() + self.WLRE.write_wang_landau_checkpoint() + new_checkpoint = np.loadtxt(filename) + np.testing.assert_almost_equal(new_checkpoint, modified_checkpoint) + + def test_02_checkpointing(self): + self.WLRE.add_collective_variable_potential_energy( + filename=self.file_input, delta=0.05) + + # run MC for long enough to sample a non-trivial histogram + for _ in range(1000): + try: + self.WLRE.reaction() + for _ in range(2): + self.WLRE.displacement_mc_move_for_particles_of_type(3) + except espressomd.reaction_ensemble.WangLandauHasConverged: + break + + filenames = ["checkpoint_wang_landau_potential_checkpoint", + "checkpoint_wang_landau_histogram_checkpoint"] + for filename in filenames: + self.check_checkpoint(filename) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/wang_landau_stats.py b/testsuite/python/wang_landau_stats.py index be4b7835aa1..68ad1c62489 100644 --- a/testsuite/python/wang_landau_stats.py +++ b/testsuite/python/wang_landau_stats.py @@ -70,8 +70,8 @@ class WangLandauReactionEnsembleTest(ut.TestCase): default_charges={0: 0, 1: -1, 2: +1}) system.setup_type_map([0, 1, 2, 3]) # initialize wang_landau - file_input = "energy_boundaries.dat" - file_output = "WL_potential_out.dat" + file_input = "energy_boundaries_stats.dat" + file_output = "WL_potential_out_stats.dat" # generate preliminary_energy_run_results here, this should be done in a # separate simulation without energy reweighting using the update energy # functions @@ -85,16 +85,7 @@ class WangLandauReactionEnsembleTest(ut.TestCase): do_not_sample_reaction_partition_function=True, full_path_to_output_filename=file_output) - def test_wang_landau_energy_recording(self): - self.WLRE.update_maximum_and_minimum_energies_at_current_state() - self.WLRE.write_out_preliminary_energy_run_results() - nbars, E_mins, E_maxs = np.loadtxt( - "preliminary_energy_run_results", unpack=True) - np.testing.assert_almost_equal(nbars, [0, 1]) - np.testing.assert_almost_equal(E_mins, [27.0, -10]) - np.testing.assert_almost_equal(E_maxs, [27.0, -10]) - - def test_wang_landau_output(self): + def test_potential_and_heat_capacity(self): self.WLRE.add_collective_variable_potential_energy( filename=self.file_input, delta=0.05) @@ -131,30 +122,6 @@ def calc_from_partition_function(quantity): heat_capacity, 1.5, places=1, msg="difference to analytical expected canonical configurational heat capacity too big") - def _wang_landau_output_checkpoint(self, filename): - # write first checkpoint - self.WLRE.write_wang_landau_checkpoint() - old_checkpoint = np.loadtxt(filename) - - # modify old_checkpoint in memory and in file (this destroys the - # information contained in the checkpoint, but allows for testing of - # the functions) - modified_checkpoint = old_checkpoint - modified_checkpoint[0] = 1 - np.savetxt(filename, modified_checkpoint) - - # check whether changes are carried out correctly - self.WLRE.load_wang_landau_checkpoint() - self.WLRE.write_wang_landau_checkpoint() - new_checkpoint = np.loadtxt(filename) - np.testing.assert_almost_equal(new_checkpoint, modified_checkpoint) - - def test_wang_landau_output_checkpoint(self): - filenames = ["checkpoint_wang_landau_potential_checkpoint", - "checkpoint_wang_landau_histogram_checkpoint"] - for filename in filenames: - self._wang_landau_output_checkpoint(filename) - if __name__ == "__main__": ut.main() From ce04cc6a403a8d4a2d726fac408eb5ca768b6289 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 18 Mar 2021 09:46:25 +0100 Subject: [PATCH 16/23] tests: Create a faster constant pH test Sample at a pH farther away from the pKa, shorten warmup phase, run with 1 MPI rank to remove large overhead from communication. --- testsuite/python/CMakeLists.txt | 1 + testsuite/python/constant_pH.py | 85 +++++++++++++++++++++++++++++++++ 2 files changed, 86 insertions(+) create mode 100644 testsuite/python/constant_pH.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 172525f1fd6..a94ba79bc5d 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -119,6 +119,7 @@ python_test(FILE rotational_dynamics.py MAX_NUM_PROC 1) python_test(FILE script_interface_object_params.py MAX_NUM_PROC 4) python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4) python_test(FILE widom_insertion.py MAX_NUM_PROC 1) +python_test(FILE constant_pH.py MAX_NUM_PROC 1) python_test(FILE constant_pH_stats.py MAX_NUM_PROC 4 LABELS long) python_test(FILE writevtf.py MAX_NUM_PROC 4) python_test(FILE lb_stokes_sphere.py MAX_NUM_PROC 4 LABELS gpu long) diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py new file mode 100644 index 00000000000..1f7a67c02ab --- /dev/null +++ b/testsuite/python/constant_pH.py @@ -0,0 +1,85 @@ +# +# Copyright (C) 2013-2021 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 . +# + +import unittest as ut +import numpy as np +import espressomd +import espressomd.reaction_ensemble + + +class ConstantpHTest(ut.TestCase): + + """Test the core implementation of the constant pH reaction ensemble.""" + + np.random.seed(42) + + def test_ideal_alpha(self): + N0 = 40 + c0 = 0.00028 + type_HA = 0 + type_A = 1 + type_H = 5 + pH = 2.0 + pKa = 2.5 + box_l = (N0 / c0)**(1.0 / 3.0) + system = espressomd.System(box_l=3 * [box_l]) + system.cell_system.skin = 0.4 + system.time_step = 0.01 + system.part.add(pos=np.random.random((2 * N0, 3)) * system.box_l, + type=N0 * [type_A, type_H]) + + RE = espressomd.reaction_ensemble.ConstantpHEnsemble( + temperature=1.0, + exclusion_radius=1, + seed=44) + RE.add_reaction( + gamma=10**(-pKa), + reactant_types=[type_HA], + reactant_coefficients=[1], + product_types=[type_A, type_H], + product_coefficients=[1, 1], + default_charges={type_HA: 0, type_A: -1, type_H: +1}) + RE.constant_pH = pH + + # equilibration + RE.reaction(800) + + # sampling + alphas = [] + for _ in range(80): + RE.reaction(15) + num_H = system.number_of_particles(type=type_H) + num_HA = system.number_of_particles(type=type_HA) + num_A = system.number_of_particles(type=type_A) + self.assertEqual(num_H, num_A) + self.assertEqual(num_A + num_HA, N0) + alphas.append(num_A / N0) + + # note: you cannot calculate the pH via -log10(/volume) in the + # constant pH ensemble, since the volume is totally arbitrary and does + # not influence the average number of protons + target_alpha = 1. / (1. + 10.**(pKa - pH)) + alpha_mean = np.mean(alphas) + alpha_error = np.std(alphas) / np.sqrt(len(alphas)) + self.assertAlmostEqual(alpha_mean, target_alpha, delta=0.015) + self.assertAlmostEqual(alpha_error, 0., delta=0.01) + + +if __name__ == "__main__": + ut.main() From 76332ae4c0ea0077b607fab364efb692172502f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 18 Mar 2021 13:00:50 +0100 Subject: [PATCH 17/23] tests: Remove UB in unit test --- src/core/unit_tests/reaction_ensemble_classes_test.cpp | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index 711dc9b7191..1f7729a36f0 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -421,15 +421,6 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { BOOST_TEST(bin_edges == bin_edges_ref, boost::test_tools::per_element()); BOOST_TEST(histogram == histogram_ref, boost::test_tools::per_element()); } - // add collective variable (with different energies in WL object) - r_algo.min_boundaries_energies.emplace_back(-1.); - r_algo.max_boundaries_energies.emplace_back(10.); - r_algo.add_new_CV_potential_energy(filename_input, delta_CV); - BOOST_REQUIRE_EQUAL(r_algo.collective_variables.size(), 2ul); - BOOST_REQUIRE_EQUAL(r_algo.min_boundaries_energies.size(), 5ul); - BOOST_REQUIRE_EQUAL(r_algo.max_boundaries_energies.size(), 5ul); - BOOST_CHECK_EQUAL(r_algo.collective_variables[1]->CV_minimum, -1.); - BOOST_CHECK_EQUAL(r_algo.collective_variables[1]->CV_maximum, 10.); // exception if file doesn't exist BOOST_CHECK_THROW(r_algo.add_new_CV_potential_energy("unknown", 0.), std::runtime_error); From a942d87bcbb55938240c0a5630d2a8caff1b037b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Mar 2021 14:18:03 +0100 Subject: [PATCH 18/23] core: Refactor Wang-Landau write function Separate the formatting code from the file opening code. Helps with unit testing. --- src/core/reaction_ensemble.cpp | 39 +++++++++++-------- src/core/reaction_ensemble.hpp | 1 + .../reaction_ensemble_classes_test.cpp | 11 +++--- 3 files changed, 28 insertions(+), 23 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index d7eed526a4c..09cf107cb64 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -59,8 +59,8 @@ void EnergyCollectiveVariable::load_CV_boundaries( // load energy boundaries from file std::ifstream infile; infile.open(energy_boundaries_filename); - if (infile.fail()) - throw std::runtime_error("ERROR: energy boundaries file for the specific " + if (!infile.is_open()) + throw std::runtime_error("energy boundaries file for the specific " "system could not be read.\n"); // Note that you cannot change the other collective variables in the @@ -1204,17 +1204,7 @@ bool WangLandauReactionEnsemble:: return false; } -/** - *Writes the Wang-Landau potential to file. - */ -void WangLandauReactionEnsemble::write_wang_landau_results_to_file( - const std::string &full_path_to_output_filename) { - - FILE *pFile; - pFile = fopen(full_path_to_output_filename.c_str(), "w"); - if (pFile == nullptr) { - throw std::runtime_error("ERROR: Wang-Landau file could not be written\n"); - } +void WangLandauReactionEnsemble::format_wang_landau_results(std::ostream &out) { for (std::size_t flattened_index = 0; flattened_index < wang_landau_potential.size(); flattened_index++) { // unravel index @@ -1236,13 +1226,28 @@ void WangLandauReactionEnsemble::write_wang_landau_results_to_file( auto const value = static_cast(unraveled_index[i]) * collective_variables[i]->delta_CV + collective_variables[i]->CV_minimum; - fprintf(pFile, "%f ", value); + out << value << " "; } - fprintf(pFile, "%f \n", wang_landau_potential[flattened_index]); + out << wang_landau_potential[flattened_index] << " \n"; } } - fflush(pFile); - fclose(pFile); + out.flush(); +} + +/** + *Writes the Wang-Landau potential to file. + */ +void WangLandauReactionEnsemble::write_wang_landau_results_to_file( + const std::string &full_path_to_output_filename) { + std::ofstream outfile; + + outfile.open(full_path_to_output_filename); + if (!outfile.is_open()) { + throw std::runtime_error("Wang-Landau file could not be opened\n"); + } + + format_wang_landau_results(outfile); + outfile.close(); } /** diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index ba28635bcb8..97209d64bc9 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -323,6 +323,7 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { double E_pot_new, std::map const &old_particle_numbers, int old_state_index, int new_state_index, bool only_make_configuration_changing_move) const override; + void format_wang_landau_results(std::ostream &out); private: void on_reaction_entry(int &old_state_index) override; diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index 1f7729a36f0..6bfc40692f7 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -42,8 +42,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -318,6 +320,7 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { class WangLandauReactionEnsembleTest : public WangLandauReactionEnsemble { public: using WangLandauReactionEnsemble::calculate_acceptance_probability; + using WangLandauReactionEnsemble::format_wang_landau_results; using WangLandauReactionEnsemble::WangLandauReactionEnsemble; }; constexpr double tol = 100 * std::numeric_limits::epsilon(); @@ -378,7 +381,6 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { // setup input auto const delta_CV = 0.5; auto const filename_input = std::string("wl_input.dat"); - auto const filename_output = std::string("wl_output.dat"); std::ofstream wl_file; wl_file.open(filename_input); wl_file << "# header 1.0 0.5\n"; @@ -401,13 +403,10 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { BOOST_CHECK_CLOSE(r_algo.max_boundaries_energies[1], 5., tol); // check Wang-Landau histogram { - r_algo.write_wang_landau_results_to_file(filename_output); - std::ifstream outfile; - outfile.open(filename_output); - BOOST_TEST_REQUIRE(!outfile.fail(), "output file not generated"); + std::basic_stringstream outfile; + r_algo.format_wang_landau_results(outfile); std::istream_iterator start(outfile), end; std::vector flat_array(start, end); - outfile.close(); BOOST_REQUIRE_EQUAL(flat_array.size(), 14ul); std::vector bin_edges; std::vector histogram; From ebb8ab76d16471fbb9c65c061c905b2ddcb9c636 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Mar 2021 15:58:17 +0100 Subject: [PATCH 19/23] tests: Simplify code --- .../reaction_ensemble_classes_test.cpp | 40 ++++++++++--------- src/python/espressomd/reaction_ensemble.pxd | 2 +- src/python/espressomd/reaction_ensemble.pyx | 26 ++++++------ testsuite/python/reaction_ensemble.py | 2 - 4 files changed, 36 insertions(+), 34 deletions(-) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index 6bfc40692f7..81e2c89c673 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -204,32 +204,32 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { // exception if no reaction was added BOOST_CHECK_THROW(r_algo.check_reaction_ensemble(), 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 reaction addition { - // 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 ref(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); - r_algo.add_reaction(ref.gamma, ref.reactant_types, - ref.reactant_coefficients, ref.product_types, - ref.product_coefficients); + r_algo.add_reaction(reaction.gamma, reaction.reactant_types, + reaction.reactant_coefficients, reaction.product_types, + reaction.product_coefficients); BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 1ul); - auto const &reaction = r_algo.reactions[0]; - BOOST_TEST(reaction.reactant_types == ref.reactant_types, + auto const &value = r_algo.reactions[0]; + BOOST_TEST(value.reactant_types == reaction.reactant_types, boost::test_tools::per_element()); - BOOST_TEST(reaction.reactant_coefficients == ref.reactant_coefficients, + BOOST_TEST(value.reactant_coefficients == reaction.reactant_coefficients, boost::test_tools::per_element()); - BOOST_TEST(reaction.product_types == ref.product_types, + BOOST_TEST(value.product_types == reaction.product_types, boost::test_tools::per_element()); - BOOST_TEST(reaction.product_coefficients == ref.product_coefficients, + BOOST_TEST(value.product_coefficients == reaction.product_coefficients, boost::test_tools::per_element()); - BOOST_CHECK_EQUAL(reaction.gamma, ref.gamma); + BOOST_CHECK_EQUAL(value.gamma, reaction.gamma); } // check acceptance probability { - SingleReaction const reaction(2., {0}, {1}, {1, 2}, {3, 4}); double probability = r_algo.calculate_acceptance_probability( reaction, -1., -1., {{1, 2}}, -1, -1, false); BOOST_CHECK_EQUAL(probability, -10.); @@ -256,12 +256,16 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { // check reaction removal { - r_algo.add_reaction(5., {1}, {1}, {2}, {1}); + SingleReaction const new_reaction(5., {type_B}, {1}, {type_C}, {1}); + r_algo.add_reaction(new_reaction.gamma, new_reaction.reactant_types, + new_reaction.reactant_coefficients, + new_reaction.product_types, + new_reaction.product_coefficients); BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 2ul); - BOOST_CHECK_EQUAL(r_algo.reactions[1].gamma, 5.); + BOOST_CHECK_EQUAL(r_algo.reactions[1].gamma, new_reaction.gamma); r_algo.delete_reaction(1); BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 1ul); - BOOST_CHECK_EQUAL(r_algo.reactions[0].gamma, 2.); + BOOST_CHECK_EQUAL(r_algo.reactions[0].gamma, reaction.gamma); r_algo.delete_reaction(0); BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 0ul); } diff --git a/src/python/espressomd/reaction_ensemble.pxd b/src/python/espressomd/reaction_ensemble.pxd index 0cd60388615..12c3bf3046e 100644 --- a/src/python/espressomd/reaction_ensemble.pxd +++ b/src/python/espressomd/reaction_ensemble.pxd @@ -41,7 +41,7 @@ cdef extern from "reaction_ensemble.hpp" namespace "ReactionEnsemble": int check_reaction_ensemble() except + double get_acceptance_rate_configurational_moves() int delete_particle(int p_id) - void add_reaction(double gamma, vector[int] _reactant_types, vector[int] _reactant_coefficients, vector[int] _product_types, vector[int] _product_coefficients) except + + void add_reaction(double gamma, vector[int] reactant_types, vector[int] reactant_coefficients, vector[int] product_types, vector[int] product_coefficients) except + void delete_reaction(int reaction_id) vector[SingleReaction] reactions diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index 474a03b660a..dd779648ce6 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -223,27 +223,27 @@ cdef class ReactionAlgorithm: def _set_params_in_es_core_add(self): cdef vector[int] reactant_types - for i in range(len(self._params["reactant_types"])): - reactant_types.push_back(self._params["reactant_types"][i]) cdef vector[int] reactant_coefficients - for i in range(len(self._params["reactant_coefficients"])): - reactant_coefficients.push_back( - self._params["reactant_coefficients"][i]) cdef vector[int] product_types - for i in range(len(self._params["product_types"])): - product_types.push_back(self._params["product_types"][i]) cdef vector[int] product_coefficients - for i in range(len(self._params["product_coefficients"])): - product_coefficients.push_back( - self._params["product_coefficients"][i]) + for value in self._params["reactant_types"]: + reactant_types.push_back(value) + for value in self._params["reactant_coefficients"]: + reactant_coefficients.push_back(value) + for value in self._params["product_types"]: + product_types.push_back(value) + for value in self._params["product_coefficients"]: + product_coefficients.push_back(value) + + # forward reaction deref(self.RE).add_reaction( self._params["gamma"], reactant_types, reactant_coefficients, product_types, product_coefficients) + # backward reaction deref(self.RE).add_reaction( 1.0 / self._params["gamma"], product_types, product_coefficients, reactant_types, reactant_coefficients) - for key in self._params["default_charges"]: # the keys are the types - deref(self.RE).charges_of_types[ - int(key)] = self._params["default_charges"][key] + for key, value in self._params["default_charges"].items(): + deref(self.RE).charges_of_types[int(key)] = value deref(self.RE).check_reaction_ensemble() def _validate_params_default_charge(self): diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py index 1a3c0e9eda7..ff5e68cd460 100644 --- a/testsuite/python/reaction_ensemble.py +++ b/testsuite/python/reaction_ensemble.py @@ -109,7 +109,6 @@ def test_ideal_titration_curve(self): average_NA /= num_samples average_NHA /= num_samples average_alpha = average_NA / float(N0) - print(average_alpha) # Note: with 40 particles, alpha=0.5 and 1000*10 reactions, standard # deviation of average alpha is about 0.003 (determined from 40 # repeated simulations). We set the desired accuracy to 5*std = 0.015 @@ -176,7 +175,6 @@ def test_change_reaction_constant(self): RE_status = RE.get_status() forward_reaction = RE_status["reactions"][0] backward_reaction = RE_status["reactions"][1] - print(forward_reaction) self.assertEqual( new_reaction_constant, forward_reaction["gamma"], From ead37ef002dc49433b8580d1335ba65e6f247dfa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Mar 2021 13:37:54 +0100 Subject: [PATCH 20/23] core: Simplify file I/O operations --- src/core/reaction_ensemble.cpp | 81 +++++++++++---------- src/core/reaction_ensemble.hpp | 9 ++- src/python/espressomd/reaction_ensemble.pxd | 8 +- 3 files changed, 53 insertions(+), 45 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 09cf107cb64..6ce05c42ddf 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -38,10 +38,12 @@ #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -1231,19 +1233,18 @@ void WangLandauReactionEnsemble::format_wang_landau_results(std::ostream &out) { out << wang_landau_potential[flattened_index] << " \n"; } } - out.flush(); } /** *Writes the Wang-Landau potential to file. */ void WangLandauReactionEnsemble::write_wang_landau_results_to_file( - const std::string &full_path_to_output_filename) { + const std::string &filename) { std::ofstream outfile; - outfile.open(full_path_to_output_filename); + outfile.open(filename); if (!outfile.is_open()) { - throw std::runtime_error("Wang-Landau file could not be opened\n"); + throw std::runtime_error("Cannot write to " + filename); } format_wang_landau_results(outfile); @@ -1287,14 +1288,14 @@ int WangLandauReactionEnsemble:: *preliminary energy reweighting run. */ void WangLandauReactionEnsemble::write_out_preliminary_energy_run_results( - const std::string &full_path_to_output_filename) { - FILE *pFile; - pFile = fopen(full_path_to_output_filename.c_str(), "w"); - if (pFile == nullptr) { - throw std::runtime_error("ERROR: Wang-Landau file could not be written\n"); - } - fprintf(pFile, "#nbar E_min E_max\n"); + const std::string &filename) { + std::ofstream outfile; + outfile.open(filename); + if (!outfile.is_open()) { + throw std::runtime_error("Cannot write to " + filename); + } + outfile << "#nbar E_min E_max\n"; for (std::size_t flattened_index = 0; flattened_index < wang_landau_potential.size(); flattened_index++) { // unravel index @@ -1308,13 +1309,12 @@ void WangLandauReactionEnsemble::write_out_preliminary_energy_run_results( auto const value = static_cast(unraveled_index[i]) * collective_variables[i]->delta_CV + collective_variables[i]->CV_minimum; - fprintf(pFile, "%f ", value); + outfile << value << " "; } - fprintf(pFile, "%f %f \n", minimum_energies_at_flat_index[flattened_index], - maximum_energies_at_flat_index[flattened_index]); + outfile << minimum_energies_at_flat_index[flattened_index] << " " + << maximum_energies_at_flat_index[flattened_index] << " \n"; } - fflush(pFile); - fclose(pFile); + outfile.close(); } /** @@ -1376,43 +1376,56 @@ int WangLandauReactionEnsemble:: * Additionally you should store the positions of the particles. * Not storing them introduces small, small statistical errors. */ -int WangLandauReactionEnsemble::write_wang_landau_checkpoint( +void WangLandauReactionEnsemble::write_wang_landau_checkpoint( const std::string &identifier) { + std::string filename; std::ofstream outfile; // write current Wang-Landau parameters (wang_landau_parameter, // monte_carlo_trial_moves, flat_index_of_current_state) - outfile.open(std::string("checkpoint_wang_landau_parameters_") + identifier); + filename = std::string("checkpoint_wang_landau_parameters_") + identifier; + outfile.open(filename); + if (!outfile.is_open()) { + throw std::runtime_error("Cannot write to " + filename); + } outfile << wang_landau_parameter << " " << monte_carlo_trial_moves << " " << get_flattened_index_wang_landau_of_current_state() << "\n"; outfile.close(); // write histogram - outfile.open(std::string("checkpoint_wang_landau_histogram_") + identifier); + filename = std::string("checkpoint_wang_landau_histogram_") + identifier; + outfile.open(filename); + if (!outfile.is_open()) { + throw std::runtime_error("Cannot write to " + filename); + } for (int i = 0; i < wang_landau_potential.size(); i++) { outfile << histogram[i] << "\n"; } outfile.close(); // write Wang-Landau potential - outfile.open(std::string("checkpoint_wang_landau_potential_") + identifier); + filename = std::string("checkpoint_wang_landau_potential_") + identifier; + outfile.open(filename); + if (!outfile.is_open()) { + throw std::runtime_error("Cannot write to " + filename); + } for (double i : wang_landau_potential) { outfile << i << "\n"; } outfile.close(); - return 0; } /** *Loads the Wang-Landau checkpoint */ -int WangLandauReactionEnsemble::load_wang_landau_checkpoint( +void WangLandauReactionEnsemble::load_wang_landau_checkpoint( const std::string &identifier) { + std::string filename; std::ifstream infile; // restore Wang-Landau parameters - infile.open(std::string("checkpoint_wang_landau_parameters_") + identifier); + filename = std::string("checkpoint_wang_landau_parameters_") + identifier; + infile.open(filename); if (infile.is_open()) { - double wang_landau_parameter_entry; int wang_landau_monte_carlo_trial_moves_entry; int flat_index_of_state_at_checkpointing; @@ -1426,13 +1439,12 @@ int WangLandauReactionEnsemble::load_wang_landau_checkpoint( } infile.close(); } else { - throw std::runtime_error("Exception opening" + - std::string("checkpoint_wang_landau_parameters_") + - identifier); + throw std::runtime_error("Cannot read " + filename); } // restore histogram - infile.open(std::string("checkpoint_wang_landau_histogram_") + identifier); + filename = std::string("checkpoint_wang_landau_histogram_") + identifier; + infile.open(filename); if (infile.is_open()) { int hist_entry; int line = 0; @@ -1442,13 +1454,12 @@ int WangLandauReactionEnsemble::load_wang_landau_checkpoint( } infile.close(); } else { - throw std::runtime_error("Exception opening/ reading " + - std::string("checkpoint_wang_landau_histogram_") + - identifier); + throw std::runtime_error("Cannot read " + filename); } // restore Wang-Landau potential - infile.open(std::string("checkpoint_wang_landau_potential_") + identifier); + filename = std::string("checkpoint_wang_landau_potential_") + identifier; + infile.open(filename); if (infile.is_open()) { double wang_landau_potential_entry; int line = 0; @@ -1458,16 +1469,12 @@ int WangLandauReactionEnsemble::load_wang_landau_checkpoint( } infile.close(); } else { - throw std::runtime_error("Exception opening " + - std::string("checkpoint_wang_landau_potential_") + - identifier); + throw std::runtime_error("Cannot read " + filename); } // possible task: restore state in which the system was when the checkpoint // was written. However as long as checkpointing and restoring the system form // the checkpoint is rare this should not matter statistically. - - return 0; } int ConstantpHEnsemble::get_random_valid_p_id() { diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 97209d64bc9..4a33c40c999 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -27,10 +27,12 @@ #include #include +#include #include #include #include #include +#include #include #include #include @@ -312,10 +314,9 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { // checkpointing, only designed to reassign values of a previous simulation to // a new simulation with the same initialization process - int load_wang_landau_checkpoint(const std::string &identifier); - int write_wang_landau_checkpoint(const std::string &identifier); - void write_wang_landau_results_to_file( - const std::string &full_path_to_output_filename); + void load_wang_landau_checkpoint(const std::string &identifier); + void write_wang_landau_checkpoint(const std::string &identifier); + void write_wang_landau_results_to_file(const std::string &filename); protected: double calculate_acceptance_probability( diff --git a/src/python/espressomd/reaction_ensemble.pxd b/src/python/espressomd/reaction_ensemble.pxd index 12c3bf3046e..2c6b94ff1e4 100644 --- a/src/python/espressomd/reaction_ensemble.pxd +++ b/src/python/espressomd/reaction_ensemble.pxd @@ -72,10 +72,10 @@ cdef extern from "reaction_ensemble.hpp" namespace "ReactionEnsemble": void add_new_CV_degree_of_association(int associated_type, double CV_minimum, double CV_maximum, vector[int] corresponding_acid_types) void add_new_CV_potential_energy(string filename, double delta_CV) int update_maximum_and_minimum_energies_at_current_state() - void write_out_preliminary_energy_run_results(string filename) - int write_wang_landau_checkpoint(string identifier) - int load_wang_landau_checkpoint(string identifier) - void write_wang_landau_results_to_file(string full_path_to_output_filename) + void write_out_preliminary_energy_run_results(string filename) except + + void write_wang_landau_checkpoint(string identifier) except + + void load_wang_landau_checkpoint(string identifier) except + + void write_wang_landau_results_to_file(string filename) except + cdef cppclass CConstantpHEnsemble "ReactionEnsemble::ConstantpHEnsemble"(CReactionAlgorithm): CConstantpHEnsemble(int seed) From 28eea677aae1fa4cbf85921c0435a374ce7325af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Mar 2021 14:02:29 +0100 Subject: [PATCH 21/23] core: Refactor Wang-Landau Collective Variable --- src/core/reaction_ensemble.cpp | 38 +++++++++++-------- src/core/reaction_ensemble.hpp | 5 ++- .../reaction_ensemble_classes_test.cpp | 11 ++---- src/python/espressomd/reaction_ensemble.pxd | 2 +- 4 files changed, 31 insertions(+), 25 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 6ce05c42ddf..6d4377e9861 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -34,6 +34,7 @@ #include #include +#include #include #include #include @@ -51,19 +52,13 @@ namespace ReactionEnsemble { -/** Save minimum and maximum energies as a function of the other collective - * variables under min_boundaries_energies, max_boundaries_energies +/** Load minimum and maximum energies as a function of the other collective + * variables. */ void EnergyCollectiveVariable::load_CV_boundaries( - WangLandauReactionEnsemble &wl_system) { + WangLandauReactionEnsemble &wl_system, std::istream &infile) { wl_system.do_energy_reweighting = true; - // load energy boundaries from file - std::ifstream infile; - infile.open(energy_boundaries_filename); - if (!infile.is_open()) - throw std::runtime_error("energy boundaries file for the specific " - "system could not be read.\n"); // Note that you cannot change the other collective variables in the // pre-production run and the production run @@ -81,9 +76,8 @@ void EnergyCollectiveVariable::load_CV_boundaries( values.push_back(value); } assert(values.size() >= 2); - auto it = values.end(); - wl_system.max_boundaries_energies.push_back(*(--it)); - wl_system.min_boundaries_energies.push_back(*(--it)); + wl_system.max_boundaries_energies.emplace_back(values.back()); + wl_system.min_boundaries_energies.emplace_back(values[values.size() - 2]); } CV_minimum = *boost::range::min_element(wl_system.min_boundaries_energies); @@ -805,16 +799,30 @@ void WangLandauReactionEnsemble::add_new_CV_degree_of_association( * Wang-Landau sampling */ void WangLandauReactionEnsemble::add_new_CV_potential_energy( - const std::string &filename, double delta_CV) { + std::istream &infile, double delta_CV) { std::shared_ptr new_collective_variable = std::make_shared(); - new_collective_variable->energy_boundaries_filename = filename; new_collective_variable->delta_CV = delta_CV; - new_collective_variable->load_CV_boundaries(*this); + new_collective_variable->load_CV_boundaries(*this, infile); collective_variables.emplace_back(new_collective_variable); initialize_wang_landau(); } +/** + * Adds a new collective variable (CV) of the type potential energy to the + * Wang-Landau sampling + */ +void WangLandauReactionEnsemble::add_new_CV_potential_energy( + const std::string &filename, double delta_CV) { + std::ifstream infile; + infile.open(filename); + if (!infile.is_open()) { + throw std::runtime_error("Cannot read " + filename); + } + add_new_CV_potential_energy(infile, delta_CV); + infile.close(); +} + /** * Returns the flattened index of the multidimensional Wang-Landau histogram */ diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 4a33c40c999..87ee27f7926 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -94,11 +94,11 @@ struct CollectiveVariable { class WangLandauReactionEnsemble; struct EnergyCollectiveVariable : public CollectiveVariable { - std::string energy_boundaries_filename; double determine_current_state() const override { return calculate_current_potential_energy_of_system(); } - void load_CV_boundaries(WangLandauReactionEnsemble &wl_system); + void load_CV_boundaries(WangLandauReactionEnsemble &wl_system, + std::istream &infile); }; /** Measure the degree of association of a chemical species. @@ -295,6 +295,7 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { const std::vector &_corresponding_acid_types); void add_new_CV_potential_energy(const std::string &filename, double delta_CV); + void add_new_CV_potential_energy(std::istream &infile, double delta_CV); std::vector> collective_variables; std::string output_filename = ""; diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index 81e2c89c673..021f3da6a24 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -40,8 +40,7 @@ #include #include #include -#include -#include +#include #include #include #include @@ -384,16 +383,14 @@ BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { using namespace boost::adaptors; // setup input auto const delta_CV = 0.5; - auto const filename_input = std::string("wl_input.dat"); - std::ofstream wl_file; - wl_file.open(filename_input); + std::basic_stringstream wl_file; wl_file << "# header 1.0 0.5\n"; wl_file << "0 1.0 4.0\n"; wl_file << "1\t2.0\t5.0\n"; - wl_file.close(); + wl_file.flush(); r_algo.do_energy_reweighting = false; // add collective variable - r_algo.add_new_CV_potential_energy(filename_input, delta_CV); + r_algo.add_new_CV_potential_energy(wl_file, delta_CV); BOOST_REQUIRE_EQUAL(r_algo.collective_variables.size(), 1ul); BOOST_REQUIRE_EQUAL(r_algo.min_boundaries_energies.size(), 2ul); BOOST_REQUIRE_EQUAL(r_algo.max_boundaries_energies.size(), 2ul); diff --git a/src/python/espressomd/reaction_ensemble.pxd b/src/python/espressomd/reaction_ensemble.pxd index 2c6b94ff1e4..f044e137392 100644 --- a/src/python/espressomd/reaction_ensemble.pxd +++ b/src/python/espressomd/reaction_ensemble.pxd @@ -70,7 +70,7 @@ cdef extern from "reaction_ensemble.hpp" namespace "ReactionEnsemble": vector[double] maximum_energies_at_flat_index bool do_not_sample_reaction_partition_function void add_new_CV_degree_of_association(int associated_type, double CV_minimum, double CV_maximum, vector[int] corresponding_acid_types) - void add_new_CV_potential_energy(string filename, double delta_CV) + void add_new_CV_potential_energy(string filename, double delta_CV) except + int update_maximum_and_minimum_energies_at_current_state() void write_out_preliminary_energy_run_results(string filename) except + void write_wang_landau_checkpoint(string identifier) except + From 3090ff7cb2df61e8f6c65e8dadded3540337900f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Mar 2021 14:12:56 +0100 Subject: [PATCH 22/23] core: Use std::accumulate --- src/core/reaction_ensemble.cpp | 23 +++++++++++------------ src/core/reaction_ensemble.hpp | 1 - 2 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 6d4377e9861..af602c2bdb8 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -44,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -954,17 +955,6 @@ double WangLandauReactionEnsemble::calculate_delta_degree_of_association( return delta; } -long WangLandauReactionEnsemble::get_num_needed_bins() const { - long needed_bins = 1; - // add 1 for degrees of association-related part of histogram (think of only - // one acid particle) - for (const auto &cv : collective_variables) { - needed_bins *= - static_cast((cv->CV_maximum - cv->CV_minimum) / cv->delta_CV) + 1; - } - return needed_bins; -} - void WangLandauReactionEnsemble::invalidate_bins() { // make values in histogram and Wang-Landau potential negative if they are not // allowed at the given degree of association, because the energy boundaries @@ -1020,7 +1010,16 @@ int WangLandauReactionEnsemble::initialize_wang_landau() { // construct (possibly higher dimensional) histogram and potential over // gamma (the space which should be equally sampled when the Wang-Landau // algorithm has converged) - auto const needed_bins = get_num_needed_bins(); + // add 1 for degrees of association-related part of histogram (think of only + // one acid particle) + auto const needed_bins = std::accumulate( + collective_variables.begin(), collective_variables.end(), 1l, + [](long acc, std::shared_ptr const &cv) { + return acc * (static_cast((cv->CV_maximum - cv->CV_minimum) / + cv->delta_CV) + + 1l); + }); + assert(needed_bins >= 0); histogram.resize(needed_bins, 0); wang_landau_potential.resize(needed_bins, 0); diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 87ee27f7926..ba59e78f991 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -382,7 +382,6 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { int initialize_wang_landau(); double calculate_delta_degree_of_association( DegreeOfAssociationCollectiveVariable ¤t_collective_variable); - long get_num_needed_bins() const; void invalidate_bins(); void reset_histogram(); double get_minimum_CV_value_on_delta_CV_spaced_grid(double min_CV_value, From b3e6d4f5a8ef6d6d834db863c653d252d7923bf9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Mar 2021 14:14:06 +0100 Subject: [PATCH 23/23] tests: Remove WidomInsertion unit test --- .../reaction_ensemble_classes_test.cpp | 21 ------------------- 1 file changed, 21 deletions(-) diff --git a/src/core/unit_tests/reaction_ensemble_classes_test.cpp b/src/core/unit_tests/reaction_ensemble_classes_test.cpp index 021f3da6a24..b2f4d6c3eae 100644 --- a/src/core/unit_tests/reaction_ensemble_classes_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -470,27 +470,6 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { } } -BOOST_AUTO_TEST_CASE(WidomInsertion_test) { - using namespace ReactionEnsemble; - constexpr double tol = 100 * std::numeric_limits::epsilon(); - - // check acceptance rate - WidomInsertion r_algo(42); - - // 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}); - r_algo.add_reaction(reaction.gamma, reaction.reactant_types, - reaction.reactant_coefficients, reaction.product_types, - reaction.product_coefficients); - - // exception if not enough particles - BOOST_CHECK_THROW(r_algo.measure_excess_chemical_potential(0), - std::runtime_error); -} - int main(int argc, char **argv) { auto mpi_env = std::make_shared(argc, argv); Communication::init(mpi_env);