MM-GBSA and MM-PBSA with gmx_MMPBSA 1.7: a tested workflow and how to read the numbers

MM-PBSA and MM-GBSA are popular because they are cheap. Take an MD trajectory of a complex, compute molecular-mechanics energies plus an implicit-solvent estimate of solvation for complex, receptor and ligand, subtract, and average. Unfortunately they are also among the most misread numbers in computational drug discovery. This tutorial runs gmx_MMPBSA on the HIV-1 protease–indinavir trajectory from the protein–ligand MD tutorial, and it shows, with real output, how much the answer depends on choices that are often left at their defaults.

Every command and number below comes from a run on 23 October 2026 with gmx_MMPBSA 1.7.0, AmberTools 26.0 and GROMACS 2026.3, on an 8-thread CPU. The input is 50 frames (10–500 ps, every 10 ps) of a short demonstration trajectory. That is too little sampling for a real study, but it is enough to show how the method behaves.

What MM/PBSA actually computes

For each frame, the method evaluates

where each G is

Three approximations matter when you read the result. (1) In the default single-trajectory protocol, receptor and ligand are cut out of the complex frames. Bonded terms cancel exactly, and strain and reorganisation are ignored. (2) The explicit water is replaced by a continuum with an internal dielectric ε_in that you choose. (3) Entropy is usually omitted, so ΔTOTAL is closer to an enthalpy-like score than to a free energy.

Because of these approximations, the result is best used to rank related ligands on the same target, or to find which residues contribute to binding. It should not be taken as an absolute binding free energy. For reviews, see Genheden and Ryde (2015) and Wang et al. (2019).

0Install

gmx_MMPBSA needs AmberTools (for sander and cpptraj), a gmx executable on the PATH, and Python. gmx_MMPBSA 1.7.0 requires Python 3.11 or 3.12, and it is installed from PyPI:

micromamba create -n mmpbsa -c conda-forge python=3.12 ambertools=26 mpi4py openmpi
micromamba activate mmpbsa
python -m pip install gmx_MMPBSA==1.7.0
gmx_MMPBSA -v      # gmx_MMPBSA 1.7.0 based on MMPBSA version 14.0, AmberTools 26.0 ...

Watch out for the conda-forge package. In October 2026, conda install gmx_mmpbsa resolved to the old 1.5.0.3, alongside NumPy 2.5. That combination runs plain GB, but it crashes as soon as you request interaction entropy:

AttributeError: module 'numpy' has no attribute 'float'.

(np.float was removed in NumPy 1.24.) Defaults also changed between these versions: 1.7.0 uses igb = 8 with PBRadii = 4 (mbondi3) and inp = 1, whereas 1.5.0.3 used igb = 5, PBRadii = 3 and inp = 2. Always set these explicitly and state the version in your methods. Run gmx_MMPBSA --create_input gb pb decomp to see the defaults of the version you have installed.

1Inputs from the MD run

You need four things from GROMACS:

Passing -cp matters. Without it, gmx_MMPBSA rebuilds an Amber topology from the PDB using the forcefields setting, and that is not necessarily the force field you simulated with. With -cp, the energies are computed with the same parameters as the MD (here AMBER99SB-ILDN plus GAFF2/AM1-BCC from ACPYPE).

The index file from the MD tutorial already contains group 1 (Protein) and group 13 (MK1). The trajectory still contains water and ions, and gmx_MMPBSA strips them itself (solvated_trajectory = 1, the default). We used the raw md.xtc. Check the receptor and ligand masks printed at the top of the output (:1-198 and :199 here) to confirm the split is what you intended.

2A first MM-GBSA input

mmgbsa.in:

&general
  sys_name = "1HSG_MK1", startframe = 2, endframe = 51, interval = 1,
  temperature = 300, PBRadii = 3,
/
&gb
  igb = 5, saltcon = 0.150,
/

• temperature = 300 matches the simulation. The default is 298.15 K, and the value enters the entropy terms.
• igb = 5 is the OBC-II GB model, and PBRadii = 3 (mbondi2) is the radius set it was parameterised with. For igb = 8 (GBn2), use PBRadii = 4 (mbondi3).
• saltcon = 0.150 matches the 0.15 M NaCl used in the MD.
• startframe = 2 skips frame 1 (t = 0), which is just the restart structure.

Run it:

gmx_MMPBSA -O -i mmgbsa.in -cs md.tpr -ci index.ndx -cg 1 13 \
           -ct md.xtc -cp topol.top -o FINAL_RESULTS_MMGBSA.dat -nogui

-nogui skips the results viewer (gmx_MMPBSA_ana) on a headless machine. 50 frames took about 70 s.

