analysis: reject calc_rh for chains shorter than 2 beads (bug-sweep #46)

src/script_interface/analysis/Analysis.cpp

@@ -211,6 +211,18 @@ 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);
+    // The hydrodynamic radius is a sum over distinct bead pairs with a
+    // 1/(N(N-1)) prefactor; for a single bead there are no pairs and the
+    // quantity is genuinely undefined (the kernel would yield 0/0 = NaN).
+    // This guard is calc_rh-specific because chain_length == 1 is a valid
+    // input for calc_re and calc_rg (single-bead R_e and R_g are both 0).
+    context()->parallel_try_catch([&]() {
+      if (chain_length < 2) {
+        throw std::domain_error(
+            "Hydrodynamic radius is not defined 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, "not defined 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"):