capped distance fails for very small (<0.1A) boxes

[richardjgowers]

Code to reproduce the behavior

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD,  GRO, PDB, TPR, XTC, TRR,  PRMncdf, NCDF

boxsize = 0.01
box = np.array([boxsize, boxsize, boxsize, 60., 60., 60], dtype=np.float32)

u = mda.Universe(PDB)
u.dimensions = box

mda.lib.nsgrid.FastNS(boxsize / 10., u.atoms.positions, box=u.dimensions)
# this hangs the interpreter then crashes

Doesn't seem to fail when boxsize is 0.1, so not likely to actually impact anyone.

Current version of MDAnalysis

• Which version are you using? (run python -c "import MDAnalysis as mda; print(mda.__version__)")

v0.21 dev

• Which version of Python (python -V)?
• Which operating system?

[richardjgowers]

I've narrowed it down to this access (but not the why):

https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/lib/nsgrid.pyx#L691

[IAlibay]

@richardjgowers is this fixed by your latest nsgrid code?

[richardjgowers]

I’d want to write a test, but I think so...

…

On Wed, Mar 10, 2021 at 17:20, Irfan Alibay ***@***.***> wrote: @richardjgowers <https://github.com/richardjgowers> is this fixed by your latest nsgrid code? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2670 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACGSGB3JFTXILQFTLOEF6NLTC6S6LANCNFSM4MO37S7A> .

[IAlibay]

I'll add it to the list for 1.0.2 then.

[IAlibay]

Just tried this locally with the current master and it seems to be fine.

[IAlibay]

Actually, it's not clear to me that the intended results from #2657 (comment) are actually happening, particularly with regards to:

u = mda.Universe(PDB)
u.dimensions = [1e-3, 1e-3, 1e-3, 60, 60, 60]
ag = u.select_atoms('around 0.0 resid 3')
len(ag)

Here len(ag) returns 0 on master, but the expectation from those comments was 1?

Since I didn't really have much interaction with those PRs/Issues, I'm a little bit lost. Can @lilyminium or @richardjgowers pick this one up?

[richardjgowers]

@IAlibay that should return nothing. An around selection excludes the selection string it's "around"ing. I.e. if you select waters within 4,0 of a lipid, you don't get the lipid (despite it being within 4.0 of itself).

This works:

import MDAnalysis as mda
from MDAnalysisTests.datafiles import PDB

u = mda.Universe(PDB)
u.dimensions = [1e-3] * 3 + [90.] * 3

ag1 = u.select_atoms('resid 2 3')

# this is currently returning nothing, as nothing except resid 3 is within 0.0 of resid 3
ag = ag1.select_atoms('around 0.0 resid 3')

assert len(ag) == 0

# set an overlap from resid 1
u.residues[0].atoms[0].position = u.select_atoms('resid 3')[0].position

ag2 = u.select_atoms('resid 1 3')

# this gives you *only* the atom you've set to overlap with resid 3
ag = ag2.select_atoms('around 0.0 resid 3')

assert len(ag) == 1

As a side note, it's pretty slow because by setting the box so small, everything is overlapping so the search is (inefficiently) dong the full NxN search.

tl;dr bug seems fixed, could add the above as a test to close this

[IAlibay]

Thanks @richardjgowers I'm going to remove the tags though, as this is a master change (although we will need to make an equivalent one for when this code gets put into develop), and I believe we aim for all GSOC contributions to be develop.

[IAlibay]

@richardjgowers this was fixed right? Can I close it?

[richardjgowers]

Yeah there's a corner case in #3183 where the system is unrealistically large, but the issue here is fixed