Reaction methods bug fixes

src/core/reaction_methods/ReactionAlgorithm.cpp

@@ -50,7 +50,7 @@ namespace ReactionMethods {
 /**
  * Performs a randomly selected reaction in the reaction ensemble
  */
-int ReactionAlgorithm::do_reaction(int reaction_steps) {
+void 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)
@@ -59,7 +59,6 @@ int ReactionAlgorithm::do_reaction(int reaction_steps) {
     int reaction_id = i_random(static_cast<int>(reactions.size()));
     generic_oneway_reaction(*reactions[reaction_id], current_E_pot);
   }
-  return 0;
 }
 
 /**
@@ -548,11 +547,11 @@ void ReactionAlgorithm::move_particle(int p_id, Utils::Vector3d const &new_pos,
   set_particle_v(p_id, vel);
 }
 
-std::vector<std::pair<int, Utils::Vector3d>>
+std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>>
 ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {
 
-  std::vector<std::pair<int, Utils::Vector3d>> old_positions;
-  old_positions.reserve(n_particles);
+  std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>> old_state;
+  old_state.reserve(n_particles);
 
   // draw particle ids at random without replacement
   int p_id = -1;
@@ -564,38 +563,58 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {
       p_id = get_random_p_id(type, random_index);
     }
     drawn_pids.emplace_back(p_id);
-    // store original position
+    // store original state
     auto const &p = get_particle_data(p_id);
-    old_positions.emplace_back(std::pair<int, Utils::Vector3d>{p_id, p.pos()});
-    // write new position
+    old_state.emplace_back(std::tuple<int, Utils::Vector3d, Utils::Vector3d>{
+        p_id, p.pos(), p.v()});
+    // write new position and new velocity
     auto const prefactor = std::sqrt(kT / p.mass());
     auto const new_pos = get_random_position_in_box();
     move_particle(p_id, new_pos, prefactor);
     check_exclusion_range(p_id);
+    if (particle_inside_exclusion_range_touched) {
+      break;
+    }
   }
 
-  return old_positions;
+  return old_state;
 }
 
 /**
- * Performs a global MC move for a particle of the provided type.
+ * Performs a global MC move for particles of a given type.
+ * Particles are selected without replacement.
+ * @param type     Type of particles to move.
+ * @param n_part   Number of particles of that type to move.
+ * @returns true if all moves were accepted.
  */
-bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
-    int type, int particle_number_of_type_to_be_changed) {
+bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type,
+                                                                int n_part) {
+
+  if (type < 0) {
+    throw std::domain_error("Parameter 'type_mc' must be >= 0");
+  }
+  if (n_part < 0) {
+    throw std::domain_error(
+        "Parameter 'particle_number_to_be_changed' must be >= 0");
+  }
+
+  if (n_part == 0) {
+    // reject
+    return false;
+  }
+
   m_tried_configurational_MC_moves += 1;
   particle_inside_exclusion_range_touched = false;
 
-  auto const particle_number_of_type = number_of_particles_with_type(type);
-  if (particle_number_of_type == 0 or
-      particle_number_of_type_to_be_changed == 0) {
+  auto const n_particles_of_type = number_of_particles_with_type(type);
+  if (n_part > n_particles_of_type) {
     // reject
     return false;
   }
 
   auto const E_pot_old = mpi_calculate_potential_energy();
 
-  auto const original_positions = generate_new_particle_positions(
-      type, particle_number_of_type_to_be_changed);
+  auto const original_state = generate_new_particle_positions(type, n_part);
 
   auto const E_pot_new = (particle_inside_exclusion_range_touched)
                              ? std::numeric_limits<double>::max()
@@ -622,9 +641,11 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
     m_accepted_configurational_MC_moves += 1;
     return true;
   }
-  // reject: restore original particle positions
-  for (auto const &item : original_positions)
+  // reject: restore original particle properties
+  for (auto const &item : original_state) {
+    set_particle_v(std::get<0>(item), std::get<2>(item));
     place_particle(std::get<0>(item), std::get<1>(item));
+  }
   return false;
 }
 

src/core/reaction_methods/ReactionAlgorithm.hpp

