drug-mmpbsa-gbsa

Files (26)

• SKILL.md
• README.md · examples/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL1087650/rep1/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL1087650/rep1/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL1087650/rep2/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL1087650/rep2/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL1087650/rep3/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL1087650/rep3/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL388978/rep1/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL388978/rep1/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL388978/rep2/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL388978/rep2/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL388978/rep3/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL388978/rep3/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL3943841/rep1/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL3943841/rep1/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL3943841/rep2/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL3943841/rep2/
• mmgbsa_frames.csv · examples/cdk2-htvs/results/CHEMBL3943841/rep3/
• mmgbsa_summary.json · examples/cdk2-htvs/results/CHEMBL3943841/rep3/
• README.md · examples/cdk2-mmpbsa/
• FINAL_RESULTS_MMPBSA.dat · examples/cdk2-mmpbsa/outputs/
• mmpbsa.in · examples/cdk2-mmpbsa/outputs/
• mmpbsa_summary.json · examples/cdk2-mmpbsa/outputs/
• compute_mmgbsa.py · scripts/
• compute_mmpbsa.py · scripts/

SKILL.md (first lines)

# drug-mmpbsa-gbsa (MM-GBSA / MM-PBSA)

## Goal

To estimate relative binding free energies from MD trajectories using the single-trajectory MM-GBSA / MM-PBSA approach. For each trajectory frame, the method strips explicit solvent, evaluates potential energies of the complex, receptor, and ligand subsystems in implicit solvent, and computes:

    dG = E_complex - E_receptor - E_ligand

The entropy term (-TdS) is omitted, which is standard practice when the goal is relative ranking rather than absolute binding affinity. MM-GBSA / MM-PBSA is most useful for re-ranking docked poses after MD refinement, providing an orthogonal signal to docking scores and geometric stability metrics.

## Choosing a backend

| Path | Script | When to use | Extras |
|---|---|---|---|
| **OpenMM GBn2 (fast)** | `compute_mmgbsa.py` | Throughput rescoring of HTVS hits; everything stays inside OpenMM with the same force field as the MD | No extra dependencies; ~1-5 minutes per compound on CPU |
| **AmberTools MMPBSA.py** | `compute_mmpbsa.py` | When you need PB (not just GB), per-method decomposition (ELE, VDW, EGB / EPB, ESURF), or a setup that matches what reviewers expect from the MM-PBSA literature | Adds `MMPBSA.py`, `cpptraj`, and `parmed` to the dependency surface (`parmed` is in `cpu+openmm`; AmberTools must be installed separately and on `PATH`); ~1-3 minutes for GB, ~5-30 minutes for PB depending on system size and frame count |

Both paths give comparable GB rankings for typical drug-protein systems, but the absolute dG numbers will differ across backends because they use different GB models, radius sets, and surface-area treatments. **Don't compare numbers across the two scripts.**

## Instructions

### 1. Basic usage (protein + ligand, short HTVS-style MD)

For the 1-5 ns production runs typical in the [HTVS workflow](../../.agents/workflows/drug-hit-finding-htvs.md):

```bash
${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --ligand_resname UNL \
  --skip_ns 0.5 \
  --stride 5 \
  --output_dir md/mmgbsa/
```

