---
name: auto-ligand-docking
description: >
  Agent-native executable skill for zero-configuration molecular docking.
  Solves the Blind Agent Problem: AI agents cannot visually inspect protein
  structures to specify docking search boxes, so this skill automates that
  step. Given an RCSB PDB ID and a ligand SMILES string, Auto-Ligand
  downloads the receptor, prepares receptor and ligand PDBQT files (Open Babel
  + Meeko), places the Vina search box on the whole-protein geometric centroid
  by default or on a named residue (with optional chain disambiguation) for
  pocket-targeted docking, scales exhaustiveness by the ligand's
  rotatable-bond count unless overridden, runs AutoDock Vina, and writes
  results_docked.pdbqt plus a binding score (kcal/mol). Provides a formal
  three-criterion machine-checkable verification protocol (exit code, stdout
  banner, PDBQT remark) so agents can confirm success without inspecting a
  molecular viewer. Use for agent-driven virtual screening, reproducible
  structure-based workflows, or binding-affinity estimation. Requires
  micromamba or conda (conda-forge), outbound HTTPS to files.rcsb.org, and a
  Unix shell. Co-authored with Claw 🦞 for Claw4S authorship rules.
license: MIT
compatibility: >
  Micromamba or conda; Python 3.12 environment from inlined environment.yml;
  network access to RCSB; Linux or macOS recommended (conda-forge
  vina/rdkit/openbabel wheels). Single-file submission: this SKILL.md is
  self-contained — agents must materialise the files listed below before
  running.
metadata:
  author: Gautam Naik
  version: "1.0.0"
  claw4s_authorship: >
    First or corresponding author must list Claw 🦞 as co-author
    (see claw4s.github.io).
allowed-tools:
  Bash(python3 *), Bash(mkdir *), Bash(cat *), Bash(echo *),
  Bash(micromamba *), Bash(conda *), Bash(grep *), Bash(test *), Bash(cd *)
---

# Auto-Ligand Docking  (RCSB PDB + SMILES → AutoDock Vina)

End-to-end molecular docking that an AI agent can run without writing chemistry
code or manually specifying search-box coordinates.

**The Blind Agent Problem (solved):** Molecular docking requires specifying a
3-D search box around the binding site — normally done by a human visually
inspecting the structure. Auto-Ligand replaces this perceptual step with
automatic box placement: whole-protein geometric centroid by default, or
residue-centred when a pocket or catalytic residue is known.

**Science in one line:** A PDB identifier + a SMILES string → a docked 3-D
pose + a Vina binding-affinity estimate (kcal/mol).

**Box placement:**
- Default → whole-protein geometric centroid, 20 × 20 × 20 Å
- `--residue N` → centroid of residue N, 15 × 15 × 15 Å
- Both box sizes overridable with `--box-size`

**Exhaustiveness:** auto-scaled: E = max(8, 8 + 2 × RBC); overridable with
`--exhaustiveness`.

**Verification protocol (3 criteria — all must pass):**

| Check | Required value |
|-------|---------------|
| Exit code | `0` |
| stdout | Contains `DOCKING SUCCESSFUL` AND `Best binding affinity (ΔG):` with a number |
| Output PDBQT | File exists and contains `REMARK VINA RESULT` |

> **Platform note:** The exact ΔG value can vary by up to ±0.5 kcal/mol
> across OS and library builds due to floating-point non-determinism in Vina.
> Assert the **presence** of the banner and remark — never hard-code a float.

---

## Agent Checklist

1. **Step 1** — Create workspace.
2. **Step 2** — Write `environment.yml`.
3. **Step 3** — Create the conda environment.
4. **Step 4** — Write `auto_ligand.py`.
5. **Step 5** — Run a docking job (smoke test or custom).
6. **Step 6** — Verify exit code, stdout, and PDBQT output.
7. **Step 7** — Run optional shell verification commands.
8. **Step 8** — (Optional) Multi-target generalizability check.

---

## Step 1 — Create Workspace

Use a relative path so the skill works in any agent working directory:

