docs: Document constant pH implementation details

espressomd · May 18, 2021 · a37938e · a37938e
1 parent 4ced2c3
commit a37938e
Show file tree

Hide file tree

Showing 4 changed files with 26 additions and 7 deletions.
diff --git a/doc/doxygen/bibliography.bib b/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}},

diff --git a/src/core/reaction_methods/utils.cpp b/src/core/reaction_methods/utils.cpp
@@ -71,16 +71,16 @@ double calculate_factorial_expression_cpH(
     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]);
+  {
+    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
-  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]);
+  {
+    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);
   }

diff --git a/src/core/reaction_methods/utils.hpp b/src/core/reaction_methods/utils.hpp
@@ -42,7 +42,7 @@ calculate_factorial_expression(SingleReaction const &current_reaction,
  * Calculates the factorial expression which occurs in the constant pH method
  * with symmetric proposal probability.
  *
- * See https://arxiv.org/abs/1702.04911 equation 15
+ * See @cite landsgesell17b
  */
 double calculate_factorial_expression_cpH(
     SingleReaction const &current_reaction,

diff --git a/src/python/espressomd/reaction_ensemble.pyx b/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):