Protein-Ligand Interaction Energy in GROMACS (Rerun Guide)

Add energygrps = Protein LIG to a copy of your production .mdp, build a new tpr with gmx grompp ... -n index.ndx -o ie.tpr, recompute energies over the finished trajectory with gmx mdrun -deffnm ie -rerun md.xtc -nb cpu, then extract Coul-SR:Protein-LIG and LJ-SR:Protein-LIG with gmx energy -f ie.edr. Their sum is not a binding free energy.

The request usually arrives after the simulation is finished. Your protein-ligand trajectory is on disk, the RMSD looks stable, and your supervisor asks for “the interaction energy between the protein and the ligand, split into electrostatic and van der Waals”. GROMACS can report exactly that, but not from the energy file you already have, and not on a GPU. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling) shows the full rerun workflow, what each number means, and how to fix the errors students hit most often. It belongs to our pillar guide to learning molecular dynamics with GROMACS, and the computational biology skills roadmap shows where trajectory analysis fits in the wider skill path.

What does GROMACS mean by protein-ligand interaction energy?

It is the sum of the pairwise non-bonded potential energy between every protein atom and every ligand atom, computed with your force field. GROMACS splits it into a short-range Coulomb term (Coul-SR) and a short-range Lennard-Jones term (LJ-SR). The GROMACS mdp options reference describes energygrps as the groups for which short-ranged non-bonded potential energies are written to the energy file, and notes that the option is not supported on GPUs.

Two consequences follow from the word “short-ranged”. First, only atom pairs inside your cut-offs (rcoulomb, rvdw) contribute. Second, with PME electrostatics the reciprocal-space part (Coul. recip. in the energy menu) is never split per group, so Coul-SR:Protein-LIG is the real-space electrostatic part only.

Justin Lemkul’s widely used protein-ligand complex tutorial is blunt about the limits: “this quantity is NOT a free energy or a binding energy.” He adds that most force fields are not parametrized to make this number physically meaningful on its own, with CHARMM a partial exception because it is fitted against quantum mechanical interaction energies with water. Treat it as a relative, same-force-field comparison tool.

Why can’t you just read it from your production energy file?

Because your production run almost certainly did not define energy groups, and it should not have. Energy group decomposition blocks GPU non-bonded kernels and slows the run for no benefit to the dynamics. The GROMACS source makes the recommended route explicit: when energy groups are present in a normal run on GPU hardware, mdrun advises running on the GPU without energy groups and then using mdrun -rerun on the trajectory with an energy-group tpr.

A rerun reads the saved coordinates frame by frame and recomputes the potential energy with whatever physics the new tpr defines. The mdrun features page states that a rerun takes only positions as input, ignores any velocities present, performs neighbour searching for every frame, and reports potential energies but not kinetic energy, temperature or pressure. That is exactly what an interaction energy needs.

How do you set up energygrps in a rerun .mdp file?

You need two things: groups named in an index file, and a rerun .mdp whose non-bonded settings match your production run.

Step 1: make sure the groups exist

The default groups already include Protein and one group per non-protein residue name, so a ligand named LIG normally has its own group. If your ligand is split across residues, or you want a pocket subset, build the groups first with our gmx make_ndx index groups tutorial. Check the exact spelling in the index file:

grep "\[" index.ndx

Step 2: copy your production .mdp and add one line

cp md.mdp ie.mdp

Open ie.mdp and add:

energygrps = Protein LIG

Leave the cut-off and electrostatics settings exactly as they were in production, otherwise you are computing a different model. For reference, a CHARMM36 setup like Lemkul’s tutorial uses lines such as these; yours may differ, so keep your own values:

cutoff-scheme  = Verlet
coulombtype    = PME
rcoulomb       = 1.2
vdwtype        = Cut-off
vdw-modifier   = Force-switch
rvdw-switch    = 1.0
rvdw           = 1.2

Integrator, thermostat and barostat lines can stay; a rerun does not integrate anything. If you want a refresher on what each keyword does, see the GROMACS .mdp file explained.

Step 3: build the interaction-energy tpr

gmx grompp -f ie.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o ie.tpr

