cluster_analysis: raise on fractal dimension of degenerate cluster (bug-sweep #45)

src/core/cluster_analysis/Cluster.cpp

@@ -177,6 +177,12 @@ Cluster::fractal_dimension([[maybe_unused]] double dr) {
     log_pcounts.push_back(log(static_cast<double>(subcluster_ids.size())));
     log_diameters.push_back(log(current_rg * 2.0));
   }
+  if (log_diameters.size() < 2) {
+    throw std::runtime_error(
+        "The fractal dimension of a cluster with fewer than two distinct "
+        "radial shells (e.g. a 2-particle cluster) is undefined. Provide a "
+        "larger cluster or a smaller dr.");
+  }
   double c0, c1, cov00, cov01, cov11, sumsq;
   gsl_fit_linear(log_diameters.data(), 1, log_pcounts.data(), 1,
                  log_pcounts.size(), &c0, &c1, &cov00, &cov01, &cov11, &sumsq);

testsuite/python/cluster_analysis.py

@@ -146,6 +146,27 @@ def test_single_cluster_analysis_lees_edwards(self):
         with self.assertRaises(RuntimeError):
             cluster.radius_of_gyration()
 
+    @ut.skipUnless(espressomd.has_features("GSL"),
+                   "Fractal dimension calculation requires GSL")
+    def test_fractal_dimension_small_cluster(self):
+        # A 2-particle cluster yields fewer than two radial shells, so the
+        # fractal dimension is mathematically undefined. The previous behaviour
+        # silently returned NaN (GSL linear fit with n<2 points); the call must
+        # raise a clear RuntimeError instead.
+        self.set_two_clusters()
+        dc = espressomd.pair_criteria.DistanceCriterion(cut_off=0.12)
+        self.cs.set_params(pair_criterion=dc)
+        self.cs.run_for_all_pairs()
+
+        clusters = [c for _, c in self.cs.clusters]
+        small_cluster, = [c for c in clusters if c.size() == 2]
+
+        # dr=0. gives a single fit point, dr large enough gives zero fit points;
+        # both are degenerate and must raise rather than return NaN.
+        for dr in (0., 10.):
+            with self.assertRaises(RuntimeError):
+                small_cluster.fractal_dimension(dr=dr)
+
     def test_single_cluster_analysis(self):
         # Place particles on a line (crossing periodic boundaries)
         for x in np.arange(-0.2, 0.21, 0.01):