Per-Residue Energy Decomposition with gmx_MMPBSA

Per-residue energy decomposition splits a single MM-GBSA or MM-PBSA binding free energy into the contribution of each residue. In gmx_MMPBSA you add a &decomp namelist with idecomp=2, dec_verbose=1 and print_res=”within 4″, rerun the calculation, then read the per-residue totals from FINAL_DECOMP_MMPBSA.dat. The residues with the most negative totals are your hotspots.

Most students stop at one number. The MM-GBSA run finishes, the binding free energy comes out at some value in kcal/mol, it goes into the results table, and the discussion says the ligand binds well. Then an examiner asks which residues are actually holding it, and the answer is a guess based on a picture. Decomposition is how you replace the guess with a number per residue. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling) covers the setup, the output format and the mistakes that waste a week, with every flag and default checked against the official gmx_MMPBSA documentation on 19 September 2026. It follows directly from our tutorial on binding free energy calculation with MM-PBSA in GROMACS, sits inside our pillar guide to learning molecular dynamics with GROMACS, and the computational biology skills roadmap shows where this step belongs in the wider skill path.

What does per-residue energy decomposition actually tell you?

A standard MM-GBSA calculation gives you one binding free energy for the whole complex. Decomposition takes the same trajectory frames and reports, for every residue you ask about, how much of that energy that residue contributed. You get the breakdown by term as well: van der Waals, electrostatic, polar solvation, non-polar solvation, and the total.

The residues with large negative totals are the ones doing the binding. The idea is older than the software and comes from experimental alanine scanning. Bogan and Thorn compiled a database of 2,325 alanine mutants with measured changes in binding free energy and reported in the Journal of Molecular Biology that “the free energy of binding is not evenly distributed across interfaces; instead, there are hot spots of binding energy made up of a small subset of residues in the dimer interface.” They also found those hot spots enriched in tryptophan, tyrosine and arginine. A computational decomposition is an attempt to reproduce that pattern from a simulation rather than from mutagenesis.

Two limits are worth stating before you build a figure on this. Decomposition is a partition of a model energy, not a prediction of an experimental mutation. Deleting a residue in a real protein changes packing and solvation in ways a per-residue split does not capture. And the entropy term is not decomposed, so what you are reading is an enthalpic and solvation partition, not a per-residue free energy in the strict thermodynamic sense.

What do you need before you run decomposition?

You need exactly what a normal gmx_MMPBSA run needs, because decomposition is not a separate calculation. It is the same calculation with extra reporting switched on.

• A finished, equilibrated production run. Decomposition on a trajectory that has not settled gives you noise with three decimal places.
• A TPR file for the complex, passed with -cs.
• A trajectory with periodic boundary artefacts removed, passed with -ct. Our gmx trjconv tutorial covers the processing.
• An index file with a receptor group and a ligand group, passed with -ci and named with -cg. Build it with gmx make_ndx.
• The GROMACS topology, passed with -cp.
• gmx_MMPBSA itself. Version 1.7.0 is the current release on PyPI and it requires Python 3.11 or 3.12.

If your MM-GBSA run has not worked yet, fix that first. Decomposition inherits every input problem the base calculation has and makes the failure slower to reach.

How do you write the &decomp namelist?

You add one namelist to the input file you already have. gmx_MMPBSA will even write a starting file for you:

gmx_MMPBSA --create_input decomp

The sample decomposition input in the official input file documentation is this:

&general
startframe=5, endframe=21, interval=1,
/

&gb
igb=8, saltcon=0.150,
/

&decomp
idecomp=2, dec_verbose=3,
print_res="within 4"
/

Three variables carry the work: idecomp chooses the scheme, dec_verbose chooses how much gets printed, and print_res chooses which residues get printed.

Which idecomp value should you choose?

There are four values. Two are per-residue and two are pairwise. The difference within each pair is only where the 1-4 non-bonded terms are booked.

For a thesis figure that names hotspot residues, use idecomp=2. It is the documented default and it gives one number per residue, which is what a bar chart needs. Reach for pairwise only when your research question is specifically about which residue talks to which, for example when you are arguing that a particular salt bridge carries the interaction.

Pairwise is expensive and the documentation is blunt about it. The data scales as O(N²), and the docs warn that large print selections “demand a large amount of memory to parse the mdout files and write decomposition output file (~500 MB for just 250 residues, since that’s 62500 pairs!)”. Run pairwise on a tight residue selection or not at all.

What does dec_verbose control?

It sets how much detail reaches the decomposition output file. The default is 1.

