How to Calculate Protein-Ligand Binding Free Energy in GROMACS (MM-PBSA Tutorial)

To calculate protein-ligand binding free energy in GROMACS, finish a stable MD run of the complex, then use the free gmx_MMPBSA tool to run an MM-PBSA or MM-GBSA end-state calculation on the trajectory. It returns a binding energy in kcal/mol as the difference between the bound complex and the separated protein and ligand, plus a per-residue breakdown of who contributes.

You have run a full molecular dynamics simulation of your docked complex and it stayed stable. A stable trajectory tells you the pose survives, but it does not tell you how strongly the ligand binds. To rank your compounds and put a defensible number in a thesis, you need a binding free energy. This guide from the StemSkills Lab team (10+ years in structural bioinformatics and molecular modeling) shows the practical MM-PBSA route with the open-source gmx_MMPBSA tool.

What is binding free energy and what does MM-PBSA calculate?

The binding free energy (ΔG_bind) is the free energy released when a ligand goes from free in solution to bound in the protein pocket. A more negative value means tighter binding. MM-PBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) estimates it from an MD trajectory as the difference between three states:

ΔGbind = Gcomplex − Gprotein − Gligand

Each of those free energies is built from the molecular mechanics energy (bonded, van der Waals, and electrostatic terms), a solvation free energy split into a polar part (from the Poisson-Boltzmann equation, or the faster Generalized Born approximation in MM-GBSA) and a nonpolar part from the solvent-accessible surface area, and an optional entropy term (−TΔS) from normal-mode or interaction-entropy analysis. MM-PBSA is an “end-state” method: it only samples the bound and free states, so it is far cheaper than rigorous alchemical methods such as free-energy perturbation, at the cost of being more approximate. As Genheden and Ryde put it in their review, MM/GBSA and MM/PBSA are “among the most popular methods to estimate ligand-binding affinities” (Genheden and Ryde, 2015).

What is gmx_MMPBSA and what do you need to run it?

gmx_MMPBSA is a free, open-source program that performs MM-PBSA and MM-GBSA calculations directly on GROMACS files. It reuses the well-tested MMPBSA.py machinery from AmberTools, so it inherits AMBER’s solvation models while reading GROMACS trajectories and topologies (Valdés-Tresanco et al., 2021). You need three things before you start:

• A finished, equilibrated MD run of the protein-ligand complex: the run input file (.tpr), the trajectory (.xtc), the topology (topol.top) with the ligand parameters included, and an index file (.ndx).
• AmberTools (version 21 or newer), which provides the underlying PB and GB solvers.
• GROMACS and gmx_MMPBSA themselves. The documented install is a conda environment, for example conda install -c conda-forge "ambertools>=21" mpi4py followed by python -m pip install gmx_MMPBSA.

Getting a correct ligand topology upstream is the part that trips people up, because MM-PBSA cannot fix a badly parameterized ligand. If you have not run the MD yet, start from our GROMACS molecular dynamics pillar guide.

MM-PBSA vs MM-GBSA: which should you use?

The honest guidance is that neither is uniformly more accurate; the result depends on the target and the parameters. A common workflow is to screen with the faster MM-GBSA and confirm the top hits with MM-PBSA. You can compute both in a single gmx_MMPBSA run by including both the &gb and &pb sections in the input file.

How do you calculate binding free energy step by step?

Assume you have a completed complex simulation with md.tpr, md.xtc, topol.top, and an index file. The workflow has four steps.

1. Make an index group for the receptor and the ligand. gmx_MMPBSA needs to know which atoms are the protein and which are the ligand:
   gmx make_ndx -f md.tpr -o index.ndx
   Note the group numbers, for example group 1 for the protein and group 13 for the ligand.
2. Write the input file mmpbsa.in. A minimal MM-GBSA setup over 100 frames looks like this:
   
   &general

startframe=1, endframe=100, interval=1,

/

&gb

igb=5, saltcon=0.150,

/
   
   Here igb=5 selects an Onufriev-Bashford-Case Generalized Born model and saltcon=0.150 sets a physiological 0.15 M salt concentration. To add MM-PBSA, append a &pb section.
3. Run gmx_MMPBSA. Point it at the topology, run file, index, group numbers, and trajectory:
   
   gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ci index.ndx -cg 1 13 -ct md.xtc -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv
   
   The -cg 1 13 flag tells it group 1 is the receptor and group 13 is the ligand. The tool extracts the two components, computes each energy term per frame, and averages.
4. Read the output. Open FINAL_RESULTS_MMPBSA.dat for the averaged ΔG with its standard deviation, or load the run interactively in the analysis GUI: gmx_MMPBSA_ana -f _GMXMMPBSA_info.

