From 2fe5a9f944a35b2c35bfc816ea8f111ad3d5fb0a Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Thu, 25 Nov 2021 18:12:44 +0000
Subject: [PATCH 01/17] added minimum vectors function

---
 package/MDAnalysis/lib/c_distances.pyx | 113 +++++++++++++++++++++++++
 package/MDAnalysis/lib/distances.py    |  38 +++++++++
 2 files changed, 151 insertions(+)

diff --git a/package/MDAnalysis/lib/c_distances.pyx b/package/MDAnalysis/lib/c_distances.pyx
index a3b285f561c..3b0af6df605 100644
--- a/package/MDAnalysis/lib/c_distances.pyx
+++ b/package/MDAnalysis/lib/c_distances.pyx
@@ -33,6 +33,9 @@ cimport cython
 import numpy
 cimport numpy
 
+from libc.math cimport fabs
+from libc.float cimport FLT_MAX, DBL_MAX
+
 cdef extern from "string.h":
     void* memcpy(void* dst, void* src, int len)
 
@@ -58,6 +61,7 @@ cdef extern from "calc_distances.h":
     void _ortho_pbc(coordinate* coords, int numcoords, float* box)
     void _triclinic_pbc(coordinate* coords, int numcoords, float* box)
     void minimum_image(double* x, float* box, float* inverse_box)
+    void minimum_image_triclinic(float* x, float* box, float* inverse_box)
 
 OPENMP_ENABLED = True if USED_OPENMP else False
 
@@ -267,3 +271,112 @@ def contact_matrix_pbc(coord, sparse_contacts, box, cutoff):
             if dist < cutoff2:
                 sparse_contacts[i, j] = True
                 sparse_contacts[j, i] = True
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box):
+    cdef int i, j
+
+    for i in range(3):
+        if box[i] > 0:
+            j = <int> (fabs(dx[i]) * inverse_box[i])
+            dx[i] -= j * box[i]
+
+
+# Lifted from calc_distances.h
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box):
+    cdef cython.floating dx_min[3]
+    cdef cython.floating dsq, dsq_min, rx
+    cdef cython.floating ry[2]
+    cdef cython.floating rz[3]
+    cdef int j, ix, iy, iz
+
+    if cython.floating is float:
+        dsq_min = FLT_MAX
+    else:
+        dsq_min = DBL_MAX
+
+    dx_min[0] = 0.0
+    dx_min[1] = 0.0
+    dx_min[2] = 0.0
+
+    # first make shift only 1 cell in any direction
+    j = <int> (fabs(dx[0]) * inverse_box[0])
+    dx[0] -= j * box[0]
+    dx[1] -= j * box[1]
+    dx[2] -= j * box[2]
+    j = <int> (fabs(dx[1]) * inverse_box[1])
+    dx[0] -= j * box[3]
+    dx[1] -= j * box[4]
+    dx[2] -= j * box[5]
+    j = <int> (fabs(dx[2]) * inverse_box[2])
+    dx[0] -= j * box[6]
+    dx[1] -= j * box[7]
+    dx[2] -= j * box[8]
+
+    # then check all images to see which combination of 1 cell shifts gives the best shift
+    for ix in range(-1, 2):
+        rx = dx[0] + box[0] * ix
+        for iy in range(-1, 2):
+            ry[0] = rx + box[3] * iy
+            ry[1] = dx[1] + box[4] * iy
+            for iz in range(-1, 2):
+                rz[0] = ry[0] + box[6] * iz
+                rz[1] = ry[1] + box[7] * iz
+                rz[2] = dx[2] + box[8] * iz
+                dsq = rz[0] * rz[0] + rz[1] * rz[1] + rz[2] * rz[2]
+                if (dsq < dsq_min):
+                    dsq_min = dsq
+                    dx_min[0] = rz[0]
+                    dx_min[1] = rz[1]
+                    dx_min[2] = rz[2]
+
+    dx[0] = dx_min[0]
+    dx[1] = dx_min[1]
+    dx[2] = dx_min[2]
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+def _minimise_vectors_ortho(cython.floating[:, :] vectors not None, cython.floating[:] box not None, cython.floating[:, :] output):
+    cdef int i, n
+    cdef cython.floating box_inverse[3]
+    cdef cython.floating[:] box_inverse_view
+
+    box_inverse[0] = 1.0 / box[0]
+    box_inverse[1] = 1.0 / box[1]
+    box_inverse[2] = 1.0 / box[2]
+
+    box_inverse_view = box_inverse
+    
+    n = len(vectors)
+    for i in range(n):
+        output[i, 0] = vectors[i, 0]
+        output[i, 1] = vectors[i, 1]
+        output[i, 2] = vectors[i, 2]
+        _minimum_image_orthogonal(output[i, :], box, box_inverse_view)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def _minimise_vectors_triclinic(cython.floating[:, :] vectors not None, cython.floating[:] box not None, cython.floating[:, :] output):
+    cdef int i, n
+    cdef cython.floating box_inverse[3]
+    cdef cython.floating[:] box_inverse_view
+
+    box_inverse[0] = 1.0 / box[0]
+    box_inverse[1] = 1.0 / box[1]
+    box_inverse[2] = 1.0 / box[2]
+
+    box_inverse_view = box_inverse
+    
+    n = len(vectors)
+    for i in range(n):
+        output[i, 0] = vectors[i, 0]
+        output[i, 1] = vectors[i, 1]
+        output[i, 2] = vectors[i, 2]
+        _minimum_image_triclinic(output[i, :], box, box_inverse_view)
diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index 9b6337e26bd..cdcd14dc3ed 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -73,6 +73,8 @@
 from .mdamath import triclinic_vectors
 from ._augment import augment_coordinates, undo_augment
 from .nsgrid import FastNS