@@ -115,7 +115,7 @@ class ReactionAlgorithm {
     m_max_exclusion_range = max_exclusion_range;
   }
 
-  virtual int do_reaction(int reaction_steps);
+  void do_reaction(int reaction_steps);
   void check_reaction_method() const;
   void remove_constraint() { m_reaction_constraint = ReactionConstraint::NONE; }
   void set_cyl_constraint(double center_x, double center_y, double radius);
@@ -133,8 +133,7 @@ class ReactionAlgorithm {
     reactions.erase(reactions.begin() + reaction_id);
   }
 
-  bool do_global_mc_move_for_particles_of_type(int type,
-                                               int particle_number_of_type);
+  bool displacement_move_for_particles_of_type(int type, int n_part);
 
   bool particle_inside_exclusion_range_touched = false;
   bool neighbor_search_order_n = true;
@@ -165,7 +164,7 @@ class ReactionAlgorithm {
   std::tuple<std::vector<StoredParticleProperty>, std::vector<int>,
              std::vector<StoredParticleProperty>>
   make_reaction_attempt(SingleReaction const &current_reaction);
-  std::vector<std::pair<int, Utils::Vector3d>>
+  std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>>
   generate_new_particle_positions(int type, int n_particles);
   void
   restore_properties(std::vector<StoredParticleProperty> const &property_list);

src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp

@@ -36,9 +36,11 @@
 #include <boost/mpi.hpp>
 
 #include <algorithm>
+#include <functional>
 #include <limits>
 #include <memory>
 #include <stdexcept>
+#include <tuple>
 #include <utility>
 #include <vector>
 
@@ -164,14 +166,15 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
     BOOST_CHECK(r_algo.particle_inside_exclusion_range_touched);
     // check moves and bookkeeping
     for (auto const &item : bookkeeping) {
-      auto const pid = item.first;
+      auto const pid = std::get<0>(item);
       BOOST_REQUIRE(pid == 0 or pid == 1);
       auto const ref_old_pos = ref_positions[pid].first;
       auto const ref_old_vel = ref_positions[pid].second;
       auto const &p = get_particle_data(pid);
       auto const &new_pos = p.pos();
       auto const &new_vel = p.v();
-      BOOST_CHECK_EQUAL(item.second, ref_old_pos);
+      BOOST_CHECK_EQUAL(std::get<1>(item), ref_old_pos);
+      BOOST_CHECK_EQUAL(std::get<2>(item), ref_old_vel);
       BOOST_CHECK_GE(new_pos, Utils::Vector3d::broadcast(0.));
       BOOST_CHECK_LE(new_pos, Utils::Vector3d::broadcast(box_l));
       BOOST_CHECK_GE((new_pos - ref_old_pos).norm(), 0.1);
@@ -192,15 +195,20 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
     place_particle(1, ref_positions[1]);
     set_particle_type(0, type_A);
     set_particle_type(1, type_A);
-    // check early exit when a MC move cannot be performed
-    BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_C, 1));
-    BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_B, 2));
-    BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_A, 0));
+    // check runtime errors when a MC move cannot be physically performed
+    auto displacement_move =
+        std::bind(&ReactionAlgorithm::displacement_move_for_particles_of_type,
+                  &r_algo, std::placeholders::_1, std::placeholders::_2);
+    BOOST_REQUIRE(!displacement_move(type_C, 1));
+    BOOST_REQUIRE(!displacement_move(type_B, 2));
+    BOOST_REQUIRE(!displacement_move(type_A, 0));
+    BOOST_CHECK_THROW(displacement_move(type_A, -2), std::domain_error);
+    BOOST_CHECK_THROW(displacement_move(-2, 1), std::domain_error);
     // force all MC moves to be rejected by picking particles inside
     // their exclusion radius
     r_algo.exclusion_range = box_l;
     r_algo.particle_inside_exclusion_range_touched = false;