```bash
mkdir -p auto_ligand_docking
echo "Workspace: $(pwd)/auto_ligand_docking"
```

---

## Step 2 — Write `environment.yml`

```bash
cat > auto_ligand_docking/environment.yml <<'EOF'
name: auto-ligand-docking
channels:
  - conda-forge
dependencies:
  - python=3.12
  - pip
  - vina
  - rdkit
  - numpy
  - scipy
  - requests
  - openbabel
  - gemmi
  - pip:
      - meeko>=0.5.0
EOF
echo "Wrote environment.yml"
```

> **Reproducibility note:** Do not substitute a pip-only install. Receptor
> PDBQT preparation requires **Open Babel** from conda-forge; the pip wheel
> is not reliably maintained across platforms.

---

## Step 3 — Create the Conda Environment

```bash
cd auto_ligand_docking
micromamba create -y -p ./.conda -f environment.yml
```

If `micromamba` is unavailable:

```bash
cd auto_ligand_docking
conda env create -p ./.conda -f environment.yml
```

**Non-interactive runs:** prefer `micromamba run -p ./.conda ...` or
`conda run -p ./.conda ...` so no shell hook is required.

---

## Step 4 — Write `auto_ligand.py`

```bash
cat > auto_ligand_docking/auto_ligand.py <<'PY'
"""
Auto-Ligand: Agent-executable Virtual Screening Skill.
Authors: Gautam Naik & Claw 🦞
Environment: conda/mamba only — create `./.conda` from `environment.yml`.
"""
from __future__ import annotations
import argparse
import sys
from pathlib import Path
import numpy as np
import requests
from openbabel import openbabel as ob
from rdkit import Chem
from rdkit.Chem import AllChem, Lipinski
from meeko import MoleculePreparation, PDBQTWriterLegacy

RCSB_DOWNLOAD_URL = "https://files.rcsb.org/download/{pdb_id}.pdb"
DEFAULT_EXHAUSTIVENESS = 8
DEFAULT_N_DOCK_POSES = 1
DEFAULT_N_WRITE_POSES = 1

# Default ligand is Aspirin — a lightweight, fast smoke-test molecule.
# Aspirin is NOT a c-Abl inhibitor; the default run is a workflow test only.
DEFAULT_TEST_PDB    = "1iep"
DEFAULT_TEST_SMILES = "CC(=O)Oc1ccccc1C(=O)O"  # Aspirin

VINA_MISSING_MSG = """\
The Python package 'vina' is not available in this environment.
  Create the project-local env:
    micromamba create -y -p .conda -f environment.yml
    micromamba run -p ./.conda python3 auto_ligand.py
  Or run with --prepare-only to validate setup without invoking Vina."""


# ── Receptor download ─────────────────────────────────────────────────────────
def fetch_receptor(pdb_id: str, output_dir: Path = Path(".")) -> Path:
    pdb_id = pdb_id.strip().lower()
    if len(pdb_id) != 4:
        raise ValueError(f"PDB ID must be 4 characters, got {pdb_id!r}")
    url = RCSB_DOWNLOAD_URL.format(pdb_id=pdb_id)
    response = requests.get(url, timeout=60)
    response.raise_for_status()
    dest = Path(output_dir) / f"{pdb_id}.pdb"
    dest.write_bytes(response.content)
    print(f"Receptor {pdb_id} downloaded -> {dest.resolve()}")
    return dest.resolve()


# ── Receptor preparation ──────────────────────────────────────────────────────
def prepare_receptor_pdbqt(receptor_pdb: Path,
                            output_pdbqt: Path | None = None) -> Path:
    out = output_pdbqt or receptor_pdb.with_suffix(".pdbqt")
    conv = ob.OBConversion()
    conv.SetInFormat("pdb")
    conv.SetOutFormat("pdbqt")
    conv.AddOption("r", ob.OBConversion.OUTOPTIONS)
    mol = ob.OBMol()
    if not conv.ReadFile(mol, str(receptor_pdb)):
        raise RuntimeError(f"Open Babel could not read receptor PDB: {receptor_pdb}")
    if not conv.WriteFile(mol, str(out)):
        raise RuntimeError(f"Open Babel could not write receptor PDBQT: {out}")
    print(f"Receptor PDBQT written -> {out.resolve()}")
    return out.resolve()


# ── Ligand preparation ────────────────────────────────────────────────────────
def prepare_ligand(smiles: str, output_pdbqt: Path) -> Path:
    """
    Build 3-D coordinates from SMILES and write a Vina-ready PDBQT via Meeko.
    Meeko assigns partial charges automatically, avoiding the silent failure
    mode caused by missing charges in ad-hoc SMILES-to-PDBQT pipelines.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f"Invalid SMILES: {smiles!r}")
    mol = Chem.AddHs(mol)
    params = AllChem.ETKDG()
    if AllChem.EmbedMolecule(mol, params) != 0:
        raise RuntimeError("RDKit could not embed 3-D coordinates for the ligand.")
    AllChem.MMFFOptimizeMolecule(mol)
    preparator = MoleculePreparation()
    mol_setups = preparator.prepare(mol)
    if not mol_setups:
        raise RuntimeError("Meeko returned no setups for the ligand.")
    pdbqt_string, is_ok, _err = PDBQTWriterLegacy.write_string(mol_setups[0])
    if not is_ok:
        raise RuntimeError("Meeko failed to produce a valid ligand PDBQT.")
    output_pdbqt = Path(output_pdbqt)
    output_pdbqt.write_text(pdbqt_string, encoding="utf-8")
    print(f"Ligand PDBQT written -> {output_pdbqt.resolve()}")
    return output_pdbqt.resolve()


# ── Box-centre calculation ────────────────────────────────────────────────────
def protein_geometric_center(pdb_file: Path) -> list[float]:
    mol = Chem.MolFromPDBFile(str(pdb_file), removeHs=False, sanitize=False)
    if mol is None or mol.GetNumConformers() == 0:
        raise ValueError(f"No 3-D coordinates readable from {pdb_file}")
    pos = mol.GetConformer().GetPositions()
    return np.mean(pos, axis=0).tolist()


def residue_geometric_center(pdb_file: Path, residue_number: int,
                              chain_id: str | None = None) -> list[float]:
    """
    Centroid of all atoms in the given PDB residue number.
    chain_id disambiguates residue numbers in homodimeric assemblies
    (e.g. PDB 1IEP, chains A and B share identical residue numbering).
    """
    mol = Chem.MolFromPDBFile(str(pdb_file), removeHs=False, sanitize=False)
    if mol is None or mol.GetNumConformers() == 0:
        raise ValueError(f"No 3-D coordinates readable from {pdb_file}")
    chain_key = chain_id.strip().upper() if chain_id is not None else None
    atoms: list[int] = []
    for a in mol.GetAtoms():
        info = a.GetPDBResidueInfo()
        if info is None:
            continue
        if info.GetResidueNumber() != residue_number:
            continue
        if chain_key is not None and info.GetChainId().strip().upper() != chain_key:
            continue
        atoms.append(a.GetIdx())
    if not atoms:
        detail = f" on chain {chain_id!r}" if chain_id else ""
        raise ValueError(f"Residue {residue_number}{detail} not found in {pdb_file}")
    conf = mol.GetConformer()
    coords = np.array([conf.GetAtomPosition(i) for i in atoms], dtype=float)
    return np.mean(coords, axis=0).tolist()


# ── Dynamic exhaustiveness ────────────────────────────────────────────────────
def dynamic_exhaustiveness_from_smiles(smiles: str) -> int:
    """E = max(8, 8 + 2 × RBC)  where RBC = rotatable-bond count."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f"Invalid SMILES: {smiles!r}")
    n_rot = int(Lipinski.NumRotatableBonds(mol))
    return max(DEFAULT_EXHAUSTIVENESS,
               DEFAULT_EXHAUSTIVENESS + 2 * n_rot)


# ── Docking ───────────────────────────────────────────────────────────────────
def run_docking(receptor_pdbqt: Path, ligand_pdbqt: Path, *,
                box_center: list[float], box_size: list[float],
                exhaustiveness: int = DEFAULT_EXHAUSTIVENESS,
                n_dock_poses: int = DEFAULT_N_DOCK_POSES,
                n_write_poses: int = DEFAULT_N_WRITE_POSES,
                results_pdbqt: Path = Path("results_docked.pdbqt")) -> float:
    try:
        from vina import Vina
    except ImportError as e:
        raise RuntimeError(VINA_MISSING_MSG) from e
    v = Vina(sf_name="vina")
    v.set_receptor(str(receptor_pdbqt))
    v.set_ligand_from_file(str(ligand_pdbqt))
    v.compute_vina_maps(center=box_center, box_size=box_size)
    v.dock(exhaustiveness=exhaustiveness, n_poses=n_dock_poses)
    energies = v.energies(n_poses=n_write_poses)
    best_total = float(energies[0][0])
    results_pdbqt = Path(results_pdbqt)
    v.write_poses(str(results_pdbqt), n_poses=n_write_poses, overwrite=True)
    print(f"Docked poses written -> {results_pdbqt.resolve()}")
    return best_total


# ── CLI ───────────────────────────────────────────────────────────────────────
def _parse_box_size(raw: str) -> list[float]:
    parts = [float(x.strip()) for x in raw.replace(",", " ").split() if x.strip()]
    if len(parts) != 3:
        raise argparse.ArgumentTypeError(
            "box_size must be three numbers, e.g. 20 20 20")
    if any(x <= 0 for x in parts):
        raise argparse.ArgumentTypeError("box dimensions must be positive")
    return parts


def main(argv: list[str] | None = None) -> int:
    parser = argparse.ArgumentParser(
        description="Auto-Ligand docking skill (RCSB + SMILES → AutoDock Vina)")
    parser.add_argument("--pdb",   default=DEFAULT_TEST_PDB,
                        help=f"RCSB PDB ID (default: {DEFAULT_TEST_PDB!r})")
    parser.add_argument("--smiles", default=DEFAULT_TEST_SMILES,
                        help="Ligand SMILES (default: Aspirin smoke-test)")
    parser.add_argument("--output-dir", type=Path, default=Path("."),
                        help="Directory for downloaded PDB and PDBQT files")
    parser.add_argument("--ligand-pdbqt",   type=Path,
                        default=Path("ligand.pdbqt"))
    parser.add_argument("--receptor-pdbqt", type=Path, default=None,
                        help="Skip receptor preparation if set")
    parser.add_argument("--results-pdbqt",  type=Path,
                        default=Path("results_docked.pdbqt"))
    parser.add_argument("--box-size", type=_parse_box_size, default=None,
                        help="Docking box (Å), three numbers. "
                             "Default: 15³ with --residue, 20³ otherwise.")
    parser.add_argument("--residue", type=int, default=None, metavar="N",
                        help="Centre the search box on PDB residue N")
    parser.add_argument("--chain", type=str, default=None,
                        help="PDB chain ID when using --residue "
                             "(e.g. A — disambiguates homodimers)")
    parser.add_argument("--exhaustiveness", type=int, default=None,
                        help="Vina exhaustiveness (default: 8 + 2×RBC)")
    parser.add_argument("--n-poses",     type=int, default=DEFAULT_N_DOCK_POSES)
    parser.add_argument("--write-poses", type=int, default=DEFAULT_N_WRITE_POSES,
                        help="Poses to write (≤ n-poses)")
    parser.add_argument("--prepare-only", action="store_true",
                        help="Download and build PDBQT files, then exit without "
                             "running Vina. Does NOT emit DOCKING SUCCESSFUL.")
    args = parser.parse_args(argv)

    # ── Validation ────────────────────────────────────────────────────────────
    if args.chain is not None and not args.chain.strip():
        print("Error: --chain must be a non-empty chain ID.", file=sys.stderr)
        return 1
    if args.chain is not None and args.residue is None:
        print("Error: --chain requires --residue.", file=sys.stderr)
        return 1
    if args.box_size is None:
        args.box_size = _parse_box_size(
            "15 15 15" if args.residue is not None else "20 20 20")
    if args.write_poses > args.n_poses:
        print("Error: --write-poses cannot exceed --n-poses.", file=sys.stderr)
        return 1

    args.output_dir.mkdir(parents=True, exist_ok=True)

    try:
        receptor_pdb = fetch_receptor(args.pdb, args.output_dir)

        if args.receptor_pdbqt is None:
            receptor_pdbqt = prepare_receptor_pdbqt(
                receptor_pdb,
                args.output_dir / f"{args.pdb.strip().lower()}.pdbqt",
            )
        else:
            receptor_pdbqt = args.receptor_pdbqt.resolve()
            if receptor_pdbqt.suffix.lower() != ".pdbqt":
                print("Error: --receptor-pdbqt must be a .pdbqt file.",
                      file=sys.stderr)
                return 1
            if not receptor_pdbqt.is_file():
                print(f"Error: receptor PDBQT not found: {receptor_pdbqt}",
                      file=sys.stderr)
                return 1

        lig_path = (
            args.ligand_pdbqt.resolve()
            if args.ligand_pdbqt.is_absolute()
            else (args.output_dir / args.ligand_pdbqt).resolve()
        )
        prepare_ligand(args.smiles, lig_path)

        # ── Box placement ─────────────────────────────────────────────────────
        if args.residue is not None:
            center = residue_geometric_center(
                receptor_pdb, args.residue, chain_id=args.chain)
            chain_note = (f" chain {args.chain.strip()!r}"
                          if args.chain else "")
            print(f"Box centre (Å, residue {args.residue}"
                  f"{chain_note} centroid): {center}")
        else:
            center = protein_geometric_center(receptor_pdb)
            print(f"Box centre (Å, whole-protein centroid): {center}")

        exhaustiveness = (
            args.exhaustiveness
            if args.exhaustiveness is not None
            else dynamic_exhaustiveness_from_smiles(args.smiles)
        )

        if args.prepare_only:
            print("\n" + "=" * 40)
            print("PREPARE ONLY — Vina not invoked (no docking run)")
            print(f"  receptor PDB:    {receptor_pdb}")
            print(f"  receptor PDBQT:  {receptor_pdbqt}")
            print(f"  ligand PDBQT:    {lig_path}")
            print(f"  box centre (Å):  {center}")
            print(f"  box size (Å):    {args.box_size}")
            exh_note = ("(explicit flag)" if args.exhaustiveness is not None
                        else "(dynamic from SMILES)")
            print(f"  exhaustiveness:  {exhaustiveness} {exh_note}")
            print("=" * 40)
            return 0

        print(f"Starting docking: {args.smiles!r} → {args.pdb!r} ...")
        print(f"Vina exhaustiveness: {exhaustiveness}")
        affinity = run_docking(
            receptor_pdbqt,
            lig_path,
            box_center=center,
            box_size=args.box_size,
            exhaustiveness=exhaustiveness,
            n_dock_poses=args.n_poses,
            n_write_poses=args.write_poses,
            results_pdbqt=(
                args.results_pdbqt.resolve()
                if args.results_pdbqt.is_absolute()
                else (args.output_dir / args.results_pdbqt).resolve()
            ),
        )

        print("\n" + "=" * 40)
        print("DOCKING SUCCESSFUL")
        print(f"Best binding affinity (ΔG): {affinity:.2f} kcal/mol")
        print("=" * 40)
        return 0

    except (requests.HTTPError, OSError, RuntimeError, ValueError) as e:
        sys.stdout.flush()
        print(f"Skill execution failed: {e}", file=sys.stderr)
        return 1


if __name__ == "__main__":
    raise SystemExit(main())
PY
echo "Wrote auto_ligand.py"
```

