diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py
index 98f089f438f..923bde2f5dd 100644
--- a/samples/grand_canonical.py
+++ b/samples/grand_canonical.py
@@ -17,10 +17,19 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 """
-This sample performs a grand canonical simulation of a salt solution.
+This example script performs a grand canonical simulation of a system in contact
+with a salt reservoir and ensures constant chemical potential.
+It takes two command line arguments as input: 1) the reservoir salt concentration
+in units of 1/sigma^3 and 2) the excess chemical potential of the reservoir in
+units of kT.
+The excess chemical potential of the reservoir needs to be determined prior to
+running the grand canonical simulation using the script called widom_insertion.py
+which simulates a part of the reservoir at the prescribed salt concentration.
+Be aware that the reservoir excess chemical potential depends on all interactions
+in the reservoir system.
 """
 import numpy as np
-import sys
+import argparse
 
 import espressomd
 from espressomd import reaction_ensemble
@@ -29,18 +38,19 @@
 required_features = ["P3M", "EXTERNAL_FORCES", "LENNARD_JONES"]
 espressomd.assert_features(required_features)
 
-# print help message if proper command-line arguments are not provided
-if len(sys.argv) != 3:
-    print("\nGot ", str(len(sys.argv) - 1), " arguments, need 2\n\nusage:" +
-          sys.argv[0] + " [cs_bulk] [excess_chemical_potential/kT]\n")
-    sys.exit()
+parser = argparse.ArgumentParser(epilog=__doc__)
+parser.add_argument('cs_bulk', type=float,
+                    help="bulk salt concentration [1/sigma^3]")
+parser.add_argument('excess_chemical_potential', type=float,
+                    help="excess chemical potential [kT] "
+                         "(obtained from Widom's insertion method)")
+args = parser.parse_args()
 
 # System parameters
 #############################################################
 
-cs_bulk = float(sys.argv[1])
-excess_chemical_potential_pair = float(
-    sys.argv[2])  # from widom insertion simulation of pair insertion
+cs_bulk = args.cs_bulk
+excess_chemical_potential_pair = args.excess_chemical_potential
 box_l = 50.0
 
 # Integration parameters
@@ -54,7 +64,7 @@
 system.cell_system.skin = 0.4
 temperature = 1.0
 system.thermostat.set_langevin(kT=temperature, gamma=.5, seed=42)
-system.cell_system.max_num_cells = 2744
+system.cell_system.max_num_cells = 14**3
 
 
 #############################################################
@@ -78,8 +88,7 @@
 for type_1 in types:
     for type_2 in types:
         system.non_bonded_inter[type_1, type_2].lennard_jones.set_params(
-            epsilon=lj_eps, sigma=lj_sig,
-            cutoff=lj_cut, shift="auto")
+            epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
 
 RE = reaction_ensemble.ReactionEnsemble(
     temperature=temperature, exclusion_radius=2.0, seed=3)
@@ -96,8 +105,8 @@
 p3m = electrostatics.P3M(prefactor=2.0, accuracy=1e-3)
 system.actors.add(p3m)
 p3m_params = p3m.get_params()
-for key in list(p3m_params.keys()):
-    print("{} = {}".format(key, p3m_params[key]))
+for key, value in p3m_params.items():
+    print("{} = {}".format(key, value))
 
 
 # Warmup
@@ -105,27 +114,23 @@
 # warmup integration (with capped LJ potential)
 warm_steps = 1000
 warm_n_times = 20
-# do the warmup until the particles have at least the distance min_dist
 # set LJ cap
-lj_cap = 20
-system.force_cap = lj_cap
+system.force_cap = 20
 
 # Warmup Integration Loop
-act_min_dist = system.analysis.min_dist()
 i = 0
-while (i < warm_n_times):
+while i < warm_n_times:
     print(i, "warmup")
     RE.reaction(100)
     system.integrator.run(steps=warm_steps)
     i += 1
-    #Increase LJ cap
-    lj_cap = lj_cap + 10
-    system.force_cap = lj_cap
+    # increase LJ cap
+    system.force_cap += 10
 
 # remove force capping
 system.force_cap = 0
 
-#MC warmup
+# MC warmup
 RE.reaction(1000)
 
 n_int_cycles = 10000
diff --git a/samples/reaction_ensemble.py b/samples/reaction_ensemble.py
index 9d3a615d19b..6f8852e4ef4 100644
--- a/samples/reaction_ensemble.py
+++ b/samples/reaction_ensemble.py
@@ -17,8 +17,10 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 """
-This sample simulates the reaction ensemble. It also illustrates how the
-constant pH method can be used.
+This sample illustrates how to use either the raction ensemble or the constant pH ensemble. You can choose
+in which ensemble you want to simulate via either providing --reaction_ensemble or --constant_pH_ensemble as command line argument to the script.
+Be aware that in the case of the reaction ensemble the dissociation constant gamma is not the thermodynamic reaction constant K, but rather K*1mol/l and therefore carries a unit!.
+In the case of the of the constant pH method gamma is the thermodynamic reaction constant!
 """
 import numpy as np
 import argparse
diff --git a/samples/wang_landau_reaction_ensemble.py b/samples/wang_landau_reaction_ensemble.py
index ab2792a61c3..a0f329bfadc 100644
--- a/samples/wang_landau_reaction_ensemble.py
+++ b/samples/wang_landau_reaction_ensemble.py
@@ -17,7 +17,15 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 """
-This sample simulates the Wang-Landau Reaction Ensemble for a harmonic bond.
+This example script simulates two reacting monomers which are bonded via a harmonic potential.
+The script aborts as soon as the abortion criterion in the Wang-Landau algorithm is met.
+The Wang-Landau simulation runs until the Wang-Landau potential is converged and then raises a Warning that it has converged, effectively aborting the simulation.
+With the setup of the Wang-Landau algorithm in this script you sample the density of states of a three dimensional reacting harmonic oscillator 
+as a function of the two collective variables 1) degree of association and 2) potential energy.
+The recorded Wang-Landau potential (which is updated during the simulation) is written to the file WL_potential_out.dat 
+In this simulation setup the Wang-Landau potential is the density of states. You can view the converged Wang-Landau potential e.g. via plotting with gnuplot: splot "WL_potential_out.dat". 
+As expected the three dimensional harmonic oscilltor has a density of states which goes like sqrt(E_pot).
+For a scientific description and different ways to use the algorithm please consult https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00791 
 """
 import numpy as np
 
diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py
index df832981c80..da01deecfe2 100644
--- a/samples/widom_insertion.py
+++ b/samples/widom_insertion.py
@@ -17,11 +17,13 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 """
-This sample measures the excess chemical potential for Widom insertion of
-charged particles using the reaction ensemble method.
+This example script measures the excess chemical potential of a charged WCA
+fluid via Widom's insertion method.
+As input this script requires you to provide particle number density in units
+of 1/sigma^3.
 """
 import numpy as np
-import sys
+import argparse
 
 import espressomd
 from espressomd import code_info
@@ -33,13 +35,16 @@
 required_features = ["LENNARD_JONES", "P3M"]
 espressomd.assert_features(required_features)
 
+parser = argparse.ArgumentParser(epilog=__doc__)
+parser.add_argument('cs_bulk', type=float,
+                    help="bulk salt concentration [1/sigma^3]")
+args = parser.parse_args()
+
 # System parameters
 #############################################################
-assert len(sys.argv) == 2, "please provide a value for cs_bulk"
-cs_bulk = float(sys.argv[1])
-box_l = 50.0
-N0 = int(cs_bulk * box_l**3)
-print("actual cs_bulk", float(N0) / box_l**3)
+cs_bulk = args.cs_bulk
+N0 = 70
+box_l = (N0 / cs_bulk)**(1.0 / 3.0)
 
 # Integration parameters
 #############################################################
@@ -51,7 +56,7 @@
 system.cell_system.skin = 0.4
 temperature = 1.0
 system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42)
-system.cell_system.max_num_cells = 2744
+system.cell_system.max_num_cells = 14**3
 
 
 #############################################################
@@ -76,41 +81,39 @@
 for type_1 in types:
     for type_2 in types:
         system.non_bonded_inter[type_1, type_2].lennard_jones.set_params(
-            epsilon=lj_eps, sigma=lj_sig,
-            cutoff=lj_cut, shift="auto")
+            epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
 
-p3m = electrostatics.P3M(prefactor=0.9, accuracy=1e-3)
+p3m = electrostatics.P3M(prefactor=2.0, accuracy=1e-3)
 system.actors.add(p3m)
 p3m_params = p3m.get_params()
-for key in list(p3m_params.keys()):
-    print("{} = {}".format(key, p3m_params[key]))
+for key, value in p3m_params.items():
+    print("{} = {}".format(key, value))
 
 # Warmup
 #############################################################
 # warmup integration (with capped LJ potential)
 warm_steps = 1000
 warm_n_times = 20
-# do the warmup until the particles have at least the distance min_dist
 # set LJ cap
-lj_cap = 20
-system.force_cap = lj_cap
+system.force_cap = 20
 
 # Warmup Integration Loop
-act_min_dist = system.analysis.min_dist()
 i = 0
-while (i < warm_n_times):
+while i < warm_n_times:
     print(i, "warmup")
     system.integrator.run(steps=warm_steps)
     i += 1
-    # Increase LJ cap
-    lj_cap = lj_cap + 10
-    system.force_cap = lj_cap
+    # increase LJ cap
+    system.force_cap += 10
 
 # remove force capping
 system.force_cap = 0
 
 RE = reaction_ensemble.WidomInsertion(
     temperature=temperature, seed=77)
+
+# add insertion reaction
+insertion_reaction_id = 0
 RE.add_reaction(reactant_types=[],
                 reactant_coefficients=[], product_types=[1, 2],
                 product_coefficients=[1, 1], default_charges={1: -1, 2: +1})
@@ -119,16 +122,15 @@
 
 n_iterations = 100
 for i in range(n_iterations):
-    for j in range(30):
-        RE.measure_excess_chemical_potential(0)  # 0 for insertion reaction
-        system.integrator.run(steps=2)
+    for j in range(50):
+        RE.measure_excess_chemical_potential(insertion_reaction_id)
     system.integrator.run(steps=500)
     if i % 20 == 0:
         print("mu_ex_pair ({:.4f}, +/- {:.4f})".format(
-            *RE.measure_excess_chemical_potential(0)  # 0 for insertion reaction
-        ))
+            *RE.measure_excess_chemical_potential(insertion_reaction_id)))
         print("HA", system.number_of_particles(type=0), "A-",
               system.number_of_particles(type=1), "H+",
               system.number_of_particles(type=2))
 
-print(RE.measure_excess_chemical_potential(0))  # 0 for insertion reaction
+print("excess chemical potential for an ion pair ",
+      RE.measure_excess_chemical_potential(insertion_reaction_id))
diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp
index 84d7c0210c9..ec63e5ac772 100644
--- a/src/core/reaction_ensemble.cpp
+++ b/src/core/reaction_ensemble.cpp
@@ -16,14 +16,6 @@ GNU General Public License for more details.
 You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
-/** reaction ensemble method according to smith94x for the reaction ensemble at
- *constant volume and temperature, for the reaction ensemble at constant
- *pressure additionally employ a barostat! NOTE: a chemical reaction consists of
- *a forward and backward reaction. Here both reactions have to be defined
- *separately. The extent of the reaction is here chosen to be +1. If the
- *reaction trial move for a dissociation of HA is accepted then there is one
- *more dissociated ion pair H+ and A-
- */
 
 /** @file */
 
