From e58a4f0a4ec6a13e036aa7cd8a747f9e05bf00f2 Mon Sep 17 00:00:00 2001
From: Brady Johnston <brady.johnston@me.com>
Date: Fri, 9 Oct 2026 12:01:43 +0100
Subject: [PATCH 1/7] fix triclinic pkdtree distances

---
 package/CHANGELOG                             |  6 ++
 .../MDAnalysis/lib/include/calc_distances.h   |  4 +-
 package/MDAnalysis/lib/nsgrid.pyx             | 37 +++++------
 .../MDAnalysisTests/lib/test_distances.py     | 62 +++++++++++++++++++
 testsuite/MDAnalysisTests/lib/test_nsgrid.py  | 33 ++++++++++
 5 files changed, 118 insertions(+), 24 deletions(-)

diff --git a/package/CHANGELOG b/package/CHANGELOG
index dede22e74ab..3bb64be96f8 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -22,6 +22,12 @@ The rules for this file:
  * 2.11.0
 
 Fixes
+ * Fixed `capped_distance` and `self_capped_distance` (and hence bond
+   guessing, selections and other users of the `nsgrid` and `pkdtree`
+   methods) missing pairs within the cutoff in triclinic boxes where both
+   b_x and c_y are non-zero, e.g. truncated octahedra; `apply_PBC` also
+   wrapped some coordinates in such boxes to just outside the primary
+   image (Issue #4906)
  * `MDAnalysis.analysis.nucleicacids.WatsonCrickDist`, `MinorPairDist`,
    and `MajorPairDist` now match residue names against the full resname
    instead of only the first character, fixing incorrect behaviour with
diff --git a/package/MDAnalysis/lib/include/calc_distances.h b/package/MDAnalysis/lib/include/calc_distances.h
index ff130d2c049..106996740ac 100644
--- a/package/MDAnalysis/lib/include/calc_distances.h
+++ b/package/MDAnalysis/lib/include/calc_distances.h
@@ -294,7 +294,7 @@ static void _triclinic_pbc(coordinate* coords, uint64_t numcoords, float* box)
             }
             if (!msr) {
                 // translate remainder of crd to central cell along a-axis
-                lbound = crd[1] * a_ax_yfactor + crd[2] * a_ax_zfactor;
+                lbound = (crd[1] - crd[2] * b_ax_zfactor) * a_ax_yfactor + crd[2] * a_ax_zfactor;
                 ubound = lbound + box[0];
                 if (crd[0] < lbound) {
                     crd[0] += box[0];
@@ -356,7 +356,7 @@ static void _triclinic_pbc(coordinate* coords, uint64_t numcoords, float* box)
                 crd[1] -= box[4];
             }
             // translate remainder of crd to central cell along a-axis
-            lbound = crd[1] * a_ax_yfactor + crd[2] * a_ax_zfactor;
+            lbound = (crd[1] - crd[2] * b_ax_zfactor) * a_ax_yfactor + crd[2] * a_ax_zfactor;
             ubound = lbound + box[0];
             if (crd[0] < lbound) {
                 crd[0] += box[0];
diff --git a/package/MDAnalysis/lib/nsgrid.pyx b/package/MDAnalysis/lib/nsgrid.pyx
index a4e245e884e..46c259bb384 100644
--- a/package/MDAnalysis/lib/nsgrid.pyx
+++ b/package/MDAnalysis/lib/nsgrid.pyx
@@ -105,13 +105,6 @@ cdef extern from "calc_distances.h" nogil:
     void _triclinic_pbc(coordinate* coords, int numcoords, float* box)
 
 
-cdef inline float fmax(float a, float b):
-    if a > b:
-        return a
-    else:
-        return b
-
-
 cdef inline float fmin(float a, float b):
     if a < b:
         return a
@@ -305,7 +298,7 @@ cdef class FastNS(object):
         max_cutoff : float
            the maximum allowable cutoff given the box shape and size
         """
-        cdef float cutoff, min_cellsize, max_cutoff, new_cellsize
+        cdef float cutoff, min_cellsize, max_cutoff
         cdef int i
 
         from MDAnalysis.lib.mdamath import triclinic_vectors
@@ -337,22 +330,20 @@ cdef class FastNS(object):
         max_cutoff = fmin(max_cutoff, self.triclinic_dimensions[YY] * degsin(self.dimensions[5]))
         max_cutoff /= 2
 
-        # for triclinic cells, we need to worry about the shortest path across the cells
-        min_cellsize = cutoff
-        if self.triclinic:
-            for i in range(3, 6):
-                # cutoff/sin(theta) to elongate the XX/YY/ZZ dimension to make smallest diagonal large enough
-                new_cellsize = cutoff / degsin(self.dimensions[i])
-                min_cellsize = fmax(new_cellsize, min_cellsize)
-
+        # Cells must be at least `cutoff` thick perpendicular to each pair of
+        # box faces, so size the grid by the box's perpendicular widths.
+        tv = np.asarray(self.triclinic_dimensions, dtype=np.float64).reshape(3, 3)
+        volume = abs(np.linalg.det(tv))
+        widths = [volume / np.linalg.norm(np.cross(tv[1], tv[2])),
+                  volume / np.linalg.norm(np.cross(tv[2], tv[0])),
+                  volume / np.linalg.norm(np.cross(tv[0], tv[1]))]
         # add 0.001 here to avoid floating point errors
         # will make cells slightly too large as a result, ah well
-        min_cellsize += 0.001
+        min_cellsize = cutoff + 0.001
         # If the cell size is too small, indexing overflow will occur. Limit the number
         # of cells in any dimension to the cube root of the maximum of 32 bit integer values.
-        self.ncells[0] = <int> min(math.floor(self.triclinic_dimensions[XX] / min_cellsize), MAX_GRID_DIM)
-        self.ncells[1] = <int> min(math.floor(self.triclinic_dimensions[YY] / min_cellsize), MAX_GRID_DIM)
-        self.ncells[2] = <int> min(math.floor(self.triclinic_dimensions[ZZ] / min_cellsize), MAX_GRID_DIM)
+        for i in range(3):
+            self.ncells[i] = <int> min(math.floor(widths[i] / min_cellsize), MAX_GRID_DIM)
 
         self.pbc = pbc
         # If there aren't enough cells in a given dimension it's equivalent to one
@@ -436,8 +427,10 @@ cdef class FastNS(object):
         """Calculate cell coordinate for coord"""
         # This assumes coordinate is inside the primary unit cell
         xyz[2] = <int> (coord[2] / self.cellsize[ZZ])
-        xyz[1] = <int> ((coord[1] - coord[2] * self.cellsize[YZ]) / self.cellsize[YY])
-        xyz[0] = <int> ((coord[0] - coord[1] * self.cellsize[XY]
+        # y with the c-vector contribution removed, needed for both y and x
+        cdef double y = coord[1] - coord[2] * self.cellsize[YZ]
+        xyz[1] = <int> (y / self.cellsize[YY])
+        xyz[0] = <int> ((coord[0] - y * self.cellsize[XY]
                          - coord[2] * self.cellsize[XZ]) / self.cellsize[XX])
         # Make sure cell coordinate indices are within the primary unit cell
         # (better safe than sorry):
diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index 06abe0a99c7..853f8674afb 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -84,6 +84,26 @@ def test_transform_StoR_pass(coord_dtype):
     assert_allclose(original_r, test_r)
 
 
+# Triclinic boxes with both b_x and c_y non-zero, see Issue #4906
+SKEWED_TRICLINIC_BOXES = (
+    # truncated octahedron
+    np.array([97.16, 97.16, 97.16, 109.47, 109.47, 109.47], dtype=np.float32),
+    # LAMMPS cell from Issue #4906
+    np.array(
+        [15.173525, 14.923327, 14.475218, 109.888, 109.3786, 115.684875],
+        dtype=np.float32,
+    ),
+    np.array([26.24, 27.86, 34.66, 106.56, 77.20, 79.13], dtype=np.float32),
+)
+
+
+def _skewed_box_coords(box, n, seed):
+    """Random coordinates spread over the box and its neighbouring images"""
+    rng = np.random.default_rng(seed)
+    frac = rng.uniform(-1.5, 2.5, size=(n, 3))
+    return (frac @ mdamath.triclinic_vectors(box)).astype(np.float32)
+
+
 class TestCappedDistances(object):
 
     npoints_1 = (1, 100)
@@ -122,6 +142,38 @@ def query_2_atomgroup(self, query_2):
 
     min_cutoff_1 = (None, 0.1)
 
+    @pytest.mark.parametrize("box", SKEWED_TRICLINIC_BOXES)
+    @pytest.mark.parametrize("method", ("pkdtree", "nsgrid"))
+    def test_capped_distance_skewed_triclinic(self, box, method):
+        # Issue #4906: pkdtree and nsgrid missed pairs in skewed boxes
+        width = mdamath.triclinic_vectors(box).diagonal().min()
+        cutoff = 0.1 * width
+        a = _skewed_box_coords(box, 1000, seed=1)
+        b = _skewed_box_coords(box, 1000, seed=2)
+
+        def pairs(m):
+            p = distances.capped_distance(
+                a, b, cutoff, box=box, method=m, return_distances=False
+            )
+            return {tuple(x) for x in p}
+
+        assert pairs(method) == pairs("bruteforce")
+
+    @pytest.mark.parametrize("box", SKEWED_TRICLINIC_BOXES)
+    @pytest.mark.parametrize("method", ("pkdtree", "nsgrid"))
+    def test_self_capped_distance_skewed_triclinic(self, box, method):
+        width = mdamath.triclinic_vectors(box).diagonal().min()
+        cutoff = 0.1 * width
+        a = _skewed_box_coords(box, 1000, seed=3)
+
+        def pairs(m):
+            p = distances.self_capped_distance(
+                a, cutoff, box=box, method=m, return_distances=False
+            )
+            return {tuple(sorted(x)) for x in p}
+
+        assert pairs(method) == pairs("bruteforce")
+
     def test_capped_distance_noresults(self):
         point1 = np.array([0.1, 0.1, 0.1], dtype=np.float32)
         point2 = np.array([0.95, 0.1, 0.1], dtype=np.float32)
@@ -1787,6 +1839,16 @@ def test_coords_in_central_image_tric(self, backend):
         assert np.all(relres >= 0.0)
         assert np.all(relres < 1.0)
 
+    @pytest.mark.parametrize("box", SKEWED_TRICLINIC_BOXES)
+    def test_coords_in_central_image_skewed_tric(self, backend, box):
+        # Issue #4906: boxes where both b_x and c_y are non-zero were wrapped
+        # to just outside the primary image along a.
+        coords = _skewed_box_coords(box, 1000, seed=4906)
+        res = distances.apply_PBC(coords, box, backend=backend)
+        relres = distances.transform_RtoS(res, box)
+        assert np.all(relres >= -1e-6)
+        assert np.all(relres < 1.0 + 1e-6)
+
 
 @pytest.mark.parametrize("backend", ["serial", "openmp"])
 class TestPeriodicAngles(object):
diff --git a/testsuite/MDAnalysisTests/lib/test_nsgrid.py b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
index 582e780172b..de3536bdd3d 100644
--- a/testsuite/MDAnalysisTests/lib/test_nsgrid.py
+++ b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
@@ -505,6 +505,39 @@ def test_issue_2670():
     assert len(ag2.select_atoms("around 0.0 resid 3")) == 1
 
 
+@pytest.mark.parametrize(
+    "box",
+    [
+        # truncated octahedron
+        [97.16, 97.16, 97.16, 109.47, 109.47, 109.47],
+        # LAMMPS cell from Issue #4906
+        [15.173525, 14.923327, 14.475218, 109.888, 109.3786, 115.684875],
+    ],
+)
+def test_issue_4906(box):
+    # In boxes with both b_x and c_y non-zero, atoms were binned into the
+    # wrong cells and cells could be thinner than the cutoff, so neighbours
+    # within the cutoff were missed.
+    box = np.array(box, dtype=np.float32)
+    rng = np.random.default_rng(4906)
+    frac = rng.uniform(0, 1, size=(1000, 3))
+    coords = (frac @ mda.lib.mdamath.triclinic_vectors(box)).astype(np.float32)
+    cutoff = 0.1 * box[:3].min()
+
+    dists = mda.lib.distances.distance_array(coords, coords, box=box)
+
+    searcher = nsgrid.FastNS(cutoff, coords, box)
+    pairs = searcher.self_search().get_pairs()
+    i, j = np.triu_indices(len(coords), k=1)
+    within = dists[i, j] <= cutoff
+    expected = set(zip(i[within], j[within]))
+    assert {tuple(sorted(p)) for p in pairs} == expected
+
+    pairs = searcher.search(coords).get_pairs()
+    i, j = np.nonzero(dists <= cutoff)
+    assert {tuple(p) for p in pairs} == set(zip(i, j))
+
+
 def high_mem_tests_enabled():
     """Returns true if ENABLE_HIGH_MEM_UNIT_TESTS is set to true."""
     env = os.getenv("ENABLE_HIGH_MEM_UNIT_TESTS", default="false").lower()

From a3e37a85dcf5b0af782be36b94ff02b693d9bcd5 Mon Sep 17 00:00:00 2001
From: Brady Johnston <brady.johnston@me.com>
Date: Fri, 9 Oct 2026 12:23:01 +0100
Subject: [PATCH 2/7] update incorrect expectations

---
 testsuite/MDAnalysisTests/analysis/test_dielectric.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/testsuite/MDAnalysisTests/analysis/test_dielectric.py b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
index 45d8e722cdb..3d22d364184 100644
--- a/testsuite/MDAnalysisTests/analysis/test_dielectric.py
+++ b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
@@ -41,7 +41,8 @@ def test_broken_molecules(self, ag):
             ag.wrap()
 
         eps = DielectricConstant(ag, make_whole=False).run()
-        assert_allclose(eps.results["eps_mean"], 721.711, rtol=1e-03)
+        # reference values changed after the apply_PBC fix for Issue #4906
+        assert_allclose(eps.results["eps_mean"], 630.326, rtol=1e-03)
 
     def test_broken_repaired_molecules(self, ag):
         # cut molecules apart
@@ -50,7 +51,8 @@ def test_broken_repaired_molecules(self, ag):
             ag.wrap()
 
         eps = DielectricConstant(ag, make_whole=True).run()
-        assert_allclose(eps.results["eps_mean"], 5.088, rtol=1e-03)
+        # reference values changed after the apply_PBC fix for Issue #4906
+        assert_allclose(eps.results["eps_mean"], 3.874, rtol=1e-03)
 
     def test_temperature(self, ag):
         eps = DielectricConstant(ag, temperature=100).run()

From 4940789fedd9ffbba600a1c07fbc5830388e973b Mon Sep 17 00:00:00 2001
From: Brady Johnston <brady.johnston@me.com>
Date: Thu, 15 Oct 2026 11:01:52 +0100
Subject: [PATCH 3/7] fix box widths slowdown; cleanup tests and add
 versionchanged to docs

---
 package/MDAnalysis/lib/distances.py           |  4 +++
 .../MDAnalysis/lib/include/calc_distances.h   |  6 ++---
 package/MDAnalysis/lib/nsgrid.pyx             | 27 +++++++++++++++----
 .../analysis/test_dielectric.py               |  6 ++++-
 .../MDAnalysisTests/lib/test_distances.py     | 25 ++++++++---------
 testsuite/MDAnalysisTests/lib/test_nsgrid.py  | 17 +++++++-----
 6 files changed, 57 insertions(+), 28 deletions(-)

diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index b3f862e005b..5e7ce4aec04 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -2086,6 +2086,10 @@ def apply_PBC(
     .. versionchanged:: 2.3.0
        Can now accept an :class:`~MDAnalysis.core.groups.AtomGroup` as an
        argument in any position and checks inputs using type hinting.
+    .. versionchanged:: 2.11.0
+       Fixed wrapping in triclinic boxes where both the *x* component of the
+       second box vector and the *y* component of the third are non-zero
+       (Issue #4906).
     """
     # coords is an array, the check_coords decorator made sure of that.
     # Mypy, however, is not aware of that so we have to tell it explicitly.
diff --git a/package/MDAnalysis/lib/include/calc_distances.h b/package/MDAnalysis/lib/include/calc_distances.h
index 106996740ac..1b9a6b339e6 100644
--- a/package/MDAnalysis/lib/include/calc_distances.h
+++ b/package/MDAnalysis/lib/include/calc_distances.h
@@ -225,8 +225,6 @@ static void _triclinic_pbc(coordinate* coords, uint64_t numcoords, float* box)
     const double bi6 = (-bi3 * box[7] - box[6] * bi0) * bi8;
     const double bi7 = -box[7] * bi4 * bi8;
     // variables and constants for single box shifts:
-    const double a_ax_yfactor = (double) box[3] * bi4;;
-    const double a_ax_zfactor = (double) box[6] * bi8;
     const double b_ax_zfactor = (double) box[7] * bi8;
 
 
@@ -294,7 +292,7 @@ static void _triclinic_pbc(coordinate* coords, uint64_t numcoords, float* box)
             }
             if (!msr) {
                 // translate remainder of crd to central cell along a-axis
-                lbound = (crd[1] - crd[2] * b_ax_zfactor) * a_ax_yfactor + crd[2] * a_ax_zfactor;
+                lbound = -(crd[1] * bi3 + crd[2] * bi6) * box[0];
                 ubound = lbound + box[0];
                 if (crd[0] < lbound) {
                     crd[0] += box[0];
@@ -356,7 +354,7 @@ static void _triclinic_pbc(coordinate* coords, uint64_t numcoords, float* box)
                 crd[1] -= box[4];
             }
             // translate remainder of crd to central cell along a-axis
-            lbound = (crd[1] - crd[2] * b_ax_zfactor) * a_ax_yfactor + crd[2] * a_ax_zfactor;
+            lbound = -(crd[1] * bi3 + crd[2] * bi6) * box[0];
             ubound = lbound + box[0];
             if (crd[0] < lbound) {
                 crd[0] += box[0];
diff --git a/package/MDAnalysis/lib/nsgrid.pyx b/package/MDAnalysis/lib/nsgrid.pyx
index 46c259bb384..afd3c0e06f0 100644
--- a/package/MDAnalysis/lib/nsgrid.pyx
+++ b/package/MDAnalysis/lib/nsgrid.pyx
@@ -193,6 +193,10 @@ cdef class FastNS(object):
 
     .. versionchanged::  1.0.2
        Rewrote to fix bugs with triclinic boxes
+    .. versionchanged:: 2.11.0
+       Fixed missed neighbours in triclinic boxes where both the *x*
+       component of the second box vector and the *y* component of the third
+       are non-zero (Issue #4906).
     """
     cdef readonly double cutoff
     cdef float[:, ::1] coords_bbox
@@ -299,6 +303,8 @@ cdef class FastNS(object):
            the maximum allowable cutoff given the box shape and size
         """
         cdef float cutoff, min_cellsize, max_cutoff
+        cdef double ax, bx, by, cx, cy, cz, volume
+        cdef double widths[3]
         cdef int i
 
         from MDAnalysis.lib.mdamath import triclinic_vectors
@@ -332,11 +338,22 @@ cdef class FastNS(object):
 
         # Cells must be at least `cutoff` thick perpendicular to each pair of
         # box faces, so size the grid by the box's perpendicular widths.
-        tv = np.asarray(self.triclinic_dimensions, dtype=np.float64).reshape(3, 3)
-        volume = abs(np.linalg.det(tv))
-        widths = [volume / np.linalg.norm(np.cross(tv[1], tv[2])),
-                  volume / np.linalg.norm(np.cross(tv[2], tv[0])),
-                  volume / np.linalg.norm(np.cross(tv[0], tv[1]))]
+        # The box matrix is lower triangular, so width = volume / face area
+        # has a closed form.
+        ax = self.triclinic_dimensions[XX]
+        bx = self.triclinic_dimensions[XY]
+        by = self.triclinic_dimensions[YY]
+        cx = self.triclinic_dimensions[XZ]
+        cy = self.triclinic_dimensions[YZ]
+        cz = self.triclinic_dimensions[ZZ]
+        volume = ax * by * cz
+        widths[0] = widths[1] = widths[2] = 0.0
+        # invalid boxes give an all-zero matrix, leave widths at zero
+        if volume > 0:
+            widths[0] = volume / math.sqrt((by * cz) ** 2 + (bx * cz) ** 2
+                                           + (bx * cy - by * cx) ** 2)
+            widths[1] = by * cz / math.sqrt(cy * cy + cz * cz)
+            widths[2] = cz
         # add 0.001 here to avoid floating point errors
         # will make cells slightly too large as a result, ah well
         min_cellsize = cutoff + 0.001
diff --git a/testsuite/MDAnalysisTests/analysis/test_dielectric.py b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
index 3d22d364184..ab7735a9159 100644
--- a/testsuite/MDAnalysisTests/analysis/test_dielectric.py
+++ b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
@@ -45,13 +45,17 @@ def test_broken_molecules(self, ag):
         assert_allclose(eps.results["eps_mean"], 630.326, rtol=1e-03)
 
     def test_broken_repaired_molecules(self, ag):
+        intact = DielectricConstant(ag, make_whole=True).run()
         # cut molecules apart
         ag.universe.transfer_to_memory()
         for ts in ag.universe.trajectory:
             ag.wrap()
 
+        # repairing the molecules should recover the intact result
         eps = DielectricConstant(ag, make_whole=True).run()
-        # reference values changed after the apply_PBC fix for Issue #4906
+        assert_allclose(
+            eps.results["eps_mean"], intact.results["eps_mean"], rtol=1e-03
+        )
         assert_allclose(eps.results["eps_mean"], 3.874, rtol=1e-03)
 
     def test_temperature(self, ag):
diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index 853f8674afb..bf95e1587f9 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -104,6 +104,11 @@ def _skewed_box_coords(box, n, seed):
     return (frac @ mdamath.triclinic_vectors(box)).astype(np.float32)
 
 
+def _away_from_cutoff(pairs, dists, cutoff, tol=1e-4):
+    """Drop pairs so close to the cutoff that rounding can decide them"""
+    return pairs[np.abs(dists - cutoff) > tol]
+
+
 class TestCappedDistances(object):
 
     npoints_1 = (1, 100)
@@ -146,31 +151,27 @@ def query_2_atomgroup(self, query_2):
     @pytest.mark.parametrize("method", ("pkdtree", "nsgrid"))
     def test_capped_distance_skewed_triclinic(self, box, method):
         # Issue #4906: pkdtree and nsgrid missed pairs in skewed boxes
-        width = mdamath.triclinic_vectors(box).diagonal().min()
-        cutoff = 0.1 * width
+        min_diag = mdamath.triclinic_vectors(box).diagonal().min()
+        cutoff = 0.1 * min_diag
         a = _skewed_box_coords(box, 1000, seed=1)
         b = _skewed_box_coords(box, 1000, seed=2)
 
         def pairs(m):
-            p = distances.capped_distance(
-                a, b, cutoff, box=box, method=m, return_distances=False
-            )
-            return {tuple(x) for x in p}
+            p, d = distances.capped_distance(a, b, cutoff, box=box, method=m)
+            return {tuple(x) for x in _away_from_cutoff(p, d, cutoff)}
 
         assert pairs(method) == pairs("bruteforce")
 
     @pytest.mark.parametrize("box", SKEWED_TRICLINIC_BOXES)
     @pytest.mark.parametrize("method", ("pkdtree", "nsgrid"))
     def test_self_capped_distance_skewed_triclinic(self, box, method):
-        width = mdamath.triclinic_vectors(box).diagonal().min()
-        cutoff = 0.1 * width
+        min_diag = mdamath.triclinic_vectors(box).diagonal().min()
+        cutoff = 0.1 * min_diag
         a = _skewed_box_coords(box, 1000, seed=3)
 
         def pairs(m):
-            p = distances.self_capped_distance(
-                a, cutoff, box=box, method=m, return_distances=False
-            )
-            return {tuple(sorted(x)) for x in p}
+            p, d = distances.self_capped_distance(a, cutoff, box=box, method=m)
+            return {tuple(sorted(x)) for x in _away_from_cutoff(p, d, cutoff)}
 
         assert pairs(method) == pairs("bruteforce")
 
diff --git a/testsuite/MDAnalysisTests/lib/test_nsgrid.py b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
index de3536bdd3d..7c46029dd59 100644
--- a/testsuite/MDAnalysisTests/lib/test_nsgrid.py
+++ b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
@@ -525,17 +525,22 @@ def test_issue_4906(box):
     cutoff = 0.1 * box[:3].min()
 
     dists = mda.lib.distances.distance_array(coords, coords, box=box)
+    # skip pairs so close to the cutoff that rounding can decide them
+    ambiguous = np.abs(dists - cutoff) <= 1e-4
+
+    def found(results):
+        pairs = results.get_pairs()
+        return {tuple(p) for p in pairs if not ambiguous[p[0], p[1]]}
 
     searcher = nsgrid.FastNS(cutoff, coords, box)
-    pairs = searcher.self_search().get_pairs()
     i, j = np.triu_indices(len(coords), k=1)
-    within = dists[i, j] <= cutoff
+    within = (dists[i, j] <= cutoff) & ~ambiguous[i, j]
     expected = set(zip(i[within], j[within]))
-    assert {tuple(sorted(p)) for p in pairs} == expected
+    pairs = {tuple(sorted(p)) for p in found(searcher.self_search())}
+    assert pairs == expected
 
-    pairs = searcher.search(coords).get_pairs()
-    i, j = np.nonzero(dists <= cutoff)
-    assert {tuple(p) for p in pairs} == set(zip(i, j))
+    i, j = np.nonzero((dists <= cutoff) & ~ambiguous)
+    assert found(searcher.search(coords)) == set(zip(i, j))
 
 
 def high_mem_tests_enabled():

From a7b88e128efa350544776008a0933ba54d637618 Mon Sep 17 00:00:00 2001
From: Oliver Beckstein <orbeckst@gmail.com>
Date: Mon, 19 Oct 2026 01:21:28 -0700
Subject: [PATCH 4/7] additional docs for PBC fix

- extensive CHANGELOG entry
- versionchanged (because the fix is important)
- code docs/comment updates
---
 package/CHANGELOG                             | 24 ++++++++---
 .../MDAnalysis/lib/include/calc_distances.h   | 16 +++----
 package/MDAnalysis/lib/nsgrid.pyx             | 43 ++++++++++++++++++-
 .../analysis/test_dielectric.py               |  3 +-
 4 files changed, 69 insertions(+), 17 deletions(-)

diff --git a/package/CHANGELOG b/package/CHANGELOG
index 078ca6651b2..f27374573dd 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -23,12 +23,24 @@ The rules for this file:
  * 2.11.0
 
 Fixes
- * Fixed `capped_distance` and `self_capped_distance` (and hence bond
-   guessing, selections and other users of the `nsgrid` and `pkdtree`
-   methods) missing pairs within the cutoff in triclinic boxes where both
-   b_x and c_y are non-zero, e.g. truncated octahedra; `apply_PBC` also
-   wrapped some coordinates in such boxes to just outside the primary
-   image (Issue #4906)
+ * Fixed PBC handling in triclinic boxes in `_triclinic_pbc()` (used by
+   `distances.apply_PBC()`) and `lib.nsgrid.NSGrid` (used by `capped_distance`
+   and `self_capped_distance`):
+   - `apply_PBC()`, `AtomGroup.wrap()`/`pack_into_box()`, `center*(wrap=True)`,
+     `AtomGroup.unwrap()` with `reference="com"`/`"cog"` and the
+     `transformations.wrap` transformation could leave atoms (or compound
+     centers) outside the primary unit cell in boxes with both b_x and c_y
+     non-zero (e.g. truncated octahedra); they were shifted by whole box
+     vectors, so periodic distances were unaffected, but results that depend on
+     absolute wrapped positions (e.g. dipole moments, densities, written
+     coordinates) may change.
+   - `capped_distance()` and `self_capped_distance()` with the `nsgrid` or
+     `pkdtree` methods (and hence bond guessing, distance-based selections)
+     could miss pairs within the cutoff; `nsgrid` was affected in triclinic
+     boxes with two or more non-90º angles (e.g. rhombic dodecahedron),
+     `pkdtree` in boxes with both b_x and c_y non-zero (e.g. truncated
+     octahedron).
+   (Issue #4906, PR #5496)
  * Fixed `HydrogenBondAutoCorrel.run()` returning NaN when a run starts
    without hydrogen bonds. Such runs are now skipped with a RuntimeWarning,
    and a ValueError is raised if all runs are skipped (Issue #5483, PR #5492)
diff --git a/package/MDAnalysis/lib/include/calc_distances.h b/package/MDAnalysis/lib/include/calc_distances.h
index 1b9a6b339e6..41677854172 100644
--- a/package/MDAnalysis/lib/include/calc_distances.h
+++ b/package/MDAnalysis/lib/include/calc_distances.h
@@ -218,14 +218,14 @@ static void _triclinic_pbc(coordinate* coords, uint64_t numcoords, float* box)
     }
 
     // constants for multi-box shifts:
-    const double bi0 = 1.0 / (double) box[0];
-    const double bi4 = 1.0 / (double) box[4];
-    const double bi8 = 1.0 / (double) box[8];
-    const double bi3 = -box[3] * bi0 * bi4;
-    const double bi6 = (-bi3 * box[7] - box[6] * bi0) * bi8;
-    const double bi7 = -box[7] * bi4 * bi8;
+    const double bi0 = 1.0 / (double) box[0];                  // 1/a_x
+    const double bi4 = 1.0 / (double) box[4];                  // 1/b_y
+    const double bi8 = 1.0 / (double) box[8];                  // 1/c_z
+    const double bi3 = -box[3] * bi0 * bi4;                    // -b_x / (a_x * b_y)
+    const double bi6 = (-bi3 * box[7] - box[6] * bi0) * bi8;   // (-b_x * c_y / (a_x * c_z) - c_x / a_x) / c_z
+    const double bi7 = -box[7] * bi4 * bi8;                    // -c_y / (b_y * c_z)
     // variables and constants for single box shifts:
-    const double b_ax_zfactor = (double) box[7] * bi8;
+    const double b_ax_zfactor = (double) box[7] * bi8;         // c_y / c_z
 
 
    /*
@@ -245,7 +245,7 @@ static void _triclinic_pbc(coordinate* coords, uint64_t numcoords, float* box)
         int msr = 0;
         float crd[3];
         double lbound, ubound;
-        
+
         crd[0] = coords[i][0];
         crd[1] = coords[i][1];
         crd[2] = coords[i][2];
diff --git a/package/MDAnalysis/lib/nsgrid.pyx b/package/MDAnalysis/lib/nsgrid.pyx
index afd3c0e06f0..c538095f936 100644
--- a/package/MDAnalysis/lib/nsgrid.pyx
+++ b/package/MDAnalysis/lib/nsgrid.pyx
@@ -196,7 +196,7 @@ cdef class FastNS(object):
     .. versionchanged:: 2.11.0
        Fixed missed neighbours in triclinic boxes where both the *x*
        component of the second box vector and the *y* component of the third
-       are non-zero (Issue #4906).
+       are non-zero or more than one angle is non-90º (Issue #4906).
     """
     cdef readonly double cutoff
     cdef float[:, ::1] coords_bbox
@@ -290,6 +290,8 @@ cdef class FastNS(object):
 
     cdef float _prepare_box(self, box, bint pbc):
         """
+        Set up the cell grid and calculate the maximum allowable cutoff.
+
         Parameters
         ----------
         box : numpy ndarray shape=(6,)
@@ -301,6 +303,11 @@ cdef class FastNS(object):
         -------
         max_cutoff : float
            the maximum allowable cutoff given the box shape and size
+
+
+        .. versionchanged:: 2.11.0
+           Fixed calculation of the cell size for triclinic boxes where
+           more than one angle is non-90º (Issue #4906).
         """
         cdef float cutoff, min_cellsize, max_cutoff
         cdef double ax, bx, by, cx, cy, cz, volume
@@ -441,7 +448,39 @@ cdef class FastNS(object):
         return xyz[0] + xyz[1] * self.cell_offsets[1] + xyz[2] * self.cell_offsets[2]
 
     cdef void coord2cellxyz(self, const float* coord, int* xyz) nogil:
-        """Calculate cell coordinate for coord"""
+        r"""Calculate cell coordinate for coord.
+
+        Parameters
+        ----------
+        coord : numpy.ndarray float of shape (3,)
+            Coordinate to calculate the cell coordinate for
+        xyz : numpy.ndarray int of shape (3,)
+            Cell coordinate for the given coordinate
+
+        Returns
+        -------
+        xyz : numpy.ndarray int of shape (3,)
+            Cell coordinate for the given coordinate
+
+        Background
+        ----------
+        The fractional coordinates :math:`\mathbf{s} = (s_a, s_b, s_c)` for a
+        point :math:`\mathbf{r} = (x, y, z)` in a general triclinic box with
+        box vectors :math:`\mathbf{a} = (a_x, 0, 0)`, :math:`\mathbf{b} = (b_x, b_y, 0)`,
+        and :math:`\mathbf{c} = (c_x, c_y, c_z)` are given by:
+
+        .. math::
+           \begin{align*}
+           s_c &= \frac{z}{c_z}\\
+           s_b &= \frac{y - s_c c_y}{b_y}\\
+           s_a &= \frac{x - s_b b_x - s_c c_x}{a_x}
+           \end{align*}
+
+
+        .. versionchanged:: 2.11.0
+           Fixed calculation of cell coordinate for triclinic boxes with both
+           :math:`b_x` and :math:`c_y` non-zero (Issue #4906).
+        """
         # This assumes coordinate is inside the primary unit cell
         xyz[2] = <int> (coord[2] / self.cellsize[ZZ])
         # y with the c-vector contribution removed, needed for both y and x
diff --git a/testsuite/MDAnalysisTests/analysis/test_dielectric.py b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
index ab7735a9159..5461e8ba1e1 100644
--- a/testsuite/MDAnalysisTests/analysis/test_dielectric.py
+++ b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
@@ -41,7 +41,8 @@ def test_broken_molecules(self, ag):
             ag.wrap()
 
         eps = DielectricConstant(ag, make_whole=False).run()
-        # reference values changed after the apply_PBC fix for Issue #4906
+        # reference values depend on correct PBC handling
+        # (reference values changed after the apply_PBC fix for Issue #4906)
         assert_allclose(eps.results["eps_mean"], 630.326, rtol=1e-03)
 
     def test_broken_repaired_molecules(self, ag):

From 955d89c0abcf6e355a8b9f6409ec015bfded2e66 Mon Sep 17 00:00:00 2001
From: Oliver Beckstein <orbeckst@gmail.com>
Date: Mon, 19 Oct 2026 01:59:14 -0700
Subject: [PATCH 5/7] add rhombic dodecahedron to tests (and made flake8 happy)

---
 .../MDAnalysisTests/lib/test_distances.py     | 63 ++++++++++---------
 testsuite/MDAnalysisTests/lib/test_nsgrid.py  | 14 +++--
 2 files changed, 44 insertions(+), 33 deletions(-)

diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index bf95e1587f9..015bfbf3cda 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -94,6 +94,8 @@ def test_transform_StoR_pass(coord_dtype):
         dtype=np.float32,
     ),
     np.array([26.24, 27.86, 34.66, 106.56, 77.20, 79.13], dtype=np.float32),
+    # rhombic dodecahedron
+    np.array([50.0, 50.0, 50.0, 60.0, 60.0, 90.0], dtype=np.float32),
 )
 
 
@@ -943,8 +945,10 @@ def S_mol_single(TRIC):
     @pytest.mark.parametrize("backend", ["serial", "openmp"])
     @pytest.mark.parametrize("S_mol", [S_mol, S_mol_single], indirect=True)
     def test_transforms(self, S_mol, tri_vec_box, box, backend):
-        # To check the cython coordinate transform, the same operation is done in numpy
-        # Is a matrix multiplication of Coords x tri_vec_box = NewCoords, so can use np.dot
+        # To check the cython coordinate transform, the same operation is done
+        # in numpy.
+        # Is a matrix multiplication of Coords x tri_vec_box = NewCoords, so
+        # can use np.dot
         S_mol1, S_mol2 = S_mol
         # Test transformation
         R_mol1 = distances.transform_StoR(S_mol1, box, backend=backend)
@@ -1168,8 +1172,8 @@ def test_HAS_DISTOPIA_distopia_too_old():
 
 class TestCythonFunctions(object):
     # Unit tests for calc_bonds calc_angles and calc_dihedrals in lib.distances
-    # Tests both numerical results as well as input types as Cython will silently
-    # produce nonsensical results if given wrong data types otherwise.
+    # Tests both numerical results as well as input types as Cython will
+    # silently produce nonsensical results if given wrong data types otherwise.
     prec = 5
 
     @staticmethod
@@ -1394,8 +1398,6 @@ def test_angles(self, backend, dtype, pos, request):
         assert_equal(
             len(angles), 4, err_msg="calc_angles results have wrong length"
         )
-        #        assert_almost_equal(angles[0], 0.0, self.prec,
-        #                           err_msg="Zero length angle calculation failed") # What should this be?
         assert_almost_equal(
             angles[1],
             np.pi,
@@ -1639,7 +1641,8 @@ def manual_dihedral(a, b, c, d):
     @pytest.mark.parametrize("backend", distopia_conditional_backend())
     def test_numpy_compliance_bonds(self, positions, backend):
         a, b, c, d = positions
-        # Checks that the cython functions give identical results to the numpy versions
+        # Checks that the cython functions give identical results to the numpy
+        # versions
         bonds = distances.calc_bonds(a, b, backend=backend)
         bonds_numpy = np.array([mdamath.norm(y - x) for x, y in zip(a, b)])
 
@@ -1653,7 +1656,8 @@ def test_numpy_compliance_bonds(self, positions, backend):
     @pytest.mark.parametrize("backend", distopia_conditional_backend())
     def test_numpy_compliance_angles(self, positions, backend):
         a, b, c, d = positions
-        # Checks that the cython functions give identical results to the numpy versions
+        # Checks that the cython functions give identical results to the numpy
+        # versions
         angles = distances.calc_angles(a, b, c, backend=backend)
         vec1 = a - b
         vec2 = c - b
@@ -1671,7 +1675,8 @@ def test_numpy_compliance_angles(self, positions, backend):
     @pytest.mark.parametrize("backend", distopia_conditional_backend())
     def test_numpy_compliance_dihedrals(self, positions, backend):
         a, b, c, d = positions
-        # Checks that the cython functions give identical results to the numpy versions
+        # Checks that the cython functions give identical results to the numpy
+        # versions
         dihedrals = distances.calc_dihedrals(a, b, c, d, backend=backend)
         ab = a - b
         bc = b - c
@@ -1853,8 +1858,8 @@ def test_coords_in_central_image_skewed_tric(self, backend, box):
 
 @pytest.mark.parametrize("backend", ["serial", "openmp"])
 class TestPeriodicAngles(object):
-    """Test case for properly considering minimum image convention when calculating angles and dihedrals
-    (Issue 172)
+    """Test case for properly considering minimum image convention when
+    calculating angles and dihedrals(Issue 172)
     """
 
     @staticmethod
@@ -1870,7 +1875,8 @@ def positions():
     prec = 5
 
     def test_angles(self, positions, backend):
-        # Shift atom coordinates a few box lengths in random directions and see if we still get same results
+        # Shift atom coordinates a few box lengths in random directions and
+        # see if we still get same results
         a, b, c, d, box = positions
         a2 = a + box * (-1, 0, 0)
         b2 = b + box * (1, 0, 1)
@@ -2320,7 +2326,7 @@ class TestOutputTypes(object):
     @pytest.mark.parametrize("backend", distopia_conditional_backend())
     def test_output_type_distance_array(self, incoords, box, backend):
         res = distances.distance_array(*incoords, box=box, backend=backend)
-        assert type(res) == np.ndarray
+        assert isinstance(res, np.ndarray)
         assert res.shape == (
             incoords[0].shape[0] % 2,
             incoords[1].shape[0] % 2,
@@ -2332,7 +2338,7 @@ def test_output_type_distance_array(self, incoords, box, backend):
     @pytest.mark.parametrize("backend", distopia_conditional_backend())
     def test_output_type_self_distance_array(self, incoords, box, backend):
         res = distances.self_distance_array(incoords, box=box, backend=backend)
-        assert type(res) == np.ndarray
+        assert isinstance(res, np.ndarray)
         assert res.shape == (0,)
         assert res.dtype.type == np.float64
 
@@ -2358,12 +2364,12 @@ def test_output_type_capped_distance(
             pairs, dist = res
         else:
             pairs = res
-        assert type(pairs) == np.ndarray
+        assert isinstance(pairs, np.ndarray)
         assert pairs.dtype.type == np.intp
         assert pairs.ndim == 2
         assert pairs.shape[1] == 2
         if ret_dist:
-            assert type(dist) == np.ndarray
+            assert isinstance(dist, np.ndarray)
             assert dist.dtype.type == np.float64
             assert dist.shape == (pairs.shape[0],)
 
@@ -2389,12 +2395,12 @@ def test_output_type_self_capped_distance(
             pairs, dist = res
         else:
             pairs = res
-        assert type(pairs) == np.ndarray
+        assert isinstance(pairs, np.ndarray)
         assert pairs.dtype.type == np.intp
         assert pairs.ndim == 2
         assert pairs.shape[1] == 2
         if ret_dist:
-            assert type(dist) == np.ndarray
+            assert isinstance(dist, np.ndarray)
             assert dist.dtype.type == np.float64
             assert dist.shape == (pairs.shape[0],)
 
@@ -2403,7 +2409,7 @@ def test_output_type_self_capped_distance(
     @pytest.mark.parametrize("backend", ["serial", "openmp"])
     def test_output_dtype_transform_RtoS(self, incoords, box, backend):
         res = distances.transform_RtoS(incoords, box, backend=backend)
-        assert type(res) == np.ndarray
+        assert isinstance(res, np.ndarray)
         assert res.dtype.type == np.float32
         assert res.shape == incoords.shape
 
@@ -2412,7 +2418,7 @@ def test_output_dtype_transform_RtoS(self, incoords, box, backend):
     @pytest.mark.parametrize("backend", ["serial", "openmp"])
     def test_output_dtype_transform_RtoS(self, incoords, box, backend):
         res = distances.transform_RtoS(incoords, box, backend=backend)
-        assert type(res) == np.ndarray
+        assert isinstance(res, np.ndarray)
         assert res.dtype.type == np.float32
         assert res.shape == incoords.shape
 
@@ -2425,9 +2431,9 @@ def test_output_type_calc_bonds(self, incoords, box, backend):
         res = distances.calc_bonds(*incoords, box=box, backend=backend)
         maxdim = max([crd.ndim for crd in incoords])
         if maxdim == 1:
-            assert type(res) == np.float64
+            assert isinstance(res, np.float64)
         else:
-            assert type(res) == np.ndarray
+            assert isinstance(res, np.ndarray)
             assert res.dtype.type == np.float64
             coord = [crd for crd in incoords if crd.ndim == maxdim][0]
             assert res.shape == (coord.shape[0],)
@@ -2441,9 +2447,9 @@ def test_output_type_calc_angles(self, incoords, box, backend):
         res = distances.calc_angles(*incoords, box=box, backend=backend)
         maxdim = max([crd.ndim for crd in incoords])
         if maxdim == 1:
-            assert type(res) == np.float64
+            assert isinstance(res, np.float64)
         else:
-            assert type(res) == np.ndarray
+            assert isinstance(res, np.ndarray)
             assert res.dtype.type == np.float64
             coord = [crd for crd in incoords if crd.ndim == maxdim][0]
             assert res.shape == (coord.shape[0],)
@@ -2457,9 +2463,9 @@ def test_output_type_calc_dihedrals(self, incoords, box, backend):
         res = distances.calc_dihedrals(*incoords, box=box, backend=backend)
         maxdim = max([crd.ndim for crd in incoords])
         if maxdim == 1:
-            assert type(res) == np.float64
+            assert isinstance(res, np.float64)
         else:
-            assert type(res) == np.ndarray
+            assert isinstance(res, np.ndarray)
             assert res.dtype.type == np.float64
             coord = [crd for crd in incoords if crd.ndim == maxdim][0]
             assert res.shape == (coord.shape[0],)
@@ -2469,7 +2475,7 @@ def test_output_type_calc_dihedrals(self, incoords, box, backend):
     @pytest.mark.parametrize("backend", ["serial", "openmp"])
     def test_output_type_apply_PBC(self, incoords, box, backend):
         res = distances.apply_PBC(incoords, box, backend=backend)
-        assert type(res) == np.ndarray
+        assert isinstance(res, np.ndarray)
         assert res.dtype.type == np.float32
         assert res.shape == incoords.shape
 
@@ -2533,7 +2539,8 @@ def test_used_openmpflag():
 @pytest.mark.parametrize("dtype", (np.float32, np.float64))
 def test_minimize_vectors(box, shift, dtype):
     # test vectors pointing in all directions
-    # these currently all obey minimum convention as they're much smaller than the box
+    # these currently all obey minimum convention as they're much smaller
+    # than the box
     vec = np.array(
         list(itertools.product(range(-1, 2), range(-1, 2), range(-1, 2))),
         dtype=dtype,
diff --git a/testsuite/MDAnalysisTests/lib/test_nsgrid.py b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
index 7c46029dd59..7addc2a1855 100644
--- a/testsuite/MDAnalysisTests/lib/test_nsgrid.py
+++ b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
@@ -94,7 +94,8 @@ def test_nsgrid_badcutoff(universe, cutoff, match):
 
 
 def test_ns_grid_noneighbor(universe):
-    """Check that grid search returns empty lists/arrays when there is no neighbors"""
+    """Check that grid search returns empty lists/arrays when there is no
+    neighbors"""
     ref_id = 0
     cutoff = 0.5
 
@@ -367,7 +368,8 @@ def test_nsgrid_probe_close_to_box_boundary():
 
 def test_zero_max_dist():
     # see issue #2656
-    # searching with max_dist = 0.0 shouldn't cause segfault (and infinite subboxes)
+    # searching with max_dist = 0.0 shouldn't cause segfault (and infinite
+    # subboxes)
     ref = np.array([1.0, 1.0, 1.0], dtype=np.float32)
     conf = np.array([2.0, 1.0, 1.0], dtype=np.float32)
 
@@ -512,12 +514,14 @@ def test_issue_2670():
         [97.16, 97.16, 97.16, 109.47, 109.47, 109.47],
         # LAMMPS cell from Issue #4906
         [15.173525, 14.923327, 14.475218, 109.888, 109.3786, 115.684875],
+        # rhombic dodecahedron
+        [50.0, 50.0, 50.0, 60.0, 60.0, 90.0],
     ],
 )
 def test_issue_4906(box):
-    # In boxes with both b_x and c_y non-zero, atoms were binned into the
-    # wrong cells and cells could be thinner than the cutoff, so neighbours
-    # within the cutoff were missed.
+    # In boxes with both b_x and c_y non-zero, or more than one non-90º angle,
+    # atoms were binned into the wrong cells and cells could be thinner than
+    # the cutoff, so neighbours were missed.
     box = np.array(box, dtype=np.float32)
     rng = np.random.default_rng(4906)
     frac = rng.uniform(0, 1, size=(1000, 3))

From 3ad297eda87832668f7bef1d4b7d58e5274ec905 Mon Sep 17 00:00:00 2001
From: Brady Johnston <brady.johnston@me.com>
Date: Mon, 19 Oct 2026 10:33:25 +0100
Subject: [PATCH 6/7] update test to properly fail on develop

---
 testsuite/MDAnalysisTests/lib/test_nsgrid.py | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/testsuite/MDAnalysisTests/lib/test_nsgrid.py b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
index 7addc2a1855..8806e474a20 100644
--- a/testsuite/MDAnalysisTests/lib/test_nsgrid.py
+++ b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
@@ -508,17 +508,21 @@ def test_issue_2670():
 
 
 @pytest.mark.parametrize(
-    "box",
+    "box, cutoff_frac",
     [
         # truncated octahedron
-        [97.16, 97.16, 97.16, 109.47, 109.47, 109.47],
+        ([97.16, 97.16, 97.16, 109.47, 109.47, 109.47], 0.1),
         # LAMMPS cell from Issue #4906
-        [15.173525, 14.923327, 14.475218, 109.888, 109.3786, 115.684875],
-        # rhombic dodecahedron
-        [50.0, 50.0, 50.0, 60.0, 60.0, 90.0],
+        (
+            [15.173525, 14.923327, 14.475218, 109.888, 109.3786, 115.684875],
+            0.1,
+        ),
+        # rhombic dodecahedron: the old cutoff / sin(angle) cell sizing only
+        # gave cells thinner than the cutoff for some cutoffs, e.g. 10.8
+        ([50.0, 50.0, 50.0, 60.0, 60.0, 90.0], 0.216),
     ],
 )
-def test_issue_4906(box):
+def test_issue_4906(box, cutoff_frac):
     # In boxes with both b_x and c_y non-zero, or more than one non-90º angle,
     # atoms were binned into the wrong cells and cells could be thinner than
     # the cutoff, so neighbours were missed.
@@ -526,7 +530,7 @@ def test_issue_4906(box):
     rng = np.random.default_rng(4906)
     frac = rng.uniform(0, 1, size=(1000, 3))
     coords = (frac @ mda.lib.mdamath.triclinic_vectors(box)).astype(np.float32)
-    cutoff = 0.1 * box[:3].min()
+    cutoff = cutoff_frac * box[:3].min()
 
     dists = mda.lib.distances.distance_array(coords, coords, box=box)
     # skip pairs so close to the cutoff that rounding can decide them

From 5635e459cd6f7f9c112738a4e7ea40bbf8872cf4 Mon Sep 17 00:00:00 2001
From: Oliver Beckstein <orbeckst@gmail.com>
Date: Mon, 19 Oct 2026 02:56:59 -0700
Subject: [PATCH 7/7] updated tests for too thin cells

---
 testsuite/MDAnalysisTests/lib/test_nsgrid.py | 55 ++++++++++++++++----
 1 file changed, 46 insertions(+), 9 deletions(-)

diff --git a/testsuite/MDAnalysisTests/lib/test_nsgrid.py b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
index 7addc2a1855..6a5c08d6ada 100644
--- a/testsuite/MDAnalysisTests/lib/test_nsgrid.py
+++ b/testsuite/MDAnalysisTests/lib/test_nsgrid.py
@@ -508,17 +508,32 @@ def test_issue_2670():
 
 
 @pytest.mark.parametrize(
-    "box",
+    "box, cutoff",
     [
-        # truncated octahedron
-        [97.16, 97.16, 97.16, 109.47, 109.47, 109.47],
-        # LAMMPS cell from Issue #4906
-        [15.173525, 14.923327, 14.475218, 109.888, 109.3786, 115.684875],
-        # rhombic dodecahedron
-        [50.0, 50.0, 50.0, 60.0, 60.0, 90.0],
+        # Each cutoff lies where the old cell sizing made cells thinner than
+        # the cutoff (perpendicular thickness given for the old code).
+        pytest.param(
+            [97.16, 97.16, 97.16, 109.47, 109.47, 109.47],
+            10.0,
+            id="truncated-octahedron",
+        ),
+        pytest.param(
+            [15.173525, 14.923327, 14.475218, 109.888, 109.3786, 115.684875],
+            1.6,
+            id="lammps-issue-4906",
+        ),
+        pytest.param(
+            # with random coordinates the missed pairs are too rare to
+            # fail reliably with the old code (although it does w/ this seed)
+            # see test_issue_4906_cell_thickness for a test that is guaranteed
+            # to fail.
+            [50.0, 50.0, 50.0, 60.0, 60.0, 90.0],
+            6.0,
+            id="rhombic-dodecahedron",
+        ),
     ],
 )
-def test_issue_4906(box):
+def test_issue_4906(box, cutoff):
     # In boxes with both b_x and c_y non-zero, or more than one non-90º angle,
     # atoms were binned into the wrong cells and cells could be thinner than
     # the cutoff, so neighbours were missed.
@@ -526,7 +541,6 @@ def test_issue_4906(box):
     rng = np.random.default_rng(4906)
     frac = rng.uniform(0, 1, size=(1000, 3))
     coords = (frac @ mda.lib.mdamath.triclinic_vectors(box)).astype(np.float32)
-    cutoff = 0.1 * box[:3].min()
 
     dists = mda.lib.distances.distance_array(coords, coords, box=box)
     # skip pairs so close to the cutoff that rounding can decide them
@@ -547,6 +561,29 @@ def found(results):
     assert found(searcher.search(coords)) == set(zip(i, j))
 
 
+def test_issue_4906_cell_thickness():
+    # Check the fix for placing atoms in the wrong cells; b_x = 0, so only the
+    # cell sizing matters. This test engineers pairs just inside the cutoff
+    # along the face normals, with cutoff so chose that the incorrect code
+    # gives cells that are too thin.
+    box = np.array([50.0, 50.0, 50.0, 60.0, 60.0, 90.0], dtype=np.float32)
+    cutoff = 6.18
+    tv = mda.lib.mdamath.triclinic_vectors(box).astype(np.float64)
+    rng = np.random.default_rng(4906)
+    n = 2000
+    p = rng.uniform(0, 1, size=(n, 3)) @ tv
+    # unit normals of the bc and ca faces
+    normals = np.array([np.cross(tv[1], tv[2]), np.cross(tv[2], tv[0])])
+    normals /= np.linalg.norm(normals, axis=1)[:, None]
+    q = p + 0.99 * cutoff * normals[rng.integers(0, 2, n)]
+    coords = np.concatenate([p, q]).astype(np.float32)
+
+    pairs = nsgrid.FastNS(cutoff, coords, box).self_search().get_pairs()
+    found = {tuple(sorted(pair)) for pair in pairs}
+    missing = [k for k in range(n) if (k, k + n) not in found]
+    assert not missing
+
+
 def high_mem_tests_enabled():
     """Returns true if ENABLE_HIGH_MEM_UNIT_TESTS is set to true."""
     env = os.getenv("ENABLE_HIGH_MEM_UNIT_TESTS", default="false").lower()