Skip to content

H-bond selection discrepency from version 2.8+ #5010

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

Closed
cara-a-k opened this issue Apr 4, 2025 · 1 comment · Fixed by #5018
Closed

H-bond selection discrepency from version 2.8+ #5010

cara-a-k opened this issue Apr 4, 2025 · 1 comment · Fixed by #5018

Comments

@cara-a-k
Copy link

cara-a-k commented Apr 4, 2025

I tested a script developed with MDAnalysis 2.5 with version 2.9 and noticed it was reporting fewer hydrogen bonds than previously. After a series of tests to narrow down the issue, I found the difference started in version 2.8 and is related to how the donors and acceptors are set.

Expected behavior

hbonds = HydrogenBondAnalysis(
    universe=u,
    between=["resname DOPC", "resname G12C"],
    d_a_cutoff=3.25,
    d_h_a_angle_cutoff=130,
    update_selections=False
)
hbonds.donors_sel="name O*"
hbonds.acceptors_sel="name O*"

should set the donor and acceptor selections in the same way as:

hbonds = HydrogenBondAnalysis(
    universe=u,
    between=["resname DOPC", "resname G12C"],
    d_a_cutoff=3.25,
    d_h_a_angle_cutoff=130,
    donors_sel="name O*",
    acceptors_sel="name O*",
    update_selections=False
)

and therefore detect the same number and type of H-bonds.

Actual behavior

From version 2.8, only the second method works to set the donor and acceptor selections. The first method wasn't working and wasn't reporting any errors. Instead it just reverted to default behaviour. In my case this resulted in some acceptors being ignored. I used one of the test systems to reproduce the problem below. In this case, it fails to exclude N atoms from the analysis.

Code to reproduce the behavior

import MDAnalysis as mda
import pandas as pd
from MDAnalysis.tests.datafiles import TPR_xvf, XTC_sub_sol
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis

u = mda.Universe(TPR_xvf, XTC_sub_sol)

hbonds = HydrogenBondAnalysis(
    universe=u,
    between=["protein", "resname SOL"],
    d_a_cutoff=3.25,
    d_h_a_angle_cutoff=130,
    update_selections=False
)
hbonds.donors_sel="name O*"
hbonds.acceptors_sel="name O*"

hbonds.run(
    start=0,
    stop=1,
    step=None,
    verbose=True
)
donor_atom = []
acceptor_atom = []
df = pd.DataFrame(hbonds.results["hbonds"])
df.columns = ["frame", "donor_index", "donor_h_index", "acceptor_index", "distance", "angle"]
print(f"{len(df)} h-bonds identified")

for i in range(len(df)):
    donor_atom.append(u.select_atoms("index {:.0f}".format(df["donor_index"][i])).names[0])
    acceptor_atom.append(u.select_atoms("index {:.0f}".format(df["acceptor_index"][i])).names[0])

atom_df = pd.DataFrame({"donor_atom":pd.Series(donor_atom), "acceptor_atom":pd.Series(acceptor_atom)})
print(f"Unique donors: {atom_df.donor_atom.unique()}")
print(f"Unique acceptors: {atom_df.acceptor_atom.unique()}")
@yuxuanzhuang
Copy link
Contributor

Thanks for reporting this!

The bug was introduced in 474be5b#diff-fe4f2ba60589629c5242f1bf6f65809b6fef8f89b82452a35afa7b2d581cae76 when we moved self._donors, self._hydrogens = self._get_dh_pairs() from _prepare to __init__ and thus they were not updated when donors_sel and acceptors_sel were set later on.

I’ll start working on a fix shortly.

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