@@ -66,7 +58,7 @@ double average_list_of_allowed_entries(const std::vector<T> &vector) {
 }
 
 /**
- * Checks whether a number is in a std:vector of numbers.
+ * Checks whether a number is in a std::vector of numbers.
  */
 template <typename T> bool is_in_list(T value, const std::vector<T> &list) {
   for (int i = 0; i < list.size(); i++) {
@@ -76,10 +68,11 @@ template <typename T> bool is_in_list(T value, const std::vector<T> &list) {
   return false;
 }
 
+/** Save minimum and maximum energies as a function of the other collective
+ *  variables under min_boundaries_energies, max_boundaries_energies
+ */
 void EnergyCollectiveVariable::load_CV_boundaries(
     WangLandauReactionEnsemble &m_current_wang_landau_system) {
-  /**save minimum and maximum energies as a function of the other collective
-   * variables under min_boundaries_energies, max_boundaries_energies **/
 
   m_current_wang_landau_system.do_energy_reweighting = true;
   // load energy boundaries from file
@@ -164,7 +157,6 @@ void ReactionAlgorithm::add_reaction(
  * set.
  */
 void ReactionAlgorithm::check_reaction_ensemble() {
-  /**checks the reaction_ensemble struct for valid parameters */
   if (reactions.empty()) {
     throw std::runtime_error("Reaction system not initialized");
   }
@@ -369,7 +361,6 @@ double ReactionEnsemble::calculate_acceptance_probability(
     std::map<int, int> &old_particle_numbers, int dummy_old_state_index,
     int dummy_new_state_index,
     bool dummy_only_make_configuration_changing_move) {
-  /**calculate the acceptance probability in the reaction ensemble */
   const double factorial_expr =
       calculate_factorial_expression(current_reaction, old_particle_numbers);
 
@@ -418,20 +409,19 @@ void WangLandauReactionEnsemble::on_end_reaction(int &accepted_state) {
   update_wang_landau_potential_and_histogram(accepted_state);
 }
 
+/**
+ * Generic one way reaction
+ * A+B+...+G +... --> K+...X + Z +...
+ * You need to use 2A --> B instead of A+A --> B since in the last case you
+ * assume distinctness of the particles, however both ways to describe the
+ * reaction are equivalent in the thermodynamic limit (large particle numbers).
+ * Furthermore, it is crucial for the function in which order you provide the
+ * reactant and product types since particles will be replaced correspondingly!
+ * If there are less reactants than products, new product particles are created
+ * randomly in the box. Matching particles simply change the types. If there
+ * are more reactants than products, old reactant particles are deleted.
+ */
 bool ReactionAlgorithm::generic_oneway_reaction(int reaction_id) {
-  /**
-   *generic one way reaction
-   *A+B+...+G +... --> K+...X + Z +...
-   *you need to use 2A --> B instead of A+A --> B since in the last case you
-   *assume
-   *distinctness of the particles, however both ways to describe the reaction
-   *are equivalent in the thermodynamic limit (large particle numbers) further
-   *it is crucial for the function in which order you provide the reactant and
-   *product types since particles will be replaced correspondingly! If there are
-   *less reactants than products, new product particles are created randomly in
-   *the box. Matching particles simply change the types. If there are more
-   *reactants than products, old reactant particles are deleted.
-   */
 
   SingleReaction &current_reaction = reactions[reaction_id];
   current_reaction.tried_moves += 1;
@@ -527,9 +517,9 @@ bool ReactionAlgorithm::generic_oneway_reaction(int reaction_id) {
     for (int p_ids_created_particle : p_ids_created_particles) {
       delete_particle(p_ids_created_particle);
     }
-    // 2)restore previously hidden reactant particles
+    // 2) restore previously hidden reactant particles
     restore_properties(hidden_particles_properties, number_of_saved_properties);
-    // 2)restore previously changed reactant particles
+    // 3) restore previously changed reactant particles
     restore_properties(changed_particles_properties,
                        number_of_saved_properties);
     reaction_is_accepted = false;
@@ -570,19 +560,19 @@ void ReactionAlgorithm::replace_particle(int p_id, int desired_type) {
  * Hides a particle from short ranged interactions and from the electrostatic
  * interaction. Additional hiding from interactions would need to be implemented
  * here.
+ *
+ * Removing the particle charge and changing its type to a non existing one
+ * deactivates all interactions with other particles, as if the particle was
+ * inexistent (currently only type-based interactions are switched off, as well
+ * as the electrostatic interaction).
+ * This function does not break bonds for simple reactions, as long as there
+ * are no reactions like 2A -->B where one of the reacting A particles occurs
+ * in the polymer (think of bond breakages if the monomer in the polymer gets
+ * deleted in the reaction). This constraint is not of fundamental reason, but
+ * there would be a need for a rule for such "collision" reactions (a reaction
+ * like the one above).
  */
 void ReactionAlgorithm::hide_particle(int p_id, int previous_type) {
-  /**
-   *remove_charge and put type to a non existing one --> no interactions anymore
-   *it is as if the particle was non existing (currently only type-based
-   *interactions are switched off, as well as the electrostatic interaction)
-   *hide_particle() does not break bonds for simple reactions. as long as there
-   *are no reactions like 2A -->B where one of the reacting A particles occurs
-   *in the polymer (think of bond breakages if the monomer in the polymer gets
-   *deleted in the reaction). This constraint is not of fundamental reason, but
-   *there would be a need for a rule for such "collision" reactions (a reaction
-   *like the one above).
-   */
 
   auto part = get_particle_data(p_id);
   double d_min = distto(partCfg(), part.r.p, p_id);
@@ -597,17 +587,13 @@ void ReactionAlgorithm::hide_particle(int p_id, int previous_type) {
   set_particle_type(p_id, non_interacting_type);
 }
 
+/**
+ * Deletes the particle with the given p_id and stores the id if the deletion
+ * created a hole in the particle id range. This method is intended to only
+ * delete unbonded particles since bonds are coupled to ids. This is used to
+ * avoid the id range becoming excessively huge.
+ */
 int ReactionAlgorithm::delete_particle(int p_id) {
-  /**
-   * Deletes the particle with the given p_id and stores if it created a hole
-   * at that position in the particle id range. This method is intended to
-   * only
-   * delete unbonded particles since bonds are coupled to ids. This is used to
-   * avoid the id
-   * range becoming excessively huge.
-   */
-
-  /**deletes the particle with the provided id  */
   int old_max_seen_id = max_seen_particle;
   if (p_id == old_max_seen_id) {
     // last particle, just delete
@@ -742,19 +728,13 @@ int ReactionAlgorithm::create_particle(int desired_type) {
   // set velocities
   set_particle_v(p_id, vel);
   double d_min = distto(partCfg(), pos_vec, p_id);
-  if (d_min < exclusion_radius)
-    particle_inside_exclusion_radius_touched =
-        true; // setting of a minimal
-              // distance is allowed to
-              // avoid overlapping
-              // configurations if there is
-              // a repulsive potential.
-              // States with very high
-              // energies have a probability
-              // of almost zero and
-              // therefore do not contribute
-              // to ensemble averages.
-
+  if (d_min < exclusion_radius) {
+    // setting of a minimal distance is allowed to avoid overlapping
+    // configurations if there is a repulsive potential. States with
+    // very high energies have a probability of almost zero and
+    // therefore do not contribute to ensemble averages.
+    particle_inside_exclusion_radius_touched = true;
+  }
   return p_id;
 }
 
@@ -876,16 +856,17 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
                                                                    // symmetric
   }
 
-  //	//correct for enhanced proposal of small radii by using the Metropolis
-  // hastings algorithm for asymmetric proposal densities.
-  //	double
-  // old_radius=std::sqrt(std::pow(particle_positions[0]-cyl_x,2)+std::pow(particle_positions[1]-cyl_y,2));
-  //	double
-  // new_radius=std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2));
-  //	bf=std::min(1.0,
-  // bf*exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius);
-  ////Metropolis-Hastings Algorithm for asymmetric proposal density
+  // // correct for enhanced proposal of small radii by using the
+  // // Metropolis-Hastings algorithm for asymmetric proposal densities
+  // double old_radius =
+  //     std::sqrt(std::pow(particle_positions[0]-cyl_x,2) +
+  //               std::pow(particle_positions[1]-cyl_y,2));
+  // double new_radius =
+  //     std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2));
+  // bf = std::min(1.0,
+  //     bf*exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius);
 
+  // Metropolis-Hastings Algorithm for asymmetric proposal density
   if (m_uniform_real_distribution(m_generator) < bf) {
     // accept
     m_accepted_configurational_MC_moves += 1;
@@ -957,7 +938,7 @@ int WangLandauReactionEnsemble::get_flattened_index_wang_landau(
   int index = -10; // negative number is not allowed as index and therefore
                    // indicates error
   std::vector<int> individual_indices(nr_collective_variables); // pre result
-  //	individual_indices.resize(nr_collective_variables,-1); //initialize
+  // individual_indices.resize(nr_collective_variables,-1); //initialize
   // individual_indices to -1
 
   // check for the current state to be an allowed state in the [range
@@ -1197,16 +1178,13 @@ int WangLandauReactionEnsemble::initialize_wang_landau() {
 }
 
 /**
- * Calculates the expression which occurs in the Wang-Landau acceptance
- * probability.
+ * Calculates the expression in the acceptance probability of the Wang-Landau
+ * reaction ensemble
  */
 double WangLandauReactionEnsemble::calculate_acceptance_probability(
     SingleReaction &current_reaction, double E_pot_old, double E_pot_new,
     std::map<int, int> &old_particle_numbers, int old_state_index,
     int new_state_index, bool only_make_configuration_changing_move) {
-  /**determine the acceptance probabilities of the reaction move
-   * in Wang-Landau reaction ensemble
-   */
   double beta = 1.0 / temperature;
   double bf;
   if (do_not_sample_reaction_partition_function ||
@@ -1325,9 +1303,9 @@ int WangLandauReactionEnsemble::do_reaction(int reaction_steps) {
 
 // boring helper functions
 
+/** Increase the Wang-Landau potential and histogram at the current nbar */
 void WangLandauReactionEnsemble::update_wang_landau_potential_and_histogram(
     int index_of_state_after_acceptance_or_rejection) {
-  /**increase the Wang-Landau potential and histogram at the current nbar */
   if (index_of_state_after_acceptance_or_rejection >= 0) {
     if (histogram[index_of_state_after_acceptance_or_rejection] >= 0) {
       histogram[index_of_state_after_acceptance_or_rejection] += 1;
@@ -1695,20 +1673,6 @@ int ConstantpHEnsemble::get_random_valid_p_id() {
   return random_p_id;
 }
 
-/**
- * Constant-pH Ensemble, for derivation see Reed and Reed 1992
- * For the constant pH reactions you need to provide the deprotonation and
- * afterwards the corresponding protonation reaction (in this order). If you
- * want to deal with multiple reactions do it multiple times. Note that there is
- * a difference in the usecase of the constant pH reactions and the above
- * reaction ensemble. For the constant pH simulation directily the
- * **apparent equilibrium constant which carries a unit** needs to be provided
- * -- this is equivalent to the gamma of the reaction ensemble above, where the
- * dimensionless reaction constant needs to be provided. Again: For the
- * constant-pH algorithm not the dimensionless reaction constant needs to be
- * provided here, but the apparent reaction constant.
- */
-
 /**
  *Performs a reaction in the constant pH ensemble
  */
@@ -1736,10 +1700,10 @@ int ConstantpHEnsemble::do_reaction(int reaction_steps) {
       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
+             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);
@@ -1756,15 +1720,15 @@ int ConstantpHEnsemble::do_reaction(int reaction_steps) {
   return 0;
 }
 
+/**
+ * Calculates the expression in the acceptance probability of the constant pH
+ * method.
+ */
 double ConstantpHEnsemble::calculate_acceptance_probability(
     SingleReaction &current_reaction, double E_pot_old, double E_pot_new,
     std::map<int, int> &dummy_old_particle_numbers, int dummy_old_state_index,
     int dummy_new_state_index,
     bool dummy_only_make_configuration_changing_move) {
-  /**
-   *Calculates the expression in the acceptance probability of the constant pH
-   *method.
-   */
   double ln_bf;
   double pKa;
   const double beta = 1.0 / temperature;
@@ -1796,9 +1760,9 @@ WidomInsertion::measure_excess_chemical_potential(int reaction_id) {
   for (int p_ids_created_particle : p_ids_created_particles) {
     delete_particle(p_ids_created_particle);
   }
-  // 2)restore previously hidden reactant particles
+  // 2) restore previously hidden reactant particles
   restore_properties(hidden_particles_properties, number_of_saved_properties);
-  // 2)restore previously changed reactant particles
+  // 3) restore previously changed reactant particles
   restore_properties(changed_particles_properties, number_of_saved_properties);
   std::vector<double> exponential = {
       exp(-1.0 / temperature * (E_pot_new - E_pot_old))};
diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp
index 99bd5f46957..c9f50dae8bd 100644
--- a/src/core/reaction_ensemble.hpp
+++ b/src/core/reaction_ensemble.hpp
@@ -106,6 +106,7 @@ struct DegreeOfAssociationCollectiveVariable : public CollectiveVariable {
   }
 };
 
+/** Base class for reaction ensemble methods */
 class ReactionAlgorithm {
 
 public:
@@ -223,9 +224,17 @@ class ReactionAlgorithm {
   get_random_position_in_box_enhanced_proposal_of_small_radii();
 };
 
-////////////////////////////////////////////////////////////////actual
-/// declaration of specific reaction algorithms
-
+///////////////////////////// actual declaration of specific reaction algorithms
+
+/** Reaction ensemble method according to smith94x.
+ *  Works for the reaction ensemble at constant volume and temperature. For the
+ *  reaction ensemble at constant pressure additionally employ a barostat!
+ *  NOTE: a chemical reaction consists of a forward and backward reaction.
+ *  Here both reactions have to be defined separately. The extent of the
+ *  reaction is here chosen to be +1. If the reaction trial move for a
+ *  dissociation of HA is accepted then there is one more dissociated ion
+ *  pair H+ and A-
+ */
 class ReactionEnsemble : public ReactionAlgorithm {
 public:
   ReactionEnsemble(int seed) : ReactionAlgorithm(seed) {}
@@ -238,6 +247,7 @@ class ReactionEnsemble : public ReactionAlgorithm {
       bool dummy_only_make_configuration_changing_move) override;
 };
 
+/** Wang-Landau reaction ensemble method */
 class WangLandauReactionEnsemble : public ReactionAlgorithm {
 public:
   WangLandauReactionEnsemble(int seed) : ReactionAlgorithm(seed) {}
@@ -341,6 +351,19 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm {
                                                       double delta_CV);
 };
 
+/**
+ * Constant-pH Ensemble, for derivation see Reed and Reed 1992.
+ * For the constant pH reactions you need to provide the deprotonation and
+ * afterwards the corresponding protonation reaction (in this order). If you
+ * want to deal with multiple reactions do it multiple times. Note that there is
+ * a difference in the usecase of the constant pH reactions and the above
+ * reaction ensemble. For the constant pH simulation directily the
+ * **apparent equilibrium constant which carries a unit** needs to be provided
+ * -- this is equivalent to the gamma of the reaction ensemble above, where the
+ * dimensionless reaction constant needs to be provided. Again: For the
+ * constant-pH algorithm not the dimensionless reaction constant needs to be
+ * provided here, but the apparent reaction constant.
+ */
 class ConstantpHEnsemble : public ReactionAlgorithm {
 public:
   ConstantpHEnsemble(int seed) : ReactionAlgorithm(seed) {}
@@ -356,6 +379,7 @@ class ConstantpHEnsemble : public ReactionAlgorithm {
   int get_random_valid_p_id();
 };
 
+/** Widom insertion method */
 class WidomInsertion : public ReactionAlgorithm {
 public:
   WidomInsertion(int seed) : ReactionAlgorithm(seed) {}
diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx
index be2034a8edb..6eb1bafbc95 100644
--- a/src/python/espressomd/reaction_ensemble.pyx
+++ b/src/python/espressomd/reaction_ensemble.pyx
@@ -13,10 +13,10 @@ class WangLandauHasConverged(Exception):
 cdef class ReactionAlgorithm:
     """
 
-    This class provides the base class for Reaction Algorithms like the Reaction Ensemble algorithm, the Wang-Landau
-    Reaction Ensemble algorithm and the constant pH method. Initialize the
-    reaction algorithm by setting the standard pressure, temperature, and the
-    exclusion radius.
+    This class provides the base class for Reaction Algorithms like the Reaction
+    Ensemble algorithm, the Wang-Landau Reaction Ensemble algorithm and the
+    constant pH method. Initialize the reaction algorithm by setting the
+    standard pressure, temperature, and the exclusion radius.
 
     Note: When creating particles the velocities of the new particles are set
     according the Maxwell-Boltzmann distribution. In this step the mass of the
@@ -26,16 +26,15 @@ cdef class ReactionAlgorithm:
     Parameters
     ----------
     temperature : :obj:`float`
-                  The temperature at which the reaction is performed.
+        The temperature at which the reaction is performed.
     exclusion_radius : :obj:`float`
-                       Minimal distance from any particle, within which new
-                       particle will not be inserted. This is useful to avoid
-                       integrator failures if particles are too close and there
-                       is a diverging repulsive interaction, or to prevent two
-                       oppositely charged particles from being placed on top of
-                       each other.  The Boltzmann factor :math:`\exp(-\\beta
-                       E)` gives these configurations a small contribution to
-                       the partition function, therefore they can be neglected.
+        Minimal distance from any particle, within which new particle will not
+        be inserted. This is useful to avoid integrator failures if particles
+        are too close and there is a diverging repulsive interaction, or to
+        prevent two oppositely charged particles from being placed on top of
+        each other. The Boltzmann factor :math:`\\exp(-\\beta E)` gives these
+        configurations a small contribution to the partition function,
+        therefore they can be neglected.
     """
     cdef object _params
     cdef CReactionAlgorithm * RE
@@ -47,8 +46,7 @@ cdef class ReactionAlgorithm:
         return "temperature", "exclusion_radius", "seed"
 
     def _set_params_in_es_core(self):
-        deref(self.RE).temperature = self._params[
-            "temperature"]
+        deref(self.RE).temperature = self._params["temperature"]
         # setting a volume is a side effect, sets the default volume of the
         # reaction ensemble as the volume of the cuboid simulation box. this
         # volume can be altered by the command "reaction ensemble volume
@@ -57,8 +55,7 @@ cdef class ReactionAlgorithm:
         # only contained in the volume of the cylinder)
         if(deref(self.RE).volume < 0):
             deref(self.RE).set_cuboid_reaction_ensemble_volume()
-        deref(self.RE).exclusion_radius = self._params[
-            "exclusion_radius"]
+        deref(self.RE).exclusion_radius = self._params["exclusion_radius"]
 
     def set_cylindrical_constraint_in_z_direction(self, center_x, center_y,
                                                   radius_of_cylinder):
@@ -70,11 +67,11 @@ cdef class ReactionAlgorithm:
         Parameters
         ----------
         center_x : :obj:`float`
-                   x coordinate of center of the cylinder.
+            x coordinate of center of the cylinder.
         center_y : :obj:`float`
-                   y coordinate of center of the cylinder.
+            y coordinate of center of the cylinder.
         radius_of_cylinder : :obj:`float`
-                   radius of the cylinder
+            radius of the cylinder
 
         """
         deref(self.RE).cyl_x = center_x
@@ -135,14 +132,14 @@ cdef class ReactionAlgorithm:
         """
         Sets the particle type for non-interacting particles.
         Default value: 100.
-        This is used to temporarily hide
-        particles during a reaction trial move, if they are to be deleted after
-        the move is accepted.
-        Please change this value if you intend to use the type 100 for some other
-        particle types with interactions. Please also note that particles in the current
-        implementation of the Reaction Ensemble are only hidden with respect to
-        Lennard-Jones and Coulomb interactions. Hiding of other interactions,
-        for example a magnetic, needs to be implemented in the code.
+        This is used to temporarily hide particles during a reaction trial
+        move, if they are to be deleted after the move is accepted. Please
+        change this value if you intend to use the type 100 for some other
+        particle types with interactions. Please also note that particles
+        in the current implementation of the Reaction Ensemble are only
+        hidden with respect to Lennard-Jones and Coulomb interactions. Hiding
+        of other interactions, for example a magnetic, needs to be implemented
+        in the code.
         """
         deref(self.RE).non_interacting_type = non_interacting_type
 
@@ -160,24 +157,22 @@ cdef class ReactionAlgorithm:
         Parameters
         ----------
         gamma : :obj:`float`
-                               Equilibrium constant of the reaction, :math:`\gamma`
-                               (see the User guide, section 6.6 for the definition and further details).
+            Equilibrium constant of the reaction, :math:`\\gamma` (see the User
+            guide, section 6.6 for the definition and further details).
         reactant_types : list of :obj:`int`
-                                List of particle types of reactants in the reaction.
-        reactant_coefficients : list
-                                List of stoichiometric coefficients of the
-                                reactants in the same order as the list of
-                                their types.
-        product_types : list
-                               List of particle types of products in the reaction.
-        product_coefficients : list
-                               List of stoichiometric coefficients of
-                               products of the reaction in the same order as
-                               the list of their types
+            List of particle types of reactants in the reaction.
+        reactant_coefficients : list of :obj:`int`
+            List of stoichiometric coefficients of the reactants in the same
+            order as the list of their types.
+        product_types : list of :obj:`int`
+            List of particle types of products in the reaction.
+        product_coefficients : list of :obj:`int`
+            List of stoichiometric coefficients of products of the reaction in
+            the same order as the list of their types
         default_charges : dictionary
-                        A dictionary of default charges for types that occur in the provided reaction.
-        check_for_electroneutrality : Boolean
-                                      Check for electroneutrality of the given reaction if True.
+            A dictionary of default charges for types that occur in the provided reaction.
+        check_for_electroneutrality : :obj:`bool`
+            Check for electroneutrality of the given reaction if True.
 
         """
         self._params["check_for_electroneutrality"] = True
@@ -267,7 +262,7 @@ cdef class ReactionAlgorithm:
         Parameters
         ----------
         reaction_steps : :obj:`int`, optional
-                          The number of reactions to be performed at once, defaults to 1.
+            The number of reactions to be performed at once, defaults to 1.
 
         """
         deref(self.RE).do_reaction(int(reaction_steps))
@@ -275,10 +270,11 @@ cdef class ReactionAlgorithm:
     def displacement_mc_move_for_particles_of_type(self, type_mc,
                                                    particle_number_to_be_changed=1):
         """
-        Performs a displacement Monte Carlo move for particles of given type. New positions
-        of the displaced particles are chosen from the whole box with a uniform probability distribution.
-        If there are multiple types, that are being moved in a simulation, they should be moved in a
-        random order to avoid artefacts.
+        Performs a displacement Monte Carlo move for particles of given type.
+        New positions of the displaced particles are chosen from the whole box
+        with a uniform probability distribution. If there are multiple types,
+        that are being moved in a simulation, they should be moved in a random
+        order to avoid artefacts.
 
         Parameters
         ----------
@@ -317,8 +313,12 @@ cdef class ReactionAlgorithm:
             for i in range(deref(self.RE).reactions[single_reaction_i].product_types.size()):
                 product_coefficients.append(
                     deref(self.RE).reactions[single_reaction_i].product_coefficients[i])
-            reaction = {"reactant_coefficients": reactant_coefficients, "reactant_types": reactant_types, "product_types": product_types, "product_coefficients":
-                        product_coefficients, "reactant_types": reactant_types, "gamma": deref(self.RE).reactions[single_reaction_i].gamma}
+            reaction = {"reactant_coefficients": reactant_coefficients,
+                        "reactant_types": reactant_types,
+                        "product_types": product_types,
+                        "product_coefficients": product_coefficients,
+                        "reactant_types": reactant_types,
+                        "gamma": deref(self.RE).reactions[single_reaction_i].gamma}
             reactions.append(reaction)
 
         return {"reactions": reactions, "temperature": deref(self.RE).temperature, "exclusion_radius": deref(self.RE).exclusion_radius}
@@ -356,7 +356,7 @@ cdef class ReactionEnsemble(ReactionAlgorithm):
             if k in self._valid_keys():
                 self._params[k] = kwargs[k]
             else:
-                raise KeyError("%s is not a vaild key" % k)
+                raise KeyError("%s is not a valid key" % k)
 
         self._set_params_in_es_core()
 
@@ -379,7 +379,7 @@ cdef class ConstantpHEnsemble(ReactionAlgorithm):
             if k in self._valid_keys():
                 self._params[k] = kwargs[k]
             else:
-                raise KeyError("%s is not a vaild key" % k)
+                raise KeyError("%s is not a valid key" % k)
 
         self._set_params_in_es_core()
 
@@ -425,7 +425,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
             if k in self._valid_keys():
                 self._params[k] = kwargs[k]
             else:
-                raise KeyError("%s is not a vaild key" % k)
+                raise KeyError("%s is not a valid key" % k)
 
         self.WLRptr.reset(new CWangLandauReactionEnsemble(int(self._params["seed"])))
         self.RE = <CReactionAlgorithm * > self.WLRptr.get()
@@ -434,11 +434,12 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
 
     def reaction(self, reaction_steps=1):
         """
-        Performs reaction_steps reactions. Sets the number of reaction steps which are
-        performed at once. Do not use too many reaction steps
-        steps consecutively without having conformation
-        changing steps in between (especially important for the Wang-Landau reaction ensemble). Providing a number for the parameter reaction steps reduces the need for the interpreter to be
-        called between consecutive reactions.
+        Performs reaction_steps reactions. Sets the number of reaction steps
+        which are performed at once. Do not use too many reaction steps
+        consecutively without having conformation-changing steps in-between
+        (especially important for the Wang-Landau reaction ensemble). Providing
+        a number for the parameter reaction steps reduces the need for the
+        interpreter to be called between consecutive reactions.
 
         """
         status_wang_landau = deref(
@@ -455,14 +456,13 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
         Parameters
         ----------
         associated_type : :obj:`int`
-                          Particle type of the associated state of the reacting species.
+            Particle type of the associated state of the reacting species.
         min : :obj:`float`
-              Minimum value of the collective variable.
+            Minimum value of the collective variable.
         max : :obj:`float`
-              Maximum value of the collective variable.
-        corresponding_acid_types : list
-                                   List of the types of the version of the
-                                   species.
+            Maximum value of the collective variable.
+        corresponding_acid_types : list of :obj:`int`
+            List of the types of the version of the species.
 
         """
         for k in kwargs:
@@ -498,23 +498,23 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
         Parameters
         ----------
         filename : :obj:`str`
-                   Filename of the energy boundary file which provides the
-                   potential energy boundaries (min E_pot, max E_pot) tabulated
-                   for all degrees of association. Make sure to only list the
-                   degrees of association which are used by the degree of
-                   association collective variable within this file. The energy
-                   boundary file can be created in a preliminary energy run. By
-                   the help of the functions
-                   :meth:`update_maximum_and_minimum_energies_at_current_state`
-                   and :meth:`write_out_preliminary_energy_run_results`. This
-                   file has to be obtained before being able to run a
-                   simulation with the energy as collective variable.
+            Filename of the energy boundary file which provides the
+            potential energy boundaries (min E_pot, max E_pot) tabulated
+            for all degrees of association. Make sure to only list the
+            degrees of association which are used by the degree of
+            association collective variable within this file. The energy
+            boundary file can be created in a preliminary energy run. By
+            the help of the functions
+            :meth:`update_maximum_and_minimum_energies_at_current_state`
+            and :meth:`write_out_preliminary_energy_run_results`. This
+            file has to be obtained before being able to run a
+            simulation with the energy as collective variable.
         delta : :obj:`float`
-                Provides the discretization of the potential energy range. Only
-                for small enough delta the results of the energy reweighted
-                averages are correct. If delta is chosen too big there are
-                discretization errors in the numerical integration which occurs
-                during the energy reweighting process.
+            Provides the discretization of the potential energy range. Only
+            for small enough delta the results of the energy reweighted
+            averages are correct. If delta is chosen too big there are
+            discretization errors in the numerical integration which occurs
+            during the energy reweighting process.
 
         """
         for k in kwargs:
@@ -546,25 +546,17 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
         Parameters
         ----------
         final_wang_landau_parameter : :obj:`float`
-                                      Sets the final Wang-Landau parameter, which is the Wang-Landau parameter after which the simulation should stop.).
+            Sets the final Wang-Landau parameter, which is the Wang-Landau
+            parameter after which the simulation should stop.
         full_path_to_output_filename : :obj:`str`
-                                       Sets the path to the output file of the
-                                       Wang-Landau algorithm which contains the
-                                       Wang-Landau potential
+            Sets the path to the output file of the Wang-Landau algorithm which
+            contains the Wang-Landau potential
         do_not_sample_reaction_partition_function : :obj:`bool`
-                                                    Avoids sampling the
-                                                    Reaction ensemble partition
-                                                    function in the Wang-Landau
-                                                    algorithm. Therefore this
-                                                    option makes all degrees of
-                                                    association equally
-                                                    probable. This option may
-                                                    be used in the sweeping
-                                                    mode of the reaction
-                                                    ensemble, since the
-                                                    reaction ensemble partition
-                                                    function can be later added
-                                                    analytically.
+            Avoids sampling the Reaction ensemble partition function in the
+            Wang-Landau algorithm. Therefore this option makes all degrees of
+            association equally probable. This option may be used in the
+            sweeping mode of the reaction ensemble, since the reaction ensemble
+            partition function can be later added analytically.
 
         """
         for k in kwargs:
@@ -606,9 +598,9 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
         Records the minimum and maximum potential energy as a function of the
         degree of association in a preliminary Wang-Landau reaction ensemble
         simulation where the acceptance probability includes the factor
-        :math:`\exp(-\\beta \\Delta E_{pot})`. The minimal and maximal
+        :math:`\\exp(-\\beta \\Delta E_{pot})`. The minimal and maximal
         potential energies which occur in the system are needed for the energy
-        reweighting simulations where the factor :math:`\exp(-\\beta \\Delta E_{pot})`
+        reweighting simulations where the factor :math:`\\exp(-\\beta \\Delta E_{pot})`
         is not included in the acceptance probability in
         order to avoid choosing the wrong potential energy boundaries.
 
@@ -638,10 +630,18 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
     def displacement_mc_move_for_particles_of_type(self, type_mc,
                                                    particle_number_to_be_changed=1):
         """
-        Performs an MC (Monte Carlo) move for particle_number_to_be_changed particle of type type_mc. Positions for the particles are drawn uniformly random within the box. The command takes into account the Wang-Landau terms in the acceptance probability.
-        If there are multiple types, that
-        need to be moved, make sure to move them in a random order to avoid
-        artefacts. For the Wang-Landau algorithm in the case of energy reweighting you would also need to move the monomers of the polymer with special moves for the MC part. Those polymer configuration changing moves need to be implemented in the case of using Wang-Landau with energy reweighting and a polymer in the system. Polymer configuration changing moves had been implemented before but were removed from espresso.
+        Performs an MC (Monte Carlo) move for particle_number_to_be_changed
+        particle of type type_mc. Positions for the particles are drawn
+        uniformly random within the box. The command takes into account the
+        Wang-Landau terms in the acceptance probability.
+        If there are multiple types, that need to be moved, make sure to move
+        them in a random order to avoid artefacts. For the Wang-Landau algorithm
+        in the case of energy reweighting you would also need to move the
+        monomers of the polymer with special moves for the MC part. Those
+        polymer configuration changing moves need to be implemented in the
+        case of using Wang-Landau with energy reweighting and a polymer in the
+        system. Polymer configuration changing moves had been implemented
+        before but were removed from espresso.
 
         """
         use_wang_landau = True
@@ -651,7 +651,9 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
 
 cdef class WidomInsertion(ReactionAlgorithm):
     """
-    This class implements the Widom insertion method in the canonical ensemble for homogeneous systems, where the excess chemical potential is not depending on the location.
+    This class implements the Widom insertion method in the canonical ensemble
+    for homogeneous systems, where the excess chemical potential is not
+    depending on the location.
 
     """
 
@@ -687,13 +689,16 @@ cdef class WidomInsertion(ReactionAlgorithm):
             if k in self._valid_keys():
                 self._params[k] = kwargs[k]
             else:
-                raise KeyError("%s is not a vaild key" % k)
+                raise KeyError("%s is not a valid key" % k)
 
         self._set_params_in_es_core()
 
     def measure_excess_chemical_potential(self, reaction_id=0):
         """
-        Measures the excess chemical potential in a homogeneous system. Returns the excess chemical potential and the standard error for the excess chemical potential. It assumes that your samples are uncorrelated in estimating the standard error.
+        Measures the excess chemical potential in a homogeneous system.
+        Returns the excess chemical potential and the standard error for the
+        excess chemical potential. It assumes that your samples are
+        uncorrelated in estimating the standard error.
 
         """
         return self.WidomInsertionPtr.get().measure_excess_chemical_potential(int(reaction_id))
diff --git a/testsuite/scripts/samples/test_grand_canonical.py b/testsuite/scripts/samples/test_grand_canonical.py
index af46713722e..e42a03a0ff4 100644
--- a/testsuite/scripts/samples/test_grand_canonical.py
+++ b/testsuite/scripts/samples/test_grand_canonical.py
@@ -20,7 +20,11 @@
 
 sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
     "@SAMPLES_DIR@/grand_canonical.py", n_int_cycles=51, n_int_steps=5,
-    warm_n_times=10, warm_steps=50, cmd_arguments=[0.01, 0.01])
+    warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659e-04, -0.336])
+# Note: the command line arguments which are submitted are 1) the salt
+# concentration in 1/sigma^3 and 2) the excess chemical potential in kT.
+# The excess chemical potential needs to be determined using the Widom
+# insertion method, e.g. with the example script widom_insertion.py
 
 
 @skipIfMissingFeatures
@@ -28,8 +32,8 @@ class Sample(ut.TestCase):
     system = sample.system
 
     def test_deviation_to_target_concentration(self):
-        # deviation < 10%
-        self.assertLess(abs(sample.deviation), 10.0)
+        # deviation < 15%
+        self.assertLess(abs(sample.deviation), 15.0)
 
 
 if __name__ == "__main__":