diff --git a/src/core/reaction_methods/ConstantpHEnsemble.cpp b/src/core/reaction_methods/ConstantpHEnsemble.cpp index 21427196ce1..4a496962a78 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.cpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.cpp @@ -21,75 +21,21 @@ #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, + 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; @@ -97,7 +43,9 @@ double ConstantpHEnsemble::calculate_acceptance_probability( 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/utils.cpp b/src/core/reaction_methods/utils.cpp index 145f90e57fe..9955940d750 100644 --- a/src/core/reaction_methods/utils.cpp +++ b/src/core/reaction_methods/utils.cpp @@ -66,4 +66,26 @@ 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 + for (int i = 0; i < 1; i++) { + int nu_i = -1 * current_reaction.reactant_coefficients[i]; + int N_i0 = old_particle_numbers.at(current_reaction.reactant_types[i]); + factorial_expr *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i( + N_i0, nu_i); + } + // factorial contribution of products + for (int i = 0; i < 1; i++) { + int nu_i = current_reaction.product_coefficients[i]; + int N_i0 = old_particle_numbers.at(current_reaction.product_types[i]); + factorial_expr *= factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i( + N_i0, nu_i); + } + return factorial_expr; +} + + } // namespace ReactionMethods diff --git a/src/core/reaction_methods/utils.hpp b/src/core/reaction_methods/utils.hpp index 4147ef8b549..b3ab6df828f 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 https://arxiv.org/abs/1702.04911 equation 15 + */ +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