Change Reaction Methods interface (#4305)

Fixes minor issues in #4251

Description of changes:
- `ReactionAlgorithm::generic_oneway_reaction()`: changed the name of the accumulators to `accumulator_potential_energy_difference_exponential` for clarity
- The name of the argument `temperature` has been replaced by `kT`. The related tests have been modified accordingly.
- The arguments `reactant_coefficients` and `product_coefficients` are now optional parameters for the constant pH ensemble, with default values of `[1]` resp. `[1, 1]`. If the user provides these parameters he gets either (i) an error message if values other than the default ones are supplied or (ii) a one-time warning stating that these arguments are only kept as input for backward compatibility and are deprecated.

doc/tutorials/constant_pH/constant_pH.ipynb

@@ -488,7 +488,7 @@
     "```python\n",
     "exclusion_radius = 1.0 if USE_WCA else 0.0\n",
     "RE = espressomd.reaction_ensemble.ConstantpHEnsemble(\n",
-    "    temperature=KT.to(\"sim_energy\").magnitude,\n",
+    "    kT=KT.to(\"sim_energy\").magnitude,\n",
     "    exclusion_radius=exclusion_radius,\n",
     "    seed=77\n",
     ")\n",

samples/grand_canonical.py

@@ -89,7 +89,7 @@
             epsilon=wca_eps, sigma=wca_sig)
 
 RE = espressomd.reaction_ensemble.ReactionEnsemble(
-    temperature=temperature, exclusion_radius=wca_sig, seed=3)
+    kT=temperature, exclusion_radius=wca_sig, seed=3)
 RE.add_reaction(
     gamma=cs_bulk**2 * np.exp(excess_chemical_potential_pair / temperature),
     reactant_types=[], reactant_coefficients=[], product_types=[1, 2],

samples/reaction_ensemble.py

@@ -72,19 +72,23 @@
 RE = None
 if args.mode == "reaction_ensemble":
     RE = espressomd.reaction_ensemble.ReactionEnsemble(
-        temperature=1,
+        kT=1,
         exclusion_radius=1,
         seed=77)
+    RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1],
+                    product_types=[1, 2], product_coefficients=[1, 1],
+                    default_charges={0: 0, 1: -1, 2: +1})
 elif args.mode == "constant_pH_ensemble":
     RE = espressomd.reaction_ensemble.ConstantpHEnsemble(
-        temperature=1, exclusion_radius=1, seed=77)
+        kT=1, exclusion_radius=1, seed=77)
     RE.constant_pH = 2
+    RE.add_reaction(gamma=K_diss, reactant_types=[0],
+                    product_types=[1, 2],
+                    default_charges={0: 0, 1: -1, 2: +1})
 else:
     raise RuntimeError(
         "Please provide either --reaction_ensemble or --constant_pH_ensemble as argument ")
-RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1],
-                product_types=[1, 2], product_coefficients=[1, 1],
-                default_charges={0: 0, 1: -1, 2: +1})
+
 print(RE.get_status())
 system.setup_type_map([0, 1, 2])
 

samples/reaction_ensemble_complex_reaction.py

@@ -82,7 +82,7 @@
 
 # use an exclusion radius of 0 to simulate an ideal gas
 RE = espressomd.reaction_ensemble.ReactionEnsemble(
-    temperature=1, exclusion_radius=0, seed=4)
+    kT=1, exclusion_radius=0, seed=4)
 
 
 RE.add_reaction(

samples/widom_insertion.py

@@ -103,7 +103,7 @@
 system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42)
 
 widom = espressomd.reaction_ensemble.WidomInsertion(
-    temperature=temperature, seed=77)
+    kT=temperature, seed=77)
 
 # add insertion reaction
 insertion_reaction_id = 0

src/core/reaction_methods/ConstantpHEnsemble.cpp

