introduce symmetric proposal probability in cpH method (#4207)

Fixes #4197

Description of changes:
- Cleanup constant pH code via removal of asymmetric proposal probabilities, now implements constant pH method with symmetric proposal probability in Landsgesell et al. 2017 (doi:10.1140/epjst/e2016-60324-3), equation 15
- Fix infinite loop when the H+ and A- types are accidentally switched when defining the acid-base reaction

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}},

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/tests/ConstantpHEnsemble_test.cpp

@@ -71,11 +71,14 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) {
         auto const p_numbers =
             std::map<int, int>{{type_A, i}, {type_B, j}, {type_C, k}};
         auto const energy = static_cast<double>(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);

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
+  {
+    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

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 @cite landsgesell17b
+ */
+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

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):