introduce symmetric proposal probability

src/core/reaction_methods/ConstantpHEnsemble.cpp

@@ -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,
-    int dummy_old_state_index, int dummy_new_state_index,
+    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;
 }
 

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 &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

src/core/reaction_methods/utils.cpp

@@ -66,4 +66,25 @@ 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

src/core/reaction_methods/utils.hpp

@@ -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