Amber Mmgbsa

End-to-end guide for running MM/GBSA binding free energy calculations using AmberTools/Amber on a pre-existing receptor–ligand complex without molecular dynamics sampling. Covers structure preparation, GAFF2 ligand parameterisation, topology construction, single-frame trajectory generation, MMPBSA.py execution, and result interpretation with explicit uncertainty boundaries. Suitable for crystal-pose scoring, post-docking rapid screening, and mechanistic energy decomposition.

ning-kun@ning-kun

Install

openclaw skills install @ning-kun/amber-mmgbsa

Amber MM/GBSA — Single-Structure Workflow

Overview

This skill describes a complete, production-ready workflow for computing MM/GBSA (Molecular Mechanics / Generalized Born Surface Area) binding free energies using AmberTools / Amber on a single static structure — no molecular dynamics sampling required.

It is designed for two primary use cases:

Rapid pose scoring — evaluate the binding energy of a crystal structure or a top docking pose without running an MD simulation.
Mechanistic decomposition — obtain a standard breakdown VDWAALS / EEL / EGB / ESURF / ΔG_bind to interpret which interactions drive or oppose binding.

Scope note: This workflow computes a single-structure enthalpy-like MM/GBSA estimate. It does not include conformational entropy, water-mediated interactions, or ensemble averaging. For quantitative ranking of multiple ligands or for results intended to compare with experimental ΔG values, a full MD trajectory + MMPBSA.py ensemble average is strongly recommended.

Prerequisites

Software

Tool	Version	Purpose
`tleap`	AmberTools 20+	System building & topology
`antechamber`	AmberTools 20+	Ligand parameterisation
`parmchk2`	AmberTools 20+	Missing GAFF2 parameters
`cpptraj`	AmberTools 20+	Trajectory extraction
`MMPBSA.py`	AmberTools 20+	MM/GBSA calculation
`parmed`	any recent	Topology inspection & merging

Set your environment:

bash

export AMBERHOME=/path/to/amber24   # adjust to your installation
export PATH=$AMBERHOME/bin:$PATH
PY=$AMBERHOME/miniconda/bin/python  # preferred Python

Input files required

File	Description
`receptor.pdb`	Receptor (protein ± DNA ± metals) with correct residue names
`ligand.pdb`	Ligand with 3D coordinates, correct bond orders, and known net charge

Recommended Directory Layout

text

project_mmgbsa/
├── input/
│   ├── receptor.pdb
│   └── ligand.pdb
├── prep/          # cleaned structures
├── build/         # topologies and coordinates
├── traj/          # single-frame trajectories
├── mmgbsa/        # MMPBSA.py inputs and outputs
└── logs/          # all log files

Step-by-Step Protocol

Step 0 — Verify Input Structures

Before any calculation, confirm:

Receptor: Standard Amber residue names (protein: standard 20 aa; DNA: DC/DA/DG/DT; water: WAT/HOH). Remove缓冲 salts and small molecules with no parameters.
Ligand: 3D coordinates present; bond orders chemically reasonable; net charge known.
Complex geometry: Receptor and ligand must already be in the correct binding pose (crystal structure or top-ranked docking pose). They cannot be randomly docked in separate files without a reference complex.

Step 1 — Clean the Receptor PDB

Rename crystallographic water to Amber-compatible WAT:

bash

mkdir -p prep logs
awk '
BEGIN{OFS=""}
/^ATOM|^HETATM/ {
  res=substr($0,18,3)
  if (res=="HOH") {
    print substr($0,1,17),"WAT",substr($0,21)
  } else {
    print $0
  }
  next
}
{print}
' input/receptor.pdb > prep/receptor_clean.pdb

Should I keep crystallographic waters? Usually remove them for MM/GBSA unless a specific water is structurally validated as a key bridging molecule. If retained, document this decision and apply it consistently across all compared systems.

Step 2 — Parameterise the Ligand (GAFF2)