The -n index.ndx flag matters. According to the grompp error text in the GROMACS source, group names in the .mdp must match either [moleculetype] names or custom index groups, and custom groups require the index file. Group name matching is case-insensitive.

How do you run mdrun -rerun over the existing trajectory?

gmx mdrun -deffnm ie -rerun md.xtc -nb cpu

-deffnm ie reads ie.tpr and writes ie.edr, ie.log and friends. -rerun md.xtc names the trajectory to re-evaluate. -nb cpu keeps non-bonded work on the CPU. If you leave -nb on its default, recent GROMACS versions print that multiple energy groups are not implemented for GPUs and fall back to the CPU. If you force -nb gpu, mdrun stops with a fatal error instead.

Three practical points:

• Use the original trajectory. Rerun the .xtc written by mdrun, not one you rotated with trjconv -fit rot+trans. Fitting rotates the coordinates without rotating the box, which breaks the periodic images the energy calculation relies on. Keep the processed files from our gmx trjconv tutorial for visualisation and RMSD.
• Skip equilibration later, not now. Rerun the whole trajectory and trim the start with -b in gmx energy, so you can test several cut-off times without rerunning.
• Expect an informational note. If your topology has constraints, the log says that rerun does not recalculate constraints. That is normal, because rerun only evaluates positions.

Use -ntomp to set CPU threads if the node is shared. Because neighbour searching runs on every frame, rerun time scales with the number of saved frames, not with the original number of MD steps.

How do you extract Coul-SR and LJ-SR with gmx energy?

gmx energy -f ie.edr -o interaction_energy.xvg

The menu now lists group-pair terms named in the pattern term:group1-group2. With energygrps = Protein LIG in a solvated system you will see entries such as:

• Coul-SR:Protein-LIG and LJ-SR:Protein-LIG: the protein-ligand interaction you want.
• Coul-14:Protein-LIG and LJ-14:Protein-LIG: 1-4 pair terms, which are zero for a non-covalent ligand because no bonded 1-4 pairs cross the two molecules.
• Protein-Protein, LIG-LIG, Protein-rest, LIG-rest and rest-rest terms. GROMACS places every atom you did not list (water, ions) in an automatic group called rest, so LIG-rest is the ligand-solvent interaction.

Type the two names (or their numbers), then press Enter on an empty line. To script it:

printf "Coul-SR:Protein-LIG\nLJ-SR:Protein-LIG\n\n" | gmx energy -f ie.edr -b 10000 -o interaction_energy.xvg

Here -b 10000 starts the analysis at 10,000 ps; choose your own start from the equilibration seen in your RMSD and RMSF analysis. Add -sum if you only want the total in the xvg file.

At the end, gmx energy prints an average, an error estimate, an RMSD and a drift for each term. The gmx energy documentation explains that the error estimate comes from block averaging over 5 blocks by default (adjustable with -nbmin and -nbmax), and that these printed statistics use full-precision data, which makes them more accurate than averaging the xvg yourself. Report the printed averages.

What numbers should you expect, and how do you report them?

Values are in kJ/mol and negative for favourable contacts. As a worked reference, Lemkul’s tutorial on the JZ4 ligand bound to T4 lysozyme reports a short-range Coulomb term of -16.5 ± 10.0 kJ/mol and a Lennard-Jones term of -96.8 ± 7.4 kJ/mol. Adding the two and propagating the errors as the square root of the sum of squares gives -113.3 ± 12.4 kJ/mol. For a mostly hydrophobic ligand in a hydrophobic pocket, a dominant LJ term is what you would expect.

Lemkul also cautions that the split into Coulomb and LJ components cannot be verified experimentally, so avoid strong claims such as “binding is driven by van der Waals forces” from these two numbers alone. A defensible Results sentence reads: “The average short-range protein-ligand interaction energy over the last 90 ns was X ± Y kJ/mol (CHARMM36m, block-averaged error).” Our guide to writing an MD methods section shows how to describe the rerun itself.

To plot the time series in Python:

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt("interaction_energy.xvg", comments=["#", "@"])
t_ns = data[:, 0] / 1000
coul, lj = data[:, 1], data[:, 2]