-    BOOST_REQUIRE(!r_algo.do_global_mc_move_for_particles_of_type(type_A, 2));
+    BOOST_REQUIRE(!r_algo.displacement_move_for_particles_of_type(type_A, 2));
     // check none of the particles moved
     for (auto const pid : {0, 1}) {
       auto const ref_old_pos = ref_positions[pid];
@@ -211,7 +219,7 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
     // force a MC move to be accepted by using a constant Hamiltonian
     r_algo.exclusion_range = 0.;
     r_algo.particle_inside_exclusion_range_touched = false;
-    BOOST_REQUIRE(r_algo.do_global_mc_move_for_particles_of_type(type_A, 1));
+    BOOST_REQUIRE(r_algo.displacement_move_for_particles_of_type(type_A, 1));
     std::vector<double> distances(2);
     // check that only one particle moved
     for (auto const pid : {0, 1}) {

src/python/espressomd/reaction_methods.py

@@ -210,23 +210,28 @@ class ReactionAlgorithm(ScriptInterfaceHelper):
             The number of reactions to be performed at once, defaults to 1.
 
     displacement_mc_move_for_particles_of_type()
-        Performs a displacement Monte Carlo move for particles of given type.
+        Performs displacement Monte Carlo moves for particles of a 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.
+        with a uniform probability distribution and new velocities are
+        sampled from the Maxwell-Boltzmann distribution.
+
+        The sequence of moves is only accepted if each individual move in
+        the sequence was accepted. Particles are sampled without replacement.
+        Therefore, calling this method once for 10 particles is not equivalent
+        to calling this method 10 times for 1 particle.
 
         Parameters
         ----------
         type_mc : :obj:`int`
             Particle type which should be moved
         particle_number_to_be_changed : :obj:`int`
             Number of particles to move, defaults to 1.
+            Particles are selected without replacement.
 
         Returns
         -------
         :obj:`bool`
-            Whether the move was accepted.
+            Whether all moves were accepted.
 
     delete_particle()
         Deletes the particle of the given p_id and makes sure that the particle

src/script_interface/particle_data/ParticleHandle.cpp

@@ -483,6 +483,9 @@ void ParticleHandle::do_construct(VariantMap const &params) {
         do_set_parameter(kv.first, kv.second);
       }
     }
+    if (params.count("type") == 0) {
+      do_set_parameter("type", 0);
+    }
   }
 }
 

src/script_interface/reaction_methods/ReactionAlgorithm.hpp

