Contact Map Analysis in GROMACS with gmx mdmat

To make a residue contact map from a GROMACS trajectory, make the protein whole with gmx trjconv, then run gmx mdmat on a protein index group. It writes dm.xpm, the mean smallest distance between every residue pair, truncated at 1.5 nm by default. Plot it with gmx xpm2ps, or compute a per-pair contact frequency in Python with MDAnalysis.

A contact map turns a whole trajectory into one picture: which residues touch, how often, and what changes between two simulations. It is the natural next step after RMSD and RMSF, and it pairs well with dynamic cross-correlation (DCCM), which shows correlated motion rather than proximity. This guide uses GROMACS 2026 option names from the official manual and a Python script we tested on MDAnalysis 2.7. It was written by the StemSkills Lab team, who have spent more than ten years on structural bioinformatics, drug design and molecular modeling.

What is a residue contact map and what does it show?

A contact map is an N × N matrix for a protein of N residues. Each cell holds a distance, or a yes/no contact, or a contact frequency for one residue pair. The diagonal is always “in contact” (a residue with itself), helices show up as thick bands just off the diagonal, β-sheets as lines parallel or perpendicular to it, and tertiary contacts as patches far from it.

Over a trajectory, the map answers questions RMSD cannot:

• Which contacts hold the fold together for the whole run, and which ones break?
• Does a mutation or a bound ligand create or remove specific contacts?
• Do two domains, or a protein and its partner, stay in contact?

The size grows with the square of the protein. For the 214-residue adenylate kinase test system shipped with MDAnalysis, the matrix has 214 × 214 = 45,796 cells. In our test run, each residue had on average about 9.9 heavy-atom contacts within 4.5 Å in the first frame, so most of the map is empty and the information sits in a small number of cells.

How do you prepare a GROMACS trajectory for contact analysis?

Two steps, both covered in detail in our guides to gmx trjconv and gmx make_ndx index groups.

1. Make the protein whole and centred. A protein split across the periodic box gives distances that are not physical. Fix it once and use the corrected file for every analysis:

gmx trjconv -s md.tpr -f md.xtc -o md_noPBC.xtc -pbc mol -center

Choose Protein for centring and Protein (or System) for output.

2. Make an index group if you need part of the protein. For example, residues 10 to 120 of the protein:

gmx make_ndx -f md.tpr -o index.ndx
> 1 & r 10-120
> q

Group 1 is Protein in the default groups GROMACS builds for a protein system. The new group is added at the end of the list, named after your selection; check its number before you continue.

How do you run gmx mdmat?

The gmx mdmat manual page describes it as making “distance matrices consisting of the smallest distance between residue pairs”. A typical command:

gmx mdmat -f md_noPBC.xtc -s md.tpr -n index.ndx -mean dm.xpm -no num.xvg -t 1.5 -b 10000 -dt 100

In the example, -b 10000 skips the first 10 ns (treat that as equilibration only if your RMSD says so), and -dt 100 reads one frame every 100 ps.

The group you pick changes the definition of distance. mdmat takes the smallest distance between any two atoms of the two residues within your selection. Pick Protein and you get an all-atom minimum distance. Pick C-alpha and you get Cα to Cα distances, which are larger and smoother. State which one you used.

dm.xpm holds distances, not contact frequencies. A cell with a mean distance of 0.5 nm could be a pair in contact half the time and far apart the rest. If you need “how often”, use the Python method below.

How do you turn the .xpm file into a figure?

GROMACS ships gmx xpm2ps for this:

gmx xpm2ps -f dm.xpm -o dm.eps -rainbow blue
ps2pdf -dEPSCrop dm.eps dm.pdf

-rainbow accepts no, blue or red and turns the default grey scale into colour. ps2pdf comes with Ghostscript and gives you a PDF you can place in a thesis.

For two simulations of the same protein, xpm2ps can put one map in each half of a single figure. With -f2, the manual says the upper left half of the first matrix is plotted with the lower right half of the second:

gmx xpm2ps -f wt_dm.xpm -f2 mut_dm.xpm -o wt_vs_mut.eps -rainbow blue

The -combine option offers halves (the default), add, sub, mult and div, so -combine sub gives a difference map. Both matrices need to cover the same residues for either plot to mean anything.

How do you make a contact-frequency heatmap in Python with MDAnalysis?

This script computes, for every residue pair, the fraction of frames in which the minimum heavy-atom distance is below 4.5 Å. It uses MDAnalysis distance_array with the box, so periodic images are handled. We checked its per-residue minimum distances against a brute-force loop over 300 random residue pairs on the MDAnalysis adenylate kinase test trajectory, with zero mismatches.

import numpy as np
import MDAnalysis as mda
from MDAnalysis.analysis.distances import distance_array

u = mda.Universe("md.tpr", "md_noPBC.xtc")
heavy = u.select_atoms("protein and not name H*")
res = heavy.residues
# index of the first heavy atom of each residue, for reduceat
starts = np.searchsorted(heavy.resindices, res.resindices)

cutoff = 4.5          # Å, minimum heavy-atom distance that counts as a contact
stride = 1            # use every frame; raise this for long trajectories
counts = np.zeros((res.n_residues, res.n_residues))
n_frames = 0

for ts in u.trajectory[::stride]:
    d = distance_array(heavy.positions, heavy.positions, box=u.dimensions)
    d = np.minimum.reduceat(d, starts, axis=0)    # min over atoms of residue i
    d = np.minimum.reduceat(d, starts, axis=1)    # min over atoms of residue j
    counts += d < cutoff
    n_frames += 1

freq = counts / n_frames      # contact frequency, 0 to 1
np.save("contact_freq.npy", freq)
np.save("resids.npy", res.resids)
print(n_frames, "frames,", res.n_residues, "residues")