+from .c_distances import _minimise_vectors_ortho, _minimise_vectors_triclinic
+
 
 # hack to select backend with backend=<backend> kwarg. Note that
 # the cython parallel code (prange) in parallel.distances is
@@ -1549,3 +1551,39 @@ def apply_PBC(coords, box, backend="serial"):
         _run("triclinic_pbc", args=(coords, box), backend=backend)
 
     return coords
+
+
+@check_coords('vectors')
+def minimise_vectors(vectors, box):
+    """Apply minimum image convention to an array of vectors
+
+    Parameters
+    ----------
+    vectors : numpy.ndarray
+        Vector array of shape ``(n, 3)``
+    box : numpy.ndarray
+        The unitcell dimensions of the system, which can be orthogonal or
+        triclinic and must be provided in the same format as returned by
+        :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:
+        ``[lx, ly, lz, alpha, beta, gamma]``.
+
+    Returns
+    -------
+    minimised_vectors : numpy.ndarray
+        Same shape as input.  The vectors from the input, but minimised
+        according to the size of the box.
+
+    .. versionadded:: 2.1.0
+    """
+    boxtype, box = check_box(box)
+    output = np.empty_like(vectors)
+
+    # use box which is same precision as input vectors
+    box = box.astype(vectors.dtype)
+
+    if boxtype == 'ortho':
+        _minimise_vectors_ortho(vectors, box, output)
+    else:
+        _minimise_vectors_triclinic(vectors, box, output)
+
+    return output

From 0df30697a53abe1b259799f17b6784aa73c5656c Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Fri, 26 Nov 2021 13:33:00 +0000
Subject: [PATCH 02/17] better docstrings

---
 package/MDAnalysis/lib/distances.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index cdcd14dc3ed..44f4310d2c5 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -1560,7 +1560,8 @@ def minimise_vectors(vectors, box):
     Parameters
     ----------
     vectors : numpy.ndarray
-        Vector array of shape ``(n, 3)``
+        Vector array of shape ``(n, 3)``, either float32 or float64.  These
+        represent many vectors (such as between two particles).
     box : numpy.ndarray
         The unitcell dimensions of the system, which can be orthogonal or
         triclinic and must be provided in the same format as returned by
@@ -1570,7 +1571,7 @@ def minimise_vectors(vectors, box):
     Returns
     -------
     minimised_vectors : numpy.ndarray
-        Same shape as input.  The vectors from the input, but minimised
+        Same shape and dtype as input.  The vectors from the input, but minimised
         according to the size of the box.
 
     .. versionadded:: 2.1.0

From 68c1320b556dfdfe24df11e089aa582a6a311a69 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Sat, 27 Nov 2021 17:40:43 +0000
Subject: [PATCH 03/17] match calc_distances.h implementation of minimum image

---
 package/MDAnalysis/lib/c_distances.pyx | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/package/MDAnalysis/lib/c_distances.pyx b/package/MDAnalysis/lib/c_distances.pyx
index 3b0af6df605..96f6d035e16 100644
--- a/package/MDAnalysis/lib/c_distances.pyx
+++ b/package/MDAnalysis/lib/c_distances.pyx
@@ -33,7 +33,7 @@ cimport cython
 import numpy
 cimport numpy
 
