Interface with mdtrajΒΆ

try pytraj online:

http://mybinder.org/badge.svg

Iterface with mdtraj

In [21]:
# It's not uncommon to use a similiar package to load data to pytraj.
# In this tutorial, we will show how to convert mdtraj's Trajectory to pytraj's Trajectory
# Why? for example, you want to load '.xtc' file extension that is not supported in pytraj/cpptraj

# load pytraj
import pytraj as pt

# load mdtraj
import mdtraj as md
from mdtraj.testing import get_fn
In [22]:
# load h5 file 
filename, topology_name = get_fn('monolayer.xtc'), get_fn('monolayer.pdb')
In [23]:
pt.load(filename, topology_name)
Out[23]:
<pytraj.Trajectory, 0 frames, include:
<Topology: 3600 atoms, 1 residues, 0 mols, PBC with box type = ortho>>
In [24]:
# if you can see, load xtc file extension results 0 frame because pytraj/cpptraj does not know how to read it.
# try mdtraj
m_traj = md.load(filename, top=topology_name)
m_traj
Out[24]:
<mdtraj.Trajectory with 101 frames, 3600 atoms, 1 residues, and unitcells at 0x2aaad3727b00>
In [25]:
# convert mdtraj.Trajectory to pytraj.Trajectory
# need to provide topology_name (a pdb file)
# and convert m_traj's dtype of float32 to float64

traj = pt.Trajectory(xyz=m_traj.xyz.astype('f8'), top=topology_name)
traj
Out[25]:
<pytraj.Trajectory, 101 frames, include:
<Topology: 3600 atoms, 1 residues, 0 mols, PBC with box type = ortho>>
In [26]:
# compare mdtraj's m_traj and pytraj' s traj
(m_traj.xyz - traj.xyz)[:2]
Out[26]:
array([[[ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        ..., 
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.]],

       [[ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        ..., 
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.]]])
In [28]:
# perfom some calculation with pytraj
pt.rmsd(traj, ref=0)
Out[28]:
array([ 0.        ,  0.78563455,  0.19675571, ...,  0.78462811,
        0.56851667,  0.57288392])
In [35]:
# mdtraj
print(set(atom.name for atom in m_traj.top.atoms))

# pytraj
print(set(atom.name for atom in traj.top.atoms))

{'Si', 'H', 'C', 'O'}
{'Si', 'H', 'C', 'O'}