Subsampling with interval matters: you rarely need every frame, and using every tenth frame over a converged part of the trajectory gives a stable average far faster than processing thousands of frames.

How do you read the results and find the key residues?

The .dat file reports the total binding energy and its components: the van der Waals and electrostatic contributions from molecular mechanics, the polar and nonpolar solvation terms, and the resulting ΔG_bind. The most useful diagnostic is often the sign and size of each term, because it tells you whether binding is driven by shape complementarity (van der Waals) or by specific polar contacts (electrostatics minus polar solvation).

To see which residues matter, turn on per-residue decomposition by adding a &decomp section with idecomp=1 to the input file. gmx_MMPBSA then reports each residue’s contribution, and the analysis GUI plots them so the hotspot residues stand out. That decomposition is exactly the evidence a thesis examiner wants, because it connects a single binding number to the specific interactions you also show in a 2D interaction diagram. For the full path from sequence to a simulated, scored complex, see our computational biology skills roadmap.

How accurate is MM-PBSA, and what are its limits?

Treat MM-PBSA numbers as relative rankings, not absolute affinities. The method is good at ordering a congeneric series of similar ligands and weak at predicting an experimental ΔG on the nose. Three limits are worth stating in any thesis:

• Entropy is usually omitted. The configurational entropy term is expensive and noisy, so many studies leave it out, which means the reported “ΔG” is really an effective enthalpy-like quantity. Say so explicitly.
• Results are sensitive to parameters. The choice of GB model, interior dielectric constant, and salt concentration all shift the number, so keep them fixed across the ligands you compare.
• Single-trajectory bias. The standard single-trajectory protocol reuses the complex trajectory for the free protein and ligand, which cancels some errors but assumes little conformational change on binding.

Used with those caveats, MM-PBSA is a widely accepted way to add a quantitative binding result on top of docking and MD.

Troubleshooting: common gmx_MMPBSA errors

• “GROMACS is not installed” or gmx not found. gmx_MMPBSA calls GROMACS internally, so gmx must be on your PATH in the same environment. Activate the conda environment where both are installed before running.
• The trajectory and tpr do not match. The atom count in the .tpr must correspond to the trajectory you pass. Feed the tool the same complex run and topology you simulated, not a stripped or re-solvated variant.
• Ligand parameters missing from the topology. If topol.top does not include the ligand’s force-field parameters (for example a GAFF topology from ACPYPE), the AMBER conversion fails. Fix the ligand parameterization before running MM-PBSA.
• Non-integer total charge warning. A fractional system charge usually means a partial-charge or protonation error in the ligand. Recheck the ligand charges rather than ignoring the warning.
• The run is extremely slow. PB solves are grid-based and slow per frame. Switch to MM-GBSA for screening, and raise interval so you process fewer frames.

Frequently asked questions

Can GROMACS calculate binding free energy on its own?

Not the MM-PBSA type directly. Native GROMACS supports rigorous alchemical free energy (for example thermodynamic integration), but for the popular end-state MM-PBSA and MM-GBSA estimates you use an add-on such as gmx_MMPBSA, which reads GROMACS files and calls the AmberTools solvers.

Is gmx_MMPBSA free?

Yes. gmx_MMPBSA is open-source software. It does require a working AmberTools installation, which is also free for academic use.

How many frames should I use?

Use frames from a converged, equilibrated portion of the trajectory and subsample with the interval setting. A few tens to a few hundred well-spaced frames usually give a stable average; adding more frames from the same equilibrium region rarely changes the answer.

Should I report MM-GBSA or MM-PBSA in my thesis?

Report whichever you validated against your data, and state the model and parameters. Many students screen with MM-GBSA and confirm the top candidates with MM-PBSA, presenting both.

Why is my binding energy positive?

A positive value suggests the pose is not stably bound, the parameters or charges are wrong, or entropy and solvation are dominating an unfavourable pose. Recheck that the MD complex was stable and that the ligand topology is correct before trusting the number.

Sources: Valdés-Tresanco MS, Valdés-Tresanco ME, Valiente PA, Moreno E. gmx_MMPBSA: a new tool to perform end-state free energy calculations with GROMACS. J Chem Theory Comput. 2021;17(10):6281-6291. doi:10.1021/acs.jctc.1c00645. Miller BR 3rd, McGee TD Jr, Swails JM, Homeyer N, Gohlke H, Roitberg AE. MMPBSA.py: an efficient program for end-state free energy calculations. J Chem Theory Comput. 2012;8(9):3314-3321. doi:10.1021/ct300418h. Genheden S, Ryde U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discov. 2015;10(5):449-461. doi:10.1517/17460441.2015.1032936.