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 3e783fed4b4..af602c2bdb8 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,36 +31,42 @@ #include #include +#include + +#include +#include #include #include #include #include +#include #include #include +#include +#include +#include +#include +#include +#include +#include #include 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 &m_current_wang_landau_system) { - - m_current_wang_landau_system.do_energy_reweighting = true; - // load energy boundaries from file - std::ifstream infile; + WangLandauReactionEnsemble &wl_system, std::istream &infile) { - infile.open( - energy_boundaries_filename); // file containing numbers in 3 columns - if (infile.fail()) - throw std::runtime_error("ERROR: energy boundaries file for the specific " - "system could not be read.\n"); + wl_system.do_energy_reweighting = true; // 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)) { @@ -69,18 +76,13 @@ 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); + 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); + CV_maximum = *boost::range::max_element(wl_system.max_boundaries_energies); } /** @@ -88,7 +90,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; @@ -98,20 +100,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; - - 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); + 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) @@ -128,7 +122,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"); } @@ -200,7 +194,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 = @@ -230,7 +224,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) { @@ -275,8 +269,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++) { @@ -304,7 +298,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) @@ -326,10 +320,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); @@ -364,10 +358,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) { @@ -463,13 +456,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++) { @@ -493,23 +484,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. @@ -629,8 +603,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); @@ -684,7 +657,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) { @@ -827,36 +800,43 @@ 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; - collective_variables.push_back(new_collective_variable); - new_collective_variable->load_CV_boundaries(*this); - collective_variables[collective_variables.size() - 1] = - new_collective_variable; + 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 */ 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 - 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] + @@ -868,26 +848,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++) { @@ -904,7 +886,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++) @@ -942,7 +925,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; @@ -972,50 +955,34 @@ double WangLandauReactionEnsemble::calculate_delta_degree_of_association( return delta; } -/** - * Initializes the Wang-Landau histogram. - */ -int WangLandauReactionEnsemble::get_num_needed_bins() { - int needed_bins = 1; - 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) - } - 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 // 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( + 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( 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 > @@ -1024,37 +991,37 @@ 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; } -/** - * 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; - 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 - - // construct (possibly higher dimensional) histogram over gamma (the room - // 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 + auto const &cv = collective_variables.back(); + // add 1 for collective variables which are of type degree of association + 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) + // 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); used_bins = needed_bins; // initialize for 1/t wang_landau algorithm @@ -1069,15 +1036,15 @@ 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 || - 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 * @@ -1086,9 +1053,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 @@ -1136,7 +1102,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 @@ -1181,7 +1147,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 @@ -1220,8 +1186,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; @@ -1239,7 +1205,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; @@ -1247,44 +1213,49 @@ 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"); - } - for (int flattened_index = 0; flattened_index < wang_landau_potential.size(); - flattened_index++) { +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 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; + out << value << " "; } - fprintf(pFile, "%f \n", wang_landau_potential[flattened_index]); + out << wang_landau_potential[flattened_index] << " \n"; } } - fflush(pFile); - fclose(pFile); +} + +/** + *Writes the Wang-Landau potential to file. + */ +void WangLandauReactionEnsemble::write_wang_landau_results_to_file( + const std::string &filename) { + std::ofstream outfile; + + outfile.open(filename); + if (!outfile.is_open()) { + throw std::runtime_error("Cannot write to " + filename); + } + + format_wang_landau_results(outfile); + outfile.close(); } /** @@ -1324,33 +1295,33 @@ 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) { - for (int flattened_index = 0; flattened_index < wang_landau_potential.size(); - flattened_index++) { + 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 - 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; + 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(); } /** @@ -1361,18 +1332,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; } @@ -1410,43 +1383,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; @@ -1460,13 +1446,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; @@ -1476,13 +1461,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; @@ -1492,16 +1476,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() { @@ -1552,8 +1532,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; @@ -1564,10 +1545,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 * @@ -1610,14 +1591,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; } @@ -1628,13 +1607,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 @@ -1643,7 +1622,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 35c1ad5818c..ba59e78f991 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -27,8 +27,12 @@ #include #include +#include +#include #include #include +#include +#include #include #include #include @@ -38,6 +42,22 @@ 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) { + 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) - + std::accumulate(reactant_coefficients.begin(), + reactant_coefficients.end(), 0); + } + // strict input to the algorithm std::vector reactant_types; std::vector reactant_coefficients; @@ -45,14 +65,14 @@ 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; - double get_acceptance_rate() { + double get_acceptance_rate() const { return static_cast(accepted_moves) / static_cast(tried_moves); - }; + } }; struct StoredParticleProperty { @@ -65,27 +85,34 @@ 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; }; 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 - load_CV_boundaries(WangLandauReactionEnsemble &m_current_wang_landau_system); + void load_CV_boundaries(WangLandauReactionEnsemble &wl_system, + std::istream &infile); }; +/** 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; - double determine_current_state() override { + /** Calculate the degree of association of the reference species */ + double determine_current_state() const override { return calculate_degree_of_association(); } @@ -95,7 +122,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 = @@ -110,8 +137,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; } }; @@ -148,20 +174,20 @@ 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, - 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); } @@ -184,15 +210,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 @@ -204,7 +231,16 @@ 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; + +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; @@ -213,10 +249,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); @@ -224,14 +256,6 @@ class ReactionAlgorithm { void append_particle_property_of_random_particle( 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) { - return -10; - }; - - void add_types_to_index(std::vector &type_list); Utils::Vector3d get_random_position_in_box(); }; @@ -250,12 +274,12 @@ class ReactionEnsemble : public ReactionAlgorithm { public: ReactionEnsemble(int seed) : ReactionAlgorithm(seed) {} -private: +protected: 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 */ @@ -271,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 = ""; @@ -278,25 +303,29 @@ 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); // 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( + 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 format_wang_landau_results(std::ostream &out); private: void on_reaction_entry(int &old_state_index) override; @@ -304,28 +333,26 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { 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 ¤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; 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; 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; - 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( @@ -333,30 +360,32 @@ 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 + /** + * 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(); 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; }; /** @@ -378,12 +407,14 @@ class ConstantpHEnsemble : public ReactionAlgorithm { double m_constant_pH = -10; int do_reaction(int reaction_steps) override; -private: +protected: 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; + +private: int get_random_valid_p_id(); }; @@ -398,8 +429,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); @@ -425,7 +457,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; } 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..b2f4d6c3eae --- /dev/null +++ b/src/core/unit_tests/reaction_ensemble_classes_test.cpp @@ -0,0 +1,478 @@ +/* + * 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 + +#include +#include +#include +#include +#include +#include +#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; +}; + +// 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(); + + // 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); +} + +// 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; + 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); +} + +// 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(); + + // 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 + 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); + } + } + + // 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) { + // 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); + } + } + } +} + +// Check the base class for all Monte Carlo algorithms. +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 + 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; + 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); + + // 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 + { + 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 &value = r_algo.reactions[0]; + BOOST_TEST(value.reactant_types == reaction.reactant_types, + boost::test_tools::per_element()); + BOOST_TEST(value.reactant_coefficients == reaction.reactant_coefficients, + boost::test_tools::per_element()); + BOOST_TEST(value.product_types == reaction.product_types, + boost::test_tools::per_element()); + BOOST_TEST(value.product_coefficients == reaction.product_coefficients, + boost::test_tools::per_element()); + BOOST_CHECK_EQUAL(value.gamma, reaction.gamma); + } + + // check acceptance probability + { + 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); + + // 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 + { + 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, new_reaction.gamma); + r_algo.delete_reaction(1); + BOOST_REQUIRE_EQUAL(r_algo.reactions.size(), 1ul); + BOOST_CHECK_EQUAL(r_algo.reactions[0].gamma, reaction.gamma); + r_algo.delete_reaction(0); + 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: + 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); + + // 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; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // system contains i x A, j x B, and k x C + auto const p_numbers = + std::map{{type_A, i}, {type_B, j}, {type_C, k}}; + auto const energy = static_cast(i + 1); + auto const f_expr = calculate_factorial_expression(reaction, p_numbers); + 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(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 { + public: + using WangLandauReactionEnsemble::calculate_acceptance_probability; + using WangLandauReactionEnsemble::format_wang_landau_results; + 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); + + // 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; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // system contains i x A, j x B, and k x C + auto const p_numbers = + std::map{{type_A, i}, {type_B, j}, {type_C, k}}; + auto const energy = static_cast(i + 1); + auto const f_expr = calculate_factorial_expression(reaction, p_numbers); + auto const prefactor = 2e6 * f_expr; + auto const boltzmann = std::exp(energy / 20.); + double acceptance; + r_algo.do_energy_reweighting = false; + 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; + 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); + } + } + } + + // check collective variables + { + using namespace boost::adaptors; + // setup input + auto const delta_CV = 0.5; + 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.flush(); + r_algo.do_energy_reweighting = false; + // add collective variable + 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); + 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 + { + std::basic_stringstream outfile; + r_algo.format_wang_landau_results(outfile); + std::istream_iterator start(outfile), end; + std::vector flat_array(start, end); + 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()); + } + // exception if file doesn't exist + BOOST_CHECK_THROW(r_algo.add_new_CV_potential_energy("unknown", 0.), + std::runtime_error); + } +} + +// 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 { + 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); + + // 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; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + // system contains i x A, j x B, and k x C + auto const p_numbers = + std::map{{type_A, i}, {type_B, j}, {type_C, k}}; + auto const energy = static_cast(i + 1); + auto const acceptance_ref = + std::exp(energy / 20. + std::log(10.) * (1. + std::log10(2.))); + auto const acceptance = r_algo.calculate_acceptance_probability( + reaction, energy, 0., p_numbers, -1, -1, false); + BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); + } + } + } +} + +int main(int argc, char **argv) { + auto mpi_env = std::make_shared(argc, argv); + Communication::init(mpi_env); + + return boost::unit_test::unit_test_main(init_unit_test, argc, argv); +} diff --git a/src/core/unit_tests/reaction_ensemble_utils_test.cpp b/src/core/unit_tests/reaction_ensemble_utils_test.cpp index 3f78f7691ef..a234bfab15c 100644 --- a/src/core/unit_tests/reaction_ensemble_utils_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_utils_test.cpp @@ -17,12 +17,15 @@ * 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 #include "reaction_ensemble.hpp" +#include #include #include #include @@ -57,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); + } + } } diff --git a/src/python/espressomd/reaction_ensemble.pxd b/src/python/espressomd/reaction_ensemble.pxd index 0cd60388615..f044e137392 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 @@ -70,12 +70,12 @@ 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) - 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) 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/CMakeLists.txt b/testsuite/python/CMakeLists.txt index f13ae95562d..a94ba79bc5d 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -119,7 +119,8 @@ 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.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) python_test(FILE lb_pressure_tensor.py MAX_NUM_PROC 1 LABELS gpu long) @@ -171,7 +172,8 @@ 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.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) python_test(FILE observable_profile.py MAX_NUM_PROC 4) diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py index 878f5f83141..1f7a67c02ab 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH.py @@ -1,5 +1,5 @@ # -# Copyright (C) 2013-2019 The ESPResSo project +# Copyright (C) 2013-2021 The ESPResSo project # # This file is part of ESPResSo. # @@ -17,105 +17,68 @@ # 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): +class ConstantpHTest(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 - type_HA = 0 - 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 - pH = 2 - pKa = pKa_minus_pH + pH - Ka = 10**(-pKa) - box_l = (N0 / c0)**(1.0 / 3.0) - system = espressomd.System(box_l=[box_l, box_l, box_l]) - np.random.seed(69) # make reaction code fully deterministic - system.cell_system.skin = 0.4 - system.time_step = 0.01 - RE = reaction_ensemble.ConstantpHEnsemble( - temperature=1.0, exclusion_radius=1, seed=44) + np.random.seed(42) - @classmethod - def setUpClass(cls): - for i in range(0, 2 * cls.N0, 2): - cls.system.part.add(id=i, pos=np.random.random(3) * - cls.system.box_l, type=cls.type_A) - cls.system.part.add(id=i + 1, pos=np.random.random(3) * - cls.system.box_l, type=cls.type_H) + 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]) - cls.RE.add_reaction( - gamma=cls.Ka, - reactant_types=[cls.type_HA], + 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=[cls.type_A, cls.type_H], + product_types=[type_A, type_H], product_coefficients=[1, 1], - default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1}) - cls.RE.constant_pH = cls.pH + default_charges={type_HA: 0, type_A: -1, type_H: +1}) + RE.constant_pH = pH - @classmethod - def ideal_alpha(cls, pH): - return 1.0 / (1 + 10**(cls.pKa - pH)) + # equilibration + RE.reaction(800) - def test_ideal_titration_curve(self): - N0 = ReactionEnsembleTest.N0 - type_A = ReactionEnsembleTest.type_A - type_H = ReactionEnsembleTest.type_H - type_HA = ReactionEnsembleTest.type_HA - system = ReactionEnsembleTest.system - RE = ReactionEnsembleTest.RE - # chemical warmup - get close to chemical equilibrium before we start # sampling - RE.reaction(40 * N0) + 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) - average_NH = 0.0 - average_NHA = 0.0 - average_NA = 0.0 - num_samples = 1000 - for _ in range(num_samples): - RE.reaction(10) - average_NH += system.number_of_particles(type=type_H) - average_NHA += system.number_of_particles(type=type_HA) - average_NA += system.number_of_particles(type=type_A) - average_NH /= float(num_samples) - average_NA /= float(num_samples) - average_NHA /= float(num_samples) - average_alpha = average_NA / float(N0) - # note you cannot calculate the pH via -log10(/volume) in the + # 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 - pH = ReactionEnsembleTest.pH - pKa = ReactionEnsembleTest.pKa - target_alpha = ReactionEnsembleTest.ideal_alpha(pH) - rel_error_alpha = abs(average_alpha - target_alpha) / target_alpha - # relative error - self.assertLess( - rel_error_alpha, - 0.015, - msg="\nDeviation from ideal titration curve is too big for the given input parameters.\n" - + f" pH: {pH:.2f}" - + f" pKa: {pKa:.2f}" - + f" average NH: {average_NH:.1f}" - + f" average NA: {average_NA:.1f}" - + f" average NHA: {average_NHA:.1f}" - + f" average alpha: {average_alpha:.3f}" - + f" target alpha: {target_alpha:.3f}" - + f" rel_error: {rel_error_alpha * 100:.1f}%" - ) + 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__": diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py new file mode 100644 index 00000000000..5ed6c79bf26 --- /dev/null +++ b/testsuite/python/constant_pH_stats.py @@ -0,0 +1,119 @@ +# +# 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 . +# + +import unittest as ut +import numpy as np +import espressomd +import espressomd.reaction_ensemble + + +class ReactionEnsembleTest(ut.TestCase): + + """Test the core implementation of the constant pH reaction ensemble.""" + + N0 = 40 + c0 = 0.00028 + type_HA = 0 + type_A = 1 + type_H = 5 + temperature = 1.0 + # choose target alpha not too far from 0.5 to get good statistics in a + # small number of steps + pKa_minus_pH = -0.2 + pH = 2 + pKa = pKa_minus_pH + pH + Ka = 10**(-pKa) + box_l = (N0 / c0)**(1.0 / 3.0) + system = espressomd.System(box_l=[box_l, box_l, box_l]) + np.random.seed(69) # make reaction code fully deterministic + system.cell_system.skin = 0.4 + system.time_step = 0.01 + RE = espressomd.reaction_ensemble.ConstantpHEnsemble( + temperature=1.0, exclusion_radius=1, seed=44) + + @classmethod + def setUpClass(cls): + for i in range(0, 2 * cls.N0, 2): + cls.system.part.add(id=i, pos=np.random.random(3) * + cls.system.box_l, type=cls.type_A) + cls.system.part.add(id=i + 1, pos=np.random.random(3) * + cls.system.box_l, type=cls.type_H) + + cls.RE.add_reaction( + gamma=cls.Ka, + reactant_types=[cls.type_HA], + reactant_coefficients=[1], + product_types=[cls.type_A, cls.type_H], + product_coefficients=[1, 1], + default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1}) + cls.RE.constant_pH = cls.pH + + @classmethod + def ideal_alpha(cls, pH): + return 1.0 / (1 + 10**(cls.pKa - pH)) + + def test_ideal_titration_curve(self): + N0 = ReactionEnsembleTest.N0 + type_A = ReactionEnsembleTest.type_A + type_H = ReactionEnsembleTest.type_H + type_HA = ReactionEnsembleTest.type_HA + system = ReactionEnsembleTest.system + RE = ReactionEnsembleTest.RE + # chemical warmup - get close to chemical equilibrium before we start + # sampling + RE.reaction(40 * N0) + + average_NH = 0.0 + average_NHA = 0.0 + average_NA = 0.0 + num_samples = 1000 + for _ in range(num_samples): + RE.reaction(10) + average_NH += system.number_of_particles(type=type_H) + average_NHA += system.number_of_particles(type=type_HA) + average_NA += system.number_of_particles(type=type_A) + average_NH /= float(num_samples) + average_NA /= float(num_samples) + average_NHA /= float(num_samples) + average_alpha = average_NA / float(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 + pH = ReactionEnsembleTest.pH + pKa = ReactionEnsembleTest.pKa + target_alpha = ReactionEnsembleTest.ideal_alpha(pH) + rel_error_alpha = abs(average_alpha - target_alpha) / target_alpha + # relative error + self.assertLess( + rel_error_alpha, + 0.015, + msg="\nDeviation from ideal titration curve is too big for the given input parameters.\n" + + f" pH: {pH:.2f}" + + f" pKa: {pKa:.2f}" + + f" average NH: {average_NH:.1f}" + + f" average NA: {average_NA:.1f}" + + f" average NHA: {average_NHA:.1f}" + + f" average alpha: {average_alpha:.3f}" + + f" target alpha: {target_alpha:.3f}" + + f" rel_error: {rel_error_alpha * 100:.1f}%" + ) + + +if __name__ == "__main__": + ut.main() 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"], diff --git a/testsuite/python/wang_landau_reaction_ensemble.py b/testsuite/python/wang_landau.py similarity index 61% rename from testsuite/python/wang_landau_reaction_ensemble.py rename to testsuite/python/wang_landau.py index f534517230e..a3e85b335ac 100644 --- a/testsuite/python/wang_landau_reaction_ensemble.py +++ b/testsuite/python/wang_landau.py @@ -22,19 +22,13 @@ 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. - - 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 - """ + """Test the interface of the Wang-Landau reaction ensemble.""" # System parameters # @@ -53,18 +47,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], @@ -87,7 +80,7 @@ class ReactionEnsembleTest(ut.TestCase): do_not_sample_reaction_partition_function=True, full_path_to_output_filename=file_output) - def test_wang_landau_energy_recording(self): + 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( @@ -96,44 +89,7 @@ def test_wang_landau_energy_recording(self): 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): - self.WLRE.add_collective_variable_potential_energy( - filename=self.file_input, delta=0.05) - - # run MC until convergence - while True: - try: - self.WLRE.reaction() - for _ in range(2): - self.WLRE.displacement_mc_move_for_particles_of_type(3) - except reaction_ensemble.WangLandauHasConverged: - break - - nbars, Epots, WL_potentials = np.loadtxt(self.file_output, unpack=True) - mask_nbar_0 = np.where(np.abs(nbars - 1.0) < 0.0001) - Epots = Epots[mask_nbar_0][1:] - WL_potentials = WL_potentials[mask_nbar_0][1:] - - def calc_from_partition_function(quantity): - probability = np.exp(WL_potentials - Epots / self.temperature) - return np.sum(quantity * probability) / np.sum(probability) - - # calculate the canonical potential energy - pot_energy = calc_from_partition_function(Epots) - # calculate the canonical configurational heat capacity - pot_energy_sq = calc_from_partition_function(Epots**2) - heat_capacity = pot_energy_sq - pot_energy**2 - - # for the calculation regarding the analytical results which are - # compared here, see Master Thesis Jonas Landsgesell p. 72 - self.assertAlmostEqual( - pot_energy, 1.5, places=1, - msg="difference to analytical expected canonical potential energy too big") - self.assertAlmostEqual( - heat_capacity, 1.5, places=1, - msg="difference to analytical expected canonical configurational heat capacity too big") - - def _wang_landau_output_checkpoint(self, filename): + def check_checkpoint(self, filename): # write first checkpoint self.WLRE.write_wang_landau_checkpoint() old_checkpoint = np.loadtxt(filename) @@ -151,11 +107,23 @@ def _wang_landau_output_checkpoint(self, filename): new_checkpoint = np.loadtxt(filename) np.testing.assert_almost_equal(new_checkpoint, modified_checkpoint) - def test_wang_landau_output_checkpoint(self): + 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._wang_landau_output_checkpoint(filename) + self.check_checkpoint(filename) if __name__ == "__main__": diff --git a/testsuite/python/wang_landau_stats.py b/testsuite/python/wang_landau_stats.py new file mode 100644 index 00000000000..68ad1c62489 --- /dev/null +++ b/testsuite/python/wang_landau_stats.py @@ -0,0 +1,127 @@ +# +# 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 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 + the density of states is known as the one of the harmonic oscillator. + """ + + # 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_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 + 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_potential_and_heat_capacity(self): + self.WLRE.add_collective_variable_potential_energy( + filename=self.file_input, delta=0.05) + + # run MC until convergence + while True: + try: + self.WLRE.reaction() + for _ in range(2): + self.WLRE.displacement_mc_move_for_particles_of_type(3) + except espressomd.reaction_ensemble.WangLandauHasConverged: + break + + nbars, Epots, WL_potentials = np.loadtxt(self.file_output, unpack=True) + mask_nbar_0 = np.where(np.abs(nbars - 1.0) < 0.0001) + Epots = Epots[mask_nbar_0][1:] + WL_potentials = WL_potentials[mask_nbar_0][1:] + + def calc_from_partition_function(quantity): + probability = np.exp(WL_potentials - Epots / self.temperature) + return np.sum(quantity * probability) / np.sum(probability) + + # calculate the canonical potential energy + pot_energy = calc_from_partition_function(Epots) + # calculate the canonical configurational heat capacity + pot_energy_sq = calc_from_partition_function(Epots**2) + heat_capacity = pot_energy_sq - pot_energy**2 + + # for the calculation regarding the analytical results which are + # compared here, see Master Thesis Jonas Landsgesell p. 72 + self.assertAlmostEqual( + pot_energy, 1.5, places=1, + msg="difference to analytical expected canonical potential energy too big") + self.assertAlmostEqual( + heat_capacity, 1.5, places=1, + msg="difference to analytical expected canonical configurational heat capacity too big") + + +if __name__ == "__main__": + ut.main()