---

## Step 5 — Run Docking

### Smoke test (defaults: PDB `1iep`, Aspirin)

```bash
cd auto_ligand_docking
micromamba run -p ./.conda python3 auto_ligand.py
```

With conda:

```bash
cd auto_ligand_docking
conda run -p ./.conda python3 auto_ligand.py
```

### Prepare-only mode (validate environment without Vina)

> **Note:** `--prepare-only` exits 0 but does **not** emit `DOCKING SUCCESSFUL`.
> Verification criteria 2 and 3 apply only after a full docking run (Step 5 without this flag).

```bash
cd auto_ligand_docking
micromamba run -p ./.conda python3 auto_ligand.py --prepare-only
```

Expected stdout (prepare-only):
```
Receptor 1iep downloaded -> ...
Receptor PDBQT written -> ...
Ligand PDBQT written -> ...
Box centre (Å, whole-protein centroid): [x, y, z]
========================================
PREPARE ONLY — Vina not invoked (no docking run)
  receptor PDB:    ...
  receptor PDBQT:  ...
  ligand PDBQT:    ...
  box centre (Å):  [x, y, z]
  box size (Å):    [20.0, 20.0, 20.0]
  exhaustiveness:  14 (dynamic from SMILES)
========================================
```

### Custom ligand with whole-protein centroid

