Wrong results from capped_distance() with "nsgrid" search and atom in the center of the box

[a-anik]

Expected behavior

capped_distance() function should always give the correct result regardless of the employed search method and atomic positions.
In the example below the expected result for cutoff=3.2 is:

• cutoff=3.2
  bruteforce: 2497 pairs
  pkdtree: 2497 pairs
  nsgrid: 2497 pairs

Actual behavior

When one of the atoms of the reference group is positioned exactly in the center of the box, capped_distance() sometimes returns erroneous results for "nsgrid" search method.

• cutoff=2.8
  bruteforce: 1115 pairs
  pkdtree: 1115 pairs
  nsgrid: 1115 pairs
• cutoff=3.2
  bruteforce: 2497 pairs
  pkdtree: 2497 pairs
  nsgrid: 2510 pairs

Code to reproduce the behavior

import MDAnalysis as mda
from MDAnalysis.lib.distances import capped_distance
from MDAnalysis.transformations.translate import center_in_box
from MDAnalysis.tests.datafiles import PDB_xvf

u = mda.Universe(PDB_xvf)
ag = u.select_atoms('index 0')
u.trajectory.ts = center_in_box(ag)(u.trajectory.ts)

box = u.dimensions
reference = u.select_atoms('protein')
configuration = u.select_atoms('not protein')

for cutoff in [2.8, 3.2]:
    print(f"* cutoff={cutoff}")
    for method in ['bruteforce', 'pkdtree', 'nsgrid']:
        pairs, distances = capped_distance(
            reference.positions,
            configuration.positions,
            max_cutoff=cutoff,
            box=box,
            method=method,
            return_distances=True,
            )
        print(f"{method}: {len(pairs)} pairs")

Current version of MDAnalysis

• MDAnalysis 1.0.0 (conda-forge)
• Python 3.7.7
• Ubuntu 18.04.5 LTS

[zemanj]

Hi @a-anik,

thank you for your detailed bug report. It looks like this is a duplicate of #2345 (see especially my comment therein).

[zemanj]

Ah, sorry, I posted the wrong bug report. This one sounds rather like a duplicate of #2229 (see my comment for an explanation). Let me know if you think that this is not the case!

[a-anik]

Yes, it looks like a part of that ongoing issue with the grid-based neighbor searches. There was a discussion in #2345 to make it opt-in. I wonder why it made the way into the latest release as a default choice for some systems. Periodic KD-trees would seem to be a better option. I found this bug while using the new hbond_analysis code. There is no way to control search parameters (and PBC behavior for that matter) in that module, so I'm stuck with off-tree patches for now.

[zemanj]

You are right, this shouldn't even be part of MDAnalysis at the moment. Taken seriously, the bugs in NSGrid have the potential to make almost any result unreliable as NSGrid is used by capped_distance, which, in turn, is used all over the place in the entire code base...
@MDAnalysis/coredevs Could you please comment as this seems serious.

[richardjgowers]

Yeah I'm looking at this

[IAlibay]

Looks like we'll want to get 1.0.1 done soon-ish (see #2798), would it make sense to temporarily take the performance hit and just disable nsgrid (or at least make _determine_method return pkdtree isntead of nsgrid)?

[zemanj]

IMHO the question is whether anyone can find the time to re-implement nsgrid. If not, it *really* should be disabled. I think that performance is secondary, and reliability is much more important. After all, MDAnalysis is scientific Software...

…

On September 2, 2020 8:21:36 PM GMT+02:00, Irfan Alibay ***@***.***> wrote: Looks like we'll want to get 1.0.1 done soon-ish (see #2798), would it make sense to temporarily take the performance hit and just disable nsgrid (or at least make `_determine_method` return pkdtree isntead of nsgrid)? -- You are receiving this because you commented. Reply to this email directly or view it on GitHub: #2919 (comment)

[orbeckst]

I have to agree with @zemanj – scientific software has to be correct first and anything else comes second.

[orbeckst]

I raised #2930 — if @richardjgowers is faster then we use the correct solution, otherwise we should disable nsgrid-use for right now, at a minimum for 1.0.1.

[richardjgowers]

Give me 48h, I think there’s a fix

…

On Wed, Sep 2, 2020 at 22:54, Oliver Beckstein ***@***.***> wrote: I raised #2930 <#2930> — if @richardjgowers <https://github.com/richardjgowers> is faster then we use the correct solution, otherwise we should disable nsgrid-use for right now, at a minimum for 1.0.1. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2919 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACGSGB6MYNNPBENJZ7WVGG3SD25ITANCNFSM4QLXDQHA> .

[orbeckst]

For 1.0.1 we will disable nsgrid (see #2930 and PR #3004 ).

[orbeckst]

Give me 48h, I think there’s a fix.

These were long 2 days ;-).

I am glad that this one is done!