2a. Generate GAFF2 mol2 with AM1-BCC charges

bash

cd prep
antechamber \
  -i ../input/ligand.pdb \
  -fi pdb \
  -o ligand.mol2 \
  -fo mol2 \
  -c bcc \
  -nc 0 \
  -at gaff2 \
  -j 4 \
  > ../logs/antechamber.log 2>&1

# Verify mol2 was produced
ls -lh ligand.mol2

Critical: The -nc value must match the true formal charge of the ligand. Common examples:

Neutral organic molecule → -nc 0
Carboxylate → -nc -1
Phosphate → -nc -2 or -3
Metal complex → confirm charge separately

2b. Check for missing force-field parameters

bash

parmchk2 \
  -i ligand.mol2 \
  -f mol2 \
  -o ligand.frcmod \
  > ../logs/parmchk2.log 2>&1

# If "frcmod" is empty or missing params are critical, either:
#   a) re-draw the problematic substructure in the PDB
#   b) manually add parameters to the frcmod file

Step 3 — Build the Receptor Topology

Protein only

bash

mkdir -p build
cat > build/tleap_rec.in <<'EOF'
source leaprc.protein.ff14SB
source leaprc.water.tip3p

rec = loadpdb ../prep/receptor_clean.pdb
desc rec
saveamberparm rec receptor.prmtop receptor.inpcrd
quit
EOF

tleap -f build/tleap_rec.in > logs/tleap_rec.log 2>&1

Protein + DNA

bash

cat > build/tleap_rec.in <<'EOF'
source leaprc.protein.ff14SB
source leaprc.DNA.OL15
source leaprc.water.tip3p

rec = loadpdb ../prep/receptor_clean.pdb
desc rec
saveamberparm rec receptor.prmtop receptor.inpcrd
quit
EOF

tleap -f build/tleap_rec.in > logs/tleap_rec.log 2>&1

Protein + DNA + metal ions or special groups

bash

cat > build/tleap_rec.in <<'EOF'
source leaprc.protein.ff14SB
source leaprc.DNA.OL15
loadoff $AMBERHOME/dat/leap/lib/terminal_monophosphate.lib
source leaprc.water.tip3p

rec = loadpdb ../prep/receptor_clean.pdb
desc rec
saveamberparm rec receptor.prmtop receptor.inpcrd
quit
EOF

tleap -f build/tleap_rec.in > logs/tleap_rec.log 2>&1

Verify success:

bash

ls -lh build/receptor.prmtop build/receptor.inpcrd
# receptor.prmtop should be > 10 KB for a realistic protein

Step 4 — Build the Ligand Topology

bash

cat > build/tleap_lig.in <<'EOF'
source leaprc.gaff2
lig = loadmol2 ../prep/ligand.mol2
loadAmberParams ../prep/ligand.frcmod
desc lig
saveamberparm lig ligand.prmtop ligand.inpcrd
quit
EOF

tleap -f build/tleap_lig.in > logs/tleap_lig.log 2>&1

# Verify and inspect charge
$PY -c "import parmed as p; r=p.load_file('build/ligand.prmtop'); print(f'Ligand charge: {r.charge:.4f}')"

Step 5 — Build the Complex Topology (Consistent Atom Types)

The key rule: receptor and ligand topologies must come from the same tleap session to ensure atom type definitions are mutually consistent.

bash

cat > build/tleap_complex.in <<'EOF'
source leaprc.protein.ff14SB
source leaprc.DNA.OL15
source leaprc.gaff2
source leaprc.water.tip3p

lig = loadmol2 ../prep/ligand.mol2
loadAmberParams ../prep/ligand.frcmod

rec = loadpdb ../prep/receptor_clean.pdb

desc rec
desc lig

# Save all three topologies from the same session
saveamberparm rec receptor_final.prmtop receptor_final.inpcrd
saveamberparm lig ligand_final.prmtop ligand_final.inpcrd
quit
EOF

