How to Build a Free Energy Landscape from a GROMACS MD Trajectory

To build a free energy landscape from a GROMACS trajectory, project the fitted coordinates onto the first two principal components with gmx covar and gmx anaeig, then feed that 2D projection to gmx sham, which converts the sampled population density into a Gibbs free energy surface in kJ/mol. Each basin on that surface is a stable conformational state of your protein.

A free energy landscape (FEL) is the plot reviewers ask for when they want proof that your simulation actually sampled distinct conformational states, not just noise around one structure. Most beginners can run principal component analysis and get a 2D scatter of points, then stall at the real question: how do you convert that scatter into a proper energy surface with labelled minima? This guide walks the full path, with the exact commands, the flags that matter, and a troubleshooting section for the flat and single-well plots that trip people up.

What is a free energy landscape in molecular dynamics?

A free energy landscape maps how the free energy of your system changes across one or two chosen coordinates. In practice you pick two reaction coordinates (most often the first two principal components), bin the trajectory into a 2D histogram of how often each region is visited, and convert that population P into a Gibbs free energy with the Boltzmann relation:

G(x, y) = -kBT ln[ P(x, y) / Pmax ]

The most populated bin becomes the global minimum at 0 kJ/mol, and every other region is reported relative to it. Deep wells are stable states the protein spends a lot of time in; ridges between them are transition regions. The idea of reducing a protein’s motion to a few dominant coordinates comes from the essential dynamics method of Amadei, Linssen and Berendsen (1993), who showed you can separate the motion into what they called “an essential subspace” of a handful of large-amplitude modes (Amadei et al., Proteins, 1993).

One number worth remembering: at 300 K, k_BT (equivalently RT per mole) is about 2.5 kJ/mol. A free energy difference of that size corresponds to roughly a factor of e, about 2.7, in relative population. So a basin that sits 5 kJ/mol above the global minimum is visited only about one-seventh as often. That is why shallow bumps on an FEL rarely mean much, and why a noisy landscape usually means undersampling rather than real states.

What do you need before you start?

You need a finished, equilibrated production trajectory and its run input file. Concretely: a .xtc (or .trr) trajectory, the matching .tpr, and a working GROMACS install (this guide uses commands stable across GROMACS 2020 and later). You also need to have already checked that your simulation is equilibrated, because an FEL built from a drifting trajectory is meaningless. If you are unsure your run converged, sort that out first; the rest of the GROMACS molecular dynamics track covers RMSD-based equilibration checks.

Step 1: Remove overall rotation and translation

Principal component analysis is meant to capture internal motion, so you must first strip out the whole-molecule tumbling and drift. Do a least-squares fit before anything else:

gmx trjconv -s md.tpr -f md.xtc -o md_fit.xtc -fit rot+trans

GROMACS will ask for a group to fit on (choose Backbone or C-alpha) and a group to write out (usually Protein). Skipping this step is the single most common reason a landscape comes out wrong, because the PCA then captures the protein spinning in the box rather than its conformational changes.

Step 2: Build the covariance matrix with gmx covar

Now compute the covariance matrix of atomic fluctuations and diagonalise it to get the eigenvectors (the principal components) and eigenvalues:

gmx covar -s md.tpr -f md_fit.xtc -o eigenval.xvg -v eigenvec.trr -av average.pdb

Choose a consistent atom selection, typically C-alpha or Backbone, for both the fit and the analysis. This writes the eigenvalues (eigenval.xvg), the eigenvectors (eigenvec.trr) and the average structure (average.pdb). The eigenvalues, sorted largest first, tell you how much positional variance each component explains; the first two usually dominate, which is why a 2D landscape on PC1 and PC2 is a reasonable summary. See the official reference for every option: gmx covar documentation.

Step 3: Project the trajectory onto PC1 and PC2 with gmx anaeig

Next, project every frame onto the first two eigenvectors to get the 2D coordinates you will build the landscape from:

gmx anaeig -v eigenvec.trr -f md_fit.xtc -eig eigenval.xvg -s md.tpr -first 1 -last 2 -2d 2dproj.xvg

The -2d output combined with -first 1 -last 2 gives a two-column file: PC1 in column one, PC2 in column two, one row per frame, and importantly no time column. That detail matters for the next step. Full option list: gmx anaeig documentation.

Step 4: Compute the landscape with gmx sham

Feed the 2D projection to gmx sham, which bins it and applies the Boltzmann inversion above:

gmx sham -f 2dproj.xvg -notime -tsham 300 -nlevels 100 -ls FES.xpm -g sham.log

The flags that actually matter here:

• -notime: tells sham the input has no time column. Because 2dproj.xvg is pure (PC1, PC2) data, forgetting this flag makes sham misread your columns and is the top cause of a garbage plot.
• -tsham: the temperature used in the -kT ln P conversion. Set it to your simulation temperature (300 K here). The default is 298.15 K.
• -nlevels: the number of levels written into the output map (default 100). More levels give a smoother contour plot.
• -ls: the Gibbs free energy landscape, written as an .xpm matrix (FES.xpm).

