-
Notifications
You must be signed in to change notification settings - Fork 823
Closed
Description
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
Reactions are currently unavailable