diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index d6aa626084..db6dab8418 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/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(parameters, "chain_length"); auto const n_chains = get_value(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(result.begin(), result.end()); } diff --git a/testsuite/python/analyze_chains.py b/testsuite/python/analyze_chains.py index 09f3a37201..eb704bd73d 100644 --- a/testsuite/python/analyze_chains.py +++ b/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"):