Strange Hydrogen Bonds from a GROMACS Trajectory.

[hariseldon99]

Hi all. So I'm getting strange hydrogen bonds when I run HydrogenBondAnalysis() from MDAnalysis on a ligand-residue selection from a GROMACS xtc trajectory file.

After the run, I execute the count_by_ids() method and tabulate the names of the donor, hydrogen and acceptor. This yields counts for bonds like N2::H7::O, but also bonds like H7::H7::O, H8::H8::HB1 etc.

What are they? These do not seem to have any electronegative donors! In any case, my code is below:

import MDAnalysis as mda
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA
from tabulate import tabulate

print(mda.__version__)

#Giving full tpr file to include bonding info
u = mda.Universe("prod_complex.tpr", "prod_complex_compact_compact.xtc")

hb = HBA(universe=u,
                donors_sel='resname PHZ',
                acceptors_sel='resnum 93',
                update_selections=False
                )
hb.run(verbose=True)

hbond_data = []
counts = hb.count_by_ids()

for donor, hydrogen, acceptor, count in counts[:10]:
        d, h, a = u.atoms[donor], u.atoms[hydrogen], u.atoms[acceptor]
        hbond_data.append([d.name,h.name,a.name, int(count)])

print(tabulate(hbond_data, headers=["Donor", "Hydrogen", "Acceptor", "Count"]))

Output:

'2.8.0-dev0'

Donor    Hydrogen    Acceptor      Count
-------  ----------  ----------  -------
H7       H7          O               353
H8       H8          HB1             217
H8       H8          O               144
H7       H7          HB1              99
H8       H8          HB2              92
N2       H7          O                90
H8       H8          HD22             64
H7       H7          HD23             63
H6       H6          HD21             49
H8       H8          HD23             45

Update

I ran

gmx check -f prod_complex.xtc -s1 prod_complex.tpr

and got distance mismatches of the order of a few picometres between the xtc file and tpr file.

...

# Atoms  32758
Precision 0.001 (nm)

...
Distance between atoms 338 and 344 is 0.155, should be 0.152
Distance between atoms 344 and 345 is 0.125, should be 0.123
Distance between atoms 344 and 346 is 0.138, should be 0.133
Distance between atoms 346 and 348 is 0.149, should be 0.145
Distance between atoms 348 and 350 is 0.147, should be 0.153
Distance between atoms 348 and 356 is 0.157, should be 0.152
Distance between atoms 350 and 353 is 0.150, should be 0.152
Distance between atoms 353 and 354 is 0.126, should be 0.125
Distance between atoms 353 and 355 is 0.126, should be 0.125
Distance between atoms 356 and 357 is 0.121, should be 0.123
Distance between atoms 356 and 358 is 0.133, should be 0.133
Distance between atoms 358 and 360 is 0.152, should be 0.145
Distance between atoms 360 and 362 is 0.156, should be 0.153
Distance between atoms 360 and 371 is 0.153, should be 0.152
Distance between atoms 362 and 365 is 0.154, should be 0.153
Distance between atoms 365 and 368 is 0.151, should be 0.152
Distance between atoms 368 and 369 is 0.120, should be 0.125
Distance between atoms 368 and 370 is 0.124, should be 0.125
...
...

Is that significant? Could that be affecting the detection of the hydrogen bonds? Seems small...

[hariseldon99]

Update: This seems to work

Apparently, I have to be super precise about my selection rules for donors, hydrogens and acceptors, or else it'll come up with these false positives...

This seems to work okay, but the bond counts were distributed over multiple entries. I artificially combined them for now...

from tabulate import tabulate
import MDAnalysis as mda
u = mda.Universe("prod_complex.tpr", "prod_complex_traj.xtc")
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA
hb = HBA(universe=u, update_selections=False)
hb.donors_sel = hb.guess_donors("resname PHZ")
hb.hydrogens_sel = hb.guess_hydrogens("resname PHZ")
hb.acceptors_sel = hb.guess_acceptors("resnum 93")
hb.run(verbose=True)

hbond_data = []
counts = hb.count_by_ids()

for donor, hydrogen, acceptor, count in counts[:10]:
        d, h, a = u.atoms[donor], u.atoms[hydrogen], u.atoms[acceptor]
        hbond_data.append([d.name,h.name,a.name, int(count)])
        
        
print(tabulate(hbond_data, headers=["Donor", "Hydrogen", "Acceptor", "Count"]))

Donor    Hydrogen    Acceptor      Count
-------  ----------  ----------  -------
N2       H7          O               216
N2       H7          O                90
N2       H8          O                34
N2       H8          O                15

Dunno why the counts for both N2::H7::0 and N2::H8::0 are distributed among two separate entries. Anyhow, I combined them

combined_hbonds = {}

for donor, hydrogen, acceptor, count in hbond_data:
    key = (donor, hydrogen, acceptor)
    if key in combined_hbonds:
        combined_hbonds[key] += int(count)
    else:
        combined_hbonds[key] = int(count)

combined_hbond_data = [[donor, hydrogen, acceptor, count] 
                for (donor, hydrogen, acceptor), count in combined_hbonds.items()]
print(tabulate(combined_hbond_data, headers=["Donor", "Hydrogen", "Acceptor", "Count"]))

Donor    Hydrogen    Acceptor      Count
-------  ----------  ----------  -------
N2       H7          O               306
N2       H8          O                49