Widom excess chemical potential calculation and RxMC performance improvement

samples/widom_insertion.py

@@ -118,18 +118,21 @@
 # up the simulation
 widom.set_non_interacting_type(max(types) + 1)
 
-n_iterations = 100
+particle_insertion_potential_energy_samples = []
+
+n_iterations = 500
 for i in range(n_iterations):
-    for j in range(50):
-        widom.measure_excess_chemical_potential(insertion_reaction_id)
+    particle_insertion_potential_energy_samples.append(
+        widom.calculate_particle_insertion_potential_energy(0))
     system.integrator.run(steps=500)
+
     if i % 20 == 0:
-        print("mu_ex_pair ({:.4f}, +/- {:.4f})".format(
-            *widom.measure_excess_chemical_potential(insertion_reaction_id)))
         print(f"HA {system.number_of_particles(type=0)}",
               f"A- {system.number_of_particles(type=1)}",
               f"H+ {system.number_of_particles(type=2)}")
 
+mu_ex_mean, mu_ex_Delta = widom.calculate_excess_chemical_potential(
+    particle_insertion_potential_energy_samples=particle_insertion_potential_energy_samples)
 
-print("excess chemical potential for an ion pair ",
-      widom.measure_excess_chemical_potential(insertion_reaction_id))
+print(
+    f"excess chemical potential for an ion pair {mu_ex_mean:.4g} +/- {mu_ex_Delta:.4g}")

src/core/reaction_methods/ReactionAlgorithm.cpp