-from libc.math cimport fabs
+from libc.math cimport fabs, round as cround
 from libc.float cimport FLT_MAX, DBL_MAX
 
 cdef extern from "string.h":
@@ -277,11 +277,12 @@ def contact_matrix_pbc(coord, sparse_contacts, box, cutoff):
 @cython.wraparound(False)
 cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box):
     cdef int i, j
+    cdef cython.floating s
 
     for i in range(3):
         if box[i] > 0:
-            j = <int> (fabs(dx[i]) * inverse_box[i])
-            dx[i] -= j * box[i]
+            s = inverse_box[i] * dx[i]
+            dx[i] = box[i] * (s - cround(s))
 
 
 # Lifted from calc_distances.h

From fb80670ca122e45b198cc98dd8905359b82f7c6b Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Sat, 27 Nov 2021 17:55:40 +0000
Subject: [PATCH 04/17] added enforce_dtype option to check_coords

---
 package/MDAnalysis/lib/util.py | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/package/MDAnalysis/lib/util.py b/package/MDAnalysis/lib/util.py
index 4d83b14004c..f94b68b5b99 100644
--- a/package/MDAnalysis/lib/util.py
+++ b/package/MDAnalysis/lib/util.py
@@ -1961,6 +1961,8 @@ def check_coords(*coord_names, **options):
         * **enforce_copy** (:class:`bool`, optional) -- Enforce working on a
           copy of the coordinate arrays. This is useful to ensure that the input
           arrays are left unchanged. Default: ``True``
+        * **enforce_dtype** (:class:`bool`, optional) -- Enforce a conversion
+          to float32.  Default: ``True``
         * **allow_single** (:class:`bool`, optional) -- Allow the input
           coordinate array to be a single coordinate with shape ``(3,)``.
         * **convert_single** (:class:`bool`, optional) -- If ``True``, single
@@ -2016,6 +2018,7 @@ def check_coords(*coord_names, **options):
     .. versionadded:: 0.19.0
     """
     enforce_copy = options.get('enforce_copy', True)
+    enforce_dtype = optional.get('enforce_dtype', True)
     allow_single = options.get('allow_single', True)
     convert_single = options.get('convert_single', True)
     reduce_result_if_single = options.get('reduce_result_if_single', True)
@@ -2061,13 +2064,14 @@ def _check_coords(coords, argname):
                 if (coords.ndim != 2) or (coords.shape[1] != 3):
                     raise ValueError("{}(): {}.shape must be (n, 3), got {}."
                                      "".format(fname, argname, coords.shape))
-            try:
-                coords = coords.astype(
-                    np.float32, order='C', copy=enforce_copy)
-            except ValueError:
-                errmsg = (f"{fname}(): {argname}.dtype must be convertible to "
-                          f"float32, got {coords.dtype}.")
-                raise TypeError(errmsg) from None
+            if enforce_dtype:
+                try:
+                    coords = coords.astype(
+                        np.float32, order='C', copy=enforce_copy)
+                except ValueError:
+                    errmsg = (f"{fname}(): {argname}.dtype must be convertible to "
+                              f"float32, got {coords.dtype}.")
+                    raise TypeError(errmsg) from None
             return coords, is_single
 
         @wraps(func)

From 563942c017c206ce76a688ec6eff8b857b206fd9 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Mon, 29 Nov 2021 17:53:31 +0000
Subject: [PATCH 05/17] fixup

---
 package/MDAnalysis/lib/util.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/package/MDAnalysis/lib/util.py b/package/MDAnalysis/lib/util.py
index f94b68b5b99..c1498513785 100644
--- a/package/MDAnalysis/lib/util.py
+++ b/package/MDAnalysis/lib/util.py
@@ -2018,7 +2018,7 @@ def check_coords(*coord_names, **options):
     .. versionadded:: 0.19.0
     """
     enforce_copy = options.get('enforce_copy', True)
-    enforce_dtype = optional.get('enforce_dtype', True)
+    enforce_dtype = options.get('enforce_dtype', True)
     allow_single = options.get('allow_single', True)
     convert_single = options.get('convert_single', True)
     reduce_result_if_single = options.get('reduce_result_if_single', True)

From ada60779738ec8ee56cfacff235269f94aa76e07 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Mon, 29 Nov 2021 17:54:49 +0000
Subject: [PATCH 06/17] use check_coord on minimise vectors

---
 package/MDAnalysis/lib/distances.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index 44f4310d2c5..d76e4c7f884 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -1553,7 +1553,7 @@ def apply_PBC(coords, box, backend="serial"):
     return coords
 
 
