Skip to content

Commit

Permalink
introduce symmetric proposal probability
Browse files Browse the repository at this point in the history
  • Loading branch information
Landsgesell Jonas (XC-DA/EPB4) authored and jonaslandsgesell committed Apr 6, 2021
1 parent f1f6b69 commit 2604465
Show file tree
Hide file tree
Showing 4 changed files with 39 additions and 63 deletions.
64 changes: 6 additions & 58 deletions src/core/reaction_methods/ConstantpHEnsemble.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -21,83 +21,31 @@

#include "particle_data.hpp"

#include "utils.hpp"

#include <cmath>
#include <map>
#include <vector>

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<int> 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 &current_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<int>(
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 &current_reaction, double E_pot_old, double E_pot_new,
std::map<int, int> const &dummy_old_particle_numbers,
std::map<int, int> 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;
}

Expand Down
6 changes: 1 addition & 5 deletions src/core/reaction_methods/ConstantpHEnsemble.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 &current_reaction, double E_pot_old,
double E_pot_new, std::map<int, int> const &dummy_old_particle_numbers,
double E_pot_new, std::map<int, int> 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
Expand Down
22 changes: 22 additions & 0 deletions src/core/reaction_methods/utils.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,4 +66,26 @@ calculate_factorial_expression(SingleReaction const &current_reaction,
return factorial_expr;
}

double
calculate_factorial_expression_cpH(SingleReaction const &current_reaction,
std::map<int, int> 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
10 changes: 10 additions & 0 deletions src/core/reaction_methods/utils.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,16 @@ double
calculate_factorial_expression(SingleReaction const &current_reaction,
std::map<int, int> 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 &current_reaction,
std::map<int, int> const &old_particle_numbers);

/**
* Calculates the factorial expression which occurs in the reaction ensemble
* acceptance probability
Expand Down

0 comments on commit 2604465

Please sign in to comment.