From 4d83948648c5c534f50fe7c07d1996b5b3672ba9 Mon Sep 17 00:00:00 2001
From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com>
Date: Mon, 24 Aug 2026 09:11:57 +0200
Subject: [PATCH] fix: size legacy RDF histogram from smallest box edge

The legacy bin setup derived the bin count from the largest box
edge while both histogram kernels discard every pair beyond half
the smallest edge. For non-cubic boxes this emitted a tail of bins
that could never be filled and normalized the curve over a range
the kernels never count, fabricating zero g(r) rows and large
negative pair-count residuals. Deriving the count from the global
minimum edge matches the kernel cutoff and the general path; cubic
parity fixtures are unaffected.
---
 PQAnalysis/analysis/rdf/rdf.py    |  4 +--
 tests/analysis/rdf/test_parity.py | 56 +++++++++++++++++++++++++++++++
 2 files changed, 58 insertions(+), 2 deletions(-)

diff --git a/PQAnalysis/analysis/rdf/rdf.py b/PQAnalysis/analysis/rdf/rdf.py
index ff8e37d1..dd6c3e1e 100644
--- a/PQAnalysis/analysis/rdf/rdf.py
+++ b/PQAnalysis/analysis/rdf/rdf.py
@@ -587,8 +587,8 @@ def _setup_legacy_bins(self, delta_r: PositiveReal):
         """Sets bins with the scalar semantics of ``thh_tools/RDF``."""
         self.r_min = 0.0
         self.delta_r = float(np.float32(delta_r))
-        self.n_bins = max(
-            int(max(cell.box_lengths) / 2.0 / self.delta_r)
+        self.n_bins = min(
+            int(min(cell.box_lengths) / 2.0 / self.delta_r)
             for cell in self.cells
         )
         self.r_max = self.delta_r * self.n_bins
diff --git a/tests/analysis/rdf/test_parity.py b/tests/analysis/rdf/test_parity.py
index 227777b2..b52802ed 100644
--- a/tests/analysis/rdf/test_parity.py
+++ b/tests/analysis/rdf/test_parity.py
@@ -9,6 +9,7 @@
 
 import hashlib
 import importlib
+import math
 from pathlib import Path
 
 import numpy as np
@@ -129,6 +130,61 @@ def test_float64_coordinate_controls_boundary_bin(
     assert np.array_equal(analysis.bins, np.array([1.0, 0.0, 0.0, 0.0]))
 
 
+@pytest.mark.parametrize(("frame_kernel", "batch_kernel"), KERNELS)
+def test_noncubic_box_bins_span_only_the_kernel_cutoff(
+    tmp_path,
+    monkeypatch,
+    frame_kernel,
+    batch_kernel,
+):
+    monkeypatch.setattr(
+        rdf_module, "legacy_rdf_frame_histogram", frame_kernel
+    )
+    monkeypatch.setattr(
+        rdf_module, "legacy_rdf_batch_histogram", batch_kernel
+    )
+
+    trajectory = tmp_path / "noncubic.xyz"
+    trajectory.write_text(
+        "2 20.0 20.0 60.0 90.0 90.0 90.0\n\n"
+        "X 0.0 0.0 0.0\n"
+        "Y 3.0 0.0 0.0\n",
+        encoding="utf-8",
+    )
+
+    analysis = RDF(
+        TrajectoryReader(str(trajectory)),
+        "X",
+        "Y",
+        delta_r=0.1,
+    )
+    result = analysis.run()
+
+    # The kernel drops every pair beyond min(20, 20, 60) / 2 = 10.0, so
+    # the histogram must cover int(10.0 / float32(0.1)) = 99 bins. Sizing
+    # it from the largest edge (299 bins) would normalize over a tail the
+    # kernel never counts.
+    delta_r = float(np.float32(0.1))
+
+    assert analysis._legacy_rdf  # pylint: disable=protected-access
+    assert analysis.n_bins == 99
+    assert analysis.r_max == 99 * delta_r
+
+    # The single pair sits at distance 3.0 -> bin floor(3.0 / delta_r) = 29.
+    expected_bins = np.zeros(99)
+    expected_bins[29] = 1.0
+    assert np.array_equal(analysis.bins, expected_bins)
+
+    # Hand-computed: one pair, one reference, one frame, number density
+    # 1 / (20 * 20 * 60), ideal-gas shell between 29 and 30 bin edges.
+    density = 1.0 / (20.0 * 20.0 * 60.0)
+    shell = 4.0 / 3.0 * math.pi * (
+        (30.0 * delta_r) ** 3 - (29.0 * delta_r) ** 3
+    )
+    assert result[1][29] == pytest.approx(1.0 / (density * shell), rel=1e-12)
+    assert result[2][-1] == 1.0
+
+
 @pytest.mark.parametrize("example_dir", ["msd"], indirect=False)
 def test_api_output_round_trips_exact_values(test_with_data_dir):
     Path("rdf.in").write_text(