analysis: raise instead of NaN for empty particle types (#5372)

Mass-related analysis functions now throw an exception when no particle is found.

Co-authored-by: Jean-Noël Grad <jgrad@icp.uni-stuttgart.de>
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>

Author: 3 people

src/core/analysis/statistics.cpp

@@ -47,6 +47,7 @@
 #include <functional>
 #include <limits>
 #include <numbers>
+#include <stdexcept>
 #include <tuple>
 #include <utility>
 #include <vector>
@@ -173,6 +174,13 @@ Utils::Vector3d center_of_mass(System::System const &system, int p_type) {
   double mass = 1.; // placeholder value to avoid division by zero
   boost::mpi::reduce(::comm_cart, local_com, com, std::plus<>(), 0);
   boost::mpi::reduce(::comm_cart, local_mass, mass, std::plus<>(), 0);
+  auto invalid_mass = (mass == 0.);
+  boost::mpi::broadcast(::comm_cart, invalid_mass, 0);
+  if (invalid_mass) {
+    throw std::runtime_error(
+        "Cannot calculate the center of mass: no particle with non-zero mass "
+        "of the given type(s) was found");
+  }
   return com / mass;
 }
 
@@ -200,7 +208,11 @@ Utils::Vector9d gyration_tensor(System::System const &system,
 
   Utils::Vector9d mat{};
   if (::comm_cart.rank() == 0) {
-    assert(not buf_pos.empty());
+    if (buf_pos.empty()) {
+      throw std::runtime_error(
+          "Cannot calculate the gyration tensor: no particle of the given "
+          "type(s) was found");
+    }
     auto const center =
         std::accumulate(buf_pos.begin(), buf_pos.end(), Utils::Vector3d{}) /
         static_cast<double>(buf_pos.size());

src/script_interface/analysis/Analysis.cpp

@@ -167,9 +167,13 @@ Variant Analysis::do_call_method(std::string const &name,
   }
   if (name == "center_of_mass") {
     auto const p_type = get_value<int>(parameters, "p_type");
-    context()->parallel_try_catch([&]() { check_particle_type(p_type); });
-    auto const local = center_of_mass(get_system(), p_type);
-    return mpi_reduce_sum(context()->get_comm(), local).as_vector();
+    Variant result;
+    context()->parallel_try_catch([&]() {
+      check_particle_type(p_type);
+      auto const local = center_of_mass(get_system(), p_type);
+      result = mpi_reduce_sum(context()->get_comm(), local).as_vector();
+    });
+    return result;
   }
   if (name == "angular_momentum") {
     auto const p_type = get_value<int>(parameters, "p_type");
@@ -216,17 +220,25 @@ Variant Analysis::do_call_method(std::string const &name,
   }
   if (name == "gyration_tensor") {
     auto const p_types = get_value<std::vector<int>>(parameters, "p_types");
-    for (auto const p_type : p_types) {
-      context()->parallel_try_catch([&]() { check_particle_type(p_type); });
-    }
-    auto const mat = gyration_tensor(get_system(), p_types);
-    return std::vector<double>(mat.begin(), mat.end());
+    Variant result;
+    context()->parallel_try_catch([&]() {
+      for (auto const p_type : p_types) {
+        check_particle_type(p_type);
+      }
+      auto const mat = gyration_tensor(get_system(), p_types);
+      result = std::vector<double>(mat.begin(), mat.end());
+    });
+    return result;
   }
   if (name == "moment_of_inertia_matrix") {
     auto const p_type = get_value<int>(parameters, "p_type");
-    context()->parallel_try_catch([&]() { check_particle_type(p_type); });
-    auto const local = moment_of_inertia_matrix(get_system(), p_type);
-    return mpi_reduce_sum(context()->get_comm(), local).as_vector();
+    Variant result;
+    context()->parallel_try_catch([&]() {
+      check_particle_type(p_type);
+      auto const local = moment_of_inertia_matrix(get_system(), p_type);
+      result = mpi_reduce_sum(context()->get_comm(), local).as_vector();
+    });
+    return result;
   }
   if (name == "structure_factor") {
     auto const order = get_value<int>(parameters, "sf_order");

testsuite/python/analyze_mass_related.py

@@ -148,6 +148,24 @@ def test_gyration_radius(self):
                 self.system.analysis.calc_rg(chain_start=0, number_of_chains=1,
                                              chain_length=len(self.system.part))
 
+    def test_empty_but_seen_type(self):
+        """
+        A particle type that was seen (so it passes the type-range check)
+        but currently has no particles must not divide by zero in mass-related
+        analysis routines.
+        """
+        # create a new type, make it "seen", then make it empty
+        empty_type = max(self.system.part.all().type) + 17
+        p = self.system.part.add(pos=[1., 1., 1.], type=empty_type)
+        p.remove()
+
+        with self.assertRaisesRegex(Exception, "no particle with non-zero mass of the given type"):
+            self.system.analysis.center_of_mass(p_type=empty_type)
+        with self.assertRaisesRegex(Exception, "no particle of the given type"):
+            self.system.analysis.gyration_tensor(p_type=empty_type)
+        with self.assertRaisesRegex(Exception, "no particle with non-zero mass of the given type"):
+            self.system.analysis.moment_of_inertia_matrix(p_type=empty_type)
+
 
 if __name__ == "__main__":
     ut.main()