mdanalysis
mdanalysis copied to clipboard
Bugs for distance calculation functions
Expected behavior
I found two distance calculation functions cannot generate the same results, which results that one bond shorter than the cutoff cannot be found when I use guess_bond()
Actual behavior
Code to reproduce the behavior
Datafile: struc.txt
from MDAnalysis.lib import distances
tmp_lammps_filename="struc.txt"
u = mda.Universe(tmp_lammps_filename, format="DATA")
a,b=distances._nsgrid_capped_self(u.atoms.positions,max_cutoff=3.4,box=u.dimensions)
# Find the row in 'a' where the pair [145,80] appears
row_indices = np.where((a == [145,80]).all(axis=1))[0]
if len(row_indices) > 0:
print(f"Found pair [145,80] at row {row_indices[0]}")
print(f"Distance: {b[row_indices[0]]} Å")
else:
# Also check for [80,145] since order might be reversed
row_indices = np.where((a == [80,145]).all(axis=1))[0]
if len(row_indices) > 0:
print(f"Found pair [80,145] at row {row_indices[0]}")
print(f"Distance: {b[row_indices[0]]} Å")
else:
print("Pair not found")
u = mda.Universe(tmp_lammps_filename, format="DATA")
a,b=distances._bruteforce_capped_self(u.atoms.positions,max_cutoff=3.4,box=u.dimensions)
# Find the row in 'a' where the pair [145,80] appears
row_indices = np.where((a == [145,80]).all(axis=1))[0]
if len(row_indices) > 0:
print(f"Found pair [145,80] at row {row_indices[0]}")
print(f"Distance: {b[row_indices[0]]} Å")
else:
# Also check for [80,145] since order might be reversed
row_indices = np.where((a == [80,145]).all(axis=1))[0]
if len(row_indices) > 0:
print(f"Found pair [80,145] at row {row_indices[0]}")
print(f"Distance: {b[row_indices[0]]} Å")
else:
print("Pair not found")
Here _nsgrid_capped_self does not work, which prints: Pair not found
Here _bruteforce_capped_self works which prints: Found pair [80,145] at row 829 Distance: 1.3481745221090606 Å ....
## Current version of MDAnalysis ##
- Which version are you using? 2.8.0
- Which version of Python (`python -V`)?: Python 3.12.0
- Which operating system?centos:9
Hi @zbaibg, thanks for reporting this issue!
It appears to be related to how NSGrid calculates images in a triclinic box, similar to #3133. Richard @richardjgowers, would you mind taking a look at this when you get a chance?
Hi @zbaibg, thanks for reporting this issue!
It appears to be related to how NSGrid calculates images in a triclinic box, similar to #3133. Richard @richardjgowers, would you mind taking a look at this when you get a chance?
It seems grid dividing is unsuitable: https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/lib/nsgrid.pyx#L350
The example incidates that the number of grid is too big results in some bonds missing, should fix the bug as soon as possible because many module use nsgrid method :
@yuxuanzhuang can you confirm that we need to label this issue a defect/bug?