capped_distance array: kernel exploding for a specific frame in this trajectory

[micaela-matta]

Expected behavior

capped_distance should return something for each step of the trajectory

Actual behavior

when iterating over the trajectory, the kernel dies as the 8th frame is being processed - so initially I suspected a memory issue.
however, when running only the 8th frame, the same thing (see below) happens.
I have used other analysis functions on this trajectory including distance_array, so it seems specific to the recent implementation of capped_distance.

Code to reproduce the behavior

import MDAnalysis as mda
from MDAnalysis.lib import distances

u = mda.Universe(TOP, TRJ)
u.trajectory[7]
A=u.atoms.select_atoms("type os oh")
H=u.atoms.select_atoms("type HW")
#this works
distances.distance_array(H.positions, A.positions, box=u.dimensions)
#this doesn't
capped_distance(H.positions, A.positions, max_cutoff=3.0, box=u.dimensions)

Currently version of MDAnalysis

• 0.19.0
• Python 3.6.6 :: Anaconda custom (64-bit)
• Mac OS X Mojave

TOP
TRJ

[richardjgowers]

This is a tricky one... to narrow it down further:

from MDAnalysis.lib.nsgrid import FastNS

fstns = FastNS(3.0, A.positions, u.dimensions)
# this position causes the crash
fstns.search(H.positions[5991][None, :])

Looking at what happens, we generate a probe position of

probe = [52.619850158691406, 41.050682067871094, 57.17596435546875]
# compared to the box dimensions
u.dimensions = [53.785854, 43.951054, 57.17597 , 90.      , 90.      , 90.      ]

So the probe is ever so slightly outside of the box, after they're meant to be packed inside. We then turn this coordinate into a cell id, (ie which mini-box inside the overall simulation box we search).

This then tries to access cell 4759 (out of 4522), and out of bounds hilarity ensues.

ping @zemanj who has been playing with precision things

[zemanj]

I'll do some time travelling tomorrow and see whether I have broken anything. Might be because of the slightly increased cell size I introduced.

[richardjgowers]

@zemanj I did play with reducing everything to single precision etc, I don't think it's your fault (so far :))

Ideally we don't introduce bound checks, I'm just thinking how we can do the coordinate to cell address transformation in a safe way...

[orbeckst]

Being a clueless by-stander, where is

So the probe is ever so slightly outside of the box,

in #2132 (comment) ? It looks to me that the probe coordinates are still inside the u.dimensions.

In [2]: import numpy as np

In [3]: probe = np.array([52.619850158691406, 41.050682067871094, 57.17596435546
   ...: 875])

In [4]: dimensions = np.array([53.785854, 43.951054, 57.17597 , 90.      , 90.
   ...:     , 90.      ])

In [5]: probe < dimensions[:3]
Out[5]: array([ True,  True,  True], dtype=bool)

[zemanj]

I reverted to the last commit before I started fiddling with lib.nsgrid but the problem persists. On the one hand, I'm really glad it wasn't me who broke the code, but, on the other hand, that doesn't solve the issue. I'm a bit short on time today, hope I can squeeze this in tomorrow.

[zemanj]

@orbeckst You're right, but you're wrong: Still it is this very tiny distance from the box boundary that causes the overflow. The reason for the failure is that _NSGrid.coord2cellid() uses single precision floating point division. In this particular case, that leads to
<int> (probe[2] / cellsize[2]) == <int> (dimensions[2] / cellsize[2])
which should never happen. Luckily, the fix for this is literally the smallest fix ever:
cdef real[DIM] cellsize ➡️ cdef dreal[DIM] cellsize 😆

I'll open a PR in a second.