1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
|
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')
|