Use 1 for a normal hotspot analysis. The sidechain and backbone split is the part that makes a claim defensible: if a residue’s contribution is mostly backbone, calling it a sidechain hotspot and proposing a mutation to alanine does not follow. Use 3 only when you need the complex, receptor and ligand terms separately, and accept the larger files.

How do you select which residues get printed?

The print_res default is "within 6", which prints every residue within 6 Angstrom across the interface. Most published protein-ligand figures use a tighter shell:

print_res="within 4"

The distance form has a useful safety behaviour. If your cutoff is so small that fewer than two residues would be selected, gmx_MMPBSA increases the cutoff by 0.1 until at least two are selected, so a typo in the distance produces a small selection rather than an empty one.

You can also name residues explicitly with a chain and number mask:

print_res="A/1,3-10,15,100 B/25"

If you use the explicit form, the documentation carries a warning in red: include at least one residue from both the receptor and the ligand in the mask. The distance form guarantees this automatically, which is another reason to prefer it. Insertion codes attach directly to the number, so A/27B is valid and A/27:B is not, and numeric ranges pick up only blank-code residues, so print_res="A/5-6B,6D-7" is an error.

Avoid print_res="all". The docs call it “generally not recommended because most residue contributions are zero”, and with pairwise decomposition it is the fastest route to a run that fills your disk.

What is the full gmx_MMPBSA command for a decomposition run?

The only change from a normal run is that you ask for the two decomposition output files. The documented command for a decomposition example is:

gmx_MMPBSA -O \
  -i mmpbsa.in \
  -cs com.tpr \
  -ct com_traj.xtc \
  -ci index.ndx \
  -cg SOLU_chain1 SOLU_chain2 \
  -cp topol.top \
  -o FINAL_RESULTS_MMPBSA.dat \
  -eo FINAL_RESULTS_MMPBSA.csv \
  -do FINAL_DECOMP_MMPBSA.dat \
  -deo FINAL_DECOMP_MMPBSA.csv

Two flags are the ones people miss. -do writes the decomposition statistics summary and -deo writes the per-residue energy terms in CSV form. Without them the decomposition is computed and then dropped.

Two more flags are worth knowing. -nogui stops the graphical analyser from opening automatically when the run finishes, which matters on a cluster or over SSH. And -cr passes a reference structure; the documentation recommends it to guarantee that the residue you selected in print_res is the residue the topology numbering points at.

How do you read FINAL_DECOMP_MMPBSA.dat?

The file is organised in sections. With dec_verbose=1 you get a Delta section, and inside it three blocks: Total Energy Decomposition (TDC), Sidechain Energy Decomposition (SDC) and Backbone Energy Decomposition (BDC). With dec_verbose=3 the same three blocks appear first for the complex, then the receptor, then the ligand, then the delta.

Each residue row carries six energy terms: internal, van der Waals, electrostatic, polar solvation, non-polar solvation, and the total. Each term is reported as an average, a standard deviation and a standard error of the mean across the frames you used. The csv_format default of 1 is what puts the standard error in there, and the CSV opens directly in a spreadsheet.

Read it in this order:

1. Go to the Delta section, TDC block. This is the binding contribution, not the absolute energy inside the complex.
2. Sort by the TOTAL column, most negative first. The top few residues are your candidate hotspots.
3. Check the standard error next to each total. A residue at -3.1 kcal/mol with a standard error of 0.2 is a result. The same -3.1 with a standard error of 2.4 is a residue that matters in some frames and not others, which is a different and more interesting claim.
4. Look at the term split. A van der Waals dominated hotspot is a packing contact. An electrostatic dominated one is a salt bridge or hydrogen bond, and it should be visible in your interaction analysis too.
5. Compare TDC against SDC. If the total is carried by the backbone, the residue identity is less important than its position.

Cross-check the top residues against a contact analysis of the same trajectory. Our guide to protein-ligand interaction analysis with ProLIF gives an occupancy number per interaction, and a residue that shows up as both a high-occupancy contact and a large negative decomposition total is a far stronger claim than either alone.

How do you turn the numbers into a thesis figure?

gmx_MMPBSA ships with a graphical analyser. It opens automatically after a run unless you passed -nogui, and you can start it later on a finished run:

gmx_MMPBSA_ana -f .

The -f flag takes info files, a containing folder, or a list of them, and defaults to the current working directory. Add -r to search recursively one level down, which is how you load several systems at once for comparison.

For decomposition it produces bar plots with error bars, heatmaps of per-residue contributions by frame, and per-wise heatmaps for pairwise runs. It also writes PyMOL representations that colour the complex by per-residue energy, which is the figure most examiners respond to. The loading dialog has a “Remove non-contributing residues” option with a default per-residue threshold of 0.5 kcal/mol, which clears the flat residues out of a crowded bar chart.

