Issue 2656 nsgrid segfault

package/CHANGELOG

@@ -23,6 +23,9 @@ mm/dd/yy richardjgowers, kain88-de, lilyminium, p-j-smith, bdice, joaomcteixeira
   * 0.21.0
 
 Fixes
+  * Fixed select_atoms("around 0.0 ...") selections and capped_distance 
+    causing a segfault
+    (Issue #2656 PR #2665)
   * `PDBWriter` writes unitary `CRYST1` record (cubic box with sides of 1 Å)
     when `u.dimensions` is `None` or `np.zeros(6)` (Issue #2679, PR #2685)
   * Ensures principal_axes() returns axes with the right hand convention (Issue #2637)

package/MDAnalysis/lib/nsgrid.pyx

@@ -535,7 +535,7 @@ cdef class _NSGrid(object):
 
     """
 
-    cdef readonly dreal cutoff  # cutoff
+    cdef readonly dreal used_cutoff  # cutoff used for cell creation
     cdef ns_int size  # total cells
     cdef ns_int ncoords  # number of coordinates
     cdef ns_int[DIM] ncells  # individual cells in every dimension
@@ -570,27 +570,30 @@ cdef class _NSGrid(object):
         cdef ns_int xi, yi, zi
         cdef real bbox_vol
         cdef dreal relative_cutoff_margin
+        cdef dreal original_cutoff
 
         self.ncoords = ncoords
 
-        # Calculate best cutoff
-        self.cutoff = cutoff
         if not force:
+            # Calculate best cutoff, with 0.01A minimum
+            cutoff = max(cutoff, 0.01)
+            original_cutoff = cutoff
             # First, we add a small margin to the cell size so that we can safely
             # use the condition d <= cutoff (instead of d < cutoff) for neighbor
             # search.
             relative_cutoff_margin = 1.0e-8
-            while self.cutoff == cutoff:
-                self.cutoff = cutoff * (1.0 + relative_cutoff_margin)
+            while cutoff == original_cutoff:
+                cutoff = cutoff * (1.0 + relative_cutoff_margin)
                 relative_cutoff_margin *= 10.0
             bbox_vol = box.c_pbcbox.box[XX][XX] * box.c_pbcbox.box[YY][YY] * box.c_pbcbox.box[YY][YY]
-            while bbox_vol/self.cutoff**3 > max_size:
-                self.cutoff *= 1.2
+            while bbox_vol / cutoff**3 > max_size:
+                cutoff *= 1.2
 
         for i in range(DIM):
-            self.ncells[i] = <ns_int> (box.c_pbcbox.box[i][i] / self.cutoff)
+            self.ncells[i] = <ns_int> (box.c_pbcbox.box[i][i] / cutoff)
             self.cellsize[i] = box.c_pbcbox.box[i][i] / self.ncells[i]
         self.size = self.ncells[XX] * self.ncells[YY] * self.ncells[ZZ]
+        self.used_cutoff = cutoff
 
         self.cell_offsets[XX] = 0
         self.cell_offsets[YY] = self.ncells[XX]
@@ -706,6 +709,7 @@ cdef class FastNS(object):
     cdef real[:, ::1] coords
     cdef real[:, ::1] coords_bbox
     cdef readonly dreal cutoff
+
     cdef _NSGrid grid
     cdef ns_int max_gridsize
     cdef bint periodic
@@ -863,8 +867,7 @@ cdef class FastNS(object):
         # Generate another grid to search
         searchcoords = search_coords.astype(np.float32, order='C', copy=False)
         searchcoords_bbox = self.box.fast_put_atoms_in_bbox(searchcoords)
-        searchgrid = _NSGrid(searchcoords_bbox.shape[0], self.grid.cutoff, self.box, self.max_gridsize, force=True)
-        searchgrid.fill_grid(searchcoords_bbox)
+        searchgrid = _NSGrid(searchcoords_bbox.shape[0], self.grid.used_cutoff, self.box, self.max_gridsize, force=True)
 
         size_search = searchcoords.shape[0]
 

testsuite/MDAnalysisTests/lib/test_nsgrid.py

@@ -29,7 +29,7 @@
 import numpy as np
 
 import MDAnalysis as mda
-from MDAnalysisTests.datafiles import GRO, Martini_membrane_gro
+from MDAnalysisTests.datafiles import GRO, Martini_membrane_gro, PDB
 from MDAnalysis.lib import nsgrid
 
 
@@ -231,3 +231,46 @@ def test_nsgrid_probe_close_to_box_boundary():
     expected_dists = np.array([2.3689647], dtype=np.float64)
     assert_equal(results.get_pairs(), expected_pairs)
     assert_allclose(results.get_pair_distances(), expected_dists, rtol=1.e-6)
+
+
+def test_zero_max_dist():
+    # see issue #2656
+    # 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)
+
+    box = np.array([10., 10., 10., 90., 90., 90.], dtype=np.float32)
+
+    res = mda.lib.distances._nsgrid_capped(ref, conf, box=box, max_cutoff=0.0)
+
+
+@pytest.fixture()
+def u_pbc_triclinic():
+    u = mda.Universe(PDB)
+    u.dimensions = [10, 10, 10, 60, 60, 60]
+    return u
+
+
+def test_around_res(u_pbc_triclinic):
+    # sanity check for issue 2656, shouldn't segfault (obviously)
+    ag = u_pbc_triclinic.select_atoms('around 0.0 resid 3')
+    assert len(ag) == 0
+
+
+def test_around_overlapping():
+    # check that around 0.0 catches when atoms *are* superimposed
+    u = mda.Universe.empty(60, trajectory=True)
+    xyz = np.zeros((60, 3))
+    x = np.tile(np.arange(12), (5,))+np.repeat(np.arange(5)*100, 12)
+    # x is 5 images of 12 atoms
+
+    xyz[:, 0] = x  # y and z are 0
+    u.load_new(xyz)
+
+    u.dimensions = [100, 100, 100, 60, 60, 60]
+    # Technically true but not what we're testing:
+    # dist = mda.lib.distances.distance_array(u.atoms[:12].positions,
+    #                                         u.atoms[12:].positions,
+    #                                         box=u.dimensions)
+    # assert np.count_nonzero(np.any(dist <= 0.0, axis=0)) == 48
+    assert u.select_atoms('around 0.0 index 0:11').n_atoms == 48