gmx_MMPBSA

What is gmx_MMPBSA?

gmx_MMPBSA calculates binding free energies from molecular dynamics trajectories using MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) and MM/GBSA (Molecular Mechanics/Generalized Born Surface Area) methods. Developed by Valdés-Tresanco and colleagues, it bridges GROMACS simulations with AMBER's MMPBSA.py engine, enabling end-state free energy calculations for protein-ligand, protein-protein, and protein-nucleic acid complexes.

These methods estimate binding affinity by calculating the free energy difference between bound and unbound states. The calculation combines molecular mechanics energies (van der Waals, electrostatic) with implicit solvation models, offering a balance between the speed of empirical scoring functions and the accuracy of rigorous alchemical methods.

How MM/PBSA and MM/GBSA work

Binding free energy ($\Delta G_{\text{bind}}$) is computed from snapshots extracted from an MD trajectory:

$\Delta G_{\text{bind}} = G_{\text{complex}} - G_{\text{receptor}} - G_{\text{ligand}}$

Each term comprises several energy components:

Component	Description
$\Delta E_{\text{vdW}}$	Van der Waals interactions (Lennard-Jones potential)
$\Delta E_{\text{elec}}$	Coulombic electrostatic interactions
$\Delta G_{\text{polar}}$	Polar solvation energy (PB or GB model)
$\Delta G_{\text{nonpolar}}$	Nonpolar solvation, proportional to solvent-accessible surface area
$-T\Delta S$	Conformational entropy (optional, computationally expensive)

MM/PBSA solves the Poisson-Boltzmann equation for electrostatic solvation, providing higher accuracy but at greater computational cost. MM/GBSA uses analytical Generalized Born approximations, running faster while maintaining reasonable accuracy for ranking compounds.

Generalized Born models

Model	Description
GB-HCT (igb=1)	Hawkins-Cramer-Truhlar model, fast but less accurate
GB-OBC1 (igb=2)	Onufriev-Bashford-Case model, variant 1
GB-OBC2 (igb=5)	OBC variant 2, recommended for most applications
GB-Neck (igb=7)	Improved treatment of interstitial regions
GB-Neck2 (igb=8)	Further refinements to neck region calculations

The OBC2 model (igb=5) provides a good balance of accuracy and speed for protein-ligand systems.

How to use gmx_MMPBSA online

ProteinIQ runs gmx_MMPBSA 1.6.5 with a tested GROMACS 2025.4 and AmberTools 23 environment, eliminating the need to configure these dependencies locally.

Inputs

Input	Description
Trajectory	One to four GROMACS trajectories (.xtc, .trr) or multi-model PDB files. gmx_MMPBSA concatenates the submitted trajectories in order.
Complex structure	Bound complex structure (.tpr recommended, .pdb supported) passed to gmx_MMPBSA as -cs.
Complex topology file	Optional GROMACS topology (.top) passed to gmx_MMPBSA as -cp.
Reference structure	Optional but recommended PDB used for reliable chain assignment, passed as -cr.
Ligand MOL2 file	Required for a protein-small-molecule complex when neither a complex nor ligand topology supplies the ligand parameters. It must be an Antechamber MOL2 file.
Index file	Required GROMACS index file (.ndx). Set the receptor and ligand group IDs in the settings.
Input parameter file	Optional gmx_MMPBSA input file (.in or .txt). Overrides calculation settings with custom parameters. Receptor and ligand group IDs still come from the form.

Settings

Calculation settings

Setting	Description
Calculation method	MM/GBSA (faster, suitable for screening), MM/PBSA (more accurate), or Both for comparison
Energy decomposition	When enabled, calculates per-residue energy contributions to identify key binding site residues

Trajectory settings

Setting	Description
Start frame	First frame to analyze (1-indexed)
End frame	Last frame to analyze. Use 0 to keep the gmx_MMPBSA source default
Frame interval	Analyze every nth frame. Set to 2 or higher to reduce computation for long trajectories

Solvation settings

Setting	Description
Salt concentration	Ionic strength in molar. The default keeps the gmx_MMPBSA source value (0.0 M)
GB model	Generalized Born variant for MM/GBSA. GB-OBC2 (igb=5) recommended

Complex index groups

Setting	Description
Receptor group	Receptor group ID from the complex index file, passed as the first -cg value
Ligand group	Ligand group ID from the complex index file, passed as the second -cg value

Output

Results include a separate total binding free energy and component breakdown for every calculation method requested:

Column	Description
DELTA TOTAL	Net binding free energy (kcal/mol). More negative indicates stronger binding.
VDWAALS	Van der Waals contribution
EEL	Electrostatic energy (gas phase)
EGB or EPB	Polar solvation (GB or PB method)
ESURF or ENPOLAR	Nonpolar solvation (surface area term)

When energy decomposition is enabled, per-residue contributions identify which amino acids drive binding.

The Files tab includes the source result report, per-frame energy CSV, decomposition files when requested, bounded execution logs, the effective input file, and COMPACT_MMXSA_RESULTS.mmxsa for analysis with the gmx_MMPBSA Python API or desktop analyzer. Completed jobs stream these artifacts directly to persistent storage. Callback-free runs use a bounded inline fallback and fail explicitly instead of discarding an oversized result.

The separate run.log records the software versions, energy units, calculation methods, input source and outcome. Full native energy statistics, mutant comparisons and residue-pair decomposition remain available when produced by the calculation.

When you provide a custom input parameter file, ProteinIQ reports the calculation methods detected from the source output rather than claiming that the separate form defaults were used.

Interpreting binding energies

Interpret these estimates in the context of the same molecular system, sampling and calculation settings. A numerical MM/PBSA or MM/GBSA threshold does not by itself establish an experimental affinity range. Absolute values carry significant uncertainty, including contributions omitted from the calculation.

Applications

• Lead optimization: Compare binding energies across congeneric compound series to guide medicinal chemistry
• Binding hotspot identification: Energy decomposition reveals which residues contribute most to binding, informing mutagenesis studies
• Selectivity analysis: Compare binding to on-target versus off-target proteins
• Protein-protein interface analysis: Identify key residues stabilizing macromolecular complexes

Limitations

The single-trajectory approach used here assumes identical conformations for complex, receptor, and ligand. This reduces noise but misses conformational reorganization effects. Results depend heavily on the quality and length of the MD simulation—short simulations with poor sampling produce unreliable energies.

Entropy contributions are not calculated by default due to computational cost and large uncertainties. For charged ligands, electrostatic and polar solvation terms often exhibit significant cancellation, amplifying small errors. Results should guide compound prioritization rather than predict absolute binding constants.