Plotting a Ramachandran map with matplotlib

ramachandran-plot_evaluated
In [1]: 
%pylab inline
import mdtraj

# Lets load up the trajectory that we simulated in a previous example
traj = mdtraj.load('ala2.h5')
atoms, bonds = traj.topology.to_dataframe()
print atoms

Populating the interactive namespace from numpy and matplotlib
       name element  resSeq resName  chainID
serial                                      
0        H1       H       0     ACE        0
1       CH3       C       0     ACE        0
2        H2       H       0     ACE        0
3        H3       H       0     ACE        0
4         C       C       0     ACE        0
5         O       O       0     ACE        0
6         N       N       1     ALA        0
7         H       H       1     ALA        0
8        CA       C       1     ALA        0
9        HA       H       1     ALA        0
10       CB       C       1     ALA        0
11      HB1       H       1     ALA        0
12      HB2       H       1     ALA        0
13      HB3       H       1     ALA        0
14        C       C       1     ALA        0
15        O       O       1     ALA        0
16        N       N       2     NME        0
17        H       H       2     NME        0
18        C       C       2     NME        0
19       H1       H       2     NME        0
20       H2       H       2     NME        0
21       H3       H       2     NME        0


/Users/rmcgibbo/miniconda/envs/2.7.9/lib/python2.7/site-packages/mdtraj-0.8.0-py2.7-macosx-10.5-x86_64.egg/mdtraj/formats/hdf5.py:330: UserWarning: No resSeq information found in HDF file, defaulting to zero-based indices
  warnings.warn('No resSeq information found in HDF file, defaulting to zero-based indices')

Because alanine dipeptide is a little nonstandard in the sense that it's basically dominated by the ACE and NME capping residues, we need to find the indicies of the atoms involved in the phi and psi angles somewhat manually. For standard cases, see compute_phi() and compute_psi() for easier solutions that don't require you to manually find the indices of each dihedral angle.

In this case, we're just specifying the four atoms that together parameterize the phi and psi dihedral angles.

In [2]: 
psi_indices, phi_indices = [6, 8, 14, 16], [4, 6, 8, 14]
angles = mdtraj.geometry.compute_dihedrals(traj, [phi_indices, psi_indices])

Lets plot our dihedral angles in a scatter plot using matplotlib. What conformational states of Alanine dipeptide did we sample?

In [3]: 
from math import pi
figure()
title('Dihedral Map: Alanine dipeptide')
scatter(angles[:, 0], angles[:, 1], marker='x', c=traj.time)
cbar = pp.colorbar()
cbar.set_label('Time [ps]')
xlabel(r'$\Phi$ Angle [radians]')
xlim(-pi, pi)
ylabel(r'$\Psi$ Angle [radians]')
ylim(-pi, pi)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-3-4cbc29ccc8f2> in <module>()
      3 title('Dihedral Map: Alanine dipeptide')
      4 scatter(angles[:, 0], angles[:, 1], marker='x', c=traj.time)
----> 5 cbar = pp.colorbar()
      6 cbar.set_label('Time [ps]')
      7 xlabel(r'$\Phi$ Angle [radians]')

NameError: name 'pp' is not defined

(ramachandran-plot.ipynb; ramachandran-plot_evaluated.ipynb; ramachandran-plot.py)

Versions