-@check_coords('vectors')
+@check_coords('vectors', enforce_copy=False, enforce_dtype=False)
 def minimise_vectors(vectors, box):
     """Apply minimum image convention to an array of vectors
 
@@ -1585,6 +1585,6 @@ def minimise_vectors(vectors, box):
     if boxtype == 'ortho':
         _minimise_vectors_ortho(vectors, box, output)
     else:
-        _minimise_vectors_triclinic(vectors, box, output)
+        _minimise_vectors_triclinic(vectors, box.ravel(), output)
 
     return output

From 2db0cb59574248ec7303184628b9c619b5e9bf51 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Mon, 29 Nov 2021 18:23:15 +0000
Subject: [PATCH 07/17] fix triclinic minimise_vectors

---
 package/MDAnalysis/lib/c_distances.pyx | 41 +++++++++++++-------------
 1 file changed, 20 insertions(+), 21 deletions(-)

diff --git a/package/MDAnalysis/lib/c_distances.pyx b/package/MDAnalysis/lib/c_distances.pyx
index 96f6d035e16..8ae5307376c 100644
--- a/package/MDAnalysis/lib/c_distances.pyx
+++ b/package/MDAnalysis/lib/c_distances.pyx
@@ -286,38 +286,35 @@ cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floatin
 
 
 # Lifted from calc_distances.h
+# This however assumes that vectors are at most a single box length, so this is modified to fulfill that
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box):
     cdef cython.floating dx_min[3]
-    cdef cython.floating dsq, dsq_min, rx
+    cdef cython.floating s, dsq, dsq_min, rx
     cdef cython.floating ry[2]
     cdef cython.floating rz[3]
-    cdef int j, ix, iy, iz
+    cdef int ix, iy, iz
+
+    # first make shift only 1 cell in any direction
+    s = cround(inverse_box[2] * dx[2])
+    dx[0] -= s * box[6]
+    dx[1] -= s * box[7]
+    dx[2] -= s * box[8]
+    s = cround(inverse_box[1] * dx[1])
+    dx[0] -= s * box[3]
+    dx[1] -= s * box[4]
+    s = cround(inverse_box[0] * dx[0])
+    dx[0] -= s * box[0]
 
     if cython.floating is float:
         dsq_min = FLT_MAX
     else:
         dsq_min = DBL_MAX
-
     dx_min[0] = 0.0
     dx_min[1] = 0.0
     dx_min[2] = 0.0
 
-    # first make shift only 1 cell in any direction
-    j = <int> (fabs(dx[0]) * inverse_box[0])
-    dx[0] -= j * box[0]
-    dx[1] -= j * box[1]
-    dx[2] -= j * box[2]
-    j = <int> (fabs(dx[1]) * inverse_box[1])
-    dx[0] -= j * box[3]
-    dx[1] -= j * box[4]
-    dx[2] -= j * box[5]
-    j = <int> (fabs(dx[2]) * inverse_box[2])
-    dx[0] -= j * box[6]
-    dx[1] -= j * box[7]
-    dx[2] -= j * box[8]
-
     # then check all images to see which combination of 1 cell shifts gives the best shift
     for ix in range(-1, 2):
         rx = dx[0] + box[0] * ix
@@ -343,7 +340,8 @@ cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.cdivision(True)
-def _minimise_vectors_ortho(cython.floating[:, :] vectors not None, cython.floating[:] box not None, cython.floating[:, :] output):
+def _minimise_vectors_ortho(cython.floating[:, :] vectors not None, cython.floating[:] box not None,
+                            cython.floating[:, :] output not None):
     cdef int i, n
     cdef cython.floating box_inverse[3]
     cdef cython.floating[:] box_inverse_view
@@ -364,14 +362,15 @@ def _minimise_vectors_ortho(cython.floating[:, :] vectors not None, cython.float
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def _minimise_vectors_triclinic(cython.floating[:, :] vectors not None, cython.floating[:] box not None, cython.floating[:, :] output):
+def _minimise_vectors_triclinic(cython.floating[:, :] vectors not None, cython.floating[:] box not None,
+                                cython.floating[:, :] output not None):
     cdef int i, n
     cdef cython.floating box_inverse[3]
     cdef cython.floating[:] box_inverse_view
 
     box_inverse[0] = 1.0 / box[0]
-    box_inverse[1] = 1.0 / box[1]
-    box_inverse[2] = 1.0 / box[2]
+    box_inverse[1] = 1.0 / box[4]
+    box_inverse[2] = 1.0 / box[8]
 
     box_inverse_view = box_inverse
     

From 04b75d0c93f789e154e9035efd4fdceb4bde435f Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Mon, 29 Nov 2021 18:43:16 +0000
Subject: [PATCH 08/17] tests for minimise_vectors

---
 .../MDAnalysisTests/lib/test_distances.py     | 25 +++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index 7df0727dffe..e2f791e85e5 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -23,6 +23,7 @@
 import pytest
 import numpy as np
 from numpy.testing import assert_equal, assert_almost_equal
+import itertools
 from itertools import combinations_with_replacement as comb
 
 import MDAnalysis
@@ -1376,3 +1377,27 @@ def test_wront_backend(self, backend_selection_pos):
 
 def test_used_openmpflag():
     assert isinstance(distances.USED_OPENMP, bool)
+
+
+# test both orthognal and triclinic boxes
+@pytest.mark.parametrize('box', (np.eye(3) * 10, np.array([[10, 0, 0], [2, 10, 0], [2, 2, 10]])))
+# try shifts of -2 to +2 in each dimension, and all combinations of shifts
+@pytest.mark.parametrize('shift', itertools.product(range(-2, 3), range(-2, 3), range(-2, 3)))
+@pytest.mark.parametrize('dtype', (np.float32, np.float64))
+def test_minimise_vectors(box, shift, dtype):
+    # test vectors pointing in all directions
+    # 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)
+    box = box.astype(dtype)
+
+    # box is 3x3 representation
+    # multiply by shift, then sum xyz components then add these to the vector
+    # this technically doesn't alter the vector because of periodic boundaries
+    shifted_vec = (vec + (box.T * shift).sum(axis=1)).astype(dtype)
+
+    box2 = mda.lib.mdamath.triclinic_box(*box).astype(dtype)
+
+    res = mda.lib.distances.minimise_vectors(shifted_vec, box2)
+
+    assert_almost_equal(res, vec, decimal=5)
+    assert res.dtype == dtype

From f4d3fa5bf7b06accda7f7f7d10febe390250d5b3 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Tue, 30 Nov 2021 10:28:22 +0000
Subject: [PATCH 09/17] improved docs

---
 package/MDAnalysis/lib/distances.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index d76e4c7f884..13c2af4023d 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -1557,6 +1557,10 @@ def apply_PBC(coords, box, backend="serial"):
 def minimise_vectors(vectors, box):
     """Apply minimum image convention to an array of vectors
 