tleap -f build/tleap_complex.in > logs/tleap_complex.log 2>&1

ls -lh build/receptor_final.prmtop build/ligand_final.prmtop

Optional — ParmEd merge (if you need a single complex topology):

python

# build/merge_complex.py
import parmed as p, os

rec = p.load_file('build/receptor_final.prmtop', 'build/receptor_final.inpcrd')
lig = p.load_file('build/ligand_final.prmtop', 'build/ligand_final.inpcrd')

complex_sys = rec + lig
complex_sys.save('build/complex.prmtop', overwrite=True)
complex_sys.save('build/complex.inpcrd', overwrite=True)

print(f"Receptor:   {rec.atoms} atoms,  charge={rec.charge:.4f}")
print(f"Ligand:     {lig.atoms} atoms,  charge={lig.charge:.4f}")
print(f"Complex:    {complex_sys.atoms} atoms, charge={complex_sys.charge:.4f}")

Step 6 — Generate Single-Frame NetCDF Trajectories

MMPBSA.py requires trajectory inputs. For a single-structure calculation, convert the inpcrd restart file into a single-frame NetCDF:

bash

mkdir -p traj logs

# Receptor
cpptraj <<'EOF' > logs/cpptraj_rec.log 2>&1
parm ../build/receptor.prmtop
trajin ../build/receptor.inpcrd
trajout receptor_traj.nc netcdf
run
EOF

# Ligand
cpptraj <<'EOF' > logs/cpptraj_lig.log 2>&1
parm ../build/ligand.prmtop
trajin ../build/ligand.inpcrd
trajout ligand_traj.nc netcdf
run
EOF

# Complex
cpptraj <<'EOF' > logs/cpptraj_com.log 2>&1
parm ../build/complex.prmtop
trajin ../build/complex.inpcrd
trajout complex_traj.nc netcdf
run
EOF