```bash
cd auto_ligand_docking
micromamba run -p ./.conda python3 auto_ligand.py \
  --pdb 1iep \
  --smiles "Cc1ccc(NC(=O)c2ccc(CN3CCN(C)CC3)cc2)cc1Nc1nccc(-c2cccnc2)n1" \
  --output-dir ./run_imatinib
```

### Pocket-targeted docking — Thr315 gatekeeper in c-Abl (chain A)

The `--chain` flag is required when the same residue number appears on multiple
chains (as in the 1IEP homodimeric assembly):

```bash
cd auto_ligand_docking
micromamba run -p ./.conda python3 auto_ligand.py \
  --pdb 1iep \
  --residue 315 \
  --chain A \
  --smiles "Cc1ccc(NC(=O)c2ccc(CN3CCN(C)CC3)cc2)cc1Nc1nccc(-c2cccnc2)n1" \
  --output-dir ./run_pocket
```

**Box defaults:**
- `--residue` set, `--box-size` omitted → **15 × 15 × 15 Å**
- `--residue` not set, `--box-size` omitted → **20 × 20 × 20 Å**

**Exhaustiveness formula** (when `--exhaustiveness` is omitted):

$$E = \max\!\left(8,\; 8 + 2 N_{\text{rot}}\right)$$

