DCCM Analysis of a GROMACS Trajectory, Step by Step

Short answer: A dynamic cross-correlation matrix (DCCM) shows which residues move together during an MD run. To build one from GROMACS, remove periodic jumps and fit the trajectory on C-alpha atoms with gmx trjconv, then compute normalized covariances of C-alpha displacements with Bio3D dccm(), ProDy, or a short NumPy script, and plot the values from -1 to +1 as a heatmap.

The red and blue correlation map turns up in the results sections of many protein-ligand MD papers, yet most tutorials stop at RMSD and RMSF. This guide, written by the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling), walks through the full route: trajectory preparation, three ways to compute the matrix, plotting, interpretation, and the errors that produce misleading maps. It sits in our GROMACS molecular dynamics learning path, after RMSD and RMSF analysis.

What does a dynamic cross-correlation matrix actually measure?

For every pair of atoms i and j (usually C-alpha atoms, one per residue), the DCCM stores a single number:

C(i,j) = <Δri · Δrj> / sqrt( <Δri2> · <Δrj2> )

Here Δr_i is the displacement of atom i from its average position in one frame, and the angle brackets mean an average over all frames. The numerator is the covariance of the two displacement vectors; the denominator normalizes it, so every value falls between -1 and +1. The formulation goes back to Ichiye and Karplus, “Collective motions in proteins: A covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations” (Proteins 11:205-217, 1991), which is still the standard citation for the method in a methods section.

• C(i,j) close to +1: the two residues tend to move in the same direction at the same time (correlated).
• C(i,j) close to -1: they move in opposite directions (anticorrelated), for example two lobes closing on a ligand.
• C(i,j) close to 0: no linear relationship between their motions.
• The diagonal is always +1, because each residue is perfectly correlated with itself.

The matrix is N x N for N residues. Hen egg-white lysozyme (PDB 1AKI, 129 residues), the system used in many beginner GROMACS tutorials, gives a 129 x 129 matrix with 16,641 entries, symmetric about the diagonal.

How do I prepare a GROMACS trajectory for DCCM analysis?

Preparation matters more than the choice of software. A DCCM computed on raw md.xtc output is dominated by whole-molecule rotation and periodic boundary jumps, not internal motion. Two gmx trjconv passes fix this (full option list in the official gmx trjconv documentation; our trjconv tutorial explains each flag).

Step 1. Make molecules whole and centre the protein.

printf "Protein\nSystem\n" | gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -pbc mol -center

The first group answers the centering prompt, the second picks what to write out. If your protein is a multimer that drifts apart across the box, try -pbc cluster instead of -pbc mol.

Step 2. Fit on C-alpha atoms and keep only C-alpha atoms.

printf "C-alpha\nC-alpha\n" | gmx trjconv -s md.tpr -f md_center.xtc -o md_fit_ca.xtc -fit rot+trans -b 10000

The first C-alpha is the fit group and the second is the output group. -b 10000 drops the first 10,000 ps (10 ns); set it to wherever your RMSD plot shows the system has settled, because including the equilibration drift inflates correlations. Add -dt 10 if you want one frame every 10 ps to keep the file small.

Step 3. Write a matching C-alpha reference structure.

echo "C-alpha" | gmx trjconv -s md.tpr -f md_fit_ca.xtc -o ca_ref.pdb -dump 0

This PDB has exactly the same atoms, in the same order, as the fitted trajectory. Every tool below needs that one-to-one match. If C-alpha is not in your default groups, create it with gmx make_ndx and pass -n index.ndx.

Step 4 (Bio3D and ProDy only). Convert XTC to DCD. Bio3D’s read.dcd() and ProDy’s parseDCD() read DCD files, not XTC. MDAnalysis converts in a few lines:

import MDAnalysis as mda

u = mda.Universe("ca_ref.pdb", "md_fit_ca.xtc")
with mda.Writer("md_fit_ca.dcd", u.atoms.n_atoms) as w:
    for ts in u.trajectory:
        w.write(u.atoms)

How do I calculate the DCCM with Bio3D in R?

Bio3D (Grant et al., Bioinformatics 22:2695-2696, 2006) is the most common R route to a DCCM. The code below follows the example in the official dccm.xyz reference page:

library(bio3d)

pdb  <- read.pdb("ca_ref.pdb")
trj  <- read.dcd("md_fit_ca.dcd")
inds <- atom.select(pdb, elety = "CA")

