core: Add sanity checks in MC displacement moves

espressomd · Oct 13, 2022 · 7893235 · 7893235
1 parent 292bf32
commit 7893235
Show file tree

Hide file tree

Showing 4 changed files with 31 additions and 9 deletions.
diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/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,14 +589,23 @@ 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 (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 n_particles_of_type = number_of_particles_with_type(type);
   if (n_part > n_particles_of_type) {
     // reject

diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/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);

diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/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,15 @@ 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_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;

diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py
@@ -103,6 +103,8 @@ def check_reaction_parameters(reactions, parameters):
             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))
+        self.assertFalse(method.displacement_mc_move_for_particles_of_type(
+            type_mc=0, particle_number_to_be_changed=100000))
 
         # check constraints
         method.set_wall_constraints_in_z_direction(
@@ -291,6 +293,12 @@ def test_exceptions(self):
         with self.assertRaisesRegex(ValueError, "Invalid value for 'kT'"):
             espressomd.reaction_methods.ReactionEnsemble(
                 kT=-1., exclusion_range=1., seed=12)
+        with self.assertRaisesRegex(ValueError, "Parameter 'particle_number_to_be_changed' must be >= 0"):
+            method.displacement_mc_move_for_particles_of_type(
+                type_mc=0, particle_number_to_be_changed=-1)
+        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'"):