@@ -128,7 +128,7 @@ class ReactionAlgorithm : public AutoParameters<ReactionAlgorithm> {
     } else if (name == "reaction") {
       RE()->do_reaction(get_value_or<int>(parameters, "reaction_steps", 1));
     } else if (name == "displacement_mc_move_for_particles_of_type") {
-      return RE()->do_global_mc_move_for_particles_of_type(
+      return RE()->displacement_move_for_particles_of_type(
           get_value<int>(parameters, "type_mc"),
           get_value_or<int>(parameters, "particle_number_to_be_changed", 1));
     } else if (name == "check_reaction_method") {

testsuite/python/CMakeLists.txt

@@ -177,6 +177,7 @@ python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4)
 python_test(FILE widom_insertion.py MAX_NUM_PROC 1)
 python_test(FILE constant_pH.py MAX_NUM_PROC 1)
 python_test(FILE constant_pH_stats.py MAX_NUM_PROC 4 LABELS long)
+python_test(FILE canonical_ensemble.py MAX_NUM_PROC 2)
 python_test(FILE writevtf.py MAX_NUM_PROC 4)
 python_test(FILE lb_stokes_sphere.py MAX_NUM_PROC 4 LABELS gpu long)
 python_test(FILE lb_pressure_tensor.py MAX_NUM_PROC 1 LABELS gpu long)

testsuite/python/canonical_ensemble.py

@@ -0,0 +1,119 @@
+#
+# Copyright (C) 2022 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# 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/>.
+#
+
+import espressomd
+import espressomd.constraints
+import espressomd.observables
+import espressomd.accumulators
+import espressomd.reaction_methods
+
+import numpy as np
+import scipy.optimize
+import unittest as ut
+import unittest_decorators as utx
+
+
+class Test(ut.TestCase):
+
+    system = espressomd.System(box_l=[1., 1., 1.])
+    system.cell_system.skin = 0.4
+    system.setup_type_map(type_list=[0])
+
+    def tearDown(self):
+        self.system.part.clear()
+        self.system.constraints.clear()
+
+    @utx.skipIfMissingFeatures("ELECTROSTATICS")
+    def test_linear_potential(self):
+        """
+        Set a particle in a box with a linear potential along the x-axis.
+        The particle distribution resulting from accepted Monte Carlo moves
+        should follow a Maxwell-Boltzmann distribution.
+        """
+
+        method = espressomd.reaction_methods.ReactionEnsemble(
+            kT=0.2, seed=42, exclusion_range=0., search_algorithm="order_n")
+        method.set_non_interacting_type(type=1)
+
+        p = self.system.part.add(pos=[0., 0., 0.], q=1, type=0)
+        obs_pos = espressomd.observables.ParticlePositions(ids=(p.id,))
+        obs_vel = espressomd.observables.ParticleVelocities(ids=(p.id,))
+        acc_pos = espressomd.accumulators.TimeSeries(obs=obs_pos)
+        acc_vel = espressomd.accumulators.TimeSeries(obs=obs_vel)
+
+        E = np.array([-1., 0., 0.])
+        field = espressomd.constraints.LinearElectricPotential(E=E, phi0=0.)
+        self.system.constraints.add(field)
+
+        for _ in range(5000):
+            accepted = method.displacement_mc_move_for_particles_of_type(
+                type_mc=0, particle_number_to_be_changed=1)
+            if accepted:
+                acc_pos.update()
+                acc_vel.update()
+                p.pos = [0., 0., 0.]
+                p.v = [0., 0., 0.]
+            else:
+                self.assertAlmostEqual(np.linalg.norm(p.v), 0., delta=1e-12)
+
+        # the x-position should follow an exponential distribution
+        # -> mean = kT, median = kT x ln(2), variance = kT^2
+        series = acc_pos.time_series()[:, p.id, 0]
+        ydata, xbins = np.histogram(series, bins=15, range=[0., 1.])
+        xdata = (xbins[1:] + xbins[:-1]) / 2.
+        ydata = ydata / float(ydata[0])
+        (a, b, c), _ = scipy.optimize.curve_fit(
+            lambda x, a, b, c: a * np.exp(-b * x) + c, xdata, ydata)
+        # check histogram profile is roughly exponential
+        self.assertAlmostEqual(a, 1., delta=0.2)
+        self.assertAlmostEqual(b, 1. / method.kT, delta=0.3)
+        self.assertAlmostEqual(c, 0., delta=0.01)
+        # check distribution parameters with high accuracy
+        ln2 = np.log(2)
+        self.assertAlmostEqual(np.mean(series), method.kT, delta=0.02)
+        self.assertAlmostEqual(np.median(series) / ln2, method.kT, delta=0.02)
+        self.assertAlmostEqual(np.sqrt(np.var(series)), method.kT, delta=0.02)
+
+        # the y- and z-position should follow a uniform distribution
+        for axis in (1, 2):
+            series = acc_pos.time_series()[:, p.id, axis]
+            ydata, _ = np.histogram(series, bins=10, range=[0., 1.])
+            ydata = ydata / np.mean(ydata)
+            np.testing.assert_allclose(ydata, 1., atol=0.25)
+
+        # the velocity vector should follow a normal distribution
+        # -> mean = 0, median = 0, variance = kT
+        for axis in (0, 1, 2):
+            series = acc_vel.time_series()[:, p.id, axis]
+            ydata, xbins = np.histogram(series, bins=25, range=[-1.5, 1.5])
+            xdata = (xbins[1:] + xbins[:-1]) / 2.
+            ydata = ydata / len(series)
+            (_, b, c), _ = scipy.optimize.curve_fit(
+                lambda x, a, b, c: a * np.exp(-b * x**2) + c, xdata, ydata)
+            # check histogram profile is roughly gaussian
+            self.assertAlmostEqual(b, 0.5 / method.kT, delta=0.45)
+            self.assertAlmostEqual(c, 0., delta=0.002)
+            # check distribution parameters with high accuracy
+            self.assertAlmostEqual(np.mean(series), 0., delta=0.05)
+            self.assertAlmostEqual(np.median(series), 0., delta=0.025)
+            self.assertAlmostEqual(np.var(series), method.kT, delta=0.025)
+
+
+if __name__ == "__main__":
+    ut.main()