How to do covariance analysis

In Gromacs

gmx covar -f md.trr -s md.tpr -n index.ndx
gmx anaeig -v eigenvec.trr -s ca.pdb -first 1 -last 2 -f traj.trr -2d all.xvg

From scratch

The source code can be found here

1. Get coordinates using MDTraj
2. Align the structures
3. Subtract by mean
4. Divided by SDs
5. Combine residual fluctuations (3Nx3N > NxN)

The key computation is to compute the covariances:

X = sel_traj.reshape(traj.n_frames, sel.shape[0] * 3)
X -= X.mean(axis=0)
X /= X.std(axis=0)
emp_cov = np.dot(X.T, X)/traj.n_frames

and to combine intra-residue covariances:

for ii in range(0, sel.shape[0]):
   for jj in range(0, sel.shape[0]):
       emp_cov_residue[ii][jj] = ( emp_cov[ii*3][jj*3] + emp_cov[ii*3+1][jj*3+1] + emp_cov[ii*3+2][jj*3+2] ) / 3