Introduction to MDTraj¶
Start by loading up a trajectory from disk. MDTraj will automatically parse the file extension and use the appropriate loader.
>>> import mdtraj as md >>> t = md.load('trajectory.xtc', top='trajectory.pdb') >>> print t <mdtraj.Trajectory with 100 frames, 22 atoms at 0x109f0a3d0>
To load files that don’t contain topology information, like Gromacs XTC files,
we need to supply something with the
top keyword argument that describes
the topology, for example a PDB file.
If I’m interested in only a subset of the frames of the trajectory, I can slice it
>>> # lets take a look at the first ten frames >>> print t[1:10] <mdtraj.Trajectory with 9 frames, 22 atoms at 0x109b91310> >>> # or maybe the last frame? >>> print t[-1] <mdtraj.Trajectory with 1 frames, 22 atoms at 0x109b97810>
There’s a lot of information in the trajectory object. The most obvious is the
Cartesian coordinates. They’re stored as a numpy array under
xyz. All of
the distances in the Trajectory are stored in nanometers. The time unit
is picoseconds. Angles are stored in degrees (not radians).
>>> print t.xyz.shape (100, 22, 3) >>> print np.mean(t.xyz) 0.89365752249053032 >>> # the simulation time (in picoseconds) of th first 10 frames >>> print t.time[0:10] array([ 0.002, 0.004, 0.006, 0.008, 0.01 , 0.012, 0.014, 0.016, 0.018, 0.02 ], dtype=float32) >>> # or the unitcell lengths in the last frame? (in nanometers of course) >>> t.unitcell_lengths[-1] array([ 2., 2., 2.], dtype=float32)
Saving the trajectory back to disk is easy.
>>> # the hdf5 format stores the topology inside the file for convenience >>> t[::2].save('halftraj.h5') >>> # the format will be parsed based on the extension, or you can call the >>> # format-specific save methods >>> t[0:10].save_dcd('first-ten-frames.dcd')
The trajectory contains a reference to a topology object, which can come in handy. For example, if you want to save a copy of your trajectory with only alpha carbons present, you can do that pretty easily.
>>> atoms_to_keep = [a.index for a in t.topology.atoms if a.name == 'CA'] >>> t.restrict_atoms(atoms_to_keep) # this acts inplace on the trajectory >>> t.save('CA-only.h5')