@@ -49,9 +49,13 @@ namespace ReactionMethods {
  * Performs a randomly selected reaction in the reaction ensemble
  */
 int ReactionAlgorithm::do_reaction(int reaction_steps) {
+  // calculate potential energy; only consider potential energy since we
+  // assume that the kinetic part drops out in the process of calculating
+  // ensemble averages (kinetic part may be separated and crossed out)
+  auto current_E_pot = calculate_current_potential_energy_of_system();
   for (int i = 0; i < reaction_steps; i++) {
     int reaction_id = i_random(static_cast<int>(reactions.size()));
-    generic_oneway_reaction(reactions[reaction_id]);
+    generic_oneway_reaction(reactions[reaction_id], current_E_pot);
   }
   return 0;
 }
@@ -278,20 +282,8 @@ std::map<int, int> ReactionAlgorithm::save_old_particle_numbers(
   return old_particle_numbers;
 }
 
-/**
- * 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.
- */
 void ReactionAlgorithm::generic_oneway_reaction(
-    SingleReaction &current_reaction) {
+    SingleReaction &current_reaction, double &E_pot_old) {
 
   current_reaction.tried_moves += 1;
   particle_inside_exclusion_radius_touched = false;
@@ -301,18 +293,6 @@ void ReactionAlgorithm::generic_oneway_reaction(
     return;
   }
 
-  // calculate potential energy
-  const double E_pot_old =
-      calculate_current_potential_energy_of_system(); // only consider potential
-                                                      // energy since we assume
-                                                      // that the kinetic part
-                                                      // drops out in the
-                                                      // process of calculating
-                                                      // ensemble averages
-                                                      // (kinetic part may be
-                                                      // separated and crossed
-                                                      // out)
-
   // find reacting molecules in reactants and save their properties for later
   // recreation if step is not accepted
   // do reaction
@@ -363,6 +343,8 @@ void ReactionAlgorithm::generic_oneway_reaction(
       delete_particle(to_be_deleted_hidden_ids[i]); // delete particle
     }
     current_reaction.accepted_moves += 1;
+    E_pot_old = E_pot_new; // Update the system energy
+
   } else {
     // reject
     // reverse reaction

src/core/reaction_methods/ReactionAlgorithm.hpp

@@ -98,7 +98,26 @@ class ReactionAlgorithm {
 
 protected:
   std::vector<int> m_empty_p_ids_smaller_than_max_seen_particle;
-  void generic_oneway_reaction(SingleReaction &current_reaction);
+  /**
+   * @brief Carry out a generic one-way chemical reaction.
+   *
+   * Generic one way reaction of the type
+   * <tt>A+B+...+G +... --> K+...X + Z +...</tt>
+   * You need to use <tt>2A --> B</tt> instead of <tt>A+A --> B</tt> since
+   * in the latter 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.
+   *
+   * @param[in,out] current_reaction  The reaction to attempt.
+   * @param[in,out] E_pot_old         The current potential energy.
+   */
+  void generic_oneway_reaction(SingleReaction &current_reaction,
+                               double &E_pot_old);
 
   std::tuple<std::vector<StoredParticleProperty>, std::vector<int>,
              std::vector<StoredParticleProperty>>

src/core/reaction_methods/WidomInsertion.cpp

@@ -29,14 +29,14 @@
 
 namespace ReactionMethods {
 
-std::pair<double, double> WidomInsertion::measure_excess_chemical_potential(
+double WidomInsertion::calculate_particle_insertion_potential_energy(
     SingleReaction &current_reaction) {
 
   if (!all_reactant_particles_exist(current_reaction))
     throw std::runtime_error("Trying to remove some non-existing particles "
                              "from the system via the inverse Widom scheme.");
 
-  const double E_pot_old = calculate_current_potential_energy_of_system();
+  auto const E_pot_old = calculate_current_potential_energy_of_system();
 
   // make reaction attempt
   std::vector<int> p_ids_created_particles;
@@ -45,13 +45,13 @@ std::pair<double, double> WidomInsertion::measure_excess_chemical_potential(
 
   // save p_id, charge and type of the reactant particle, only thing we
   // need to hide the particle and recover it
-  const int number_of_saved_properties = 3;
+  auto constexpr number_of_saved_properties = 3;
 
   std::tie(changed_particles_properties, p_ids_created_particles,
            hidden_particles_properties) =
       make_reaction_attempt(current_reaction);
 
-  const double E_pot_new = calculate_current_potential_energy_of_system();
+  auto const E_pot_new = calculate_current_potential_energy_of_system();
   // reverse reaction attempt
   // reverse reaction
   // 1) delete created product particles
@@ -62,15 +62,11 @@ 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 / 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
-  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])};
+  // calculate the particle insertion potential energy
+  auto const E_pot_insertion = E_pot_new - E_pot_old;
+
+  return E_pot_insertion;
 }
 
 } // namespace ReactionMethods

src/core/reaction_methods/WidomInsertion.hpp

@@ -29,8 +29,8 @@ namespace ReactionMethods {
 class WidomInsertion : public ReactionAlgorithm {
 public:
   WidomInsertion(int seed) : ReactionAlgorithm(seed) {}
-  std::pair<double, double>
-  measure_excess_chemical_potential(SingleReaction &current_reaction);
+  double calculate_particle_insertion_potential_energy(
+      SingleReaction &current_reaction);
 };
 
 } // namespace ReactionMethods

src/core/reaction_methods/tests/ReactionEnsemble_test.cpp

@@ -28,9 +28,12 @@
 #include "reaction_methods/ReactionEnsemble.hpp"
 #include "reaction_methods/utils.hpp"
 
+#include "EspressoSystemStandAlone.hpp"
 #include "communication.hpp"
 #include "particle_data.hpp"
 
+#include <utils/Vector.hpp>
+
 #include <boost/mpi.hpp>
 
 #include <cmath>
@@ -39,54 +42,113 @@
 #include <memory>
 #include <stdexcept>
 
+namespace espresso {
+// ESPResSo system instance
+std::unique_ptr<EspressoSystemStandAlone> system;
+} // namespace espresso
+
 // Check the Monte Carlo algorithm where moves depend on the system
 // configuration and energy.
 BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) {
   using namespace ReactionMethods;
   class ReactionEnsembleTest : public ReactionEnsemble {
   public:
     using ReactionEnsemble::calculate_acceptance_probability;
+    using ReactionEnsemble::generic_oneway_reaction;
     using ReactionEnsemble::ReactionEnsemble;
   };
   constexpr double tol = 100 * std::numeric_limits<double>::epsilon();
 
-  ReactionEnsembleTest r_algo(42);
-  r_algo.volume = 10.;
-  r_algo.kT = 20.;
-
-  // exception if no reaction was added
-  BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error);
-
-  // create a reaction A -> 3 B + 4 C
-  int const type_A = 0;
-  int const type_B = 1;
-  int const type_C = 2;
-  SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4});
-
-  // check acceptance probability
-  constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i;
-  for (int i = 0; i < 3; ++i) {
-    for (int j = 0; j < 3; ++j) {
-      for (int k = 0; k < 3; ++k) {
-        // system contains i x A, j x B, and k x C
-        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);
-        auto const f_expr = calculate_factorial_expression(reaction, p_numbers);
-        // 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.kT);
-        auto const acceptance = r_algo.calculate_acceptance_probability(
-            reaction, energy, 0., p_numbers);
-        BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol);
+  // check basic interface
+  {
+    ReactionEnsembleTest r_algo(42);
+    r_algo.volume = 10.;
+    r_algo.kT = 20.;
+
+    // exception if no reaction was added
+    BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error);
+
+    // create a reaction A -> 3 B + 4 C
+    int const type_A = 0;
+    int const type_B = 1;
+    int const type_C = 2;
+    SingleReaction const reaction(2., {type_A}, {1}, {type_B, type_C}, {3, 4});
+
+    // check acceptance probability
+    constexpr auto g = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i;
+    for (int i = 0; i < 3; ++i) {
+      for (int j = 0; j < 3; ++j) {
+        for (int k = 0; k < 3; ++k) {
+          // system contains i x A, j x B, and k x C
+          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);
+          auto const f_expr =
+              calculate_factorial_expression(reaction, p_numbers);
+          // 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.kT);
+          auto const acceptance = r_algo.calculate_acceptance_probability(
+              reaction, energy, 0., p_numbers);
+          BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol);
+        }
       }
     }
   }