@@ -36,7 +36,7 @@ namespace ReactionMethods {
 double ConstantpHEnsemble::calculate_acceptance_probability(
     SingleReaction const &current_reaction, double E_pot_old, double E_pot_new,
     std::map<int, int> const &old_particle_numbers) const {
-  auto const beta = 1.0 / temperature;
+  auto const beta = 1.0 / kT;
   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) *

src/core/reaction_methods/ReactionAlgorithm.cpp

@@ -85,13 +85,13 @@ void ReactionAlgorithm::check_reaction_method() const {
     throw std::runtime_error("Reaction system not initialized");
   }
 
-  if (temperature < 0) {
-    throw std::runtime_error("Temperatures cannot be negative. Please provide "
-                             "a temperature (in k_B T) to the simulation. "
-                             "Normally it should be 1.0. This will be used "
-                             "directly to calculate beta:=1/(k_B T) which "
+  if (kT < 0) {
+    throw std::runtime_error("kT cannot be negative,"
+                             "normally it should be 1.0. This will be used"
+                             "directly to calculate beta:=1/(k_B T) which"
                              "occurs in the exp(-beta*E)\n");
   }
+
 #ifdef ELECTROSTATICS
   // check for the existence of default charges for all types that take part in
   // the reactions
@@ -336,9 +336,9 @@ void ReactionAlgorithm::generic_oneway_reaction(
   auto const bf = calculate_acceptance_probability(
       current_reaction, E_pot_old, E_pot_new, old_particle_numbers);
 
-  std::vector<double> exponential = {
-      exp(-1.0 / temperature * (E_pot_new - E_pot_old))};
-  current_reaction.accumulator_exponentials(exponential);
+  std::vector<double> exponential = {exp(-1.0 / kT * (E_pot_new - E_pot_old))};
+  current_reaction.accumulator_potential_energy_difference_exponential(
+      exponential);
 
   if (m_uniform_real_distribution(m_generator) < bf) {
     // accept
@@ -500,7 +500,7 @@ int ReactionAlgorithm::create_particle(int desired_type) {
   }
 
   // we use mass=1 for all particles, think about adapting this
-  move_particle(p_id, get_random_position_in_box(), std::sqrt(temperature));
+  move_particle(p_id, get_random_position_in_box(), std::sqrt(kT));
   set_particle_type(p_id, desired_type);
 #ifdef ELECTROSTATICS
   set_particle_q(p_id, charges_of_types[desired_type]);
@@ -540,7 +540,7 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {
     auto const &p = get_particle_data(p_id);
     old_positions.emplace_back(std::pair<int, Utils::Vector3d>{p_id, p.r.p});
     // write new position
-    auto const prefactor = std::sqrt(temperature / p.p.mass);
+    auto const prefactor = std::sqrt(kT / p.p.mass);
     auto const new_pos = get_random_position_in_box();
     move_particle(p_id, new_pos, prefactor);
     check_exclusion_radius(p_id);
@@ -573,7 +573,7 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
                              ? std::numeric_limits<double>::max()
                              : calculate_current_potential_energy_of_system();
 
-  double beta = 1.0 / temperature;
+  double beta = 1.0 / kT;
 
   // Metropolis algorithm since proposal density is symmetric
   auto const bf = std::min(1.0, exp(-beta * (E_pot_new - E_pot_old)));

src/core/reaction_methods/ReactionAlgorithm.hpp

@@ -53,7 +53,7 @@ class ReactionAlgorithm {
 
   std::vector<SingleReaction> reactions;
   std::map<int, double> charges_of_types;
-  double temperature = -10.0;
+  double kT = -10.0;
   /**
    * Hard sphere radius. If particles are closer than this value,
    * it is assumed that their interaction energy gets approximately

src/core/reaction_methods/ReactionEnsemble.cpp

@@ -36,7 +36,7 @@ double ReactionEnsemble::calculate_acceptance_probability(
   const double factorial_expr =
       calculate_factorial_expression(current_reaction, old_particle_numbers);
 
-  const double beta = 1.0 / temperature;
+  const double beta = 1.0 / kT;
   // calculate Boltzmann factor
   return std::pow(volume, current_reaction.nu_bar) * current_reaction.gamma *
          factorial_expr * exp(-beta * (E_pot_new - E_pot_old));

src/core/reaction_methods/SingleReaction.hpp

@@ -51,7 +51,8 @@ struct SingleReaction {
   double gamma = {};
   // calculated values that are stored for performance reasons
   int nu_bar = {}; ///< change in particle numbers for the reaction
-  Utils::Accumulator accumulator_exponentials = Utils::Accumulator(1);
+  Utils::Accumulator accumulator_potential_energy_difference_exponential =
+      Utils::Accumulator(1);
   int tried_moves = 0;
   int accepted_moves = 0;
   double get_acceptance_rate() const {

src/core/reaction_methods/WidomInsertion.cpp

@@ -62,17 +62,15 @@ std::pair<double, double> WidomInsertion::measure_excess_chemical_potential(
   restore_properties(hidden_particles_properties, number_of_saved_properties);
   // 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))};
-  current_reaction.accumulator_exponentials(exponential);
+  std::vector<double> exponential = {exp(-1.0 / kT * (E_pot_new - E_pot_old))};
+  current_reaction.accumulator_potential_energy_difference_exponential(
+      exponential);
 
   // calculate mean excess chemical potential and standard error of the mean
-  std::pair<double, double> result = std::make_pair(
-      -temperature * log(current_reaction.accumulator_exponentials.mean()[0]),
-      std::abs(-temperature /
-               current_reaction.accumulator_exponentials.mean()[0] *
-               current_reaction.accumulator_exponentials.std_error()[0]));
-  return result;
+  auto const &accumulator =
+      current_reaction.accumulator_potential_energy_difference_exponential;
+  return {-kT * log(accumulator.mean()[0]),
+          std::abs(-kT / accumulator.mean()[0] * accumulator.std_error()[0])};
 }
 
 } // namespace ReactionMethods

src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp

@@ -51,7 +51,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) {
   constexpr double tol = 100 * std::numeric_limits<double>::epsilon();
 
   ConstantpHEnsembleTest r_algo(42);
-  r_algo.temperature = 20.;
+  r_algo.kT = 20.;
   r_algo.m_constant_pH = 1.;
 
   // exception if no reaction was added
@@ -76,7 +76,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) {
         // acceptance = f_expr * exp(- E / T + nu_bar * log(10) * (pH - nu_bar *
         // pKa))
         auto const acceptance_ref =
-            f_expr * std::exp(energy / r_algo.temperature +
+            f_expr * std::exp(energy / r_algo.kT +
                               std::log(10.) * (r_algo.m_constant_pH +
                                                std::log10(reaction.gamma)));
         auto const acceptance = r_algo.calculate_acceptance_probability(

src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp

@@ -107,11 +107,11 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
     BOOST_CHECK_EQUAL(probability, -10.);
   }
 
-  // exception if temperature is negative
+  // exception if kT is negative
   BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error);
 
-  // set temperature
-  r_algo.temperature = 1.;
+  // set kT
+  r_algo.kT = 1.;
 
 #ifdef ELECTROSTATICS
   // exception if reactant types have no charge information

src/core/reaction_methods/tests/ReactionEnsemble_test.cpp

@@ -52,7 +52,7 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) {
 
   ReactionEnsembleTest r_algo(42);
   r_algo.volume = 10.;
-  r_algo.temperature = 20.;
+  r_algo.kT = 20.;
 
   // exception if no reaction was added
   BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error);
@@ -76,7 +76,7 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) {
         // acceptance = V^{nu_bar} * gamma * f_expr * exp(- E / T)
         auto const acceptance_ref = std::pow(r_algo.volume, reaction.nu_bar) *
                                     reaction.gamma * f_expr *
-                                    std::exp(energy / r_algo.temperature);
+                                    std::exp(energy / r_algo.kT);
         auto const acceptance = r_algo.calculate_acceptance_probability(
             reaction, energy, 0., p_numbers);
         BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol);

src/python/espressomd/reaction_ensemble.pxd

@@ -48,7 +48,7 @@ cdef extern from "reaction_methods/ReactionAlgorithm.hpp" namespace "ReactionMet
 
         vector[SingleReaction] reactions
         map[int, double] charges_of_types
-        double temperature
+        double kT
         double exclusion_radius
         double volume
         bool box_is_cylindric_around_z_axis