Rapid calculation of the RMSD drift of a simulation.¶

In [1]: import mdtraj as md

Find two files that are distributed with MDTraj for testing purposes – we can us them to make our plot

In [2]: import mdtraj.testing

In [3]: crystal_fn = mdtraj.testing.get_fn('native.pdb')

In [4]: trajectory_fn = mdtraj.testing.get_fn('frame0.xtc')

Load up the trajectories from disk

In [5]: crystal = md.load(crystal_fn)

In [6]: trajectory = md.load(trajectory_fn, top=crystal)  # load the xtc. the crystal structure defines the topology

In [7]: print trajectory
<mdtraj.Trajectory with 501 frames, 22 atoms>

RMSD with exchangeable hydrogen atoms is generally not a good idea so let’s take a look at just the heavy atoms

In [8]: rmsds_to_crystal = md.rmsd(trajectory, crystal, 0)

In [9]: heavy_atoms = [atom.index for atom in crystal.topology.atoms if atom.element.symbol != 'H']

In [10]: heavy_rmds_to_crystal = md.rmsd(trajectory, crystal, 0, atom_indices=heavy_atoms)

Plot the results

In [11]: import matplotlib.pyplot as pp

In [12]: pp.figure();

In [13]: pp.plot(trajectory.time, rmsds_to_crystal, 'r', label='all atom');

In [14]: pp.plot(trajectory.time, heavy_rmds_to_crystal, 'b', label='heavy atom');

In [15]: pp.legend();

In [16]: pp.title('RMSDs to crystal');

In [17]: pp.xlabel('simulation time (ps)');

In [18]: pp.ylabel('RMSD (nm)');