+
+  // check that the system energy is updated after a succesful reaction
+  {
+    ReactionEnsembleTest test_reaction(42);
+    test_reaction.volume = 1.;
+    test_reaction.kT = 1.;
+    test_reaction.exclusion_radius = 0;
+
+    // create a generic identity exchange reaction D <-> E
+    int const type_D = 0;
+    int const type_E = 1;
+
+    test_reaction.charges_of_types[type_D] = 0;
+    test_reaction.charges_of_types[type_E] = 0;
+
+    // track particles
+    init_type_map(type_D);
+    init_type_map(type_E);
+
+    auto const gamma = 1e100; // reaction completly shifted to product creation
+    ReactionMethods::SingleReaction reaction(gamma, {type_D}, {1}, {type_E},
+                                             {1});
+
+    // resize system box
+    auto const box_l = Utils::Vector3d{0.5, 0.4, 0.7};
+    espresso::system->set_box_l(box_l);
+
+    // create a D particle in the system
+    auto const pid = 0;
+    auto const ref_position = Utils::Vector3d{0.1, 0.2, 0.3};
+    place_particle(pid, ref_position);
+    set_particle_type(pid, type_D);
+
+    // sentinel value to check energies get updated
+    double energy = 1.0;
+
+    // for an ideal system with gamma ~ inf, the reaction is always accepted
+    test_reaction.generic_oneway_reaction(reaction, energy);
+
+    // the potential energy of the new state has to be 0 since it is an ideal
+    // system
+    double const energy_ref = 0.0;
+    BOOST_CHECK_CLOSE(energy, energy_ref, tol);
+
+    // the reaction was updated
+    BOOST_CHECK_EQUAL(reaction.tried_moves, 1);
+    BOOST_CHECK_EQUAL(reaction.accepted_moves, 1);
+  }
 }
 
 int main(int argc, char **argv) {
   auto mpi_env = std::make_shared<boost::mpi::environment>(argc, argv);
+  espresso::system = std::make_unique<EspressoSystemStandAlone>(argc, argv);
   Communication::init(mpi_env);
 
   return boost::unit_test::unit_test_main(init_unit_test, argc, argv);

src/python/espressomd/reaction_ensemble.pxd

@@ -73,4 +73,4 @@ cdef extern from "reaction_methods/WidomInsertion.hpp" namespace "ReactionMethod
 
     cdef cppclass CWidomInsertion "ReactionMethods::WidomInsertion"(CReactionAlgorithm):
         CWidomInsertion(int seed)
-        pair[double, double] measure_excess_chemical_potential(SingleReaction & current_reaction) except +
+        double calculate_particle_insertion_potential_energy(SingleReaction & current_reaction) except +