diff --git a/doc/doxygen/bibliography.bib b/doc/doxygen/bibliography.bib index fc47479ac01..acb7f11187d 100644 --- a/doc/doxygen/bibliography.bib +++ b/doc/doxygen/bibliography.bib @@ -589,17 +589,6 @@ @Article{tironi95a doi = {10.1063/1.469273}, } -@Article{troester05a, - title = {{Wang-Landau sampling with self-adaptive range}}, - author = {Tr\"{o}ster, Andreas and Dellago, Christoph}, - journal = {Physical Review E}, - volume = {71}, - number = {6}, - pages = {066705}, - year = {2005}, - doi = {10.1103/PhysRevE.71.066705}, -} - @Article{tyagi10a, author = {Tyagi, Sandeep and S\"{u}zen, Mehmet and Sega, Marcello and Barbosa, Marcia C. and Kantorovich, Sofia S. and Holm, Christian}, title = {{An iterative, fast, linear-scaling method for computing induced charges on arbitrary dielectric boundaries}}, @@ -630,17 +619,6 @@ @Article{wang01a doi = {10.1063/1.1398588}, } -@Article{yan02b, -author = {Yan, Qiliang and Faller, Roland and {de Pablo}, Juan J.}, -title = {{Density-of-states Monte Carlo method for simulation of fluids}}, -journal = {The Journal of Chemical Physics}, -volume = {116}, -number = {20}, -pages = {8745-8749}, -year = {2002}, -doi = {10.1063/1.1463055}, -} - @Article{yeh99a, title = {{Ewald summation for systems with slab geometry}}, author = {Yeh, In-Chul and Berkowitz, Max L.}, diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index 62883d11dc4..7ba73f7ca74 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -1770,27 +1770,6 @@ be converted to :math:`K_c` as where :math:`p^{\ominus}=1\,\mathrm{atm}` is the standard pressure. Consider using the python module pint for unit conversion. -.. _Wang-Landau Reaction Ensemble: - -Wang-Landau Reaction Ensemble -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -Combination of the Reaction Ensemble with the Wang-Landau algorithm -:cite:`wang01a`. Allows for enhanced sampling of the reacting system -and for the determination of the density of states with respect -to the reaction coordinate or with respect to some other collective -variable :cite:`landsgesell17a`. Here the 1/t Wang-Landau -algorithm :cite:`belardinelli07a` is implemented since it -does not suffer from systematic errors. - -Multiple reactions and multiple collective variables can be set. - -An example script can be found here: - -* `Wang-Landau reaction ensemble `__ - -For a description of the available methods, see :class:`espressomd.reaction_ensemble.ReactionEnsemble`. - .. _Grand canonical ensemble simulation using the Reaction Ensemble: Grand canonical ensemble simulation diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib index da7f7bc2ac0..8a2cf684b95 100644 --- a/doc/sphinx/zrefs.bib +++ b/doc/sphinx/zrefs.bib @@ -101,18 +101,6 @@ @article{beenakker86a doi = {10.1063/1.451199} } -@article{belardinelli07a, - title={Fast algorithm to calculate density of states}, - author={Belardinelli, R. E. and Pereyra, V. D.}, - journal={Phys. Rev. E}, - volume={75}, - number={4}, - pages={046701}, - year={2007}, - doi={10.1103/PhysRevE.75.046701}, - publisher={APS} -} - @ARTICLE{bereau15, author = {Tristan Bereau}, title = {Multi-timestep integrator for the modified Andersen barostat}, @@ -513,17 +501,6 @@ @Book{kruger12a doi={10.1007/978-3-8348-2376-2}, } -@Article{landsgesell17a, -title = {{W}ang-{L}andau Reaction Ensemble Method: {S}imulation of Weak Polyelectrolytes and General Acid-Base Reactions}, -author = {Landsgesell, Jonas and Holm, Christian and Smiatek, Jens}, -journal = {J. Chem. Theory Comput.}, -volume = {13}, -number = {2}, -pages = {852--862}, -year = {2017}, -doi = {10.1021/acs.jctc.6b00791} -} - @Article{landsgesell17b, title = {Simulation of weak polyelectrolytes: {A} comparison between the constant p{H} and the reaction ensemble method}, author = {Landsgesell, Jonas and Holm, Christian and Smiatek, Jens}, @@ -929,17 +906,6 @@ @article{wagner02 doi = {10.1023/A:1014595628808} } -@article{wang01a, - title={Efficient, multiple-range random walk algorithm to calculate the density of states}, - author={Wang, Fugao and Landau, David P}, - journal={Phys. Rev. Lett.}, - volume={86}, - number={10}, - pages={2050}, - year={2001}, - publisher={APS} -} - @ARTICLE{wolff04a, author = {Ulli Wolff}, title = {{M}onte {C}arlo errors with less errors}, diff --git a/samples/wang_landau_reaction_ensemble.py b/samples/wang_landau_reaction_ensemble.py deleted file mode 100644 index ba597cf64d4..00000000000 --- a/samples/wang_landau_reaction_ensemble.py +++ /dev/null @@ -1,118 +0,0 @@ -# -# 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 . -# -""" -Simulate two reacting monomers which are bonded via a harmonic potential. -The script aborts as soon as the abortion criterion in the Wang-Landau -algorithm is met: the Wang-Landau simulation runs until the Wang-Landau -potential has converged and then raises a Warning that it has converged, -effectively aborting the simulation. - -With the setup of the Wang-Landau algorithm in this script, you sample the -density of states of a three-dimensional reacting harmonic oscillator as -a function of the two collective variables 1) degree of association and -2) potential energy. - -The recorded Wang-Landau potential (which is updated during the simulation) -is written to the file :file:`WL_potential_out.dat`. - -In this simulation setup the Wang-Landau potential is the density of states. -You can view the converged Wang-Landau potential e.g. via plotting with -gnuplot: ``splot "WL_potential_out.dat"``. As expected the three-dimensional -harmonic oscillator has a density of states which goes like -:math:`\\sqrt{E_{\\text{pot}}}`. - -For a scientific description and different ways to use the algorithm please -consult https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00791 -""" -import numpy as np - -import espressomd -from espressomd import reaction_ensemble -from espressomd.interactions import HarmonicBond - -# System parameters -############################################################# -box_l = 6 * np.sqrt(2) - -# Integration parameters -############################################################# -system = espressomd.System(box_l=[box_l] * 3) -np.random.seed(seed=42) -system.time_step = 0.02 -system.cell_system.skin = 0.4 - - -############################################################# -# Setup System # -############################################################# - - -# Particle setup -############################################################# -# type 0 = HA -# type 1 = A- -# type 2 = H+ - -N0 = 1 # number of titratable units -K_diss = 0.0088 - -inert_monomer = system.part.add(pos=[0, 0, 0] * system.box_l, type=3) -reactive_monomer = system.part.add( - pos=[1.0, 1.0, 1.0] * system.box_l / 2.0, type=1) -proton_1 = system.part.add(pos=np.random.random() * system.box_l, type=2) -proton_2 = system.part.add(pos=np.random.random() * system.box_l, type=2) - -# create a harmonic bond between the two reacting particles => 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 -h = HarmonicBond(r_0=0, k=1) -system.bonded_inter[0] = h -inert_monomer.add_bond((h, reactive_monomer.id)) - - -RE = reaction_ensemble.WangLandauReactionEnsemble( - temperature=1, exclusion_radius=0, seed=77) -RE.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}) -print(RE.get_status()) -system.setup_type_map([0, 1, 2, 3]) - - -# initialize wang_landau -# generate preliminary_energy_run_results here, this should be done in a -# separate simulation without energy reweighting using the update energy -# functions -np.savetxt("energy_boundaries.dat", np.c_[[0, 1], [0, 0], [9, 9]], - header="nbar E_min E_max") - -RE.add_collective_variable_degree_of_association( - associated_type=0, min=0, max=1, corresponding_acid_types=[0, 1]) -RE.add_collective_variable_potential_energy( - filename="energy_boundaries.dat", delta=0.05) -RE.set_wang_landau_parameters( - final_wang_landau_parameter=1e-3, - do_not_sample_reaction_partition_function=True, - full_path_to_output_filename="WL_potential_out.dat") - -i = 0 -while True: - RE.reaction() - RE.displacement_mc_move_for_particles_of_type(3) diff --git a/src/core/reaction_methods/CMakeLists.txt b/src/core/reaction_methods/CMakeLists.txt index 2cd4e8b9ed9..1fc7b299af4 100644 --- a/src/core/reaction_methods/CMakeLists.txt +++ b/src/core/reaction_methods/CMakeLists.txt @@ -22,7 +22,6 @@ target_sources( PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ReactionEnsemble.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ConstantpHEnsemble.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/WangLandauReactionEnsemble.cpp ${CMAKE_CURRENT_SOURCE_DIR}/WidomInsertion.cpp ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp) diff --git a/src/core/reaction_methods/ConstantpHEnsemble.cpp b/src/core/reaction_methods/ConstantpHEnsemble.cpp index 14acf2cd9d1..d071ae7a1bb 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.cpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.cpp @@ -35,7 +35,7 @@ namespace ReactionMethods { */ double ConstantpHEnsemble::calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, - std::map const &old_particle_numbers, int, int, bool) const { + std::map const &old_particle_numbers) 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 * diff --git a/src/core/reaction_methods/ConstantpHEnsemble.hpp b/src/core/reaction_methods/ConstantpHEnsemble.hpp index cd77c14ad19..b533999cfa9 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.hpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.hpp @@ -46,8 +46,8 @@ class ConstantpHEnsemble : public ReactionAlgorithm { protected: double calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, - double E_pot_new, std::map const &old_particle_numbers, int, - int, bool) const override; + double E_pot_new, + std::map const &old_particle_numbers) const override; }; } // namespace ReactionMethods diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index 7c130be42fe..af78b2b7c2d 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -293,12 +293,9 @@ void ReactionAlgorithm::generic_oneway_reaction( current_reaction.tried_moves += 1; particle_inside_exclusion_radius_touched = false; - int old_state_index = -1; // for Wang-Landau algorithm - on_reaction_entry(old_state_index); if (!all_reactant_particles_exist(current_reaction)) { // makes sure, no incomplete reaction is performed -> only need to consider // rollback of complete reactions - on_reaction_rejection_directly_after_entry(old_state_index); return; } @@ -336,14 +333,8 @@ void ReactionAlgorithm::generic_oneway_reaction( ? std::numeric_limits::max() : calculate_current_potential_energy_of_system(); - int new_state_index = -1; // save new_state_index for Wang-Landau algorithm - int accepted_state = -1; // for Wang-Landau algorithm - on_attempted_reaction(new_state_index); - - constexpr bool only_make_configuration_changing_move = false; auto const bf = calculate_acceptance_probability( - current_reaction, E_pot_old, E_pot_new, old_particle_numbers, - old_state_index, new_state_index, only_make_configuration_changing_move); + current_reaction, E_pot_old, E_pot_new, old_particle_numbers); std::vector exponential = { exp(-1.0 / temperature * (E_pot_new - E_pot_old))}; @@ -351,8 +342,6 @@ void ReactionAlgorithm::generic_oneway_reaction( if (m_uniform_real_distribution(m_generator) < bf) { // accept - accepted_state = new_state_index; - // delete hidden reactant_particles (remark: don't delete changed particles) // extract ids of to be deleted particles auto len_hidden_particles_properties = @@ -374,7 +363,6 @@ void ReactionAlgorithm::generic_oneway_reaction( current_reaction.accepted_moves += 1; } else { // reject - accepted_state = old_state_index; // reverse reaction // 1) delete created product particles for (int p_ids_created_particle : p_ids_created_particles) { @@ -386,7 +374,6 @@ void ReactionAlgorithm::generic_oneway_reaction( restore_properties(changed_particles_properties, number_of_saved_properties); } - on_end_reaction(accepted_state); } /** @@ -566,22 +553,14 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) { * 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) { + int type, int particle_number_of_type_to_be_changed) { m_tried_configurational_MC_moves += 1; particle_inside_exclusion_radius_touched = false; - int old_state_index = -1; - if (use_wang_landau) { - on_reaction_entry(old_state_index); - } - int particle_number_of_type = number_of_particles_with_type(type); if (particle_number_of_type == 0 or particle_number_of_type_to_be_changed == 0) { // reject - if (use_wang_landau) { - on_mc_rejection_directly_after_entry(old_state_index); - } return false; } @@ -596,20 +575,8 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type( double beta = 1.0 / temperature; - int new_state_index = -1; - double bf = 1.0; - std::map dummy_old_particle_numbers; - SingleReaction temp_unimportant_arbitrary_reaction; - - if (use_wang_landau) { - new_state_index = on_mc_use_WL_get_new_state(); - bf = calculate_acceptance_probability( - temp_unimportant_arbitrary_reaction, E_pot_old, E_pot_new, - dummy_old_particle_numbers, old_state_index, new_state_index, true); - } else { - // Metropolis algorithm since proposal density is symmetric - bf = std::min(1.0, bf * exp(-beta * (E_pot_new - E_pot_old))); - } + // Metropolis algorithm since proposal density is symmetric + auto const bf = std::min(1.0, exp(-beta * (E_pot_new - E_pot_old))); // // correct for enhanced proposal of small radii by using the // // Metropolis-Hastings algorithm for asymmetric proposal densities @@ -618,24 +585,16 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type( // std::pow(particle_positions[0][1]-cyl_y,2)); // double new_radius = // std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2)); - // bf = std::min(1.0, - // bf*exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius); + // auto const bf = std::min(1.0, + // exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius); // Metropolis-Hastings algorithm for asymmetric proposal density if (m_uniform_real_distribution(m_generator) < bf) { // accept m_accepted_configurational_MC_moves += 1; - if (use_wang_landau) { - on_mc_accept(new_state_index); - } return true; } - // reject - // modify wang_landau histogram and potential - if (use_wang_landau) { - on_mc_reject(old_state_index); - } - // restore original particle positions + // reject: restore original particle positions for (auto const &item : original_positions) place_particle(std::get<0>(item), std::get<1>(item)); return false; diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp index 3bf1433a2d1..6d04654d83e 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.hpp +++ b/src/core/reaction_methods/ReactionAlgorithm.hpp @@ -92,24 +92,13 @@ class ReactionAlgorithm { } bool do_global_mc_move_for_particles_of_type(int type, - int particle_number_of_type, - bool use_wang_landau); + int particle_number_of_type); - bool particle_inside_exclusion_radius_touched; + bool particle_inside_exclusion_radius_touched = false; protected: std::vector m_empty_p_ids_smaller_than_max_seen_particle; void generic_oneway_reaction(SingleReaction ¤t_reaction); - virtual void on_reaction_entry(int &old_state_index) {} - virtual void - on_reaction_rejection_directly_after_entry(int &old_state_index) {} - virtual void on_attempted_reaction(int &new_state_index) {} - virtual void on_end_reaction(int &accepted_state) {} - - 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; } std::tuple, std::vector, std::vector> @@ -136,9 +125,7 @@ class ReactionAlgorithm { 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 { + double E_pot_new, std::map const &old_particle_numbers) const { return -10; } diff --git a/src/core/reaction_methods/ReactionEnsemble.cpp b/src/core/reaction_methods/ReactionEnsemble.cpp index d97fc441ce5..09400bb1f84 100644 --- a/src/core/reaction_methods/ReactionEnsemble.cpp +++ b/src/core/reaction_methods/ReactionEnsemble.cpp @@ -32,7 +32,7 @@ namespace ReactionMethods { */ double ReactionEnsemble::calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, - std::map const &old_particle_numbers, int, int, bool) const { + std::map const &old_particle_numbers) const { const double factorial_expr = calculate_factorial_expression(current_reaction, old_particle_numbers); diff --git a/src/core/reaction_methods/ReactionEnsemble.hpp b/src/core/reaction_methods/ReactionEnsemble.hpp index 5820165fd08..c484e7176cd 100644 --- a/src/core/reaction_methods/ReactionEnsemble.hpp +++ b/src/core/reaction_methods/ReactionEnsemble.hpp @@ -41,8 +41,8 @@ class ReactionEnsemble : public ReactionAlgorithm { protected: double calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, - double E_pot_new, std::map const &old_particle_numbers, int, - int, bool) const override; + double E_pot_new, + std::map const &old_particle_numbers) const override; }; } // namespace ReactionMethods diff --git a/src/core/reaction_methods/WangLandauReactionEnsemble.cpp b/src/core/reaction_methods/WangLandauReactionEnsemble.cpp deleted file mode 100644 index 95a2e066858..00000000000 --- a/src/core/reaction_methods/WangLandauReactionEnsemble.cpp +++ /dev/null @@ -1,858 +0,0 @@ -/* - * Copyright (C) 2010-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 . - */ - -#include "reaction_methods/WangLandauReactionEnsemble.hpp" - -#include "utils.hpp" - -#include "energy.hpp" -#include "particle_data.hpp" - -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace ReactionMethods { - -double EnergyCollectiveVariable::determine_current_state() const { - return calculate_current_potential_energy_of_system(); -} - -double -DegreeOfAssociationCollectiveVariable::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 = - number_of_particles_with_type(corresponding_acid_type); - total_number_of_corresponding_acid += num_of_current_type; - } - if (total_number_of_corresponding_acid == 0) { - throw std::runtime_error("Have you forgotten to specify all " - "corresponding acid types? Total particle " - "number of corresponding acid type is zero\n"); - } - int num_of_associated_acid = number_of_particles_with_type(associated_type); - double degree_of_association = - static_cast(num_of_associated_acid) / - static_cast(total_number_of_corresponding_acid); - return degree_of_association; -} - -/** Load minimum and maximum energies as a function of the other collective - * variables. - */ -void EnergyCollectiveVariable::load_CV_boundaries( - WangLandauReactionEnsemble &wl_system, std::istream &infile) { - - wl_system.do_energy_reweighting = true; - - // Note that you cannot change the other collective variables in the - // pre-production run and the production run - - // The data is formatted as space-separated floating 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)) { - std::istringstream iss(line); - std::vector values; - double value = -1.0; - while (iss >> value) { - values.push_back(value); - } - 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); -} - -void WangLandauReactionEnsemble::on_reaction_entry(int &old_state_index) { - old_state_index = get_flattened_index_wang_landau_of_current_state(); - if (old_state_index >= 0) { - if (histogram[old_state_index] >= 0) - monte_carlo_trial_moves += 1; - } -} - -void WangLandauReactionEnsemble::on_reaction_rejection_directly_after_entry( - int &old_state_index) { - // 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) { - new_state_index = get_flattened_index_wang_landau_of_current_state(); -} - -void WangLandauReactionEnsemble::on_end_reaction(int &accepted_state) { - update_wang_landau_potential_and_histogram(accepted_state); -} - -void WangLandauReactionEnsemble::on_mc_rejection_directly_after_entry( - int &old_state_index) { - if (do_energy_reweighting) - update_wang_landau_potential_and_histogram(old_state_index); -} - -void WangLandauReactionEnsemble::on_mc_accept(int &new_state_index) { - if (do_energy_reweighting) { - // modify wang_landau histogram and potential - update_wang_landau_potential_and_histogram(new_state_index); - } -} - -void WangLandauReactionEnsemble::on_mc_reject(int &old_state_index) { - if (do_energy_reweighting) - update_wang_landau_potential_and_histogram(old_state_index); -} - -int WangLandauReactionEnsemble::on_mc_use_WL_get_new_state() { - return get_flattened_index_wang_landau_of_current_state(); -} - -/** - * Adds a new collective variable (CV) of the type degree of association to the - * Wang-Landau sampling - */ -void WangLandauReactionEnsemble::add_new_CV_degree_of_association( - int associated_type, double CV_minimum, double CV_maximum, - const std::vector &corresponding_acid_types) { - std::shared_ptr - new_collective_variable = - std::make_shared(); - new_collective_variable->associated_type = associated_type; - new_collective_variable->CV_minimum = CV_minimum; - new_collective_variable->CV_maximum = CV_maximum; - new_collective_variable->corresponding_acid_types = corresponding_acid_types; - new_collective_variable->delta_CV = - calculate_delta_degree_of_association(*new_collective_variable); - collective_variables.push_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( - std::istream &infile, double delta_CV) { - std::shared_ptr new_collective_variable = - std::make_shared(); - new_collective_variable->delta_CV = delta_CV; - 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 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) { - - // 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] + - delta_collective_variables_values[CV_i] * 0.98 || - current_state[CV_i] < - collective_variables_minimum_values[CV_i] - - delta_collective_variables_values[CV_i] * 0.01) { - return -10; - } - } - - std::vector individual_indices(nr_collective_variables); - for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - 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 - 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++) { - factor *= nr_subindices_of_collective_variable[j]; - } - index += factor * individual_indices[CV_i]; - } - return index; -} - -/** - * Returns the flattened index of the multidimensional Wang-Landau histogram for - * the current state of the simulation. - */ -int WangLandauReactionEnsemble:: - get_flattened_index_wang_landau_of_current_state() { - 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++) - current_state[CV_i] = collective_variables[CV_i]->determine_current_state(); - // get collective_variables_minimum_values - std::vector collective_variables_minimum_values( - nr_collective_variables); - for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - collective_variables_minimum_values[CV_i] = - collective_variables[CV_i]->CV_minimum; - } - // get collective_variables_maximum_values - std::vector collective_variables_maximum_values( - nr_collective_variables); - for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - collective_variables_maximum_values[CV_i] = - collective_variables[CV_i]->CV_maximum; - } - // get delta_collective_variables_values - std::vector delta_collective_variables_values( - nr_collective_variables); - for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - delta_collective_variables_values[CV_i] = - collective_variables[CV_i]->delta_CV; - } - int index = get_flattened_index_wang_landau( - current_state, collective_variables_minimum_values, - collective_variables_maximum_values, delta_collective_variables_values, - nr_collective_variables); - return index; -} - -/** - * Returns the minimum value of the collective variable on a delta_CV spaced - * grid which starts at 0 - */ -double WangLandauReactionEnsemble::get_minimum_CV_value_on_delta_CV_spaced_grid( - 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; - return minimum_CV_value_on_delta_CV_spaced_grid; -} - -/** - * Calculates the smallest difference in the degree of association which can be - * observed when changing the degree of association by one single reaction. - */ -double WangLandauReactionEnsemble::calculate_delta_degree_of_association( - DegreeOfAssociationCollectiveVariable ¤t_collective_variable) { - // calculate Delta in the degree of association so that EVERY reaction step is - // driven. - int total_number_of_corresponding_acid = 0; - for (int corresponding_acid_type : - current_collective_variable.corresponding_acid_types) { - int num_of_current_type = - number_of_particles_with_type(corresponding_acid_type); - total_number_of_corresponding_acid += num_of_current_type; - } - double delta = 1.0 / total_number_of_corresponding_acid; - // now modify the minimum value of the CV to lie on the grid - current_collective_variable.CV_minimum = - get_minimum_CV_value_on_delta_CV_spaced_grid( - current_collective_variable.CV_minimum, delta); - return delta; -} - -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 - - 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) { - // assume the energy collective variable to be the last added - // collective variable - EnergyCollectiveVariable_index = - static_cast(collective_variables.size()) - 1; - } - double current_energy = - 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( - static_cast(flattened_index), EnergyCollectiveVariable_index); - std::shared_ptr energy_CV = - collective_variables[EnergyCollectiveVariable_index]; - if (current_energy > - max_boundaries_energies[flat_index_without_energy_CV] || - current_energy < min_boundaries_energies[flat_index_without_energy_CV] - - energy_CV->delta_CV) { - histogram[flattened_index] = int_fill_value; - wang_landau_potential[flattened_index] = double_fill_value; - empty_bins_in_memory++; - } - } - - used_bins = - static_cast(wang_landau_potential.size()) - empty_bins_in_memory; -} - -int WangLandauReactionEnsemble::initialize_wang_landau() { - - nr_subindices_of_collective_variable.resize(collective_variables.size(), 0); - auto const &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 - - if (do_energy_reweighting) { - invalidate_bins(); - } - return ES_OK; -} - -/** Calculate the expression in the acceptance probability of the Wang-Landau - * reaction ensemble. - * Modify Boltzmann factor according to Wang-Landau algorithm in @cite yan02b. - */ -double WangLandauReactionEnsemble::calculate_acceptance_probability( - SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, - std::map const &old_particle_numbers, int old_state_index, - int new_state_index, bool only_make_configuration_changing_move) const { - - double beta = 1.0 / temperature; - double bf = 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 * - factorial_expr; - } - - if (!do_energy_reweighting) { - bf *= exp(-beta * (E_pot_new - E_pot_old)); - } - - // 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 - // array of the histogram and the wang_landau_potential values is "cuboid"). - if (old_state_index >= 0 && new_state_index >= 0) { - if (histogram[new_state_index] >= 0 && histogram[old_state_index] >= 0) { - // Modify Boltzmann factor (bf) according to Wang-Landau algorithm - // @cite yan02b. - // This makes the new state being accepted with the conditional - // probability bf (bf is a transition probability = conditional - // probability from the old state to move to the new state). - bf = std::min(1.0, bf * exp(wang_landau_potential[old_state_index] - - wang_landau_potential[new_state_index])); - } else { - if (histogram[old_state_index] < 0) - bf = 10; // accept the reaction if we found a state in gamma - // (histogram[new_state_index] >= 0) or to sample new configs - // which might lie in gamma(histogram[new_state_index] < 0) - else if (histogram[new_state_index] < 0) - bf = -10; // reject the reaction, since the new state - // is not in gamma while the old sate was in gamma - } - } else if (old_state_index < 0 && new_state_index >= 0) { - bf = 10; // accept the reaction if we found a new state in gamma - // (new_state_index >= 0) or to sample new configs which - // might lie in gamma (new_state_index < 0) - } else if (old_state_index >= 0 && new_state_index < 0) { - bf = -10; // reject the reaction, since the new state is - // not in gamma while the old sate was in gamma - } - return bf; -} - -/** Perform a randomly selected reaction using the Wang-Landau algorithm. - * - * Make sure to perform additional configuration changing steps, after the - * reaction step! Like in @cite yan02b. This can be done with MD in the case - * of the no-energy-reweighting case, or with the function - * @ref ReactionAlgorithm::do_global_mc_move_for_particles_of_type. - * - * Perform additional Monte Carlo moves to sample configurational - * partition function according to @cite yan02b. Do as many steps - * as needed to get to a new conformation. - */ -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(static_cast(reactions.size())); - generic_oneway_reaction(reactions[reaction_id]); - if (can_refine_wang_landau_one_over_t() && m_WL_tries % 10000 == 0) { - // check for convergence - if (achieved_desired_number_of_refinements_one_over_t()) { - write_wang_landau_results_to_file(output_filename); - return -10; // return negative value to indicate that the Wang-Landau - // algorithm has converged - } - refine_wang_landau_parameter_one_over_t(); - } - } - // shift Wang-Landau potential minimum to zero - if (m_WL_tries % (std::max(90000, 9 * reaction_steps)) == 0) { - // for numerical stability here we also subtract the minimum positive value - // of the wang_landau_potential from the wang_landau potential, allowed - // since only the difference in the wang_landau potential is of interest. - double minimum_wang_landau_potential = - find_minimum_non_negative_value(wang_landau_potential); - for (double &i : wang_landau_potential) { - if (i >= 0) // check for whether we are in the - // valid range of the collective variable - i -= minimum_wang_landau_potential; - } - // write out preliminary Wang-Landau potential results - write_wang_landau_results_to_file(output_filename); - } - return 0; -} - -/** Increase the Wang-Landau potential and histogram at the current nbar */ -void WangLandauReactionEnsemble::update_wang_landau_potential_and_histogram( - int index_of_state_after_acceptance_or_rejection) { - if (index_of_state_after_acceptance_or_rejection >= 0) { - if (histogram[index_of_state_after_acceptance_or_rejection] >= 0) { - histogram[index_of_state_after_acceptance_or_rejection] += 1; - wang_landau_potential[index_of_state_after_acceptance_or_rejection] += - wang_landau_parameter; - } - } -} - -/** - *Determines whether we can reduce the Wang-Landau parameter - */ -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 - // configuration space better. Use flatness - // criterion according to 1/t algorithm as long as - // you are not in 1/t regime. - if (do_energy_reweighting) - minimum_required_value = 20; // get faster in energy reweighting case - - return *(std::min_element(histogram.begin(), histogram.end())) > - minimum_required_value || - m_system_is_in_1_over_t_regime; -} - -/** - *Reset the Wang-Landau histogram. - */ -void WangLandauReactionEnsemble::reset_histogram() { - printf("Histogram is flat. Refining. Previous Wang-Landau modification " - "parameter was %f.\n", - wang_landau_parameter); - fflush(stdout); - - for (int i = 0; i < wang_landau_potential.size(); i++) { - if (histogram[i] >= 0) { // checks for validity of index i (think of energy - // collective variables, in a cubic memory layout - // there will be indices which are not allowed by - // the energy boundaries. These values will be - // initialized with a negative fill value) - histogram[i] = 0; - } - } -} - -/** - *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) / - 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; - if (!m_system_is_in_1_over_t_regime) { - m_system_is_in_1_over_t_regime = true; - printf("Refining: Wang-Landau parameter is now 1/t.\n"); - } - } else { - reset_histogram(); - wang_landau_parameter = wang_landau_parameter / 2.0; - } -} - -/** - *Determine whether the desired number of refinements was achieved. - */ -bool WangLandauReactionEnsemble:: - 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; - } - return false; -} - -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 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 (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 << " "; - } - out << wang_landau_potential[flattened_index] << " \n"; - } - } -} - -/** - *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(); -} - -/** - *Update the minimum and maximum observed energies using the current state. - *Needed for preliminary energy reweighting runs. - */ -int WangLandauReactionEnsemble:: - update_maximum_and_minimum_energies_at_current_state() { - if (minimum_energies_at_flat_index.empty() || - maximum_energies_at_flat_index.empty()) { - minimum_energies_at_flat_index.resize(wang_landau_potential.size(), - double_fill_value); - maximum_energies_at_flat_index.resize(wang_landau_potential.size(), - double_fill_value); - } - - const double E_pot_current = calculate_current_potential_energy_of_system(); - int index = get_flattened_index_wang_landau_of_current_state(); - - // update stored energy values - if (((E_pot_current < minimum_energies_at_flat_index[index]) || - std::abs(minimum_energies_at_flat_index[index] - double_fill_value) < - std::numeric_limits::epsilon())) { - minimum_energies_at_flat_index[index] = E_pot_current; - } - if (((E_pot_current > maximum_energies_at_flat_index[index]) || - std::abs(maximum_energies_at_flat_index[index] - double_fill_value) < - std::numeric_limits::epsilon())) { - maximum_energies_at_flat_index[index] = E_pot_current; - } - - return 0; -} - -/** - *Write out an energy boundary file using the energy boundaries observed in a - *preliminary energy reweighting run. - */ -void WangLandauReactionEnsemble::write_out_preliminary_energy_run_results( - const std::string &filename) { - - std::ofstream outfile; - outfile.open(filename); - if (!outfile.is_open()) { - throw std::runtime_error("Cannot write to " + filename); - } - outfile << "#nbar E_min E_max\n"; - for (std::size_t flattened_index = 0; - flattened_index < wang_landau_potential.size(); flattened_index++) { - // unravel index - 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 (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 << " "; - } - outfile << minimum_energies_at_flat_index[flattened_index] << " " - << maximum_energies_at_flat_index[flattened_index] << " \n"; - } - outfile.close(); -} - -/** - *Returns the flattened index of a given flattened index without the energy - *collective variable. - */ -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(), - 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] = static_cast(unraveled_index[i]) * - collective_variables[i]->delta_CV + - collective_variables[i]->CV_minimum; - } - - // get collective_variables_minimum_values - std::vector collective_variables_minimum_values( - nr_collective_variables); - for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - collective_variables_minimum_values[CV_i] = - collective_variables[CV_i]->CV_minimum; - } - // get collective_variables_maximum_values - std::vector collective_variables_maximum_values( - nr_collective_variables); - for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - collective_variables_maximum_values[CV_i] = - collective_variables[CV_i]->CV_maximum; - } - // get delta_collective_variables_values - std::vector delta_collective_variables_values( - nr_collective_variables); - for (int CV_i = 0; CV_i < nr_collective_variables; CV_i++) { - delta_collective_variables_values[CV_i] = - collective_variables[CV_i]->delta_CV; - } - int index = get_flattened_index_wang_landau( - current_state, collective_variables_minimum_values, - collective_variables_maximum_values, delta_collective_variables_values, - nr_collective_variables); - return index; -} - -/** - * Writes the Wang-Landau parameter, the histogram and the potential to a file. - * You can restart a Wang-Landau simulation using this information. - * Additionally you should store the positions of the particles. - * Not storing them introduces small, small statistical errors. - */ -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) - 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 - 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 - 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(); -} - -/** - *Loads the Wang-Landau checkpoint - */ -void WangLandauReactionEnsemble::load_wang_landau_checkpoint( - const std::string &identifier) { - std::string filename; - std::ifstream infile; - - // restore Wang-Landau parameters - 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; - int line = 0; - while (infile >> wang_landau_parameter_entry >> - wang_landau_monte_carlo_trial_moves_entry >> - flat_index_of_state_at_checkpointing) { - wang_landau_parameter = wang_landau_parameter_entry; - monte_carlo_trial_moves = wang_landau_monte_carlo_trial_moves_entry; - line += 1; - } - infile.close(); - } else { - throw std::runtime_error("Cannot read " + filename); - } - - // restore histogram - filename = std::string("checkpoint_wang_landau_histogram_") + identifier; - infile.open(filename); - if (infile.is_open()) { - int hist_entry; - int line = 0; - while (infile >> hist_entry) { - histogram[line] = hist_entry; - line += 1; - } - infile.close(); - } else { - throw std::runtime_error("Cannot read " + filename); - } - - // restore Wang-Landau potential - filename = std::string("checkpoint_wang_landau_potential_") + identifier; - infile.open(filename); - if (infile.is_open()) { - double wang_landau_potential_entry; - int line = 0; - while (infile >> wang_landau_potential_entry) { - wang_landau_potential[line] = wang_landau_potential_entry; - line += 1; - } - infile.close(); - } else { - 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. -} - -} // namespace ReactionMethods diff --git a/src/core/reaction_methods/WangLandauReactionEnsemble.hpp b/src/core/reaction_methods/WangLandauReactionEnsemble.hpp deleted file mode 100644 index 50aab96e128..00000000000 --- a/src/core/reaction_methods/WangLandauReactionEnsemble.hpp +++ /dev/null @@ -1,182 +0,0 @@ -/* - * Copyright (C) 2010-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 . - */ -#ifndef REACTION_METHODS_WANG_LANDAU_REACTION_ENSEMBLE_HPP -#define REACTION_METHODS_WANG_LANDAU_REACTION_ENSEMBLE_HPP - -#include "reaction_methods/ReactionAlgorithm.hpp" - -#include -#include -#include -#include -#include -#include - -namespace ReactionMethods { - -struct CollectiveVariable { - double CV_minimum = {}; - double CV_maximum = {}; - double delta_CV = {}; - // 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 { - double determine_current_state() const override; - 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; - /** Calculate the degree of association of the reference species */ - double determine_current_state() const override { - return calculate_degree_of_association(); - } - -private: - /** - * Return the degree of association for the current collective variable. This - * is needed since you may use multiple degrees of association as collective - * variable for the Wang-Landau algorithm. - */ - double calculate_degree_of_association() const; -}; - -/** Wang-Landau reaction ensemble method */ -class WangLandauReactionEnsemble : public ReactionAlgorithm { -public: - WangLandauReactionEnsemble(int seed) : ReactionAlgorithm(seed) {} - bool do_energy_reweighting = false; - bool do_not_sample_reaction_partition_function = false; - double final_wang_landau_parameter = 0.00001; - - void add_new_CV_degree_of_association( - int associated_type, double CV_minimum, double CV_maximum, - 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 = ""; - - std::vector min_boundaries_energies; - std::vector max_boundaries_energies; - - // 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; - - // 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 - 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; - void - on_reaction_rejection_directly_after_entry(int &old_state_index) override; - void on_attempted_reaction(int &new_state_index) override; - void on_end_reaction(int &accepted_state) override; - 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; - - // logarithm to basis e of the degeneracy of the states - std::vector wang_landau_potential; - - std::vector nr_subindices_of_collective_variable; - - // 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; - - 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( - int flattened_index_with_EnergyCollectiveVariable, - int collective_variable_index_energy_observable); // needed for energy - - int get_flattened_index_wang_landau( - 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() const; - bool m_system_is_in_1_over_t_regime = false; - bool achieved_desired_number_of_refinements_one_over_t() const; - void refine_wang_landau_parameter_one_over_t(); - - /** - * 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); - void invalidate_bins(); - void reset_histogram(); - double get_minimum_CV_value_on_delta_CV_spaced_grid(double min_CV_value, - double delta_CV) const; -}; - -} // namespace ReactionMethods -#endif diff --git a/src/core/reaction_methods/tests/CMakeLists.txt b/src/core/reaction_methods/tests/CMakeLists.txt index 0d12efaf629..fcc465c381a 100644 --- a/src/core/reaction_methods/tests/CMakeLists.txt +++ b/src/core/reaction_methods/tests/CMakeLists.txt @@ -25,9 +25,6 @@ unit_test(NAME ReactionAlgorithm_test SRC ReactionAlgorithm_test.cpp DEPENDS EspressoCore Boost::mpi MPI::MPI_CXX) unit_test(NAME ReactionEnsemble_test SRC ReactionEnsemble_test.cpp DEPENDS EspressoCore Boost::mpi MPI::MPI_CXX) -unit_test(NAME WangLandauReactionEnsemble_test SRC - WangLandauReactionEnsemble_test.cpp DEPENDS EspressoCore Boost::mpi - MPI::MPI_CXX) unit_test(NAME particle_tracking_test SRC particle_tracking_test.cpp DEPENDS EspressoCore Boost::mpi MPI::MPI_CXX) unit_test(NAME reaction_methods_utils_test SRC reaction_methods_utils_test.cpp diff --git a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp index 340d1be5ed4..09813744ce7 100644 --- a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp @@ -80,7 +80,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { std::log(10.) * (r_algo.m_constant_pH + std::log10(reaction.gamma))); auto const acceptance = r_algo.calculate_acceptance_probability( - reaction, energy, 0., p_numbers, -1, -1, false); + reaction, energy, 0., p_numbers); BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); } } diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp index cad8d8d2b48..537c2eec4c8 100644 --- a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp +++ b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp @@ -27,16 +27,23 @@ #include "reaction_methods/ReactionAlgorithm.hpp" +#include "EspressoSystemStandAlone.hpp" #include "Particle.hpp" #include "communication.hpp" #include "particle_data.hpp" #include +#include #include #include #include +namespace espresso { +// ESPResSo system instance +std::unique_ptr system; +} // namespace espresso + // Check the base class for all Monte Carlo algorithms. BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { using namespace ReactionMethods; @@ -70,6 +77,11 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { int const type_C = 2; SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); + // track particles + init_type_map(type_A); + init_type_map(type_B); + init_type_map(type_C); + // check reaction addition { r_algo.add_reaction(reaction.gamma, reaction.reactant_types, @@ -90,8 +102,8 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { // check acceptance probability { - double probability = r_algo.calculate_acceptance_probability( - reaction, -1., -1., {{1, 2}}, -1, -1, false); + double probability = + r_algo.calculate_acceptance_probability(reaction, -1., -1., {{1, 2}}); BOOST_CHECK_EQUAL(probability, -10.); } @@ -144,8 +156,8 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { place_particle(1, ref_positions[1].first); set_particle_v(0, ref_positions[0].second); set_particle_v(1, ref_positions[1].second); - // track particles - init_type_map(0); + set_particle_type(0, type_A); + set_particle_type(1, type_A); // update particle positions and velocities BOOST_CHECK(!r_algo.particle_inside_exclusion_radius_touched); r_algo.particle_inside_exclusion_radius_touched = false; @@ -167,12 +179,57 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { BOOST_CHECK_GE((new_pos - ref_old_pos).norm(), 0.1); BOOST_CHECK_GE((new_vel - ref_old_vel).norm(), 10.); } + // cleanup + remove_particle(0); + remove_particle(1); + } + + // check Monte Carlo moves + { + // set up particles + auto const box_l = 1.; + std::vector ref_positions{{0.1, 0.2, 0.3}, + {0.4, 0.5, 0.6}}; + place_particle(0, ref_positions[0]); + place_particle(1, ref_positions[1]); + set_particle_type(0, type_A); + set_particle_type(1, type_A); + // check early exit when a MC move cannot be performed + BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_C, 1)); + BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_B, 2)); + BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_A, 0)); + // force all MC moves to be rejected by picking particles inside + // their exclusion radius + r_algo.exclusion_radius = box_l; + r_algo.particle_inside_exclusion_radius_touched = false; + BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_A, 2)); + // check none of the particles moved + for (auto const pid : {0, 1}) { + auto const ref_old_pos = ref_positions[pid]; + auto const &p = get_particle_data(pid); + auto const &new_pos = p.r.p; + BOOST_CHECK_LE((new_pos - ref_old_pos).norm(), tol); + } + // force a MC move to be accepted by using a constant Hamiltonian + r_algo.exclusion_radius = 0.; + r_algo.particle_inside_exclusion_radius_touched = false; + BOOST_REQUIRE(r_algo.do_global_mc_move_for_particles_of_type(type_A, 1)); + std::vector distances(2); + // check that only one particle moved + for (auto const pid : {0, 1}) { + auto const &p = get_particle_data(pid); + distances[pid] = (ref_positions[pid] - p.r.p).norm(); + } + BOOST_CHECK_LE(std::min(distances[0], distances[1]), tol); + BOOST_CHECK_GE(std::max(distances[0], distances[1]), 0.1); + // cleanup + remove_particle(0); + remove_particle(1); } } int main(int argc, char **argv) { - auto mpi_env = std::make_shared(argc, argv); - Communication::init(mpi_env); + espresso::system = std::make_unique(argc, argv); return boost::unit_test::unit_test_main(init_unit_test, argc, argv); } diff --git a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp index ef1188551f2..6b78d0e023c 100644 --- a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp @@ -78,7 +78,7 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { reaction.gamma * f_expr * std::exp(energy / r_algo.temperature); auto const acceptance = r_algo.calculate_acceptance_probability( - reaction, energy, 0., p_numbers, -1, -1, false); + reaction, energy, 0., p_numbers); BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); } } diff --git a/src/core/reaction_methods/tests/SingleReaction_test.cpp b/src/core/reaction_methods/tests/SingleReaction_test.cpp index 42559922ba1..ea9d08e0f20 100644 --- a/src/core/reaction_methods/tests/SingleReaction_test.cpp +++ b/src/core/reaction_methods/tests/SingleReaction_test.cpp @@ -71,3 +71,17 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) { } } } + +BOOST_AUTO_TEST_CASE(SingleReaction_default_test) { + using namespace ReactionMethods; + + SingleReaction default_reaction{}; + BOOST_CHECK_EQUAL(default_reaction.gamma, 0.); + BOOST_CHECK_EQUAL(default_reaction.nu_bar, 0); + BOOST_CHECK_EQUAL(default_reaction.tried_moves, 0); + BOOST_CHECK_EQUAL(default_reaction.accepted_moves, 0); + BOOST_CHECK(default_reaction.reactant_types.empty()); + BOOST_CHECK(default_reaction.reactant_coefficients.empty()); + BOOST_CHECK(default_reaction.product_types.empty()); + BOOST_CHECK(default_reaction.product_coefficients.empty()); +} diff --git a/src/core/reaction_methods/tests/WangLandauReactionEnsemble_test.cpp b/src/core/reaction_methods/tests/WangLandauReactionEnsemble_test.cpp deleted file mode 100644 index 77c77d77286..00000000000 --- a/src/core/reaction_methods/tests/WangLandauReactionEnsemble_test.cpp +++ /dev/null @@ -1,205 +0,0 @@ -/* - * 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 Wang-Landau Reaction Ensemble method. */ - -#define BOOST_TEST_NO_MAIN -#define BOOST_TEST_MODULE Wang Landau Reaction Ensemble test -#define BOOST_TEST_ALTERNATIVE_INIT_API -#define BOOST_TEST_DYN_LINK -#include - -#include "unit_tests/ParticleFactory.hpp" - -#include "reaction_methods/WangLandauReactionEnsemble.hpp" -#include "reaction_methods/utils.hpp" - -#include "communication.hpp" -#include "particle_data.hpp" - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -// 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 ReactionMethods; - 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 the Monte Carlo algorithm where moves depend on the system -// configuration and/or energy. -BOOST_AUTO_TEST_CASE(WangLandauReactionEnsemble_test) { - using namespace ReactionMethods; - 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_method(), std::runtime_error); - - // create a reaction A -> 3 B + 4 C - int const type_A = 0; - int const type_B = 1; - int const type_C = 2; - SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4}); - - // check acceptance probability - 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); - } -} - -int main(int argc, char **argv) { - auto mpi_env = std::make_shared(argc, argv); - Communication::init(mpi_env); - - return boost::unit_test::unit_test_main(init_unit_test, argc, argv); -} diff --git a/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp b/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp index 9fff11ee9fe..135596e48c3 100644 --- a/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp +++ b/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp @@ -30,32 +30,6 @@ #include #include -BOOST_AUTO_TEST_CASE(find_minimum_non_negative_value_test) { - using namespace ReactionMethods; - - BOOST_CHECK_EQUAL(find_minimum_non_negative_value({1, 2, 3}), 1); - BOOST_CHECK_EQUAL(find_minimum_non_negative_value({3, 2, 1}), 1); - BOOST_CHECK_EQUAL(find_minimum_non_negative_value({-1, 2, 3}), 2); - BOOST_CHECK_EQUAL(find_minimum_non_negative_value({-1, -2, -3}), -3); - BOOST_CHECK_THROW(find_minimum_non_negative_value({}), std::runtime_error); -} - -BOOST_AUTO_TEST_CASE(average_list_of_allowed_entries_test) { - using namespace ReactionMethods; - constexpr double tol = 100 * std::numeric_limits::epsilon(); - - BOOST_CHECK_CLOSE(average_list_of_allowed_entries(std::vector{1, 2}), - 1.5, tol); - BOOST_CHECK_CLOSE( - average_list_of_allowed_entries(std::vector{1, 2, -2}), 1.5, tol); - BOOST_CHECK_CLOSE( - average_list_of_allowed_entries(std::vector{1.5, -3.}), 1.5, tol); - BOOST_CHECK_CLOSE( - average_list_of_allowed_entries(std::vector{-1, -2, -3}), 0.0, tol); - BOOST_CHECK_CLOSE(average_list_of_allowed_entries(std::vector{}), 0.0, - tol); -} - BOOST_AUTO_TEST_CASE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i_test) { using namespace ReactionMethods; constexpr double tol = 100 * std::numeric_limits::epsilon(); diff --git a/src/core/reaction_methods/utils.hpp b/src/core/reaction_methods/utils.hpp index 623c6126c17..5ab65babf00 100644 --- a/src/core/reaction_methods/utils.hpp +++ b/src/core/reaction_methods/utils.hpp @@ -21,10 +21,7 @@ #include "reaction_methods/ReactionAlgorithm.hpp" -#include #include -#include -#include namespace ReactionMethods { @@ -54,56 +51,5 @@ double calculate_factorial_expression_cpH( */ double factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(int Ni0, int nu_i); -/** - * Calculate the average of an array (used for the histogram of the - * Wang-Landau algorithm). It excludes values which are initialized to be - * negative. Those values indicate that the Wang-Landau algorithm should not - * sample those values. The values still occur in the list because we can only - * store "rectangular" value ranges. - */ -template -double average_list_of_allowed_entries(const std::vector &rng) { - T result = 0; - int counter_allowed_entries = 0; - for (auto &val : rng) { - if (val >= 0) { // checks for validity of index i (think of energy - // collective variables, in a cubic memory layout - // there will be indices which are not allowed by - // the energy boundaries. These values will be - // initialized with a negative fill value) - result += val; - counter_allowed_entries += 1; - } - } - if (counter_allowed_entries) { - return static_cast(result) / - static_cast(counter_allowed_entries); - } - return 0.0; -} - -/** - * Finds the minimum non negative value in the provided range and returns - * this value. - */ -inline double find_minimum_non_negative_value(std::vector const &rng) { - if (rng.empty()) - throw std::runtime_error("range is empty\n"); - // think of negative histogram values that indicate not - // allowed energies in the case of an energy observable - auto const it = std::min_element(rng.begin(), rng.end(), - [](double const &a, double const &b) { - if (a <= 0) - return false; - if (b <= 0) - return true; - return a < b; - }); - if (it == rng.end() or *it < 0) { - return rng.back(); - } - return *it; -} - } // namespace ReactionMethods #endif diff --git a/src/python/espressomd/reaction_ensemble.pxd b/src/python/espressomd/reaction_ensemble.pxd index 7a66f192b2f..0afb10082e2 100644 --- a/src/python/espressomd/reaction_ensemble.pxd +++ b/src/python/espressomd/reaction_ensemble.pxd @@ -38,7 +38,7 @@ cdef extern from "reaction_methods/ReactionAlgorithm.hpp" namespace "ReactionMet cdef cppclass CReactionAlgorithm "ReactionMethods::ReactionAlgorithm": int CReactionAlgorithm(int seed) int do_reaction(int reaction_steps) except + - bool do_global_mc_move_for_particles_of_type(int type, int particle_number_of_type, bool use_wang_landau) + bool do_global_mc_move_for_particles_of_type(int type, int particle_number_of_type) void set_cuboid_reaction_ensemble_volume() int check_reaction_method() except + double get_acceptance_rate_configurational_moves() @@ -65,24 +65,6 @@ cdef extern from "reaction_methods/ReactionEnsemble.hpp" namespace "ReactionMeth cdef cppclass CReactionEnsemble "ReactionMethods::ReactionEnsemble"(CReactionAlgorithm): CReactionEnsemble(int seed) -cdef extern from "reaction_methods/WangLandauReactionEnsemble.hpp" namespace "ReactionMethods": - - cdef cppclass CWangLandauReactionEnsemble "ReactionMethods::WangLandauReactionEnsemble"(CReactionAlgorithm): - CWangLandauReactionEnsemble(int seed) - double wang_landau_parameter - double final_wang_landau_parameter - string output_filename - vector[double] minimum_energies_at_flat_index - 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) except + - int update_maximum_and_minimum_energies_at_current_state() - void write_out_preliminary_energy_run_results(string filename) except + - void write_wang_landau_checkpoint(string identifier) except + - void load_wang_landau_checkpoint(string identifier) except + - void write_wang_landau_results_to_file(string filename) except + - cdef extern from "reaction_methods/ConstantpHEnsemble.hpp" namespace "ReactionMethods": cdef cppclass CConstantpHEnsemble "ReactionMethods::ConstantpHEnsemble"(CReactionAlgorithm): diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index fd82526ef32..30f545da207 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -21,16 +21,12 @@ from cython.operator cimport dereference as deref import numpy as np -class WangLandauHasConverged(Exception): - pass - - cdef class ReactionAlgorithm: """ - This class provides the base class for Reaction Algorithms like the Reaction - Ensemble algorithm, the Wang-Landau Reaction Ensemble algorithm and the - constant pH method. Initialize the reaction algorithm by setting the + This class provides the base class for Reaction Algorithms like + the Reaction Ensemble algorithm and the constant pH method. + Initialize the reaction algorithm by setting the standard pressure, temperature, and the exclusion radius. Note: When creating particles the velocities of the new particles are set @@ -299,9 +295,8 @@ cdef class ReactionAlgorithm: """ - use_wang_landau = False deref(self.RE).do_global_mc_move_for_particles_of_type( - type_mc, particle_number_to_be_changed, use_wang_landau) + type_mc, particle_number_to_be_changed) def get_status(self): """ @@ -484,249 +479,6 @@ cdef class ConstantpHEnsemble(ReactionAlgorithm): deref(self.constpHptr).m_constant_pH = pH -cdef class WangLandauReactionEnsemble(ReactionAlgorithm): - """ - This Class implements the Wang-Landau Reaction Ensemble. - """ - - cdef unique_ptr[CWangLandauReactionEnsemble] WLRptr - - def __init__(self, *args, **kwargs): - self._params = {"temperature": 1, - "exclusion_radius": 0} - for k in self._required_keys(): - if k not in kwargs: - raise ValueError( - "At least the following keys have to be given as keyword arguments: " + self._required_keys().__str__() + " got " + kwargs.__str__()) - self._params[k] = kwargs[k] - for k in kwargs: - if k in self._valid_keys(): - self._params[k] = kwargs[k] - else: - raise KeyError(f"{k} is not a valid key") - - self.WLRptr.reset(new CWangLandauReactionEnsemble(int(self._params["seed"]))) - self.RE = self.WLRptr.get() - - self._set_params_in_es_core() - - def reaction(self, reaction_steps=1): - """ - Performs reaction_steps reactions. Sets the number of reaction steps - which are performed at once. Do not use too many reaction steps - consecutively without having conformation-changing steps in-between - (especially important for the Wang-Landau reaction ensemble). Providing - a number for the parameter reaction steps reduces the need for the - interpreter to be called between consecutive reactions. - - """ - status_wang_landau = deref( - self.WLRptr).do_reaction(int(reaction_steps)) - if status_wang_landau < 0: - raise WangLandauHasConverged( - "The Wang-Landau algorithm has converged.") - - def add_collective_variable_degree_of_association(self, *args, **kwargs): - """ - Adds the degree of association as a collective variable (reaction coordinate) for the Wang-Landau Reaction Ensemble. - Several collective variables can be set simultaneously. - - Parameters - ---------- - associated_type : :obj:`int` - Particle type of the associated state of the reacting species. - min : :obj:`float` - Minimum value of the collective variable. - max : :obj:`float` - Maximum value of the collective variable. - corresponding_acid_types : list of :obj:`int` - List of the types of the version of the species. - - """ - for k in kwargs: - if k in self._valid_keys_add_collective_variable_degree_of_association(): - self._params[k] = kwargs[k] - else: - raise KeyError(f"{k} is not a valid key") - - for k in self._required_keys_add_collective_variable_degree_of_association(): - if k not in kwargs: - raise ValueError( - "At least the following keys have to be given as keyword arguments: " + self._required_keys_add_collective_variable_degree_of_association().__str__() + " got " + kwargs.__str__()) - self._params[k] = kwargs[k] - - cdef vector[int] _corresponding_acid_types - for i in range(len(self._params["corresponding_acid_types"])): - _corresponding_acid_types.push_back( - self._params["corresponding_acid_types"][i]) - deref(self.WLRptr).add_new_CV_degree_of_association( - self._params["associated_type"], self._params["min"], self._params["max"], _corresponding_acid_types) - - def _valid_keys_add_collective_variable_degree_of_association(self): - return "associated_type", "min", "max", "corresponding_acid_types" - - def _required_keys_add_collective_variable_degree_of_association(self): - return "associated_type", "min", "max", "corresponding_acid_types" - - def add_collective_variable_potential_energy(self, *args, **kwargs): - """ - Adds the potential energy as a collective variable (reaction coordinate) for the Wang-Landau Reaction Ensemble. - Several collective variables can be set simultaneously. - - Parameters - ---------- - filename : :obj:`str` - Filename of the energy boundary file which provides the - potential energy boundaries (min E_pot, max E_pot) tabulated - for all degrees of association. Make sure to only list the - degrees of association which are used by the degree of - association collective variable within this file. The energy - boundary file can be created in a preliminary energy run. By - the help of the functions - :meth:`update_maximum_and_minimum_energies_at_current_state` - and :meth:`write_out_preliminary_energy_run_results`. This - file has to be obtained before being able to run a - simulation with the energy as collective variable. - delta : :obj:`float` - Provides the discretization of the potential energy range. Only - for small enough delta the results of the energy reweighted - averages are correct. If delta is chosen too big there are - discretization errors in the numerical integration which occurs - during the energy reweighting process. - - """ - for k in kwargs: - if k in self._valid_keys_add_collective_variable_potential_energy(): - self._params[k] = kwargs[k] - else: - raise KeyError(f"{k} is not a valid key") - - for k in self._required_keys_add_collective_variable_potential_energy(): - if k not in kwargs: - raise ValueError( - "At least the following keys have to be given as keyword arguments: " + self._required_keys_add_collective_variable_degree_of_association().__str__() + " got " + kwargs.__str__()) - self._params[k] = kwargs[k] - filname_potential_energy_boundaries_file = self._params[ - "filename"].encode("utf-8") - deref(self.WLRptr).add_new_CV_potential_energy( - filname_potential_energy_boundaries_file, self._params["delta"]) - - def _valid_keys_add_collective_variable_potential_energy(self): - return "filename", "delta" - - def _required_keys_add_collective_variable_potential_energy(self): - return "filename", "delta" - - def set_wang_landau_parameters(self, *args, **kwargs): - """ - Sets the final Wang-Landau parameter. - - Parameters - ---------- - final_wang_landau_parameter : :obj:`float` - Sets the final Wang-Landau parameter, which is the Wang-Landau - parameter after which the simulation should stop. - full_path_to_output_filename : :obj:`str` - Sets the path to the output file of the Wang-Landau algorithm which - contains the Wang-Landau potential - do_not_sample_reaction_partition_function : :obj:`bool` - Avoids sampling the Reaction ensemble partition function in the - Wang-Landau algorithm. Therefore this option makes all degrees of - association equally probable. This option may be used in the - sweeping mode of the reaction ensemble, since the reaction ensemble - partition function can be later added analytically. - - """ - for k in kwargs: - if k in self._valid_keys_set_wang_landau_parameters(): - self._params[k] = kwargs[k] - else: - raise KeyError(f"{k} is not a valid key") - - deref(self.WLRptr).final_wang_landau_parameter = self._params[ - "final_wang_landau_parameter"] - deref(self.WLRptr).output_filename = self._params[ - "full_path_to_output_filename"].encode("utf-8") - deref(self.WLRptr).do_not_sample_reaction_partition_function = self._params[ - "do_not_sample_reaction_partition_function"] - - def _valid_keys_set_wang_landau_parameters(self): - return "final_wang_landau_parameter", "full_path_to_output_filename", "do_not_sample_reaction_partition_function" - - def load_wang_landau_checkpoint(self): - """ - Loads the dumped Wang-Landau potential file. - - """ - checkpoint_name = "checkpoint".encode("utf-8") - deref(self.WLRptr).load_wang_landau_checkpoint(checkpoint_name) - - def write_wang_landau_checkpoint(self): - """ - Dumps the Wang-Landau potential to a checkpoint file. Can be used to - checkpoint the Wang-Landau histogram, potential, parameter and the - number of executed trial moves. - - """ - checkpoint_name = "checkpoint".encode("utf-8") - deref(self.WLRptr).write_wang_landau_checkpoint(checkpoint_name) - - def update_maximum_and_minimum_energies_at_current_state(self): - """ - Records the minimum and maximum potential energy as a function of the - degree of association in a preliminary Wang-Landau reaction ensemble - simulation where the acceptance probability includes the factor - :math:`\\exp(-\\beta \\Delta E_{pot})`. The minimal and maximal - potential energies which occur in the system are needed for the energy - reweighting simulations where the factor :math:`\\exp(-\\beta \\Delta E_{pot})` - is not included in the acceptance probability in - order to avoid choosing the wrong potential energy boundaries. - - """ - self.WLRptr.get( - ).update_maximum_and_minimum_energies_at_current_state() - - def write_out_preliminary_energy_run_results(self): - """ - This writes out the minimum and maximum potential energy as a function - of the degree of association to a file. It requires that previously - :meth:`update_maximum_and_minimum_energies_at_current_state` was used. - - """ - filename = "preliminary_energy_run_results".encode("utf-8") - deref(self.WLRptr).write_out_preliminary_energy_run_results(filename) - - def write_wang_landau_results_to_file(self, filename): - """ - This writes out the Wang-Landau potential as a function of the used - collective variables. - - """ - deref(self.WLRptr).write_wang_landau_results_to_file( - filename.encode("utf-8")) - - def displacement_mc_move_for_particles_of_type(self, type_mc, - particle_number_to_be_changed=1): - """ - Performs an MC (Monte Carlo) move for ``particle_number_to_be_changed`` - particle of type ``type_mc``. Positions for the particles are drawn - uniformly and randomly within the box. The command takes into account - the Wang-Landau terms in the acceptance probability. - If there are multiple types, that need to be moved, make sure to move - them in a random order to avoid artefacts. For the Wang-Landau algorithm - in the case of energy reweighting you would also need to move the - monomers of the polymer with special moves for the MC part. Those - polymer configuration-changing moves need to be implemented in the - case of using Wang-Landau with energy reweighting and a polymer in the - system. Polymer configuration-changing moves had been implemented - before but were removed from ESPResSo. - - """ - use_wang_landau = True - deref(self.WLRptr).do_global_mc_move_for_particles_of_type( - type_mc, particle_number_to_be_changed, use_wang_landau) - - cdef class WidomInsertion(ReactionAlgorithm): """ This class implements the Widom insertion method in the canonical ensemble diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index a20c73ee4bb..09329ec2c98 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -175,8 +175,6 @@ python_test(FILE analyze_distance.py MAX_NUM_PROC 1) python_test(FILE analyze_acf.py MAX_NUM_PROC 1) python_test(FILE comfixed.py MAX_NUM_PROC 2) python_test(FILE rescale.py MAX_NUM_PROC 2) -python_test(FILE wang_landau.py MAX_NUM_PROC 1) -python_test(FILE wang_landau_stats.py MAX_NUM_PROC 1 LABELS long) python_test(FILE array_properties.py MAX_NUM_PROC 4) python_test(FILE analyze_distribution.py MAX_NUM_PROC 1) python_test(FILE observable_profile.py MAX_NUM_PROC 4) diff --git a/testsuite/python/wang_landau.py b/testsuite/python/wang_landau.py deleted file mode 100644 index b13010dc52a..00000000000 --- a/testsuite/python/wang_landau.py +++ /dev/null @@ -1,131 +0,0 @@ -# -# Copyright (C) 2013-2019 The ESPResSo project -# -# This file is part of ESPResSo. -# -# ESPResSo is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ESPResSo is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -"""Testmodule for the Wang-Landau Reaction Ensemble. -""" -import numpy as np -import unittest as ut - -import espressomd -import espressomd.interactions -import espressomd.reaction_ensemble - - -class WangLandauReactionEnsembleTest(ut.TestCase): - - """Test the interface of the Wang-Landau reaction ensemble.""" - - # System parameters - # - box_l = 6 * np.sqrt(2) - temperature = 1.0 - - # Integration parameters - # - system = espressomd.System(box_l=[box_l, box_l, box_l]) - np.random.seed(seed=42) - system.time_step = 0.01 - system.cell_system.skin = 0 - system.cell_system.set_n_square(use_verlet_lists=False) - - # - # Setup System - # - - K_diss = 0.0088 - - p1 = system.part.add(type=3, pos=[0, 0, 0]) - p2 = system.part.add(type=1, pos=system.box_l / 2.0) - system.part.add(type=2, pos=np.random.random() * system.box_l) - system.part.add(type=2, pos=np.random.random() * system.box_l) - - harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0, k=1) - system.bonded_inter[0] = harmonic_bond - p1.add_bond((harmonic_bond, p2)) - - WLRE = espressomd.reaction_ensemble.WangLandauReactionEnsemble( - temperature=temperature, exclusion_radius=0, seed=86) - WLRE.add_reaction( - gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], - product_types=[1, 2], product_coefficients=[1, 1], - default_charges={0: 0, 1: -1, 2: +1}) - system.setup_type_map([0, 1, 2, 3]) - # initialize wang_landau - file_input = "energy_boundaries.dat" - file_output = "WL_potential_out.dat" - # generate preliminary_energy_run_results here, this should be done in a - # separate simulation without energy reweighting using the update energy - # functions - np.savetxt(file_input, np.transpose([[0, 1], [0, 0], [9, 9]]), - delimiter='\t', header="nbar E_potmin E_potmax") - - WLRE.add_collective_variable_degree_of_association( - associated_type=0, min=0, max=1, corresponding_acid_types=[0, 1]) - WLRE.set_wang_landau_parameters( - final_wang_landau_parameter=0.8 * 1e-2, - do_not_sample_reaction_partition_function=True, - full_path_to_output_filename=file_output) - - def test_01_energy_recording(self): - self.WLRE.update_maximum_and_minimum_energies_at_current_state() - self.WLRE.write_out_preliminary_energy_run_results() - nbars, E_mins, E_maxs = np.loadtxt( - "preliminary_energy_run_results", unpack=True) - np.testing.assert_almost_equal(nbars, [0, 1]) - np.testing.assert_almost_equal(E_mins, [27.0, -10]) - np.testing.assert_almost_equal(E_maxs, [27.0, -10]) - - def check_checkpoint(self, filename): - # write first checkpoint - self.WLRE.write_wang_landau_checkpoint() - old_checkpoint = np.loadtxt(filename) - - # modify old_checkpoint in memory and in file (this destroys the - # information contained in the checkpoint, but allows for testing of - # the functions) - modified_checkpoint = old_checkpoint - modified_checkpoint[0] = 1 - np.savetxt(filename, modified_checkpoint) - - # check whether changes are carried out correctly - self.WLRE.load_wang_landau_checkpoint() - self.WLRE.write_wang_landau_checkpoint() - new_checkpoint = np.loadtxt(filename) - np.testing.assert_almost_equal(new_checkpoint, modified_checkpoint) - - def test_02_checkpointing(self): - self.WLRE.add_collective_variable_potential_energy( - filename=self.file_input, delta=0.05) - - # run MC for long enough to sample a non-trivial histogram - for _ in range(1000): - try: - self.WLRE.reaction() - for _ in range(2): - self.WLRE.displacement_mc_move_for_particles_of_type(3) - except espressomd.reaction_ensemble.WangLandauHasConverged: - break - - filenames = ["checkpoint_wang_landau_potential_checkpoint", - "checkpoint_wang_landau_histogram_checkpoint"] - for filename in filenames: - self.check_checkpoint(filename) - - -if __name__ == "__main__": - ut.main() diff --git a/testsuite/python/wang_landau_stats.py b/testsuite/python/wang_landau_stats.py deleted file mode 100644 index 757a04f1975..00000000000 --- a/testsuite/python/wang_landau_stats.py +++ /dev/null @@ -1,128 +0,0 @@ -# -# 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 - - p1 = system.part.add(type=3, pos=[0, 0, 0]) - p2 = system.part.add(type=1, pos=system.box_l / 2.0) - system.part.add(type=2, pos=np.random.random() * system.box_l) - system.part.add(type=2, pos=np.random.random() * system.box_l) - - harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0, k=1) - system.bonded_inter[0] = harmonic_bond - p1.add_bond((harmonic_bond, p2)) - - 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() diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt index 382e83b8007..90a7c00b95d 100644 --- a/testsuite/scripts/samples/CMakeLists.txt +++ b/testsuite/scripts/samples/CMakeLists.txt @@ -78,7 +78,6 @@ sample_test(FILE test_visualization_ljliquid.py) sample_test(FILE test_visualization_elc.py) sample_test(FILE test_visualization_npt.py) sample_test(FILE test_visualization_poiseuille.py) -sample_test(FILE test_wang_landau_reaction_ensemble.py) sample_test(FILE test_widom_insertion.py) sample_test(FILE test_object_in_fluid__motivation.py) sample_test(FILE test_immersed_boundary.py) diff --git a/testsuite/scripts/samples/test_wang_landau_reaction_ensemble.py b/testsuite/scripts/samples/test_wang_landau_reaction_ensemble.py deleted file mode 100644 index ab3704e2a8d..00000000000 --- a/testsuite/scripts/samples/test_wang_landau_reaction_ensemble.py +++ /dev/null @@ -1,41 +0,0 @@ -# Copyright (C) 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 importlib_wrapper - - -def shorten_loop(code): - # stop reaction before the exception is raised - breakpoint = "while True:" - assert breakpoint in code - code = code.replace(breakpoint, "for _ in range(6):", 1) - return code - - -sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( - "@SAMPLES_DIR@/wang_landau_reaction_ensemble.py", - substitutions=shorten_loop) - - -@skipIfMissingFeatures -class Sample(ut.TestCase): - system = sample.system - - -if __name__ == "__main__": - ut.main()