If you would rather build the figure yourself, take FINAL_DECOMP_MMPBSA.csv into Python or a spreadsheet. Plot the TOTAL column per residue with the standard error as the error bar, and label only the residues past your own threshold. State the threshold in the caption.

What goes wrong, and how do you fix it?

ValueError: could not convert string to float: ‘*************’. The documented cause is an inconsistent structure or a trajectory that still has periodic boundary artefacts, usually when the complex is longer than a box edge. The fix is upstream: regenerate the structure with gmx editconf -f md.tpr -o md.pdb, look at it, then process the trajectory with gmx trjconv -s md.tpr -f md.xtc -o md_noPBC.xtc -pbc mol -center -n -ur compact using a group that contains both molecules as the centering group. The second documented cause is a genuinely enormous energy value that will not fit in the output field, which happens when the selected group is very large.

The decomposition files are empty or were never written. You left out -do and -deo. Rerun with them, or use -rewrite-output to parse the existing output files again rather than repeating the whole calculation.

The residue numbers in the output do not match your paper. gmx_MMPBSA reports topology numbering, which is often not crystallographic numbering, particularly when a construct starts at residue 1 or a chain was extracted. Pass a reference structure with -cr and map the numbering explicitly before you write any residue name in the text.

The run takes far longer than the base MM-GBSA calculation. Check whether you set idecomp=3 or 4 with a wide print_res. Pairwise decomposition over a broad selection scales as the square of the residue count. Switch to idecomp=2, or tighten the selection.

Only one or two residues appear in the output. Your print_res distance was too small, and the automatic 0.1 increase stopped as soon as two residues qualified. Widen it to "within 4" or "within 6".

The numbers look implausible in magnitude. Check the solvation settings in the base run before blaming decomposition. With the PB model, inp=1 is the default and models the non-polar term as a single term proportional to solvent-accessible surface area; the two-term cavity plus dispersion model needs inp=2 set explicitly, and the dispersion term can be large.

How do you defend a hotspot claim in a viva?

Four sentences protect the result. State the scheme and selection you used, for example per-residue decomposition with idecomp=2 and print_res="within 4". State the frames, because a decomposition over 16 frames and one over 1,000 are different claims. Report the standard error alongside every residue you name. And say plainly that decomposition partitions a calculated energy rather than predicting a mutational effect, then point at whatever independent evidence you have, such as contact occupancy from the trajectory or a published alanine scan of the same target.

Cite the software properly as well. The method implementation is described by Valdés-Tresanco, Valdés-Tresanco, Valiente and Moreno in gmx_MMPBSA: A New Tool to Perform End-State Free Energy Calculations with GROMACS, Journal of Chemical Theory and Computation 17(10), 6281-6291, 2021. Crossref records 3,370 works citing it as of 19 September 2026, which tells you how standard the tool has become. The developers also ask you to cite the underlying MMPBSA.py paper by Miller and co-workers, Journal of Chemical Theory and Computation 8(9), 3314-3321, 2012. Our guide to writing the methods section of an MD study shows where these go in a manuscript.

Frequently asked questions

Does decomposition work with MM-PBSA as well as MM-GBSA?

Yes. The &decomp namelist is independent of whether you are running the GB or PB model, and the polar solvation column in the output will be EGB or EPB accordingly. GB is much cheaper, so most per-residue analyses are run with GB.

Can I decompose the entropy contribution?

No. Entropy methods such as normal mode analysis, interaction entropy and C2 entropy produce a value for the system, not per residue. Any per-residue number you report is an enthalpic and solvation partition, and your methods section should say so.

How many frames do I need for stable per-residue values?

There is no fixed number, and anyone quoting one is guessing. Judge it from the data: run the decomposition on the first and second halves of your production trajectory separately and compare the ranking of the top residues. If the ranking is stable and the standard errors are small relative to the differences between residues, you have enough.

What counts as a hotspot residue?

There is no universal cutoff in the computational output. Pick a threshold, state it, and apply it consistently. The gmx_MMPBSA analyser uses 0.5 kcal/mol as its default threshold for hiding non-contributing residues, which is a reasonable floor for “this residue does something”. Residues past 1 to 2 kcal/mol with small standard errors are the ones worth naming.

Why is a residue’s contribution positive?

A positive total means that residue opposes binding in the model, usually through a desolvation penalty that the electrostatic gain does not cover. Charged residues at the rim of an interface frequently come out positive. This is a normal result and often the more interesting one, because it points at where the ligand could be modified.

Do I need to rerun the whole simulation to add decomposition?

No. Decomposition reuses the same trajectory and topology. You only rerun gmx_MMPBSA with the extra namelist and the -do and -deo flags.