Plotting a Ramachandran map with matplotlib

In [1]:

from __future__ import print_function
%matplotlib inline
import mdtraj as md

Matplotlib created a temporary config/cache directory at /build/matplotlib-su8qb2pw because the default path (/homeless-shelter/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.

Lets load up the trajectory that we simulated in a previous example

In [2]:

traj = md.load('ala2.h5')
atoms, bonds = traj.topology.to_dataframe()
atoms

Out[2]:

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

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 [3]:

psi_indices, phi_indices = [6, 8, 14, 16], [4, 6, 8, 14]
angles = md.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 [4]:

from pylab import *
from math import pi

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

Out[4]:

(-3.141592653589793, 3.141592653589793)

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