In [ ]:

    

%pylab inline
import mdtraj as md


    

In [ ]:

    

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

import mdtraj.testing
crystal_fn = mdtraj.testing.get_fn('native.pdb')
trajectory_fn = mdtraj.testing.get_fn('frame0.xtc')


    

In [ ]:

    

crystal = md.load(crystal_fn)
trajectory = md.load(trajectory_fn, top=crystal)  # load the xtc. the crystal structure defines the topology
print trajectory


    

In [ ]:

    

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

rmsds_to_crystal = md.rmsd(trajectory, crystal, 0)
heavy_atoms = [atom.index for atom in crystal.topology.atoms if atom.element.symbol != 'H']
heavy_rmds_to_crystal = md.rmsd(trajectory, crystal, 0, atom_indices=heavy_atoms)


    

In [ ]:

    

figure()
plot(trajectory.time, rmsds_to_crystal, 'r', label='all atom')
plot(trajectory.time, heavy_rmds_to_crystal, 'b', label='heavy atom')
legend()
title('RMSDs to crystal')
xlabel('simulation time (ps)')
ylabel('RMSD (nm)')