---
name: amber-mmgbsa
description: |
  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.
metadata:
  {
    "openclaw":
      {
        "emoji": "🧪",
        "examples":
          [
            "run MM/GBSA on a receptor and ligand without MD sampling",
            "single-structure MMPBSA.py workflow for a crystal complex",
            "how to calculate binding energy from receptor.pdb and ligand.pdb with Amber",
            "Amber MM/GBSA tutorial no dynamics"
          ]
      }
  }
---

# 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

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

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

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

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