mdanalysis icon indicating copy to clipboard operation
mdanalysis copied to clipboard

Bugs for distance calculation functions

Open zbaibg opened this issue 9 months ago • 3 comments

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

zbaibg avatar Feb 06 '25 02:02 zbaibg

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?

yuxuanzhuang avatar Feb 08 '25 03:02 yuxuanzhuang

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 :

liuyujie714 avatar Mar 24 '25 07:03 liuyujie714

@yuxuanzhuang can you confirm that we need to label this issue a defect/bug?

orbeckst avatar Jun 12 '25 16:06 orbeckst