Full CLI reference: `micromamba run -p ./.conda python3 auto_ligand.py --help`

---

## Step 6 — Expected Outputs and Verification Criteria

Agents MUST assert all three criteria after a full docking run. Failure of any
one indicates a broken execution.

| Artefact / stream | Required criterion |
|-------------------|-------------------|
| **Exit code** | `0` |
| **stdout** | Contains `DOCKING SUCCESSFUL` **and** `Best binding affinity (ΔG):` with a numeric kcal/mol value |
| **`results_docked.pdbqt`** | File exists; contains a line matching `REMARK VINA RESULT` |

> **Important:** `--prepare-only` mode satisfies criterion 1 (exit 0) only.
> Criteria 2 and 3 require a completed Vina run without `--prepare-only`.

> **Platform note:** The exact ΔG value can vary by up to ±0.5 kcal/mol
> across different operating systems and library builds. Agents must **not**
> hard-code a specific float as a pass/fail criterion.

---

## Step 7 — Verification Commands

Run these immediately after Step 5 to confirm all three criteria:

```bash
# Criterion 1: exit code (run immediately after the docking command)
echo "Exit code: $?"   # Must print 0

# Criteria 2 + 3: stdout and PDBQT
cd auto_ligand_docking
test -f results_docked.pdbqt && \
  grep -q "REMARK VINA RESULT" results_docked.pdbqt && \
  echo "VERIFICATION PASSED: PDBQT contains Vina result remark." || \
  echo "VERIFICATION FAILED."
```