The -g sham.log file records the binning and the free energy range. The authoritative flag reference is the gmx sham documentation.

How do you turn the .xpm output into a publishable contour plot?

The FES.xpm from sham is an ASCII image matrix, not a figure. You have two honest routes to a usable plot.

Quick look: convert it to PostScript to eyeball the result immediately:

gmx xpm2ps -f FES.xpm -o FES.eps -rainbow blue

Publication contour: parse the .xpm yourself and draw a filled contour with matplotlib. The header of an .xpm lists each colour character and the free energy value (in kJ/mol) it maps to; the data lines below are the rows of the grid. Read those two pieces, build a 2D NumPy array of free energies, and plot with plt.contourf, adding a colour bar labelled “Free energy (kJ/mol)”. This gives you full control over axis labels (PC1, PC2, in nm), the colour map, and the contour spacing. If you would rather not write the parser, the same result comes from bundling gmx analysis outputs, but writing a short reader is a genuinely useful skill for the rest of your trajectory work.

How do you read minima and basins on the landscape?

Read it like a topographic map of energy. The colour scale runs from the global minimum at 0 kJ/mol (the most sampled, most stable region) up to the sparsely visited high-energy regions. Practical reading rules:

• Basins (deep wells) are metastable conformational states. A landscape with two or three clear basins is strong evidence your protein interconverts between distinct conformations.
• Barriers (ridges) between basins mark the transition regions; their height, in kJ/mol, is a rough estimate of how hard the interconversion is.
• A single well means either the protein is genuinely rigid on your timescale, or, more often, you did not sample long enough to escape the starting basin.

To find the structures sitting at the bottom of a basin, note the PC1 and PC2 values at the minimum, then use gmx anaeig with -extr or filter frames by their projection values to extract representative snapshots for the figure.

Which reaction coordinates should you use?

PC1 vs PC2 is the default, but it is not the only sensible choice. Pick the pair that best separates the states you care about.

Whatever pair you choose, the sham step is identical: give it a two-column .xvg and the -notime flag.

Troubleshooting: why is my landscape empty or a single well?

These are the failure modes that account for most beginner problems, and their fixes.

• The whole plot is one colour or looks like a single blob. You almost certainly left off -notime, so sham treated your first column as time and binned only the second. Add -notime and rerun.
• Only one deep well, no structure. Usually undersampling. A short trajectory stays trapped in its starting basin. Extend the simulation, or check that your production run was long enough for the motions you expect. It can also be a genuinely rigid system.
• The landscape captures spinning, not conformation. You skipped the -fit rot+trans step, so PCA modelled global rotation. Refit with gmx trjconv and rebuild from Step 1.
• Energies look wrong or the scale is odd. Confirm -tsham matches your simulation temperature. A mismatch rescales the whole surface.
• Axes have no physical meaning. PC projections are in nm and depend on your atom selection. Keep the same selection (for example C-alpha) across covar and anaeig so the numbers are comparable between runs.

Frequently asked questions

Is a free energy landscape from gmx sham a true free energy?

It is a relative free energy derived from the sampled populations of your unbiased trajectory, not an absolute value. It is only as reliable as your sampling. For rare events or high barriers, plain MD undersamples, and you would need enhanced sampling methods such as metadynamics or umbrella sampling to get quantitatively trustworthy barriers.

How long should my trajectory be for a meaningful landscape?

Long enough to visit the states you claim exist, more than once. There is no universal number; a small, flexible peptide may show clear basins in tens of nanoseconds, while a domain motion in a large protein may need hundreds of nanoseconds or more. If your basins appear and disappear when you split the trajectory in half, you have not sampled enough.

Can I make a 1D free energy profile instead of a 2D map?

Yes. Give gmx sham a single-column coordinate (for example just PC1, or RMSD) and it returns a 1D free energy profile. This is useful when one coordinate cleanly separates your states and a 2D map would be overkill.

Do I have to use principal components as the axes?

No. Any two order parameters work, as the comparison table shows. Principal components are popular because they are data-driven and require no prior hypothesis, but RMSD, radius of gyration, native contacts or a custom distance can be more interpretable when you already know the motion you are studying.

Why is my deepest basin not at the centre of the plot?

That is normal. The axes are the principal component values, and the most populated conformation can sit anywhere in that space. The position tells you the collective-coordinate values of your most stable state; there is nothing special about the centre of the plot.

The bottom line

A free energy landscape is four commands away once your trajectory is equilibrated: fit out global motion with gmx trjconv, build the components with gmx covar, project onto PC1 and PC2 with gmx anaeig, and Boltzmann-invert the populations with gmx sham. The craft is in the reading: distinguishing real basins from undersampling, choosing coordinates that separate your states, and reporting energies honestly in kJ/mol. Get those right and you have a figure that tells a clear structural story. For the full staged path from structure preparation to analysis, see the computational biology skills roadmap.

Written by the StemSkills Lab team, with 10+ years of combined experience in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling.