core: Prevent impossible MC displacement moves

When the number of particles is zero, negative or larger than the
number of available particles, raise an exception.

Co-authored-by: Pablo Miguel Blanco Andrés <[email protected]>

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;
 }
 
 /**
@@ -590,20 +589,24 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {
  */
 bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type,
                                                                 int n_part) {
-  m_tried_configurational_MC_moves += 1;
-  particle_inside_exclusion_range_touched = false;
 
-  if (n_part == 0) {
-    // reject
-    return false;
+  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 >= 1");
   }
-
   auto const n_particles_of_type = number_of_particles_with_type(type);
   if (n_part > n_particles_of_type) {
-    // reject
-    return false;
+    throw std::domain_error(
+        "Parameter 'particle_number_to_be_changed' must be less or equal to "
+        "the number of particles in the system");
   }
 
+  m_tried_configurational_MC_moves += 1;
+  particle_inside_exclusion_range_touched = false;
+
   auto const E_pot_old = mpi_calculate_potential_energy();
 
   auto const original_state = generate_new_particle_positions(type, n_part);

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);

src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp

@@ -36,6 +36,7 @@
 #include <boost/mpi.hpp>
 
 #include <algorithm>
+#include <functional>
 #include <limits>
 #include <memory>
 #include <stdexcept>
@@ -194,10 +195,14 @@ 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.displacement_move_for_particles_of_type(type_C, 1));
-    BOOST_REQUIRE(!r_algo.displacement_move_for_particles_of_type(type_B, 2));
-    BOOST_REQUIRE(!r_algo.displacement_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_CHECK_THROW(displacement_move(type_C, 0), std::domain_error);
+    BOOST_CHECK_THROW(displacement_move(type_B, 2), std::domain_error);
+    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;

src/python/espressomd/reaction_methods.py

@@ -218,7 +218,9 @@ class ReactionAlgorithm(ScriptInterfaceHelper):
         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.
+        to calling this method 10 times for 1 particle. A displacement move
+        cannot take place when there are not enough particles available in
+        the system; an exception will be raised.
 
         Parameters
         ----------

testsuite/python/reaction_methods_interface.py

@@ -101,8 +101,6 @@ def check_reaction_parameters(reactions, parameters):
             self.assertAlmostEqual(method.constant_pH, 10., delta=1e-10)
             method.constant_pH = 8.
             self.assertAlmostEqual(method.constant_pH, 8., delta=1e-10)
-        self.assertFalse(method.displacement_mc_move_for_particles_of_type(
-            type_mc=0, particle_number_to_be_changed=0))
 
         # check constraints
         method.set_wall_constraints_in_z_direction(
@@ -291,6 +289,16 @@ def test_exceptions(self):
         with self.assertRaisesRegex(ValueError, "Invalid value for 'kT'"):
             espressomd.reaction_methods.ReactionEnsemble(
                 kT=-1., exclusion_range=1., seed=12)
+        for i in (-2, -1, 0):
+            with self.assertRaisesRegex(ValueError, "must be >= 1"):
+                method.displacement_mc_move_for_particles_of_type(
+                    type_mc=0, particle_number_to_be_changed=i)
+        with self.assertRaisesRegex(ValueError, "must be less or equal to the number of particles in the system"):
+            method.displacement_mc_move_for_particles_of_type(
+                type_mc=0, particle_number_to_be_changed=10000000)
+        with self.assertRaisesRegex(ValueError, "Parameter 'type_mc' must be >= 0"):
+            method.displacement_mc_move_for_particles_of_type(
+                type_mc=-1, particle_number_to_be_changed=1)
 
         # check invalid exclusion ranges and radii
         with self.assertRaisesRegex(ValueError, "Invalid value for 'exclusion_range'"):