diff --git a/doc/doxygen/bibliography.bib b/doc/doxygen/bibliography.bib index fb48a3b3e8f..fc47479ac01 100644 --- a/doc/doxygen/bibliography.bib +++ b/doc/doxygen/bibliography.bib @@ -390,6 +390,17 @@ @Article{ladd94a doi = {10.1017/S0022112094001771}, } +@Article{landsgesell17b, + author = {Landsgesell, Jonas and Holm, Christian and Smiatek, Jens}, + title = {{Simulation of weak polyelectrolytes: A comparison between the constant pH and the reaction ensemble method}}, + journal = {European Physical Journal Special Topics}, + year = {2017}, + volume = {226}, + number = {4}, + pages = {725-736}, + doi = {10.1140/epjst/e2016-60324-3}, +} + @Article{marsili10a, author = {Marsili, Simone and Signorini, Giorgio Federico and Chelli, Riccardo and Marchi, Massimo and Procacci, Piero}, title = {{{ORAC}: A molecular dynamics simulation program to explore free energy surfaces in biomolecular systems at the atomistic level}}, diff --git a/src/core/reaction_methods/ConstantpHEnsemble.cpp b/src/core/reaction_methods/ConstantpHEnsemble.cpp index 21427196ce1..82fa8a9a042 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.cpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.cpp @@ -21,83 +21,31 @@ #include "particle_data.hpp" +#include "utils.hpp" + #include #include #include namespace ReactionMethods { -int ConstantpHEnsemble::get_random_valid_p_id() { - int random_p_id = i_random(get_maximal_particle_id() + 1); - // draw random p_ids till we draw a pid which exists - while (not particle_exists(random_p_id)) - random_p_id = i_random(get_maximal_particle_id() + 1); - return random_p_id; -} - -/** - *Performs a reaction in the constant pH ensemble - */ -int ConstantpHEnsemble::do_reaction(int reaction_steps) { - - for (int i = 0; i < reaction_steps; ++i) { - // get a list of reactions where a randomly selected particle type occurs in - // the reactant list. the selection probability of the particle types has to - // be proportional to the number of occurrences of the number of particles - // with - // this type - - // for optimizations this list could be determined during the initialization - std::vector list_of_reaction_ids_with_given_reactant_type; - while (list_of_reaction_ids_with_given_reactant_type - .empty()) { // avoid selecting a (e.g. salt) particle which - // does not take part in a reaction - int random_p_id = get_random_valid_p_id(); // only used to determine which - // reaction is attempted. - auto part = get_particle_data(random_p_id); - - int type_of_random_p_id = part.p.type; - - // construct list of reactions with the above reactant type - for (int reaction_i = 0; reaction_i < reactions.size(); reaction_i++) { - SingleReaction ¤t_reaction = reactions[reaction_i]; - for (int reactant_i = 0; reactant_i < 1; - reactant_i++) { // reactant_i < 1 since it is assumed in this place - // that the types A and HA occur in the first place - // only. These are the types that should be - // switched, H+ should not be switched - if (current_reaction.reactant_types[reactant_i] == - type_of_random_p_id) { - list_of_reaction_ids_with_given_reactant_type.push_back(reaction_i); - break; - } - } - } - } - // randomly select a reaction to be performed - int reaction_id = - list_of_reaction_ids_with_given_reactant_type[i_random(static_cast( - list_of_reaction_ids_with_given_reactant_type.size()))]; - generic_oneway_reaction(reaction_id); - } - return 0; -} - /** * Calculates the expression in the acceptance probability of the constant pH * method. */ double ConstantpHEnsemble::calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, - std::map const &dummy_old_particle_numbers, - int dummy_old_state_index, int dummy_new_state_index, + std::map const &old_particle_numbers, int dummy_old_state_index, + int dummy_new_state_index, bool dummy_only_make_configuration_changing_move) const { auto const beta = 1.0 / temperature; auto const pKa = -current_reaction.nu_bar * log10(current_reaction.gamma); auto const ln_bf = (E_pot_new - E_pot_old) - current_reaction.nu_bar / beta * log(10) * (m_constant_pH - pKa); - auto const bf = exp(-beta * ln_bf); + const double factorial_expr = calculate_factorial_expression_cpH( + current_reaction, old_particle_numbers); + auto const bf = factorial_expr * exp(-beta * ln_bf); return bf; } diff --git a/src/core/reaction_methods/ConstantpHEnsemble.hpp b/src/core/reaction_methods/ConstantpHEnsemble.hpp index fb088d9aa36..a2bdb7668b0 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.hpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.hpp @@ -42,17 +42,13 @@ class ConstantpHEnsemble : public ReactionAlgorithm { public: ConstantpHEnsemble(int seed) : ReactionAlgorithm(seed) {} double m_constant_pH = -10; - int do_reaction(int reaction_steps) override; protected: double calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, - double E_pot_new, std::map const &dummy_old_particle_numbers, + double E_pot_new, std::map const &old_particle_numbers, int dummy_old_state_index, int dummy_new_state_index, bool dummy_only_make_configuration_changing_move) const override; - -private: - int get_random_valid_p_id(); }; } // namespace ReactionMethods diff --git a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp index e55cae694d2..340d1be5ed4 100644 --- a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp @@ -71,11 +71,14 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { auto const p_numbers = std::map{{type_A, i}, {type_B, j}, {type_C, k}}; auto const energy = static_cast(i + 1); - // acceptance = exp(- E / T + nu_bar * log(10) * (pH - nu_bar * pKa)) + auto const f_expr = + calculate_factorial_expression_cpH(reaction, p_numbers); + // acceptance = f_expr * exp(- E / T + nu_bar * log(10) * (pH - nu_bar * + // pKa)) auto const acceptance_ref = - std::exp(energy / r_algo.temperature + - std::log(10.) * - (r_algo.m_constant_pH + std::log10(reaction.gamma))); + f_expr * std::exp(energy / r_algo.temperature + + 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); BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); diff --git a/src/core/reaction_methods/utils.cpp b/src/core/reaction_methods/utils.cpp index 145f90e57fe..65154e66d1f 100644 --- a/src/core/reaction_methods/utils.cpp +++ b/src/core/reaction_methods/utils.cpp @@ -66,4 +66,25 @@ calculate_factorial_expression(SingleReaction const ¤t_reaction, return factorial_expr; } +double calculate_factorial_expression_cpH( + SingleReaction const ¤t_reaction, + std::map const &old_particle_numbers) { + double factorial_expr = 1.0; + // factorial contribution of reactants + { + int nu_i = -1 * current_reaction.reactant_coefficients[0]; + int N_i0 = old_particle_numbers.at(current_reaction.reactant_types[0]); + factorial_expr *= + factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); + } + // factorial contribution of products + { + int nu_i = current_reaction.product_coefficients[0]; + int N_i0 = old_particle_numbers.at(current_reaction.product_types[0]); + factorial_expr *= + factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N_i0, nu_i); + } + return factorial_expr; +} + } // namespace ReactionMethods diff --git a/src/core/reaction_methods/utils.hpp b/src/core/reaction_methods/utils.hpp index 4147ef8b549..623c6126c17 100644 --- a/src/core/reaction_methods/utils.hpp +++ b/src/core/reaction_methods/utils.hpp @@ -38,6 +38,16 @@ double calculate_factorial_expression(SingleReaction const ¤t_reaction, std::map const &old_particle_numbers); +/** + * Calculates the factorial expression which occurs in the constant pH method + * with symmetric proposal probability. + * + * See @cite landsgesell17b + */ +double calculate_factorial_expression_cpH( + SingleReaction const ¤t_reaction, + std::map const &old_particle_numbers); + /** * Calculates the factorial expression which occurs in the reaction ensemble * acceptance probability diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index c2c3dbf0388..a6bb757e481 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -431,6 +431,14 @@ cdef class ReactionEnsemble(ReactionAlgorithm): self._set_params_in_es_core() cdef class ConstantpHEnsemble(ReactionAlgorithm): + """ + This class implements the constant pH Ensemble. + + When adding an acid-base reaction, the acid and base particle types + are always assumed to be at index 0 of the lists passed to arguments + ``reactant_types`` and ``product_types``. + + """ cdef unique_ptr[CConstantpHEnsemble] constpHptr def __init__(self, *args, **kwargs):