Key parameters:
- `--topology`: solvated complex PDB from [drug-complex-system-builder](../drug-complex-system-builder/SKILL.md)
- `--trajectory`: production DCD from [drug-protein-ligand-md](../drug-protein-ligand-md/SKILL.md)
- `--ligand_sdf`: ligand SDF with correct bond orders (for OpenFF parameterization and AM1-BCC charges)
- `--charge_method` (`compute_mmgbsa.py`): ligand partial charges, `am1bcc` (default; needs AmberTools' `sqm` on `PATH`, not in any uv environment) or `mmff94` / `gasteiger` (RDKit; cruder, for screening and tests). Use the same method the complex was built with.
- `--ligand_resname`: residue name of the ligand in the PDB topology (default: UNL)
- `--skip_ns`: skip the first N nanoseconds of trajectory as equilibration (default: 0.5). For longer production runs, increase proportionally: roughly 10% of production length is a reasonable rule of thumb, capped at ~5 ns for 50+ ns runs.
- `--stride`: evaluate every Nth frame after skipping (default: 5)
- `--solute_dielectric`: interior dielectric of the solute (default: 1.0). See section 3 below on when to raise this.

### 2. With a cofactor (e.g., NADPH) and longer MD

When the receptor has a bound cofactor that should be part of the receptor subsystem, and the MD production was long enough (>=10 ns) to justify a longer equilibration skip:

```bash
${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --cofactor_sdf md/nadph.sdf \
  --ligand_resname UNL \
  --cofactor_resname NDP \
  --skip_ns 2.0 \
  --output_dir md/mmgbsa/
```

### 3. Tuning the solute dielectric

The interior dielectric of the solute (`--solute_dielectric`, OpenMM's `soluteDielectric` parameter) is a major knob for MM-GBSA rankings and is often the difference between sensible and nonsense results. The default of 1.0 treats the solute interior as electronically non-polarizable; the explicit partial charges and atom geometries are still there, but no implicit electronic polarization is added on top. This is the most physically pure choice and works well for **nonpolar binding sites with few charged residues**. For **polar or highly charged pockets**, raise it:

| Binding site character | Suggested `--solute_dielectric` |
|---|---|
| Mostly hydrophobic, few ionizable residues | 1.0 (default) |
| Mixed polar/nonpolar | 2.0 |
| Many charged residues (Asp/Glu/Lys/Arg clusters, salt bridges) | 4.0 |

Higher dielectrics damp electrostatic contributions and generally improve ranking agreement with experiment for charged systems, at the cost of losing sensitivity to directional electrostatic interactions. If you are unsure, run the rescoring at 1.0 and 2.0 on a small subset with known rank-order and pick the value that tracks better. The chosen value is recorded in `mmgbsa_summary.json` for provenance.

```bash
${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --ligand_resname UNL \
  --solute_dielectric 2.0 \
  --output_dir md/mmgbsa/
```

### 4. Custom selections

For non-standard residue naming or multi-chain receptors:

```bash
${CLAUDE_SKILL_DIR}/../../venv/run cpu+openmm python ${CLAUDE_SKILL_DIR}/scripts/compute_mmgbsa.py \
  --topology md/system/complex_solvated.pdb \
  --trajectory md/run/production.dcd \
  --ligand_sdf md/ligand.sdf \
  --ligand_sel "resname UNK and not name H*" \
  --output_dir md/mmgbsa/
```

### 5. Interpret results

The script outputs:

- `mmgbsa_frames.csv`: per-frame energies (E_complex, E_receptor, E_ligand, dG)
- `mmgbsa_summary.json`: mean, std, SEM of dG, and the dielectric/frame parameters used

```json
{
    "dG_mean_kcal_mol": -42.5,
    "dG_std_kcal_mol": 8.3,
    "dG_sem_kcal_mol": 1.2,
    "n_frames_evaluated": 50,
    "solute_dielectric": 1.0,
    "solvent_dielectric": 78.5
}
```

More negative dG indicates stronger predicted binding. **Only relative ranking within a congeneric series is interpretable.** Absolute dG values from single-trajectory MM-GBSA without entropy corrections are not binding free energies in any physical sense, and specific numbers should not be reported as "predicted affinities" or compared across chemically distinct scaffolds. When ranking compounds:

- Group comparisons by formal charge state (see Constraints below).
- Compare mean dG values, but also inspect the per-frame trace in `mmgbsa_frames.csv` for stability. A compound with a noisy or drifting dG signal is less trustworthy than one with a tight distribution, even if the mean looks favorable.

From learningmatter-mit/AtomisticSkills, MIT.

More from learningmatter-mit/AtomisticSkills