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#

ToolVersionPurposetleapAmberTools 20+System building & topologyantechamberAmberTools 20+Ligand parameterisationparmchk2AmberTools 20+Missing GAFF2 parameterscpptrajAmberTools 20+Trajectory extractionMMPBSA.pyAmberTools 20+MM/GBSA calculationparmedany recentTopology 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#

FileDescriptionreceptor.pdbReceptor (protein ± DNA ± metals) with correct residue namesligand.pdbLigand 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:

ParameterRecommended valueNotesigb``5GB-neck2 (recommended for proteins)saltcon``0.150physiological ionic strength (M)surften``0.005standard surface tension termmolsurf``0use LCPO surface arearadiopt``1radii 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#

TermPhysical meaningSign conventionVDWAALSvan der Waals / hydrophobic interactionsNegative = favorableEELGas-phase electrostatic energyNegative usually favorable; large positive often indicates charge–charge repulsionEGBPolar solvation free energy (GB)Negative = favorable (desolvation gain)ESURFNonpolar solvation (surface area term)Negative = favorableΔG_bindNet binding free energyNegative = 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 ComponentValue (kcal/mol)InterpretationVDWAALS−36.99Favorable: hydrophobic / shape complementarityEEL+702.60Unfavorable: charge–charge repulsion (both partners negatively charged)EGB−662.23Favorable: polar desolvation gain upon bindingESURF−4.14Favorable: nonpolar surface burialΔG_bind****−0.76Net 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:

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#

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.