Baker-Hubbard Hydrogen Bond Identification

In [1]:
from __future__ import print_function
%matplotlib inline
import matplotlib.pyplot as plt
import itertools
import mdtraj as md
import mdtraj.testing

Load up some example data. This is a little 28 residue peptide

In [2]:
t = md.load_pdb('')
<mdtraj.Trajectory with 20 frames, 423 atoms, 28 residues, without unitcells>
/d/nyc/mcgibbon-0/garden/CentOS7/prefixes/mdtraj/1.9.4-01c7/lib/python3.7/site-packages/mdtraj/formats/pdb/ UserWarning: Unlikely unit cell vectors detected in PDB file likely resulting from a dummy CRYST1 record. Discarding unit cell vectors.
  'cell vectors.', category=UserWarning)

md.baker_hubbard identifies hydrogen bonds baced on cutoffs for the Donor-H...Acceptor distance and angle. The criterion employed is $\theta > 120$ and $r_\text{H...Acceptor} < 2.5 A$ in at least 10% of the trajectory. The return value is a list of the indices of the atoms (donor, h, acceptor) that satisfy this criteria.

In [3]:
hbonds = md.baker_hubbard(t, periodic=False)
label = lambda hbond : '%s -- %s' % (t.topology.atom(hbond[0]), t.topology.atom(hbond[2]))
for hbond in hbonds:
GLU1-N -- GLU1-OE2
GLU1-N -- GLU1-OE1
GLY6-N -- SER4-O
CYS7-N -- GLY5-O
TYR11-N -- VAL8-O
MET12-N -- LYS20-O
ARG13-NH1 -- TYR11-O
THR14-N -- ARG18-O
ASP16-N -- ASP16-OD1
GLY17-N -- THR14-O
ARG18-N -- THR14-OG1
ARG18-NE -- ASP16-OD2
LYS20-N -- MET12-O
THR22-N -- GLY10-O
THR14-OG1 -- ASP16-OD1
THR28-OG1 -- ILE27-O

Let's compute the actual distances between the donors and acceptors

In [4]:
da_distances = md.compute_distances(t, hbonds[:, [0,2]], periodic=False)

And plot a histogram for a few of them

In [5]:
color = itertools.cycle(['r', 'b', 'gold'])
for i in [2, 3, 4]:
    plt.hist(da_distances[:, i], color=next(color), label=label(hbonds[i]), alpha=0.5)
plt.xlabel('Donor-acceptor distance [nm]')
Text(0.5, 0, 'Donor-acceptor distance [nm]')

(hbonds.ipynb; hbonds_eval.ipynb;