MDAnalysis works in Å, GROMACS in nm, so 4.5 Å is 0.45 nm. The 4.5 Å value matches the default radius of the MDAnalysis contacts.Contacts class; other common choices are a 0.6 nm atom cutoff (the gmx mindist -d default) or 8 Å between Cα or Cβ atoms. Pick one, justify it, and keep it fixed across every system you compare.

Plot it, hiding the trivial contacts between a residue and its two nearest neighbours:

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

freq = np.load("contact_freq.npy")
n = freq.shape[0]
mask = np.abs(np.subtract.outer(np.arange(n), np.arange(n))) < 3
freq_plot = np.where(mask, np.nan, freq)

fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(freq_plot, origin="lower", cmap="viridis", vmin=0, vmax=1)
ax.set_xlabel("Residue index")
ax.set_ylabel("Residue index")
fig.colorbar(im, label="Contact frequency")
fig.savefig("contact_map.png", dpi=300, bbox_inches="tight")

If you are new to the library, start with our MDAnalysis tutorial for GROMACS trajectories. The original reference is Michaud-Agrawal et al., J. Comput. Chem. 32, 2319 to 2327 (2011).

How do you compare wild-type and mutant (or apo and holo) contact maps?

Run the script on both trajectories, saving contact_freq_wt.npy and contact_freq_mut.npy, then subtract. This works only if both systems have the same residues in the same order, which is true for a point mutant and false for a deletion.

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

wt = np.load("contact_freq_wt.npy")
mut = np.load("contact_freq_mut.npy")
resids = np.load("resids.npy")
diff = mut - wt                      # positive: contact gained in the mutant

fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(diff, origin="lower", cmap="RdBu_r", vmin=-1, vmax=1)
fig.colorbar(im, label="Contact frequency, mutant minus wild type")
fig.savefig("contact_diff.png", dpi=300, bbox_inches="tight")

i, j = np.triu_indices_from(diff, k=3)
order = np.argsort(-np.abs(diff[i, j]))[:10]
for k in order:
    print(f"residues {resids[i[k]]} and {resids[j[k]]}: {diff[i[k], j[k]]:+.2f}")

The printout lists the ten pairs that changed most, using the real residue numbers from your topology. A change of +0.6 means the pair is in contact in 60% more frames in the mutant. Before you interpret it, check that the change appears in independent replicas (see how many replicate MD simulations you need) and not only in one run. Our guide to wild-type vs mutant MD simulations covers the rest of that comparison.

Which contact analysis tool should you use?

Why does my contact map look wrong? Common errors and fixes

The .xpm file is hundreds of megabytes

Cause: -frames writes one matrix per frame. The manual warns that unwise options “may generate a large output file”. A 214-residue protein over 10,000 frames is 10,000 matrices of 45,796 cells. Fix: drop -frames unless you need time resolution, and use -dt to thin the trajectory.

The run is slow and the map is enormous

Cause: you selected System, so water and ions are included. Fix: choose Protein, or a custom group from index.ndx passed with -n.

Distant patches appear and vanish between frames

Cause: the protein, or a complex, is broken across the periodic boundary. Fix: rerun gmx trjconv -pbc mol -center and analyse md_noPBC.xtc. For a protein complex, check that both partners stay together in the corrected file in VMD or PyMOL.

Python runs out of memory

Cause: distance_array builds an atoms × atoms matrix. With 5,000 heavy atoms that is 25 million values, about 200 MB per frame in double precision. Fix: select one domain or chain, or use Cα atoms only ("protein and name CA") with an 8 Å cutoff.

Residue numbers in the figure do not match the paper

Cause: the matrix is indexed from 0, while your protein may start at residue 2 or have gaps. Fix: label axes and printouts with the saved resids.npy, as in the difference script.

The WT and mutant difference map is noise

Cause: frequencies from short or single runs, or a different equilibration window in each. Fix: use the same -b window and stride for both systems and average over replicas.

How do you report a contact map in a thesis or paper?

Give the tool and version, the atom selection, the cutoff, the time window and the stride, for example: “Residue contact frequencies were computed with MDAnalysis 2.7 as the fraction of frames with a minimum heavy-atom distance below 4.5 Å, over 10 to 100 ns, every 100 ps.” Show the full map, then a short table of the pairs that changed most. Link each change to structure: a broken salt bridge (check it with salt bridge analysis in GROMACS) or a loop that moves away.

For where this fits in the full analysis sequence, see the GROMACS molecular dynamics learning path and the computational biology skills roadmap.

Frequently asked questions

What does gmx mdmat calculate?

It calculates the smallest distance between every pair of residues in the selected group and, by default, writes the average over the trajectory to dm.xpm. It can also write one matrix per frame with -frames and a contact count with -no.

What cutoff should I use for a residue contact?

Common choices are 4.5 Å for a minimum heavy-atom distance, 0.6 nm for any-atom contacts (the gmx mindist default) and 8 Å between Cα or Cβ atoms. None is universally correct; choose one, state it and use it for every system.

How do I open a GROMACS .xpm file?

Convert it with gmx xpm2ps -f dm.xpm -o dm.eps, then to PDF with ps2pdf -dEPSCrop dm.eps dm.pdf. For full control over the figure, compute the matrix in Python instead.

Is a contact map the same as a DCCM?

No. A contact map shows which residues are close in space. A dynamic cross-correlation map shows which residues move together, which can include residues far apart. They answer different questions and are often shown side by side.

Can I make a contact map between a protein and a ligand?

Yes, though a residue-by-residue matrix is less useful when one partner is a small molecule. Use gmx mindist with the ligand as one group, or a per-residue contact frequency against the ligand in MDAnalysis.