3Reading the output

The key block of FINAL_RESULTS_MMGBSA.dat is the Delta (Complex - Receptor - Ligand) table (kcal/mol, 50 frames):

The bonded terms (ΔBOND, ΔANGLE, ΔDIHED) are exactly zero in the single-trajectory protocol, as expected. Electrostatics illustrates the usual pattern: a large favourable gas-phase term is largely cancelled by the desolvation penalty, so the net polar contribution is about +20 kcal/mol.

Now the comparison everyone wants to make. The experimental K_i of indinavir for HIV-1 protease is 0.52 nM (Dorsey et al., 1994). The Kd → ΔG calculator converts that to −12.7 kcal/mol at 300 K. The MM-GBSA value is more than four times larger. That is normal, and it is not something to fix by tuning parameters. ΔTOTAL has no conformational entropy, no loss of translational and rotational freedom, and no reorganisation energy. It also comes from an implicit-solvent model whose absolute scale is not calibrated. Use it to compare ligands, not to predict K_i.

4Same frames, different solvent model: GB vs PB

PB solves the Poisson–Boltzmann equation on a grid rather than using the GB approximation. It is slower, and it is often assumed to be “more accurate”. In practice the two can differ by more than the effect you care about. mmpbsa.in:

&general
  sys_name = "1HSG_MK1_PB", startframe = 2, endframe = 51, interval = 5,
  temperature = 300, PBRadii = 3,
/
&pb
  istrng = 0.150, fillratio = 4.0, radiopt = 0, inp = 1,
/

inp selects the non-polar model, and it is the most consequential line in the file. inp = 1 uses a surface-area term only. inp = 2 (the default in 1.5.x, but not in 1.7.0) splits non-polar solvation into a cavity term (ENPOLAR) and a dispersion term (EDISPER). radiopt = 0 keeps the radii from the topology (mbondi2) rather than PBSA’s optimised set, so GB and PB see the same atoms.

We ran both non-polar options, and two GB models, on the same trajectory (PB on every fifth frame, which is 10 frames, at 3–5 minutes per run):

The gas-phase terms are identical in every row: ΔVDWAALS ≈ −68 and ΔEEL ≈ −44 to −45 kcal/mol. The spread of about 75 kcal/mol comes entirely from the solvation model. One of the four settings even predicts that a sub-nanomolar inhibitor does not bind. None of them is wrong in a software sense. They are different approximations with different absolute offsets. This is why MM/PBSA numbers from different papers, or from different defaults, cannot be compared. It is also why “PB is more accurate than GB” should be tested on your own system, not assumed.

5The internal dielectric constant

Protein–ligand interfaces with charged groups are very sensitive to ε_in. Raising it (2 or 4 is common in the literature) damps both the gas-phase electrostatics and the desolvation penalty. Indinavir is neutral in this setup, but it makes several hydrogen bonds to the catalytic aspartates and the flap water region. Rerunning the igb = 5 calculation with only intdiel = 4.0 changed in &gb:

Both electrostatic terms shrink about four-fold. The net polar term falls from +19.8 to +4.3 kcal/mol, so this ligand looks better bound at higher ε_in. With a different balance of terms, for example a charged ligand, the effect can go the other way.

None of these values is “right”. If you compare a congeneric series, use the same ε_in for all of them, and report it. Studies that tune ε_in to maximise correlation with experiment on a small set are overfitting.

6Per-residue decomposition

Decomposition is often the most useful MM/GBSA output. It tells you which residues the model thinks the ligand depends on. Add a namelist to the GB input:

&decomposition
  idecomp = 2, dec_verbose = 0, print_res = "within 4",
/

idecomp = 2 adds the 1–4 terms to the internal energy of each residue (per-residue, not pairwise). print_res = "within 4" limits the output to residues within 4 Å of the ligand. The namelist is called &decomposition. Some older inputs you will find online call it &decomp, so check the name with gmx_MMPBSA --create_input decomp for your version.

The 50-frame decomposition took about 8 minutes, compared with about 70 s for plain GB. Per-residue totals (kcal/mol) for the most favourable and least favourable residues, from the decomposition file (named with -do; written as CSV by default):

In per-residue mode, each residue gets half of every pairwise interaction, and the ligand gets the other half (−32.7 kcal/mol here). The two flap isoleucines and the 27–28 and 48–50 backbone come out on top, which matches what is known about how this inhibitor class binds.

The catalytic aspartates come out unfavourable. That is a modelling artefact worth recognising. In the MD tutorial both Asp25 residues were left deprotonated, which puts two negative charges right next to indinavir’s central hydroxyl. Their large desolvation penalties outweigh the gas-phase attraction. With one aspartate protonated, as most protease–inhibitor studies do, this picture would change. A decomposition is only as good as the protonation states behind it.

7Entropy: interaction entropy and C2

gmx_MMPBSA offers two cheap entropy estimates that reuse the energies already computed. Add them to &general:

  interaction_entropy = 1, c2_entropy = 1,

Both work from the fluctuations of the protein–ligand interaction energy over the frames, ΔE_int = ΔE_vdW + ΔE_elec. Interaction entropy (IE; Duan et al., 2016) takes an exponential average, −TΔS = RT ln⟨exp(βδE_int)⟩. The C2 method is the second-order cumulant of that average, −TΔS = σ²_int/(2RT). C2 depends only on the variance, so it is less sensitive to single outlier frames.

Ekberg and Ryde (2021) showed that the exponential average is extremely poorly conditioned. When the standard deviation of the interaction energy (σ_IE) is above about 15 kJ/mol (3.6 kcal/mol), IE practically cannot be converged. Above about 25 kJ/mol (6.0 kcal/mol), C2 gives unrealistically large entropies. Check σ_IE before you report either number.

On this trajectory, σ_IE was 9.72 kcal/mol, well above both limits. gmx_MMPBSA 1.7.0 prints warnings to that effect, and links to the Ekberg and Ryde paper:

WARNING: THE INTERACTION ENERGY STANDARD DEVIATION [ σ(Int. Energy) ]
IS GREATER THAN 3.6 kcal/mol (~15 kJ/mol). THUS, THE INTERACTION ENTROPY VALUES ARE
NOT RELIABLE.

These numbers show the problem the warnings describe. The two “cheap” estimates differ by 65 kcal/mol, and the C2 value alone is larger than the whole binding energy. For a single complex, reporting ΔTOTAL without entropy and saying so is more honest than adding an unconverged −TΔS. When comparing similar ligands, the entropy is often assumed to cancel. That is an assumption too, so state it.

8Convergence and uncertainty

• Look at the per-frame values (-eo frames.csv, or gmx_MMPBSA_ana), not only the mean.
• Do not over-trust the SEM. Consecutive 10 ps frames are correlated. gmx_MMPBSA 1.7.0 reports block-based SEMs for that reason, but independent replicas (with different starting velocities) remain the honest measure of uncertainty.
• Check that the pose stays put (ligand RMSD after fitting on the protein). An average over two binding modes means nothing.

9Single vs multiple trajectories

Everything above uses the single-trajectory protocol, which ignores how the receptor and ligand reorganise on binding. The multiple-trajectory protocol uses separate apo-receptor and free-ligand runs (-rs/-rt/-rg, -ls/-lt/-lg) to capture that reorganisation. You need it for systems such as HIV-1 protease, whose flaps close over the inhibitor. The price is far more noise, because the internal energies no longer cancel.

10Reporting checklist

Report the software versions; the protocol, number of frames, simulation length and number of replicas; the GB model and radii, or the PB inp/radiopt/istrng settings; ε_in, salt and temperature; and how entropy was treated. State that the values are relative unless they have been calibrated.

References

1. 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:6281–6291. doi:10.1021/acs.jctc.1c00645
2. Miller BR III, 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:3314–3321. doi:10.1021/ct300418h
3. Genheden S, Ryde U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discov 2015;10:449–461. doi:10.1517/17460441.2015.1032936
4. Wang E, Sun H, Wang J, et al. End-Point Binding Free Energy Calculation with MM/PBSA and MM/GBSA: Strategies and Applications in Drug Design. Chem Rev 2019;119:9478–9508. doi:10.1021/acs.chemrev.9b00055
5. Onufriev A, Bashford D, Case DA. Exploring protein native states and large-scale conformational changes with a modified generalized Born model. Proteins 2004;55:383–394. doi:10.1002/prot.20033
6. Duan L, Liu X, Zhang JZH. Interaction Entropy: A New Paradigm for Highly Efficient and Reliable Computation of Protein–Ligand Binding Free Energy. J Am Chem Soc 2016;138:5722–5728. doi:10.1021/jacs.6b02682
7. Ekberg V, Ryde U. On the Use of Interaction Entropy and Related Methods to Estimate Binding Entropies. J Chem Theory Comput 2021;17:5379–5391. doi:10.1021/acs.jctc.1c00374
8. Dorsey BD, Levin RB, McDaniel SL, et al. L-735,524: The Design of a Potent and Orally Bioavailable HIV Protease Inhibitor. J Med Chem 1994;37:3443–3451. doi:10.1021/jm00047a001