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:

1. Rapid pose scoring — evaluate the binding energy of a crystal structure or a top docking pose without running an MD simulation.
2. 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

Set your environment:

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

Input files required

Recommended Directory Layout

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:

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

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

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

mkdir -p build
cat > build/tleap_rec.in \x3C\x3C'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

cat > build/tleap_rec.in \x3C\x3C'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

cat > build/tleap_rec.in \x3C\x3C'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:

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

Step 4 — Build the Ligand Topology

cat > build/tleap_lig.in \x3C\x3C'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.

cat > build/tleap_complex.in \x3C\x3C'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):

# 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:

mkdir -p traj logs

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

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

# Complex
cpptraj \x3C\x3C'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:

&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:

&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:

Step 8 — Run MMPBSA.py

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:

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

The net binding energy is:

Δ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:

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:

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:

# 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:

# Check trajectory frame count
cpptraj -p build/complex.prmtop -y traj/complex_traj.nc -e 1 \x3C\x3C'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:

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

Quick-Reference Command Summary

# 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.