From 96b3c4f6ed446b64cff2cf37a7d0e06f035952d7 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 15 Jun 2026 11:48:35 +0200 Subject: [PATCH 1/2] analysis: reject calc_rh for chains shorter than 2 beads (bug-sweep #46) The hydrodynamic radius is computed as a sum over distinct bead pairs with a 0.5*N*(N-1) prefactor. For a single-bead chain (chain_length==1) there are no pairs: the prefactor is 0 and the accumulated inverse distance ri stays 0, so calc_rh evaluated prefac/ri = 0.0/0.0 = NaN and returned a silent NaN array instead of reporting the ill-posed request. check_topology only rejects chain_length<=0 and is shared with calc_re and calc_rg, where single-bead chains are valid (R_e==0, R_g==0), so the guard cannot live there. Add a calc_rh-specific precondition in the script-interface dispatch that throws std::domain_error (surfaced as a Python ValueError, consistent with the other chain-analysis domain checks) when chain_length < 2. Co-Authored-By: Claude Opus 4.8 --- src/script_interface/analysis/Analysis.cpp | 12 ++++++++++++ testsuite/python/analyze_chains.py | 12 ++++++++++++ 2 files changed, 24 insertions(+) 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"): From 28c6d6976feb8bfa52045fc4e37039b8f3c9c49e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 15 Jul 2026 14:31:20 +0200 Subject: [PATCH 2/2] Assertion --- src/core/analysis/statistics_chain.cpp | 2 ++ src/script_interface/analysis/Analysis.cpp | 8 +------- testsuite/python/analyze_chains.py | 2 +- 3 files changed, 4 insertions(+), 8 deletions(-) diff --git a/src/core/analysis/statistics_chain.cpp b/src/core/analysis/statistics_chain.cpp index d74d4dd943..eb4820ffb6 100644 --- a/src/core/analysis/statistics_chain.cpp +++ b/src/core/analysis/statistics_chain.cpp @@ -37,6 +37,7 @@ #include #include +#include #include #include #include @@ -200,6 +201,7 @@ std::array calc_rg(System::System const &system, int chain_start, std::array 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; diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index db6dab8418..e6cc94650d 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -211,16 +211,10 @@ 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"); + "Hydrodynamic radius is undefined for chains shorter than 2 beads"); } }); auto const result = calc_rh(system, chain_start, chain_length, n_chains); diff --git a/testsuite/python/analyze_chains.py b/testsuite/python/analyze_chains.py index eb704bd73d..ab856fcba5 100644 --- a/testsuite/python/analyze_chains.py +++ b/testsuite/python/analyze_chains.py @@ -217,7 +217,7 @@ def test_exceptions(self): # 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"): + 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