diff --git a/package/CHANGELOG b/package/CHANGELOG
index 31ada12c873..72b880c8216 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -23,6 +23,24 @@ The rules for this file:
  * 2.11.0
 
 Fixes
+ * 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)
  * `InterRDF_s` with a parallel backend returned wrong counts for one
    pair of groups and failed for more than two (`nested_array_sum`
    summed exactly two pairs) (Issue #5496, PR #5497)
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 ff130d2c049..41677854172 100644
--- a/package/MDAnalysis/lib/include/calc_distances.h
+++ b/package/MDAnalysis/lib/include/calc_distances.h
@@ -218,16 +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 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;
+    const double b_ax_zfactor = (double) box[7] * bi8;         // c_y / c_z
 
 
    /*
@@ -247,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];
@@ -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] * 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] * 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 a4e245e884e..c538095f936 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
@@ -200,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 or more than one angle is non-90º (Issue #4906).
     """
     cdef readonly double cutoff
     cdef float[:, ::1] coords_bbox
@@ -293,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,)
@@ -304,8 +303,15 @@ 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, new_cellsize
+        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
@@ -337,22 +343,31 @@ 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.
+        # 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 += 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
@@ -433,11 +448,45 @@ 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])
-        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/analysis/test_dielectric.py b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
index 45d8e722cdb..5461e8ba1e1 100644
--- a/testsuite/MDAnalysisTests/analysis/test_dielectric.py
+++ b/testsuite/MDAnalysisTests/analysis/test_dielectric.py
@@ -41,16 +41,23 @@ 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 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):
+        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()
-        assert_allclose(eps.results["eps_mean"], 5.088, rtol=1e-03)
+        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):
         eps = DielectricConstant(ag, temperature=100).run()
diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index 06abe0a99c7..015bfbf3cda 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -84,6 +84,33 @@ 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),
+    # rhombic dodecahedron
+    np.array([50.0, 50.0, 50.0, 60.0, 60.0, 90.0], 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)
+
+
+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)
@@ -122,6 +149,34 @@ 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
+        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, 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):
+        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, 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")
+
     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)
@@ -890,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)
@@ -1115,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
@@ -1341,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,
@@ -1586,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)])
 
@@ -1600,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
@@ -1618,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
@@ -1787,11 +1845,21 @@ 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):
-    """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
@@ -1807,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)
@@ -2257,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,
@@ -2269,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
 
@@ -2295,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],)
 
@@ -2326,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],)
 
@@ -2340,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
 
@@ -2349,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
 
@@ -2362,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],)
@@ -2378,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],)
@@ -2394,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],)
@@ -2406,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
 
@@ -2470,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 582e780172b..39cd57df444 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)
 
@@ -505,6 +507,73 @@ def test_issue_2670():
     assert len(ag2.select_atoms("around 0.0 resid 3")) == 1
 
 
+@pytest.mark.parametrize(
+    "box, cutoff_frac",
+    [
+        # truncated octahedron
+        ([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],
+            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, 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.
+    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 = 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
+    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)
+    i, j = np.triu_indices(len(coords), k=1)
+    within = (dists[i, j] <= cutoff) & ~ambiguous[i, j]
+    expected = set(zip(i[within], j[within]))
+    pairs = {tuple(sorted(p)) for p in found(searcher.self_search())}
+    assert pairs == expected
+
+    i, j = np.nonzero((dists <= cutoff) & ~ambiguous)
+    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()