ls -lh traj/*.nc

If the ligand and receptor topologies share a consistent type table (from Step 5), they can also share a single complex_traj.nc in single-trajectory mode.

Step 7 — Write MMPBSA.py Input Files

7a. Standard MM/GBSA (total energy only)

Save as mmgbsa/mmpbsa.in:

bash

&general
  startframe=1,
  endframe=1,
  interval=1,
  verbose=1,
  keep_files=0,
/
&gb
  igb=5,
  saltcon=0.150,
  surften=0.005,
  surfoff=0.000,
  molsurf=0,
  radiopt=1,
/

7b. MM/GBSA with per-residue decomposition

Save as mmgbsa/mmpbsa_decomp.in:

bash

&general
  startframe=1,
  endframe=1,
  interval=1,
  verbose=1,
  keep_files=0,
/
&gb
  igb=5,
  saltcon=0.150,
  surften=0.005,
  surfoff=0.000,
  molsurf=0,
  radiopt=1,
/
&decomp
  idecomp=1,
  dec_verbose=0,
  print_res="within 6",
/

Key parameter reference:

Parameter	Recommended value	Notes
`igb`	`5`	GB-neck2 (recommended for proteins)
`saltcon`	`0.150`	physiological ionic strength (M)
`surften`	`0.005`	standard surface tension term
`molsurf`	`0`	use LCPO surface area
`radiopt`	`1`	radii from prmtop (matching force field)
`print_res`	`"within 6"`	only residues within 6 Å of ligand

Step 8 — Run MMPBSA.py

bash

mkdir -p mmgbsa logs

MMPBSA.py \
  -O \
  -i mmgbsa/mmpbsa.in \
  -cp build/complex.prmtop \
  -rp build/receptor.prmtop \
  -lp build/ligand.prmtop \
  -y traj/complex_traj.nc \
  -o mmgbsa/FINAL_RESULTS_MMPBSA.dat \
  > logs/mmpbsa.log 2>&1

With decomposition:

bash

MMPBSA.py \
  -O \
  -i mmgbsa/mmpbsa_decomp.in \
  -cp build/complex.prmtop \
  -rp build/receptor.prmtop \
  -lp build/ligand.prmtop \
  -y traj/complex_traj.nc \
  -o mmgbsa/FINAL_RESULTS_MMPBSA.dat \
  -do mmgbsa/FINAL_DECOMP_MMPBSA.dat \
  > logs/mmpbsa_decomp.log 2>&1

Interpreting the Results

Standard Energy Terms

Term	Physical meaning	Sign convention
`VDWAALS`	van der Waals / hydrophobic interactions	Negative = favorable
`EEL`	Gas-phase electrostatic energy	Negative usually favorable; large positive often indicates charge–charge repulsion
`EGB`	Polar solvation free energy (GB)	Negative = favorable (desolvation gain)
`ESURF`	Nonpolar solvation (surface area term)	Negative = favorable
`ΔG_bind`	Net binding free energy	Negative = favorable

The net binding energy is:

text

ΔG_bind = VDWAALS + EEL + EGB + ESURF

When Large Cancellation Occurs

In systems with highly charged components (e.g., DNA-binding proteins, phosphate-containing ligands, multi-metal active sites), you may observe:

text

EEL  = +700   (large unfavorable)
EGB  = −660   (large favorable, compensating)
ΔG_bind ≈ −0.7  (small net, from large cancellation)

This pattern is arithmetically valid but physically fragile. The final ΔG_bind value in such cases is extremely sensitive to:

Small errors in atomic charges
Conformational sampling
The GB solvent model approximation
Whether metals/waters are included

Correct interpretation: Report that the net binding is weakly favorable or inconclusive, and note that electrostatic terms largely cancel. Do not treat a near-zero ΔG_bind from large cancellation as evidence of "no binding" or as a precise quantitative affinity.

Recommended Results Table for Reports

Present results as:

Energy Component	Value (kcal/mol)	Interpretation
VDWAALS	−36.99	Favorable: hydrophobic / shape complementarity
EEL	+702.60	Unfavorable: charge–charge repulsion (both partners negatively charged)
EGB	−662.23	Favorable: polar desolvation gain upon binding
ESURF	−4.14	Favorable: nonpolar surface burial
ΔG_bind	−0.76	Net weakly favorable; dominated by electrostatic cancellation

Follow with a method note:

"ΔG_bind was calculated using the MM/GBSA method (igb=5, saltcon=0.15 M) on a single crystal/pose structure without conformational sampling or entropy correction. Absolute values are estimates and should not be directly compared to experimental binding free energies without further validation."

Common Errors and Troubleshooting

Error 1 — `Unknown residue` in tleap log

Cause: Non-standard residue name in PDB.
Fix: Remap to Amber-compatible names using pdb4amber or manually edit the PDB.

Error 2 — `antechamber` fails to produce mol2

Cause: Malformed PDB (missing bonds), wrong net charge, or unsupported elements.
Fix:

bash

# Check and fix the PDB
babel -ipdb ligand.pdb -osmi   # view SMILES
# Re-draw or re-fetch the ligand structure

Error 3 — Atom count mismatch between topologies and trajectories

Cause: Running tleap separately for receptor and ligand without a shared session.
Fix: Rebuild both from a single tleap_complex.in (Step 5).

Error 4 — Extreme EEL / EGB values (hundreds of kcal/mol)

Cause: System contains exposed charged groups, metals, or phosphate groups; GB model amplifies the gas-phase/solvation contrast.
Fix: Confirm this is expected for your system. Do not round/force these values. Report them honestly and flag the result as "large-cancellation regime."

Error 5 — `MMPBSA.py` crashes with `Segmentation fault`

Cause: Usually insufficient memory or a mismatch between trajectory frame count and topology.
Fix:

bash

# Check trajectory frame count
cpptraj -p build/complex.prmtop -y traj/complex_traj.nc -e 1 <<'EOF'
EOF
# Set startframe/endframe explicitly in mmpbsa.in

Minimum Deliverables for a Publication-Ready Result

When presenting MM/GBSA results in a paper or report, include:

Energy table — all terms with units (kcal/mol)
Method paragraph — force field, GB model, radii set, salt concentration, single vs. multi-frame
System description — protein + DNA ± metals ± waters, net charges, ligand type
Explicit limitation statement — single structure, no entropy correction, GB model limitations
Source data — enough to reproduce: input PDBs, topologies, mmpbsa.in, and raw output log

Quick-Reference Command Summary

bash

# 1. Ligand parameterisation
antechamber -i ligand.pdb -fi pdb -o ligand.mol2 -fo mol2 -c bcc -nc CHARGE -at gaff2 -j 4
parmchk2 -i ligand.mol2 -f mol2 -o ligand.frcmod

# 2. Receptor topology
tleap -f tleap_rec.in

# 3. Ligand topology
tleap -f tleap_lig.in

# 4. Single-frame trajectories
cpptraj -p complex.prmtop -y complex.inpcrd -o complex_traj.nc netcdf

# 5. MM/GBSA
MMPBSA.py -i mmpbsa.in -cp complex.prmtop -rp receptor.prmtop -lp ligand.prmtop -y complex_traj.nc -o FINAL_RESULTS.dat

When to Upgrade to MD-Based MM/GBSA

If your single-structure result is close to zero, involves large electrostatic cancellation, or will be used to rank multiple ligands, upgrade to the full MD-based workflow:

Run at least 100 ns of production MD per system
Extract 100–200 snapshots (every 0.5–1 ns)
Use three-trajectory mode (complex / receptor / ligand trajectories)
Compute block averages to estimate statistical uncertainty
Add normal mode analysis or quasi-harmonic entropy if feasible
Run at least 2–3 independent replicates per system

One-Sentence Takeaway

Amber single-structure MM/GBSA is a fast, interpretable binding energy estimator; when electrostatic terms dominate and cancel, treat the absolute ΔG_bind as a directional indicator rather than a quantitative affinity — and always upgrade to MD ensemble sampling for comparative rankings.

Amber Mmgbsa

Install

Amber MM/GBSA — Single-Structure Workflow

Overview

Prerequisites

Software

Input files required

Recommended Directory Layout

Step-by-Step Protocol

Step 0 — Verify Input Structures

Step 1 — Clean the Receptor PDB

Step 2 — Parameterise the Ligand (GAFF2)

2a. Generate GAFF2 mol2 with AM1-BCC charges

2b. Check for missing force-field parameters

Step 3 — Build the Receptor Topology

Protein only

Protein + DNA

Protein + DNA + metal ions or special groups

Step 4 — Build the Ligand Topology

Step 5 — Build the Complex Topology (Consistent Atom Types)

Step 6 — Generate Single-Frame NetCDF Trajectories

Step 7 — Write MMPBSA.py Input Files

7a. Standard MM/GBSA (total energy only)

7b. MM/GBSA with per-residue decomposition

Step 8 — Run MMPBSA.py

Interpreting the Results

Standard Energy Terms

When Large Cancellation Occurs

Recommended Results Table for Reports

Common Errors and Troubleshooting

Error 1 — Unknown residue in tleap log

Error 2 — antechamber fails to produce mol2

Error 3 — Atom count mismatch between topologies and trajectories

Error 4 — Extreme EEL / EGB values (hundreds of kcal/mol)

Error 5 — MMPBSA.py crashes with Segmentation fault

Minimum Deliverables for a Publication-Ready Result

Quick-Reference Command Summary

When to Upgrade to MD-Based MM/GBSA

One-Sentence Takeaway

Related skills

Error 1 — `Unknown residue` in tleap log

Error 2 — `antechamber` fails to produce mol2

Error 5 — `MMPBSA.py` crashes with `Segmentation fault`