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

HydrogenBondAnalysis 'between' keyword returns incorrect donor-acceptor distances #4091

Closed
p-j-smith opened this issue Mar 24, 2023 · 2 comments · Fixed by #4092
Closed

HydrogenBondAnalysis 'between' keyword returns incorrect donor-acceptor distances #4091

p-j-smith opened this issue Mar 24, 2023 · 2 comments · Fixed by #4092
Assignees
Milestone

Comments

@p-j-smith
Copy link
Member

Hi, I saw a question on the discussion forum about using the between keyword of HydrogenBondAnalysis.

After looking into it I've found an issue whereby using the between keyword leads to the incorrect donor-atom distances being returned. This is because when the hydrogen bonds are filtered to remove those formed by atoms not in the between selection, the donors, hydrogens, and acceptors are filtered but not the distances.

The donor-hydrogen-acceptor angles are not affected because these are calculated after the filtering has been done.

Code to reproduce the behavior

import numpy as np

import MDAnalysis as mda
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis
import MDAnalysisData

np.set_printoptions(suppress=True, precision=3)

data = MDAnalysisData.ifabp_water.fetch_ifabp_water()
u = mda.Universe(data.topology, data.trajectory)

hbonds = HydrogenBondAnalysis(
    universe=u,
    hydrogens_sel='index 1',
    acceptors_sel='index 2188',
    d_a_cutoff=3.30,
    d_h_a_angle_cutoff=150,
)
hbonds.run(stop=1)
print(f"{hbonds.results.hbonds=}")

hbonds_between = HydrogenBondAnalysis(
    universe=u,
    hydrogens_sel='resname TIP3 or protein',
    acceptors_sel='resname TIP3 or protein',
    d_a_cutoff=3.30,
    d_h_a_angle_cutoff=150,
    between=["index 0", "index 2188"],  # between uses donor atoms, not hydrogens
)
hbonds_between.run(stop=1)
print(f"{hbonds_between.results.hbonds=}")

Output:

hbonds.results.hbonds=array([[   0.   ,    0.   ,    1.   , 2188.   ,    2.741,  165.581]])
hbonds_between.results.hbonds=array([[   0.   ,    0.   ,    1.   , 2188.   ,    2.981,  165.581]])

The timestep, donor, hydrogen, and acceptor atoms, and the angle are all correct.
However, the donor-atom distance (5th column) is wrong.

I'll try to find some time in the next couple of days to submit a fix.

Current version of MDAnalysis

  • Which version are you using? (run python -c "import MDAnalysis as mda; print(mda.__version__)")
    2.5.0-dev0
  • Which version of Python (python -V)?
    Python 3.10.8
  • Which operating system?
    macOS 12.6.3
@orbeckst
Copy link
Member

Thank you for looking into the email report and verifying it and raising the issue, @p-j-smith !

@p-j-smith
Copy link
Member Author

Thank you for looking into the email report and verifying it and raising the issue, @p-j-smith !

No worries! It's been a while since I've had time to keep up with what's going on with MDAnalysis, but it's nice to come back to it 😄

@orbeckst orbeckst added this to the 2.5.0 milestone Apr 27, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
2 participants