MM-GBSA Tutorial: How to Calculate Binding Free Energy After Molecular Dynamics

What MM-GBSA is

MM-GBSA stands for Molecular Mechanics Generalized Born Surface Area. It is a post-processing method for estimating the binding free energy of a protein-ligand complex from a molecular dynamics trajectory. The name describes exactly what it combines: Molecular Mechanics (MM) energy terms from the force field, plus a Generalized Born (GB) implicit solvation model, plus a Surface Area (SA) term for the non-polar contribution to solvation.

In plain terms: you run a protein-ligand MD simulation, collect a set of snapshots from the trajectory, and for each snapshot calculate the energy of the complex, the energy of the protein alone, and the energy of the ligand alone. The difference gives you the binding free energy — how much energy is released when the ligand binds. Average this across hundreds of snapshots and you get a much more reliable estimate than any single-pose docking score.

How it differs from docking scores

Understanding why MM-GBSA is better than docking scores — and where it’s still limited — requires understanding what docking scores actually calculate. A docking score is a rapid approximation of binding free energy evaluated at a single static pose using an empirical or knowledge-based scoring function. It takes seconds to compute and ignores protein flexibility, entropic effects, and explicit solvation.

MM-GBSA improves on docking in three specific ways. It averages over many trajectory frames rather than evaluating one static pose, so it captures conformational flexibility. It uses the actual force field to compute energy terms rather than an empirical approximation. And it includes an implicit solvation model that accounts for the desolvation penalty of burying the ligand in the binding site.

In the hierarchy of binding free energy methods, MM-GBSA sits between raw docking scores and rigorous alchemical free energy methods like FEP:

• 1
  Docking score (AutoDock Vina, Glide)

  Single pose, empirical scoring. Seconds per compound. Useful for large-scale screening.
  Very fast · Low accuracy
• 2
  MM-GBSA / MM-PBSA

  Ensemble average over MD frames. Minutes to hours. Standard for hit validation and compound ranking.
  Hours · Better accuracy
• 3
  FEP / Thermodynamic integration

  Rigorous alchemical perturbation. Days per compound pair. Gold standard for lead optimization.
  Days · High accuracy
• 4
  Experimental (ITC, SPR, fluorescence)

  Direct measurement of ΔG, Kd, or IC50. Ground truth.
  Weeks · Ground truth

The energy equation explained

In practice, most MM-GBSA calculations omit the −TΔS term because entropy estimation with normal mode analysis is noisy and computationally expensive. The result without entropy is sometimes called ΔH_bind or “MM-GBSA score” rather than a true ΔG. This is acceptable for relative ranking of compounds against the same target — but it’s why MM-GBSA values cannot be directly compared to experimental Kd values without correction.

How to run MM-GBSA with gmx_MMPBSA

The standard tool for GROMACS-based MM-GBSA calculation is gmx_MMPBSA — an open-source Python package that interfaces GROMACS trajectories with the AMBER MMPBSA.py engine.

Installation

conda activate md
conda install -c conda-forge gmx_mmpbsa -y
# Verify installation
gmx_MMPBSA --version

Prepare input files

You need four files from your protein-ligand MD simulation:

• md.tpr — the GROMACS run input file (contains topology)
• md_centered.xtc — the processed trajectory
• md.tpr — used again as receptor/ligand topology source
• An index file defining protein, ligand, and complex groups

Create or verify the index groups needed:

gmx make_ndx -f md.tpr -o index.ndx
# Verify these groups exist in the index:
# Protein — all protein atoms
# IBP     — your ligand atoms
# Protein_IBP — the complex (protein + ligand together)

Create the MM-GBSA input file

&general
  startframe   = 1,      ! first frame to analyze
  endframe     = 1000,   ! last frame (adjust to your trajectory length)
  interval     = 5,      ! analyze every 5th frame — balance speed vs precision
  verbose      = 2,
  forcefields  = "oldff/leaprc.ff99SBildn",
/
&gb
  igb          = 5,      ! GB model 5 (OBC2) — best for drug-like molecules
  saltcon      = 0.15,   ! physiological salt concentration
/

Run gmx_MMPBSA

gmx_MMPBSA \
  -O \
  -i mmpbsa.in \
  -cs md.tpr \
  -ct md_centered.xtc \
  -ci index.ndx \
  -cg 1 13 \
  -cp topol.top \
  -o mmgbsa_results.dat \
  -do mmgbsa_decomp.dat \
  -nogui

# -cg 1 13 specifies group indices for receptor (1) and ligand (13)
# Check your index.ndx to confirm the correct group numbers
# -do enables per-residue energy decomposition

Runtime depends on frame count and system size. For 200 frames of a typical protein-ligand system, expect 10–30 minutes on a modern CPU. Use the interval parameter to control how many frames are analyzed — more frames give lower statistical uncertainty but take longer.

Interpreting the results

The output file contains average values and standard deviations for each energy component. The key number is ΔTOTAL — the overall MM-GBSA binding free energy estimate. More negative values indicate stronger predicted binding.

Beyond ΔTOTAL, the per-component breakdown tells you what’s driving binding. A large negative ΔE_vdW means hydrophobic packing dominates; a large negative ΔE_elec paired with a large positive ΔG_GB means electrostatic interactions drive binding but incur a significant desolvation penalty. This decomposition guides medicinal chemistry — if van der Waals terms dominate, growing into hydrophobic pockets is productive; if electrostatics dominate, adjusting charge state or adding H-bond donors/acceptors will be most impactful.

The per-residue decomposition (from the -do flag) shows which specific residues contribute most to binding. The top-contributing residues are the pharmacophore anchors — the ones where chemical modifications will have the largest effect on affinity.

Limitations you must know

How to report MM-GBSA in a paper

A complete MM-GBSA methods section should include: the tool name and version (gmx_MMPBSA), the GB model used (igb=5, OBC2), salt concentration, the number of frames analyzed and the interval, which portion of the trajectory was used (production phase only, after RMSD convergence), whether entropy was included or excluded and why, and the force field used for energy calculations. Report values as mean ± standard deviation across the analyzed frames.