+    This function is required for calculating the correct vectors between two points.  A naive approach of
+    ``ag1.positions - ag2.positions`` will not provide the minimum vectors between particles, even if
+    all particles are within the primary unit cell (box).
+
     Parameters
     ----------
     vectors : numpy.ndarray

From dfa0be4faad0e1b78eb8cb9f3b491b36e4bfc572 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Tue, 30 Nov 2021 11:13:36 +0000
Subject: [PATCH 10/17] fixup imports

---
 testsuite/MDAnalysisTests/lib/test_distances.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index e2f791e85e5..c8e090ec861 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -1395,9 +1395,9 @@ def test_minimise_vectors(box, shift, dtype):
     # this technically doesn't alter the vector because of periodic boundaries
     shifted_vec = (vec + (box.T * shift).sum(axis=1)).astype(dtype)
 
-    box2 = mda.lib.mdamath.triclinic_box(*box).astype(dtype)
+    box2 = mdamath.triclinic_box(*box).astype(dtype)
 
-    res = mda.lib.distances.minimise_vectors(shifted_vec, box2)
+    res = distances.minimise_vectors(shifted_vec, box2)
 
     assert_almost_equal(res, vec, decimal=5)
     assert res.dtype == dtype

From 18d85759c2e9500f4c1d13caf9ade0ca0c46f768 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Thu, 2 Dec 2021 11:28:23 +0000
Subject: [PATCH 11/17] minimise -> minimize