xyz <- fit.xyz(fixed = pdb$xyz, mobile = trj,
               fixed.inds = inds$xyz, mobile.inds = inds$xyz)

cij <- dccm(xyz, ncore = 4)
plot(cij, resno = pdb)
write.table(cij, "dccm_bio3d.tsv", sep = "\t")

The fit.xyz() call is a second superposition inside R. It is redundant after Step 2, but harmless, and it protects you if you ever skip the trjconv fit. ncore runs the calculation in parallel; ncore = NULL uses every detected core. Bio3D can also draw the strongest correlations onto the structure in PyMOL with pymol(cij, pdb), which is the fastest way to see whether a correlated block is a real structural unit.

How do I calculate it with ProDy in Python?

ProDy (Bakan, Meireles and Bahar, Bioinformatics 27:1575-1577, 2011) computes cross-correlations from modes, so the route is PCA first, then calcCrossCorr(). The setup follows the ProDy essential dynamics tutorial:

from prody import *

structure = parsePDB("ca_ref.pdb")
ensemble = parseDCD("md_fit_ca.dcd")
ensemble.setCoords(structure)
ensemble.setAtoms(structure.calpha)
ensemble.superpose()

pca = PCA("md")
pca.buildCovariance(ensemble)
pca.calcModes(n_modes=None)      # all modes, not the default 20

cc = calcCrossCorr(pca)
showCrossCorr(pca)

The n_modes argument is the detail that catches people. calcModes() computes 20 modes by default, and a cross-correlation built from 20 modes is a filtered map of the large-scale motions only. With all modes it reproduces the full DCCM. Both are useful, but label which one you are showing. parseDCD() loads the whole file into memory; for very long runs the tutorial’s DCDFile route streams frames instead. If PCA itself is new to you, read our PCA and essential dynamics guide first.

Can I compute a DCCM with plain NumPy and MDAnalysis?

Yes, and writing it yourself is the best way to understand the formula. MDAnalysis (Michaud-Agrawal et al., J. Comput. Chem. 32:2319-2327, 2011) reads the XTC directly, so no DCD conversion is needed:

import numpy as np
import MDAnalysis as mda
import matplotlib.pyplot as plt

u = mda.Universe("ca_ref.pdb", "md_fit_ca.xtc")
ca = u.select_atoms("name CA")

X = np.array([ca.positions.copy() for ts in u.trajectory])   # (frames, N, 3)
d = X - X.mean(axis=0)                                       # displacements
cov = np.einsum("tik,tjk->ij", d, d) / X.shape[0]            # <dri . drj>
norm = np.sqrt(np.diag(cov))
C = cov / np.outer(norm, norm)

np.savetxt("dccm_numpy.dat", C, fmt="%.4f")

plt.imshow(C, cmap="bwr", vmin=-1, vmax=1, origin="lower")
plt.colorbar(label="C(i,j)")
plt.xlabel("Residue index"); plt.ylabel("Residue index")
plt.savefig("dccm.png", dpi=300)

The einsum line sums the x, y and z products for every pair at once. Memory is easy to estimate: 10,000 frames of 300 C-alpha atoms is 10,000 x 300 x 3 = 9 million coordinates, about 36 MB in 32-bit floats. If you skipped the trjconv fit, align in Python first with MDAnalysis.analysis.align.AlignTraj(u, u, select="name CA", in_memory=True).run(). Our MDAnalysis tutorial covers the loading basics.

Does GROMACS have a built-in DCCM tool?

Not directly. gmx covar gets you halfway. Its documentation says the -xpma option writes the atomic covariance matrix, where for each atom pair “the sum of the xx, yy and zz covariances is written.” That is the numerator of the DCCM formula, not the normalized value, so large-amplitude loops dominate the map. To get a true DCCM from it you would export the full matrix with -ascii and divide each entry by the square roots of the matching diagonal terms. For most students, one of the three routes above is simpler.

Which tool should I use for DCCM analysis?

All correctly prepared routes give the same full-mode matrix to within rounding, which is a useful sanity check: run two of them once and compare.

How do I read a DCCM heatmap?

Start with the blocks along the diagonal. A red square there is a stretch of sequence that moves as a unit, usually a helix, a beta-sheet or a domain. Then look off the diagonal:

