analysis: reject calc_rh for chains shorter than 2 beads (espressomd#5371)

Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>

Author: 2 people

src/core/analysis/statistics_chain.cpp

@@ -37,6 +37,7 @@
 #include <boost/mpi/collectives/all_reduce.hpp>
 
 #include <array>
+#include <cassert>
 #include <cmath>
 #include <functional>
 #include <stdexcept>
@@ -200,6 +201,7 @@ std::array<double, 4> calc_rg(System::System const &system, int chain_start,
 
 std::array<double, 2> calc_rh(System::System const &system, int chain_start,
                               int chain_length, int n_chains) {
+  assert(chain_length >= 2);
   auto const &cell_structure = *system.cell_structure;
   GatherPos prefetch{*system.box_geo};
   double r_H = 0.0, r_H2 = 0.0;

src/script_interface/analysis/Analysis.cpp

@@ -221,6 +221,12 @@ Variant Analysis::do_call_method(std::string const &name,
     auto const chain_length = get_value<int>(parameters, "chain_length");
     auto const n_chains = get_value<int>(parameters, "number_of_chains");
     check_topology(*system.cell_structure, chain_start, chain_length, n_chains);
+    context()->parallel_try_catch([&]() {
+      if (chain_length < 2) {
+        throw std::domain_error(
+            "Hydrodynamic radius is undefined for chains shorter than 2 beads");
+      }
+    });
     auto const result = calc_rh(system, chain_start, chain_length, n_chains);
     return std::vector<double>(result.begin(), result.end());
   }

testsuite/python/analyze_chains.py

@@ -214,6 +214,18 @@ def test_exceptions(self):
                 method(chain_start=0, number_of_chains=0, chain_length=1)
             with self.assertRaisesRegex(ValueError, "needs at least 1 chain"):
                 method(chain_start=0, number_of_chains=-1, chain_length=1)
+        # the hydrodynamic radius is defined via a sum over distinct bead pairs;
+        # for a single-bead chain there are no pairs and R_H is undefined
+        # (0/0 = NaN), so it must be rejected
+        with self.assertRaisesRegex(ValueError, "undefined for chains shorter than 2 beads"):
+            analysis.calc_rh(chain_start=0, number_of_chains=num_poly,
+                             chain_length=1)
+        # single-bead chains are well-defined for the end-to-end distance and
+        # the radius of gyration (both are 0), so these must still succeed
+        for method in (analysis.calc_re, analysis.calc_rg):
+            result = method(chain_start=0, number_of_chains=num_poly,
+                            chain_length=1)
+            self.assertTrue(np.all(np.isfinite(result)))
         self.assertIsNone(analysis.call_method("unknown"))
         if espressomd.has_features("VIRTUAL_SITES_RELATIVE"):
             with self.assertRaisesRegex(RuntimeError, "Center of mass is not well-defined"):