Fix: Incorrect capped_distance in some triclinic box situations

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)

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.

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];

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):

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()