• Red off-diagonal patches connect distant parts of the sequence that move together, such as two strands of the same sheet or two loops that line the same pocket.
• Blue off-diagonal patches mark anticorrelated motion, the classic signature of hinge bending or of two domains closing and opening.
• Pale regions are uncorrelated. A binding-site loop that is pale in the apo run but red with a neighbouring helix in the holo run suggests the ligand couples them.

For an apo versus holo comparison, compute both matrices on the same residue set, with the same trimming and frame spacing, then plot the difference C_holo - C_apo on a -1 to +1 colour scale. Map the strongest changes onto the structure before writing a sentence about them.

Two limits belong in your discussion section. First, correlation is not causation: a DCCM cannot tell you which residue drives which. Second, the dot product misses correlated motions that are perpendicular to each other, so a zero does not always mean independence. Lange and Grubmüller addressed this with a mutual-information measure in “Generalized correlation for biomolecular dynamics” (Proteins 62:1053-1061), which is what Bio3D’s method = "lmi" option computes.

How many frames do I need for a reliable DCCM?

There is no universal number, because it depends on how slow the motions are, so test convergence instead of guessing. Split the production trajectory into two halves, compute a DCCM for each, and compare them: correlate the upper-triangle values of the two matrices with np.corrcoef, and look at the two heatmaps side by side. If the major blocks appear in one half and not the other, the run has not sampled that motion enough, and the full-length map is not yet a result. Replicate runs with different starting velocities are an even stronger check.

What goes wrong, and how do I fix it?

• Almost the whole map is deep red. The trajectory was not fitted, so overall rotation and translation make every residue look correlated. Rerun Step 2 with -fit rot+trans.
• Stripes of strong correlation through one residue range. Usually a periodic boundary jump: part of the protein crossed the box edge. Redo Step 1 with -pbc mol -center, or -pbc nojump first on the raw file for multimers.
• Atom-count mismatch errors when loading. The PDB and trajectory do not hold the same atoms. Regenerate ca_ref.pdb from md_fit_ca.xtc exactly as in Step 3.
• Bio3D or ProDy cannot open the XTC file. Neither reads XTC; convert with the MDAnalysis snippet in Step 4.
• ProDy map looks smoother than Bio3D’s. You used the default 20 modes. Call calcModes(n_modes=None) for the full matrix.
• Memory errors or very slow runs. You kept all atoms. An all-atom matrix for a 3,000-atom protein has 9 million entries. Stay with C-alpha atoms and stride frames with -dt.
• Residue numbers on the axis do not match the paper numbering. The array index starts at 0 or 1, not at your first residue number. Pass resno = pdb in Bio3D, or relabel the axis ticks from ca.resids in Python. See our guide on plotting GROMACS data in Python for figure formatting.

Where does DCCM fit in a full MD analysis?

DCCM usually comes after stability checks (RMSD, RMSF, radius of gyration) and alongside PCA, and before free-energy or interaction analyses. A typical results section reads: the system is stable, these regions fluctuate, these regions move together, and here is how the ligand changes that. The computational biology skills roadmap shows where this sits relative to docking, MD setup and free-energy methods.

Frequently asked questions

Should I use C-alpha atoms or all atoms for DCCM?

Use C-alpha atoms. They give one value per residue, match how DCCMs are reported in the literature, and keep the matrix small. All-atom matrices add side-chain noise and grow with the square of the atom count.

What is a “strong” correlation in a DCCM?

There is no fixed cutoff. Whatever threshold you choose, state it in the methods and check that the patterns you discuss are stable across trajectory halves or replicates.

Can I compute a DCCM for a protein-ligand complex that includes the ligand?

Yes. Add one representative ligand atom (or its centre of mass per frame) to the selection in the NumPy route. Keep the fit on protein C-alpha atoms so the ligand’s motion is measured relative to the protein.

Is a DCCM the same as PCA?

No, but they come from the same covariance matrix. PCA diagonalizes the covariance to find the main collective motions; DCCM normalizes it pair by pair to show which residues move together. Computing cross-correlations from all PCA modes gives back the DCCM.

How should I cite DCCM analysis in my thesis?

Cite Ichiye and Karplus (1991) for the method and the software paper for the tool you used: Grant et al. (2006) for Bio3D, Bakan et al. (2011) for ProDy, or Michaud-Agrawal et al. (2011) for MDAnalysis, plus GROMACS for the simulation.