plt.plot(t_ns, coul, label="Coul-SR")
plt.plot(t_ns, lj, label="LJ-SR")
plt.plot(t_ns, coul + lj, label="Total", linewidth=2)
plt.xlabel("Time (ns)")
plt.ylabel("Interaction energy (kJ/mol)")
plt.legend()
plt.savefig("interaction_energy.png", dpi=300)

For styling, running averages and multi-panel figures, see how to plot GROMACS xvg files in Python.

How is interaction energy different from MM-PBSA and a docking score?

All three produce a negative number in energy units, which is why students mix them up. They answer different questions.

If you need a solvation-aware estimate, move on to our MM-PBSA binding free energy guide, and for a rigorous pathway-based number see umbrella sampling in GROMACS. To see which residues make the contacts behind the energy, pair this with ProLIF interaction fingerprints.

What goes wrong, and how do you fix it?

The Protein-LIG terms are missing from the gmx energy menu

Either energygrps was not in the .mdp used for grompp, or you opened the production md.edr instead of ie.edr. Confirm the line is in ie.mdp, rebuild ie.tpr, rerun, and read ie.edr. Editing the .mdp without rerunning grompp changes nothing, because mdrun reads only the tpr.

“Group LIG referenced in the .mdp file was not found in the list of index groups”

The name in energygrps does not match any group. Check the spelling against grep "\[" index.ndx; many ligands are named UNL, MOL or a three-letter PDB code instead of LIG. If the group is custom, add -n index.ndx to grompp.

“Number of atoms in trajectory (N) does not match the run input file (M)”

Your xtc probably contains a subset of atoms, typically because production used compressed-x-grps to save only the protein and ligand. The tpr describes the full solvated system. Rerun a trajectory that contains all atoms (a full-precision .trr if you wrote one), or build a tpr for the same subset; the mdrun features page points to gmx convert-tpr for subset work. Plan ahead: if interaction energy is on your analysis list, save full-system coordinates.

“Nonbonded interactions on the GPU were required, but not supported for these simulation settings”

You passed -nb gpu (or a job script template did). Replace it with -nb cpu. Also remove GPU-only options such as -pme gpu or -update gpu from the rerun line.

“My trajectory has no velocities, will rerun fail?”

No. Rerun uses positions only and ignores velocities, so a compressed xtc is fine. What you will not get are kinetic energy, temperature and pressure; do not look for them in ie.edr.

The numbers jump wildly or look unphysical

Check that you reran the unfitted trajectory, that ie.mdp used the same cut-offs and modifiers as production, and that the ligand did not leave the pocket (plot the protein-ligand distance). A ligand drifting into solvent produces an interaction energy that climbs towards zero, which is a real result, not a bug.

Frequently asked questions

Is protein-ligand interaction energy the same as binding free energy?

No. It is the average short-range non-bonded energy between two groups. It ignores desolvation, entropy and the long-range PME reciprocal term, so it cannot be converted into a binding free energy or a dissociation constant.

Can I calculate interaction energy on a GPU in GROMACS?

Not the energy-group decomposition itself. Run production on the GPU without energygrps, then run the short mdrun -rerun on CPU with an energy-group tpr. With the default -nb auto, mdrun falls back to CPU on its own.

Why is my Coulomb term close to zero for a charged ligand?

Possible reasons include strong screening by nearby water and ions, a pose where charged groups face solvent rather than the protein, or partial charges from ligand parameterization. Check your ligand topology with our ligand parameterization guide.

Can I get per-residue interaction energies this way?

Yes, by listing more groups, for example energygrps = LIG r_45 r_89 r_120 after creating residue groups with make_ndx. Every extra group adds pair terms and slows the rerun, so keep the list short. For full per-residue decomposition, MM-PBSA tools are more practical.

Which time window should I average over?

Exclude the initial settling period visible in your ligand RMSD and energy time series, then average the stable part with gmx energy -b. Report the window explicitly, for example “the last 90 ns of a 100 ns production run”.

Where does this fit in your MD workflow?

Interaction energy is a late analysis step: after the protein-ligand MD simulation, after index groups and trajectory cleanup, and usually before MM-PBSA. For the method background behind the whole workflow, Lemkul’s peer-reviewed tutorial collection, From Proteins to Perturbed Hamiltonians (Living Journal of Computational Molecular Science, 2018), is the standard citation for these GROMACS tutorials.