---
 package/MDAnalysis/lib/c_distances.pyx          |  4 ++--
 package/MDAnalysis/lib/distances.py             | 12 ++++++------
 testsuite/MDAnalysisTests/lib/test_distances.py |  4 ++--
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/package/MDAnalysis/lib/c_distances.pyx b/package/MDAnalysis/lib/c_distances.pyx
index 8ae5307376c..488a0fc19a6 100644
--- a/package/MDAnalysis/lib/c_distances.pyx
+++ b/package/MDAnalysis/lib/c_distances.pyx
@@ -340,7 +340,7 @@ cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.cdivision(True)
-def _minimise_vectors_ortho(cython.floating[:, :] vectors not None, cython.floating[:] box not None,
+def _minimize_vectors_ortho(cython.floating[:, :] vectors not None, cython.floating[:] box not None,
                             cython.floating[:, :] output not None):
     cdef int i, n
     cdef cython.floating box_inverse[3]
@@ -362,7 +362,7 @@ def _minimise_vectors_ortho(cython.floating[:, :] vectors not None, cython.float
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def _minimise_vectors_triclinic(cython.floating[:, :] vectors not None, cython.floating[:] box not None,
+def _minimize_vectors_triclinic(cython.floating[:, :] vectors not None, cython.floating[:] box not None,
                                 cython.floating[:, :] output not None):
     cdef int i, n
     cdef cython.floating box_inverse[3]
diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index 13c2af4023d..b4d795b1ab5 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -73,7 +73,7 @@
 from .mdamath import triclinic_vectors
 from ._augment import augment_coordinates, undo_augment
 from .nsgrid import FastNS
-from .c_distances import _minimise_vectors_ortho, _minimise_vectors_triclinic
+from .c_distances import _minimize_vectors_ortho, _minimize_vectors_triclinic
 
 
 # hack to select backend with backend=<backend> kwarg. Note that
@@ -1554,7 +1554,7 @@ def apply_PBC(coords, box, backend="serial"):
 
 
 @check_coords('vectors', enforce_copy=False, enforce_dtype=False)
-def minimise_vectors(vectors, box):
+def minimize_vectors(vectors, box):
     """Apply minimum image convention to an array of vectors
 
     This function is required for calculating the correct vectors between two points.  A naive approach of
@@ -1574,8 +1574,8 @@ def minimise_vectors(vectors, box):
 
     Returns
     -------
-    minimised_vectors : numpy.ndarray
-        Same shape and dtype as input.  The vectors from the input, but minimised
+    minimized_vectors : numpy.ndarray
+        Same shape and dtype as input.  The vectors from the input, but minimized
         according to the size of the box.
 
     .. versionadded:: 2.1.0
@@ -1587,8 +1587,8 @@ def minimise_vectors(vectors, box):
     box = box.astype(vectors.dtype)
 
     if boxtype == 'ortho':
-        _minimise_vectors_ortho(vectors, box, output)
+        _minimize_vectors_ortho(vectors, box, output)
     else:
-        _minimise_vectors_triclinic(vectors, box.ravel(), output)
+        _minimize_vectors_triclinic(vectors, box.ravel(), output)
 
     return output
diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index c8e090ec861..533ef274d96 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -1384,7 +1384,7 @@ def test_used_openmpflag():
 # try shifts of -2 to +2 in each dimension, and all combinations of shifts
 @pytest.mark.parametrize('shift', itertools.product(range(-2, 3), range(-2, 3), range(-2, 3)))
 @pytest.mark.parametrize('dtype', (np.float32, np.float64))
-def test_minimise_vectors(box, shift, dtype):
+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
     vec = np.array(list(itertools.product(range(-1, 2), range(-1, 2), range(-1, 2))), dtype=dtype)
@@ -1397,7 +1397,7 @@ def test_minimise_vectors(box, shift, dtype):
 
     box2 = mdamath.triclinic_box(*box).astype(dtype)
 
-    res = distances.minimise_vectors(shifted_vec, box2)
+    res = distances.minimize_vectors(shifted_vec, box2)
 
     assert_almost_equal(res, vec, decimal=5)
     assert res.dtype == dtype

From 366e9ac3958721757c3758ec8d820cfea2c8e595 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Thu, 2 Dec 2021 13:28:27 +0000
Subject: [PATCH 12/17] pep8 docstring length fixes

---
 package/MDAnalysis/lib/distances.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index b4d795b1ab5..60f31020118 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -1557,9 +1557,10 @@ def apply_PBC(coords, box, backend="serial"):
 def minimize_vectors(vectors, box):
     """Apply minimum image convention to an array of vectors
 
-    This function is required for calculating the correct vectors between two points.  A naive approach of
-    ``ag1.positions - ag2.positions`` will not provide the minimum vectors between particles, even if
-    all particles are within the primary unit cell (box).
+    This function is required for calculating the correct vectors between two
+    points.  A naive approach of ``ag1.positions - ag2.positions`` will not
+    provide the minimum vectors between particles, even if all particles are
+    within the primary unit cell (box).
 
     Parameters
     ----------
@@ -1575,8 +1576,8 @@ def minimize_vectors(vectors, box):
     Returns
     -------
     minimized_vectors : numpy.ndarray
-        Same shape and dtype as input.  The vectors from the input, but minimized
-        according to the size of the box.
+        Same shape and dtype as input.  The vectors from the input, but
+        minimized according to the size of the box.
 
     .. versionadded:: 2.1.0
     """

From e4b8e3d0594c6033d27bad0ab3825e1d85af547a Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Thu, 2 Dec 2021 13:28:35 +0000
Subject: [PATCH 13/17] changelog changes

---
 package/CHANGELOG | 1 +
 1 file changed, 1 insertion(+)

diff --git a/package/CHANGELOG b/package/CHANGELOG
index 851b401f9f3..31c7229292a 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -29,6 +29,7 @@ Fixes
 Enhancements
   * Add option for custom compiler flags for C/C++ on build and remove march
     option from setup.cfg (Issue #3428, PR #3429).
+  * Added MDAnalysis.lib.distances.minimize_vectors
 
 Changes
   * Dropped python 3.6 support and raised minimum numpy version to 1.18.0

From 8955b2f05192c6dfaa00bd26c3f7da4217c7cdbf Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Thu, 2 Dec 2021 14:58:00 +0000
Subject: [PATCH 14/17] clarified comment

---
 package/MDAnalysis/lib/c_distances.pyx | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/package/MDAnalysis/lib/c_distances.pyx b/package/MDAnalysis/lib/c_distances.pyx
index 488a0fc19a6..0e704df32f6 100644
--- a/package/MDAnalysis/lib/c_distances.pyx
+++ b/package/MDAnalysis/lib/c_distances.pyx
@@ -286,7 +286,8 @@ cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floatin
 
 
 # Lifted from calc_distances.h
-# This however assumes that vectors are at most a single box length, so this is modified to fulfill that
+# The function in calc_distances.h however assumes that vectors are at most a
+# single box length, so this is modified to first fulfill that assumption
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box):

From 1a4cba8fce0eba13aa44e12bbc7908417435cc75 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Sun, 9 Jan 2022 16:39:00 +0000
Subject: [PATCH 15/17] added nogil sections to minimize vectors

---
 package/MDAnalysis/lib/c_distances.pyx | 26 ++++++++++++++------------
 1 file changed, 14 insertions(+), 12 deletions(-)

diff --git a/package/MDAnalysis/lib/c_distances.pyx b/package/MDAnalysis/lib/c_distances.pyx
index c113709311f..2db74ad8f6f 100644
--- a/package/MDAnalysis/lib/c_distances.pyx
+++ b/package/MDAnalysis/lib/c_distances.pyx
@@ -276,7 +276,7 @@ def contact_matrix_pbc(coord, sparse_contacts, box, cutoff):
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box):
+cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box) nogil:
     cdef int i, j
     cdef cython.floating s
 
@@ -291,7 +291,7 @@ cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floatin
 # single box length, so this is modified to first fulfill that assumption
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box):
+cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box) nogil:
     cdef cython.floating dx_min[3]
     cdef cython.floating s, dsq, dsq_min, rx
     cdef cython.floating ry[2]
@@ -355,11 +355,12 @@ def _minimize_vectors_ortho(cython.floating[:, :] vectors not None, cython.float
     box_inverse_view = box_inverse
     
     n = len(vectors)
-    for i in range(n):
-        output[i, 0] = vectors[i, 0]
-        output[i, 1] = vectors[i, 1]
-        output[i, 2] = vectors[i, 2]
-        _minimum_image_orthogonal(output[i, :], box, box_inverse_view)
+    with nogil:
+        for i in range(n):
+            output[i, 0] = vectors[i, 0]
+            output[i, 1] = vectors[i, 1]
+            output[i, 2] = vectors[i, 2]
+            _minimum_image_orthogonal(output[i, :], box, box_inverse_view)
 
 
 @cython.boundscheck(False)
@@ -377,8 +378,9 @@ def _minimize_vectors_triclinic(cython.floating[:, :] vectors not None, cython.f
     box_inverse_view = box_inverse
     
     n = len(vectors)
-    for i in range(n):
-        output[i, 0] = vectors[i, 0]
-        output[i, 1] = vectors[i, 1]
-        output[i, 2] = vectors[i, 2]
-        _minimum_image_triclinic(output[i, :], box, box_inverse_view)
+    with nogil:
+        for i in range(n):
+            output[i, 0] = vectors[i, 0]
+            output[i, 1] = vectors[i, 1]
+            output[i, 2] = vectors[i, 2]
+            _minimum_image_triclinic(output[i, :], box, box_inverse_view)

From 7477de8e05fd3b5c74a16d4e60a3b3a1b2ac76d6 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Sun, 9 Jan 2022 18:29:23 +0000
Subject: [PATCH 16/17] documentation for minimize_vectors

---
 package/MDAnalysis/lib/c_distances.pyx | 42 +++++++++++++++++++++++---
 1 file changed, 37 insertions(+), 5 deletions(-)

diff --git a/package/MDAnalysis/lib/c_distances.pyx b/package/MDAnalysis/lib/c_distances.pyx
index 2db74ad8f6f..4f2fbb31dd1 100644
--- a/package/MDAnalysis/lib/c_distances.pyx
+++ b/package/MDAnalysis/lib/c_distances.pyx
@@ -276,7 +276,22 @@ def contact_matrix_pbc(coord, sparse_contacts, box, cutoff):
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box) nogil:
+cdef inline void _minimum_image_orthogonal(cython.floating[:] dx,
+                                           cython.floating[:] box,
+                                           cython.floating[:] inverse_box) nogil:
+    """Minimize dx to be the shortest vector
+
+    Parameters
+    ----------
+    dx : numpy.array, shape (3,)
+      vector to minimize
+    box : numpy.array, shape (3,)
+      box length in each dimension
+    inverse_box : numpy.array, shape (3,)
+      inverse of box
+
+    Operates in-place on dx!
+    """
     cdef int i, j
     cdef cython.floating s
 
@@ -286,12 +301,28 @@ cdef inline void _minimum_image_orthogonal(cython.floating[:] dx, cython.floatin
             dx[i] = box[i] * (s - cround(s))
 
 
-# Lifted from calc_distances.h
-# The function in calc_distances.h however assumes that vectors are at most a
-# single box length, so this is modified to first fulfill that assumption
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cdef inline void _minimum_image_triclinic(cython.floating[:] dx, cython.floating[:] box, cython.floating[:] inverse_box) nogil:
+cdef inline void _minimum_image_triclinic(cython.floating[:] dx,
+                                          cython.floating[:] box,
+                                          cython.floating[:] inverse_box) nogil:
+    """Minimise dx to be the shortest vector
+
+    Parameters
+    ----------
+    dx : numpy.array, shape (3,)
+      the vector to apply minimum image convention to
+    box : numpy.array, shape (9,)
+      flattened 3x3 representation of the unit cell
+    inverse_box : numpy.array, shape (3,)
+      inverse of the **diagonal** of the 3x3 representation
+
+    This version is near identical to the version in calc_distances.h, with the
+    difference being that this version does not assume that coordinates are at
+    most a single box length apart.
+
+    Operates in-place on dx!
+    """
     cdef cython.floating dx_min[3]
     cdef cython.floating s, dsq, dsq_min, rx
     cdef cython.floating ry[2]
@@ -371,6 +402,7 @@ def _minimize_vectors_triclinic(cython.floating[:, :] vectors not None, cython.f
     cdef cython.floating box_inverse[3]
     cdef cython.floating[:] box_inverse_view
 
+    # this is only inverse of diagonal, used for initial shift to ensure vector is within a single shift away
     box_inverse[0] = 1.0 / box[0]
     box_inverse[1] = 1.0 / box[4]
     box_inverse[2] = 1.0 / box[8]

From eedc243be8b0daead675eaaffe61bbc8478bcb85 Mon Sep 17 00:00:00 2001
From: richard <richard@nextmovesoftware.com>
Date: Tue, 11 Jan 2022 11:02:51 +0000
Subject: [PATCH 17/17] switch tests to use allclose

---
 testsuite/MDAnalysisTests/lib/test_distances.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py
index 533ef274d96..45c9d75ada7 100644
--- a/testsuite/MDAnalysisTests/lib/test_distances.py
+++ b/testsuite/MDAnalysisTests/lib/test_distances.py
@@ -22,7 +22,7 @@
 #
 import pytest
 import numpy as np
-from numpy.testing import assert_equal, assert_almost_equal
+from numpy.testing import assert_equal, assert_almost_equal, assert_allclose
 import itertools
 from itertools import combinations_with_replacement as comb
 
@@ -1399,5 +1399,5 @@ def test_minimize_vectors(box, shift, dtype):
 
     res = distances.minimize_vectors(shifted_vec, box2)
 
-    assert_almost_equal(res, vec, decimal=5)
+    assert_allclose(res, vec, atol=0.00001)
     assert res.dtype == dtype