In this example, we're going to find a "centroid" (representitive structure) for a group of conformations. This group might potentially come from clustering, using method like Ward hierarchical clustering.
Note that there are many possible ways to define the centroids. This is just one.
from __future__ import print_function %matplotlib inline import mdtraj as md import numpy as np
Load up a trajectory to use for the example.
traj = md.load('ala2.h5') print(traj)
<mdtraj.Trajectory with 100 frames, 22 atoms, 3 residues, without unitcells>
Lets compute all pairwise rmsds between conformations.
atom_indices = [a.index for a in traj.topology.atoms if a.element.symbol != 'H'] distances = np.empty((traj.n_frames, traj.n_frames)) for i in range(traj.n_frames): distances[i] = md.rmsd(traj, traj, i, atom_indices=atom_indices)
The algorithim we're going to use is relatively simple:
Using $\beta=1$, this is implemented with the following code:
beta = 1 index = np.exp(-beta*distances / distances.std()).sum(axis=1).argmax() print(index)
centroid = traj[index] print(centroid)
<mdtraj.Trajectory with 1 frames, 22 atoms, 3 residues, without unitcells>