Neighbor search with nsgrid.FastNS yields different pairs for swapped coords and search_coords

[ezavod]

Expected behavior

Let's say one has A and B particles. It should make no difference if one constructs FastNS with coordinates of A and searches with coordinates of B or vice versa. The pairs of A-B that are found within a certain cutoff should be invariant against exchanging A around B or B around A.

Actual behavior

In some cases one gets different pairs depending on the order.

Code to reproduce the behavior

This is simple.gro (notice the extra 5 in the last decimal place of the first box vector):

Written by MDAnalysis
   2
    1RA       A    1   0.000   0.000   0.000
    2RB       B    2   5.500   0.000   0.000
   5.745585  5.00000   5.00000

This script searches pairs of A-B within a cutoff of 3. With PBC A and B are well within this distance.

import MDAnalysis
from MDAnalysis.lib.nsgrid import FastNS

u = MDAnalysis.Universe("simple.gro")

max_cutoff = 3
configuration = u.select_atoms("name A").positions
reference = u.select_atoms("name B").positions
box = u.trajectory[0].dimensions

gridsearch = FastNS(max_cutoff, configuration, box=box)
results = gridsearch.search(reference)
pairs = results.get_pairs()
print("Constructed with A; searched with B:", pairs)

gridsearch = FastNS(max_cutoff, reference, box=box)
results = gridsearch.search(configuration)
pairs = results.get_pairs()
print("Constructed with B; searched with A:", pairs)

Output:

Constructed with A; searched with B: [[0 0]]
Constructed with B; searched with A: []

If one now removes the last decimal place of the first box vector every thing works as expected. One finds the pair no matter which order. So the bug only happens under special circumstances.

Current version of MDAnalysis

• MDAnalysis version 0.19.2
• Python version 3.7.2
• openSUSE Tumbleweed 64bit

I also played around with the current develop branch (9428866) but this problem is still persistent. However one has to adjust the box size slightly, probably because of the extra margin in #2136.

Impact
I noticed this bug while using around selections. AroundSelection, SphericalLayerSelection, SphericalZoneSelection and PointSelection in core/selection.pyall use capped_distance. In there a suitable search method is chosen, which is nsgrid for all but small systems. Therefore I suspect that all geometric selections have this problem for typical system sizes.
How often this happens depends of cause on the number of particles near the box boundary. For my systems in ~10% of frames a wrong number of neighbors was found.

Edit
This also occurs without PBC. This is simple_no_pbc.gro:

Written by MDAnalysis
    2
    1RA       A    1   0.000   0.000   2.929
    2RB       B    2   0.000   0.000   2.823
   4.50000   5.50000  10.98375

Using the script from above this seems to work, but if one changes the z coordinates a bit,
i.e. by inserting configuration[0][2] += 1e-6 before the first search,
it does not work anymore. It does not matter if one uses pbc=False explicitly in the constructor. The x, y coordinates of the particles have seemingly no influence, however altering the x,y box vectors changes the behavior. Small changes in any z component also changes the outcome.
This happens in ~1/250 frames for my systems anyway.

[orbeckst]

@richardjgowers @ayushsuhane – opinions?

[orbeckst]

@ezavod did you try the same with the development branch?

[ayushsuhane]

@ezavod Thanks for posting the bug. It seems like it is one of those problematic edge cases. Although I must agree the behavior is strange.

I am not sure at this point but I suspect the problem is in one of the floating point operations/ floating to integer conversion somewhere in nsgrid.pyx. @orbeckst I will try to get to the root of the problem in a couple of days.

Meanwhile, if someone wants to fix the problem, I would suggest checking if the bug resides in capped_distances or nsgrid. One way of finding it could be to pass different arguments (method='bruteforce, 'pkdtree') in capped_distance, and check if it gives the correct solution.

[ezavod]

@orbeckst yes, both errors also occur with latest development branch (9428866).
@ayushsuhane I also tested capped_distance with all three methods. pkdtree and bruteforce work both without errors for my trajectories.
I agree that this might be somewhere in nsgrid.pyx. I don't had the time to fully understand whats going on there, but coord2cellidseems to return different indices.

Edit: clicked wrong button and accidentally closed the issue

[orbeckst]

@ayushsuhane @richardjgowers @jbarnoud @NinadBhat @tylerjereddy do we have sufficiently extensive tests for capped_distances()? Can we use hypothesis to build some smarter tests? This might also help to gain a bit more insights into where the capped_distances() machinery might be failing.

As pointed out by @ezavod we now use capped_distances() extensively and we are happy about the speed-ups but they cannot come at the cost of correctness.

[orbeckst]

Also cc @zemanj

[ezavod]

It might be a good start for testing to compare the results of capped_distances for the various (bruteforce, pkdtree, nsgrid) methods against each other. The implementation between them is quite different so this might catch some major problems. However no substitute for rigorous testing.

[zemanj]

I tracked down the error causing this behavior. It is not directly related to #2345, but it is due to the precision problems I mentioned in a comment there:

The fact that the cell grid is orthorhombic even for triclinic boxes necessitates the awkward probe shifting procedure, which is not only slow, but is also very error prone: Suppose a coordinate lies on a cell boundary in (at least) one direction. Due to round-off errors, probe coordinates shifted by one cell size can lead to the same or the second next cell being selected as a neighboring cell, depending on whether the particle lies on an upper or lower cell boundary. We have worked on reducing the probability of such cases by doing some computations in double precision. However, this just makes the problem less likely and is not a real fix.

Here, we have exactly that problem: In the case of simple.gro, one of the coordinates (0, 0, 0) lies exactly on the boundary of a grid cell (the cell with index 0).
For the given combination of box dimensions and cutoff, FastNS divides the box into a grid with 15 x 13 x 13 cells. Since the second coordinate lies close to the other side of the box in x-direction, its cell is the one with index 14 (the last cell of the first row in x-direction).
The probe-shifting algorithm now takes the first coordinate (the one at the origin) and computes 27 different probe coordinates by adding or subtracting one (or zero) cell sizes in each direction.
For each of the probe coordinates, the corresponding grid cell is determined and searched.

The precision problem arises when the probe is shifted by one cell length in negative x-direction:
The resulting probe coordinate lies outside of the box and is therefore wrapped back into the box by adding one box vector. It just happens that these floating point additions yield a coordinate just below the last cell boundary in x-direction, so that instead of the cell with index 14, it wrongly determines the cell with index 13 as the neighboring cell, which, in reality, is the second nearest neighbor cell. This is a very good example for a bug that won't go away without getting rid of the clumsy probe shifting approach (for further details see the end of my comment in #2345.)

EDIT: I also checked the second example (simple_no_pbc.gro) and the issue is the same. If unmodified, the configuration coordinate lies directly below a grid cell boundary in the z-direction.
When 1.0e-6 is added, the coordinate lies exactly on the boundary of the next cell. Even though in this case, the probe coordinate doesn't need to be shifted back into the box in z-direction, simply subtracting one cell length in the z-direction yields a wrong cell index (one cell layer too low).