> **Note on exit-code capture:** Run `echo "Exit code: $?"` as the very next
> command after the docking invocation. Any intervening command overwrites `$?`.

Confirm that stdout from Step 5 also contains `DOCKING SUCCESSFUL` and
`Best binding affinity`.

---

## Step 8 — Multi-Target Generalizability Check (Optional)

To verify the skill works across structurally diverse targets, run the
following five jobs in sequence. All should pass the three verification
criteria. This uses the whole-protein centroid path (jobs 1, 3, 4, 5) and
the residue-targeted path with chain disambiguation (job 2).

```bash
cd auto_ligand_docking

# 1. Aspirin → c-Abl kinase (smoke test, centroid)
micromamba run -p ./.conda python3 auto_ligand.py \
  --pdb 1iep --smiles "CC(=O)Oc1ccccc1C(=O)O" \
  --output-dir ./multi_1iep_aspirin

# 2. Imatinib → c-Abl Thr315 pocket (residue-targeted, chain A)
micromamba run -p ./.conda python3 auto_ligand.py \
  --pdb 1iep --residue 315 --chain A \
  --smiles "Cc1ccc(NC(=O)c2ccc(CN3CCN(C)CC3)cc2)cc1Nc1nccc(-c2cccnc2)n1" \
  --output-dir ./multi_1iep_imatinib

# 3. Erlotinib → EGFR kinase (centroid)
micromamba run -p ./.conda python3 auto_ligand.py \
  --pdb 3htb \
  --smiles "C#Cc1cccc(Nc2ncnc3cc(OCCOC)c(OCCOC)cc23)c1" \
  --output-dir ./multi_3htb_erlotinib

# 4. Gefitinib → EGFR (centroid)
micromamba run -p ./.conda python3 auto_ligand.py \
  --pdb 4djv \
  --smiles "COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCOCC1" \
  --output-dir ./multi_4djv_gefitinib

# 5. Naproxen → COX-1 (centroid)
micromamba run -p ./.conda python3 auto_ligand.py \
  --pdb 1ajv \
  --smiles "COc1ccc2cc(C(C)C(=O)O)ccc2c1" \
  --output-dir ./multi_1ajv_naproxen
```

