Improve exception handling and fix various regressions (#4479)

Description of changes:
- improve exception handling in MPI code (partial fix for #4219)
   - several runtime errors used to be queued for far too long and would only surface during integration
      - runtime errors from non-bonded interactions instantiated with a cutoff too large now immediately raise an exception, from which the user needs to recover by either reducing the cutoff or disabling the interaction
      - runtime errors from virtual sites relative applied to particles too far away now immediately raise an exception, from which the user needs to recover by increasing the minimum global cutoff
   - all automatically-generated script interface methods used to call `handle_errors("")` with an empty string; now it's called with the method name to help identify which feature threw the exception
   - virtual sites fatal errors are now runtime errors
   - quaternion fatal errors are now runtime errors
   - invalid particle ids fatal errors are now value errors
- fix bugs
   - virtual sites can no longer relate to themselves
   - particles can no longer be created with a negative id
   - LB node property `boundary` no longer returns a random integer when `LB_BOUNDARIES` is not compiled in
- improve code coverage
   - common exceptions are now properly covered by tests (e.g. error messages from the cell system and integrator)
   - rotation code is now unit tested
   - h5md exceptions are now tested
- fix regressions in the testsuite
   - checkpoint tests were only running on 1 MPI core
   - the LB VTK had a wall boundary misplaced
   - virtual sites relative test cases can now run in any order
- fix regression in the OpenGL visualizer
   - activating the `LB_draw_boundaries` option without any other `LB_draw_*` option no longer throws an exception
   - rescaling the box size in `lj-demo.py` now properly rescales the particle coordinates (partial fix for #4347)

samples/grand_canonical.py

@@ -38,7 +38,7 @@
 import espressomd.reaction_ensemble
 import espressomd.electrostatics
 
-required_features = ["P3M", "EXTERNAL_FORCES", "WCA"]
+required_features = ["P3M", "WCA"]
 espressomd.assert_features(required_features)
 
 parser = argparse.ArgumentParser(epilog=__doc__ + epilog)

samples/lj-demo.py

@@ -454,7 +454,7 @@ def main_loop():
         new_box = np.ones(3) * controls.volume**(1. / 3.)
         if np.any(np.array(system.box_l) != new_box):
             for p in system.part:
-                p.pos *= new_box / system.box_l[0]
+                p.pos = p.pos * new_box / system.box_l[0]
             print("volume changed")
             system.force_cap = lj_cap
         system.box_l = new_box

src/core/grid_based_algorithms/lb_collective_interface.cpp

@@ -109,14 +109,14 @@ auto mpi_lb_get_populations(Utils::Vector3i const &index) {
 
 REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_populations)
 
-auto mpi_lb_get_boundary_flag(Utils::Vector3i const &index) {
+boost::optional<int> mpi_lb_get_boundary_flag(Utils::Vector3i const &index) {
   return detail::lb_calc(index, [&](auto index) {
 #ifdef LB_BOUNDARIES
     auto const linear_index =
         get_linear_index(lblattice.local_index(index), lblattice.halo_grid);
     return lbfields[linear_index].boundary;
 #else
-    return false;
+    return 0;
 #endif
   });
 }

src/core/integrate.cpp

@@ -61,6 +61,7 @@
 #include <boost/range/algorithm/min_element.hpp>
 
 #include <algorithm>
+#include <cassert>
 #include <cmath>
 #include <csignal>
 #include <functional>
@@ -419,6 +420,9 @@ int integrate(int n_steps, int reuse_forces) {
 
 int python_integrate(int n_steps, bool recalc_forces_par,
                      bool reuse_forces_par) {
+
+  assert(n_steps >= 0);
+
   // Override the signal handler so that the integrator obeys Ctrl+C
   SignalHandler sa(SIGINT, [](int) { ctrl_C = 1; });
 
@@ -431,12 +435,6 @@ int python_integrate(int n_steps, bool recalc_forces_par,
     reuse_forces = -1;
   }
 
-  /* go on with integrate <n_steps> */
-  if (n_steps < 0) {
-    runtimeErrorMsg() << "illegal number of steps (must be >0)";
-    return ES_ERROR;
-  }
-
   /* if skin wasn't set, do an educated guess now */
   if (!skin_set) {
     auto const max_cut = maximal_cutoff();
@@ -552,7 +550,7 @@ REGISTER_CALLBACK(mpi_set_time_step_local)
 
 void mpi_set_time_step(double time_s) {
   if (time_s <= 0.)
-    throw std::invalid_argument("time_step must be > 0.");
+    throw std::domain_error("time_step must be > 0.");
   if (lb_lbfluid_get_lattice_switch() != ActiveLB::NONE)
     check_tau_time_step_consistency(lb_lbfluid_get_tau(), time_s);
   mpi_call_all(mpi_set_time_step_local, time_s);

src/core/io/writer/h5md_core.hpp

@@ -217,9 +217,8 @@ struct incompatible_h5mdfile : public std::exception {
 
 struct left_backupfile : public std::exception {
   const char *what() const noexcept override {
-    return "A backup of the .h5 file exists. This usually means \
-that either you forgot to call the 'close' method or your simulation \
-crashed.";
+    return "A backup of the .h5 file exists. This usually means that either "
+           "you forgot to call the 'close' method or your simulation crashed.";
   }
 };
 

src/core/particle_data.cpp

@@ -487,8 +487,9 @@ void build_particle_node() { mpi_who_has(); }
  *  @brief Get the mpi rank which owns the particle with id.
  */
 int get_particle_node(int id) {
-  if (id < 0)
-    throw std::runtime_error("Invalid particle id!");
+  if (id < 0) {
+    throw std::domain_error("Invalid particle id: " + std::to_string(id));
+  }
 
   if (particle_node.empty())
     build_particle_node();
@@ -724,6 +725,9 @@ void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) {
 }
 
 int place_particle(int p_id, Utils::Vector3d const &pos) {
+  if (p_id < 0) {
+    throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
+  }
   if (particle_exists(p_id)) {
     mpi_place_particle(get_particle_node(p_id), p_id, pos);
 

src/core/reaction_methods/ReactionAlgorithm.hpp

@@ -61,18 +61,7 @@ class ReactionAlgorithm {
     }
     this->kT = kT;
     this->exclusion_range = exclusion_range;
-    for (auto const &item : exclusion_radius_per_type) {
-      auto type = item.first;
-      auto exclusion_radius = item.second;
-      if (exclusion_radius < 0) {
-        std::stringstream ss;
-        ss << "Invalid excluded_radius value for type " << type
-           << " value: " << exclusion_radius;
-        std::string error_message = ss.str();
-        throw std::domain_error(error_message);
-      }
-    }
-    this->exclusion_radius_per_type = exclusion_radius_per_type;
+    set_exclusion_radius_per_type(exclusion_radius_per_type);
     update_volume();
   }
 
@@ -109,6 +98,20 @@ class ReactionAlgorithm {
     volume = new_volume;
   }
   void update_volume();
+  void
+  set_exclusion_radius_per_type(std::unordered_map<int, double> const &map) {
+    for (auto const &item : map) {
+      auto const type = item.first;
+      auto const exclusion_radius = item.second;
+      if (exclusion_radius < 0.) {
+        throw std::domain_error("Invalid excluded_radius value for type " +
+                                std::to_string(type) + ": radius " +
+                                std::to_string(exclusion_radius));
+      }
+    }
+    exclusion_radius_per_type = map;
+  }
+
   virtual int do_reaction(int reaction_steps);
   void check_reaction_method() const;
   void remove_constraint() { m_reaction_constraint = ReactionConstraint::NONE; }

src/core/reaction_methods/tests/particle_tracking_test.cpp

@@ -57,6 +57,9 @@ BOOST_FIXTURE_TEST_CASE(particle_type_map_test, ParticleFactory) {
   // exception for random index that exceeds the number of particles
   BOOST_CHECK_THROW(get_random_p_id(type, 10), std::runtime_error);
 
+  // exception for untracked particle types
+  BOOST_CHECK_THROW(get_random_p_id(type + 1, 0), std::runtime_error);
+
   // check particle selection
   BOOST_CHECK_EQUAL(get_random_p_id(type, 0), pid);
 }

src/core/rotation.cpp

@@ -117,7 +117,8 @@ static void define_Qdd(Particle const &p, Utils::Quaternion<double> &Qd,
  *
  *  For very high angular velocities (e.g. if the product of @p time_step
  *  with the largest component of @ref ParticleMomentum::omega "p.omega()"
- *  is superior to ~2.0), the calculation might fail.
+ *  is superior to ~2.0) and for @p time_step superior or equal to unity,
+ *  the calculation might fail.
  *
  *  \todo implement for fixed_coord_flag
  */

src/core/unit_tests/CMakeLists.txt

@@ -37,6 +37,8 @@ unit_test(NAME p3m_test SRC p3m_test.cpp DEPENDS EspressoUtils)
 unit_test(NAME link_cell_test SRC link_cell_test.cpp DEPENDS EspressoUtils)
 unit_test(NAME Particle_test SRC Particle_test.cpp DEPENDS EspressoUtils
           Boost::serialization)
+unit_test(NAME rotation_test SRC rotation_test.cpp DEPENDS EspressoUtils
+          EspressoCore)
 unit_test(NAME field_coupling_couplings SRC field_coupling_couplings_test.cpp
           DEPENDS EspressoUtils)
 unit_test(NAME field_coupling_fields SRC field_coupling_fields_test.cpp DEPENDS

src/core/unit_tests/Particle_test.cpp

@@ -243,8 +243,6 @@ BOOST_AUTO_TEST_CASE(rattle_constructors) {
 }
 #endif // BOND_CONSTRAINT
 
-#ifdef EXTERNAL_FORCES
-#ifdef ROTATION
 BOOST_AUTO_TEST_CASE(particle_bitfields) {
   auto p = Particle();
 
@@ -255,24 +253,34 @@ BOOST_AUTO_TEST_CASE(particle_bitfields) {
   BOOST_CHECK(not p.can_rotate_around(1));
 
   // check setting of one axis
+#ifdef EXTERNAL_FORCES
   p.set_fixed_along(1, true);
-  p.set_can_rotate_around(1, true);
   BOOST_CHECK(p.is_fixed_along(1));
-  BOOST_CHECK(p.can_rotate_around(1));
   BOOST_CHECK(p.has_fixed_coordinates());
+#endif
+#ifdef ROTATION
+  p.set_can_rotate_around(1, true);
+  BOOST_CHECK(p.can_rotate_around(1));
   BOOST_CHECK(p.can_rotate());
+#endif
 
   // check that unsetting is properly registered
+#ifdef EXTERNAL_FORCES
   p.set_fixed_along(1, false);
-  p.set_can_rotate_around(1, false);
   BOOST_CHECK(not p.has_fixed_coordinates());
+#endif
+#ifdef ROTATION
+  p.set_can_rotate_around(1, false);
   BOOST_CHECK(not p.can_rotate());
+#endif
 
   // check setting of all flags at once
+#ifdef ROTATION
   p.set_can_rotate_all_axes();
   BOOST_CHECK(p.can_rotate_around(0));
   BOOST_CHECK(p.can_rotate_around(1));
   BOOST_CHECK(p.can_rotate_around(2));
+  p.set_cannot_rotate_all_axes();
+  BOOST_CHECK(not p.can_rotate());
+#endif
 }
-#endif // ROTATION
-#endif // EXTERNAL_FORCES