diff --git a/package/CHANGELOG b/package/CHANGELOG index f3b47084194..8032585ef40 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -17,6 +17,8 @@ mm/dd/yy * 0.19.3 Fixes + * fixed lib.nsgrid segfault for coordinates very close to the box boundary + (Issue #2132, PR #2136) Deprecations diff --git a/package/MDAnalysis/lib/nsgrid.pyx b/package/MDAnalysis/lib/nsgrid.pyx index e2404cb43c5..94cfaea70df 100644 --- a/package/MDAnalysis/lib/nsgrid.pyx +++ b/package/MDAnalysis/lib/nsgrid.pyx @@ -539,7 +539,9 @@ cdef class _NSGrid(object): cdef ns_int ncoords # number of coordinates cdef ns_int[DIM] ncells # individual cells in every dimension cdef ns_int[DIM] cell_offsets # Cell Multipliers - cdef real[DIM] cellsize # cell size in every dimension + # cellsize MUST be double precision, otherwise coord2cellid() may fail for + # coordinates very close to the upper box boundaries! See Issue #2132 + cdef dreal[DIM] cellsize # cell size in every dimension cdef ns_int nbeads_per_cell # maximum beads cdef ns_int *nbeads # size (Number of beads in every cell) cdef ns_int *beadids # size * nbeads_per_cell (Beadids in every cell) @@ -572,14 +574,14 @@ cdef class _NSGrid(object): # Calculate best cutoff self.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) - relative_cutoff_margin *= 10.0 if not force: + # 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) + 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 diff --git a/testsuite/MDAnalysisTests/lib/test_nsgrid.py b/testsuite/MDAnalysisTests/lib/test_nsgrid.py index f745581ea8c..a6dacf0c8f7 100644 --- a/testsuite/MDAnalysisTests/lib/test_nsgrid.py +++ b/testsuite/MDAnalysisTests/lib/test_nsgrid.py @@ -47,7 +47,6 @@ def run_grid_search(u, ref_id, cutoff=3): return searcher.search(searchcoords) - def test_pbc_box(): """Check that PBC box accepts only well-formated boxes""" pbc = True @@ -214,3 +213,20 @@ def test_nsgrid_selfsearch(box, result): searchresults = searcher.self_search() pairs = searchresults.get_pairs() assert_equal(len(pairs)//2, result) + +def test_nsgrid_probe_close_to_box_boundary(): + # FastNS.search used to segfault with this box, cutoff and reference + # coordinate prior to PR #2136, so we ensure that this remains fixed. + # See Issue #2132 for further information. + ref = np.array([[55.783722, 44.190044, -54.16671]], dtype=np.float32) + box = np.array([53.785854, 43.951054, 57.17597, 90., 90., 90.], dtype=np.float32) + cutoff = 3.0 + # search within a configuration where we know the expected outcome: + conf = np.ones((1, 3), dtype=np.float32) + searcher = nsgrid.FastNS(cutoff, conf, box) + results = searcher.search(ref) + # check if results are as expected: + expected_pairs = np.zeros((1, 2), dtype=np.int64) + 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)