Drug Redocking Rmsd

drug-redocking-rmsd

Goal

To quantitatively validate a docking protocol by computing the symmetry-corrected in-place heavy-atom RMSD between docked poses and the crystallographic reference ligand. A top-scored pose (pose 1) RMSD below 2.0 A is the standard threshold for a successful self-docking control.

Self-docking is a necessary, not sufficient, check. It verifies that your receptor preparation, box definition, and scoring function can recover a known pose in its own binding site. It does not verify that the protocol will work on new compounds. For a production virtual screen, complement self-docking with cross-docking into different receptor conformations when available (see the HTVS workflow Stage 3), and pair this RMSD check with drug-pose-validation to catch poses that are geometrically near-native but physically implausible (internal clashes, strained torsions).

Instructions

1. Compute RMSD from a crystal PDB reference

When the reference ligand is extracted from a PDB (HETATM records, no bond orders), provide the SMILES so the script can assign bond orders via template matching:

${CLAUDE_SKILL_DIR}/../../venv/run cpu python ${CLAUDE_SKILL_DIR}/scripts/compute_rmsd.py \
  --docked docking/ligand_docked.pdbqt \
  --reference crystal_ligand.pdb \
  --smiles "NS(=O)(=O)c1ccc(Nc2nc3[nH]cnc3c(OCC3CCCCC3)n2)cc1" \
  --output_dir validation/

2. Compute RMSD from an SDF reference

When the reference ligand is an SDF with proper bond orders (e.g., from a database or ligand-prep), no SMILES is needed:

${CLAUDE_SKILL_DIR}/../../venv/run cpu python ${CLAUDE_SKILL_DIR}/scripts/compute_rmsd.py \
  --docked docking/ligand_docked.pdbqt \
  --reference crystal_ligand.sdf \
  --output_dir validation/

3. Tune the pass/fail threshold

The default threshold is 2.0 A, which is the classical success criterion from the original docking validation literature. Modern docking programs often do substantially better, and the threshold should scale with ligand size and flexibility:

Below about 1.5 A is typically considered "good" and below 1.0 A is "very good" for modern docking of small rigid compounds. Above roughly 3 A, numeric ordering loses meaning (a 4 A pose is not usefully "better" than a 6 A pose; both are wrong).

${CLAUDE_SKILL_DIR}/../../venv/run cpu python ${CLAUDE_SKILL_DIR}/scripts/compute_rmsd.py \
  --docked docking/ligand_docked.pdbqt \
  --reference crystal_ligand.sdf \
  --threshold 1.5 \
  --output_dir validation/

4. Interpret the output

The script writes rmsd_results.json:

{
    "reference": "crystal_ligand.sdf",
    "docked": "ligand_docked.pdbqt",
    "n_poses": 5,
    "threshold": 2.0,
    "top_pose_rmsd": 0.823,
    "gate_pass": true,
    "gate_criterion": "Top-scored docked pose (pose 1) heavy-atom RMSD below threshold. ...",
    "best_rmsd": 0.823,
    "best_pose": 1,
    "best_rmsd_note": "best_rmsd is the minimum RMSD across all poses. Use as a diagnostic only: ...",
    "poses": [
        {"pose": 1, "rmsd_heavy_atom": 0.823, "pass": true},
        {"pose": 2, "rmsd_heavy_atom": 1.451, "pass": true},
        {"pose": 3, "rmsd_heavy_atom": 4.102, "pass": false}
    ]
}

• gate_pass is the protocol-validation verdict: true iff pose 1 (the top-scored pose) is within threshold. A near-native pose further down the list is not enough; if the scoring function cannot rank it first, the protocol is not working.
• top_pose_rmsd is the pose-1 RMSD, the value gate_pass keys off.
• best_rmsd / best_pose are diagnostics only. If best_pose > 1 but best_rmsd < threshold, the sampling is finding near-native conformations but the scoring function is failing to prioritize them. This is a scoring problem, not a sampling problem, and warrants rescoring or re-ranking rather than redoing the search.
• rmsd_heavy_atom is the symmetry-corrected in-place RMSD over all heavy atoms in Angstroms. The script uses RDKit's CalcRMS which enumerates molecular automorphisms and (by default) symmetrizes conjugated terminal groups like carboxylates and nitros.

5. Use in the HTVS validation gate

This skill is designed to be called during the protocol validation gate of the HTVS workflow (Stage 3). If gate_pass is false, revisit receptor preparation, protonation states, or box placement before proceeding to the production screen. If best_pose > 1 while top_pose_rmsd > threshold, the scoring function rather than the search is the bottleneck; consider alternative scoring functions or rescoring with a more expensive method.

See examples/README.md for a worked case using real NU6102 / CDK2 self-docking data from the cdk2-htvs HTVS campaign. That example also documents a real correctness discrepancy between this version of the skill and an earlier (buggy) version, and is worth reading if you have cached validation results from a previous run.