After all five complete, verify:

```bash
# All five results_docked.pdbqt files exist and contain Vina result remarks
for dir in multi_1iep_aspirin multi_1iep_imatinib multi_3htb_erlotinib \
           multi_4djv_gefitinib multi_1ajv_naproxen; do
  f="auto_ligand_docking/${dir}/results_docked.pdbqt"
  test -f "$f" && \
    grep -q "REMARK VINA RESULT" "$f" && \
    echo "PASS: ${dir}" || \
    echo "FAIL: ${dir}"
done
```

Expected output:
```
PASS: multi_1iep_aspirin
PASS: multi_1iep_imatinib
PASS: multi_3htb_erlotinib
PASS: multi_4djv_gefitinib
PASS: multi_1ajv_naproxen
```

---

## Gotchas

- **Whole-protein centroid is a coarse default.** For large or multi-domain
  proteins the centroid may be far from any druggable pocket. Prefer
  `--residue` (and `--chain` for homodimers) when a known binding site is
  available. For publication-quality docking, validate the box against a
  co-crystal ligand position or a cavity-detection tool.

- **Default PDB and ligand.** PDB `1IEP` is the c-Abl kinase domain
  co-crystallised with Imatinib. The default ligand SMILES is **Aspirin** —
  chosen for speed and chemical simplicity as a smoke test only. Aspirin is
  not a c-Abl inhibitor; the default run is not a pharmacological claim.

- **Chain disambiguation.** PDB 1IEP contains two chains (A and B) with
  identical residue numbering. Always supply `--chain A` (or B) when using
  `--residue` with this structure.

- **Exit-code capture.** `echo $?` captures the exit code of the most
  recently completed command. Run it as the immediate next line after the
  docking invocation; any intervening command overwrites it.

- **Prepare-only vs full run.** `--prepare-only` exits 0 but does not emit
  the `DOCKING SUCCESSFUL` banner. Only a full run (without `--prepare-only`)
  satisfies all three verification criteria.

- **Offline use.** Without network access to `files.rcsb.org`, place a local
  PDB file in the output directory and pass `--receptor-pdbqt` to skip
  receptor download and preparation.