Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugs for distance calculation functions #4906

Open
zbaibg opened this issue Feb 6, 2025 · 1 comment
Open

Bugs for distance calculation functions #4906

zbaibg opened this issue Feb 6, 2025 · 1 comment

Comments

@zbaibg
Copy link

zbaibg commented Feb 6, 2025

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
@yuxuanzhuang
Copy link
Contributor

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants