introduce symmetric proposal probability in cpH method

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