Change Reaction Methods interface

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