Metadata-Version: 2.4
Name: rustina
Version: 0.7.0
Classifier: Development Status :: 5 - Production/Stable
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Rust
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Dist: rdkit>=2024.3 ; extra == 'prep'
Requires-Dist: molscrub>=0.2.2 ; extra == 'prep'
Requires-Dist: meeko>=0.7.1,<0.8 ; python_full_version < '3.10' and extra == 'prep'
Requires-Dist: meeko>=0.7.1 ; python_full_version >= '3.10' and extra == 'prep'
Requires-Dist: scipy>=1.12.0 ; extra == 'prep'
Requires-Dist: gemmi>=0.6.0 ; extra == 'prep'
Requires-Dist: joblib ; extra == 'prep'
Requires-Dist: pandas ; extra == 'prep'
Requires-Dist: numpy>=1.24 ; extra == 'ranker-dev'
Requires-Dist: scipy>=1.10 ; extra == 'ranker-dev'
Requires-Dist: torch>=2.1 ; extra == 'ranker-dev'
Requires-Dist: rdkit>=2023.9 ; extra == 'ranker-dev'
Requires-Dist: spyrmsd>=0.6 ; extra == 'ranker-dev'
Requires-Dist: pandas>=2.0 ; extra == 'ranker-dev'
Requires-Dist: pyarrow>=14 ; extra == 'ranker-dev'
Requires-Dist: posebusters>=0.6.5 ; extra == 'ranker-dev'
Requires-Dist: pytest>=7 ; extra == 'ranker-dev'
Provides-Extra: prep
Provides-Extra: ranker-dev
License-File: LICENSE
Summary: CPU molecular docking with native GNINA-compatible CNN rescoring
Author: Adrien H. Cerdan
Requires-Python: >=3.9
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Homepage, https://github.com/AdrienCerdan/rustina
Project-URL: Issues, https://github.com/AdrienCerdan/rustina/issues
Project-URL: Repository, https://github.com/AdrienCerdan/rustina

# Rustina

Rustina is a molecular-docking engine written in Rust. It provides Vina,
Vinardo and AutoDock 4.2 empirical scoring, Monte Carlo iterated local search
on CPU and GPU, native GNINA-compatible CNN rescoring, and a learned GNN pose
ranker. The supported interfaces are Python and the command-line program.

Rustina v0.6 targets Linux x86-64 and Python 3.9 or newer.

## Search modes

Docking is a Monte Carlo iterated local search. There is one docking path; the
search policy is chosen with `--search-mode` (CLI) or `search_mode=` (Python).

- **`qvina2-race`** (default) — the pinned `qvina2` budget plus a per-run
  stagnation window: a run whose incumbent score has not improved for
  `steps/8` MC steps stops early, while still-improving runs keep their full
  budget. The decision is per-run and score-based, so seeded replay is
  deterministic regardless of thread count. Requires at least 16 runs
  (`--runs` defaults to 32). The 393-target campaign measured 65.4% vs 63.9%
  top-1 at 1.23x less wall time than plain `qvina2` at `--runs 8`.
- **`qvina2`** — the pinned compatibility mode: the Vina budget with the
  official QuickVina 2 BFGS-history significance test.
- **`vina`** — the AutoDock Vina 1.2.3 Monte Carlo/BFGS policy.
- **`gpu`** — the wgpu batched lane search. See below.

The compatibility modes (`vina`, `qvina2`, and by inheritance `qvina2-race`)
reject custom steps, BFGS iteration counts, and RILC: their budgets and
validation rules are pinned so benchmark numbers stay reproducible.

`--start-seeds registration` (Python/API parity via `RUSTINA_SEED=registration`)
seeds the aligned fraction of runs with deterministic pharmacophore-triplet
registration poses instead of random pocket placements. It is opt-in and was
measured accuracy-neutral to +1 point, PoseBusters-validity identical, and
cost-neutral.

Additional research search policies exist in the library and still parse from
both interfaces so archived commands keep running, but they are hidden from
`--help`, unsupported, and must not be used for reported numbers. The
exhaustive 3D-FFT correlation engine is likewise research code: it is reachable
from `examples/correlation_dock.rs` and from `rustina::pipeline::dock`, not from
the CLI or the Python package. `CLAUDE.md` records what it has measured.

## Scoring functions

`--scoring` / `scoring=` selects the intermolecular potential. All three work
under every search mode, on CPU and on GPU.

| value | terms | notes |
|---|---|---|
| `vina` (default) | steric, hydrophobic, h-bond | no electrostatics, no desolvation |
| `vinardo` | steric, hydrophobic, h-bond | reparameterised Vina form |
| `ad4` | 12-6, directional 12-10 h-bond, Coulomb, desolvation | charge-dependent |

AD4 is the only function with an electrostatic and a desolvation term, so it is
the one to reach for on charged ligands and charged pockets — a cationic
cofactor against a carboxylate-lined site is a case where Vina and Vinardo
place the ligand essentially at random. Selecting `ad4` makes the PDBQT partial
charge column load-bearing on **both** the receptor and the ligand; Rustina
fails rather than scoring a zeroed Coulomb term in silence. AD4 also treats
polar hydrogens as real scoring atoms (the 12-10 runs from the hydrogen, not
from the donor heavy atom), and its torsional entropy is additive rather than a
divisor, so AD4 affinities are not interchangeable with Vina/Vinardo ones.

## GPU

`--search-mode gpu` / `search_mode="gpu"` runs the batched wgpu lane search. It
supports `vina`, `vinardo` and `ad4` — AD4 evaluates electrostatics,
desolvation and the torsional entropy penalty on the device, so supply
meaningful partial charges. There is no automatic CPU fallback: if no adapter
is available the run fails rather than silently changing engine.

- `--gpu-lanes` / `gpu_lanes` overrides the lane count (0 = auto).
- `--gpu-profile` / `gpu_profile` selects the search budget: `fast`
  (128 lanes x 32 steps), `balanced` (default; auto lanes x auto depth), or
  `deep` (512 x 128).
- `rustina.screen(..., scheduler="ligand")` submits valid ligands together and
  reuses the receptor maps. Screen batches run as pipelined waves: one wave is
  searched on the device while the previous wave is rescored and ranked on the
  CPU. The automatic wave size targets 2048 lane workgroups (16 ligands for the
  `fast` profile); `RUSTINA_GPU_WAVE` overrides it.

## Waters

`--waters` / `waters=` selects who owns the water and whether the ligand may
push it out.

| value | receptor waters | ligand `W` pseudo-waters | scoring |
|---|---|---|---|
| `fixed` (default) | kept, rigid, scored as ordinary receptor atoms | refused | any |
| `strip` | removed before the grids are built | refused | any |
| `displaceable` | removed | required, scored against the water map | `ad4` only |
| `toggle` | removed from the grids, restored as switchable sites | refused | any |

`fixed` is the historical behaviour and is byte-identical to the flag not
existing: a `HOH`, `WAT`, `DOD` or `H2O` record left in the receptor PDBQT is a
fixed, non-displaceable part of the site that the ligand can hydrogen bond to
and cannot move.

`displaceable` is AutoDock 4.2 hydrated docking (Forli & Olson 2012), where the
water moves onto the **ligand**. Prepare the ligand with Meeko's `--hydrate`
(or the reference `wet.py`), which attaches a rigid `W` pseudo-atom along each
hydrogen-bond vector of every polar ligand atom — 3.0 A from an acceptor, 2.0 A
from a polar hydrogen. Each `W` rides the pose and is scored against a water map
synthesised from the receptor's OA and HD maps, exactly as `mapwater.py` builds
it: a hydration site the pose keeps pays `0.6 x min(OA, HD)`, and one the pose
displaces earns a flat `-0.2`. Which waters survive is therefore an output of
the search, not an input to it. `W` atoms are otherwise structurally invisible —
they are excluded from heavy-atom counts, RMSD, TORSDOF and intramolecular
pairs.

Two combinations are hard errors rather than warnings, because each would
otherwise produce a plausible wrong number:

- `displaceable` under `vina` or `vinardo`. Neither function has a water term
  and the AD4 constants do not transfer, so the run would look like hydrated
  docking and be dry docking with some invisible atoms attached.
- a ligand carrying `W` atoms under `fixed`, `strip` or `toggle`. Under `fixed`
  and `toggle` the crystallographic waters are still modelled, so the hydration
  shell would be counted twice.
- `toggle` against a receptor carrying no water records at all, which would
  otherwise be a flag that silently did nothing.

`toggle` is GOLD's model (Verdonk et al. 2005), where the water stays on the
**receptor** and gains two degrees of freedom: an on/off switch and a rotation
about its own oxygen. Each water contributes

```text
E_w = min( 0,  min over spins s of [ sigma_p + (E_prot(s) - E_prot_ref) + E_lig(s) ] )
```

so a water is kept only when what it earns from the pose outweighs the entropic
cost `sigma_p` of ordering it, and a pose that overlaps a water simply makes the
on-state expensive and the off-state win. The outer `min(0, .)` is the water
switched off, which is why no clash cap is needed on this term — the toggle is
itself the bound. Which waters survive is again an output of the search.

- `--water-penalty` / `water_penalty=` is `sigma_p`, in kcal/mol. The default
  `2.0` is **not** GOLD's fitted constant; it is Dunitz's (1994) bound on the
  entropic cost of ordering one water at 300 K, i.e. the right order of
  magnitude and an upper bound. Treat it as uncalibrated and report the value
  used. Larger values displace more waters.
- `--water-spins` / `water_spins=` is the number of orientations sampled per
  water (default 12). It is ignored outside `--scoring ad4`: Vina and Vinardo
  are united-atom with an isotropic heavy-atom hydrogen-bond term, so a water's
  orientation is invisible to them and one orientation is exact. Under AD4 the
  12-10 runs from the polar hydrogen, so the rotation is real.
- The apo reference `E_prot_ref` is a deliberate deviation from GOLD. A buried
  water is well bonded to the protein whether or not a ligand is present, so
  charging that raw would switch every water on regardless of the pose and make
  `sigma_p` uninterpretable. Subtracting the best protein-water energy the site
  reaches on its own measures the water's *change* on binding, which puts
  `sigma_p` on the score's own kcal/mol axis.
- Waters do not see each other: each site is scored independently against the
  protein and the ligand. GOLD does not model water-water coupling either.
- `toggle` is a host-side term and is refused under `--search-mode gpu`.

Not implemented, deliberately: `dry.py`'s post-run filtering of displaced waters
from the output pose and its B-factor encoding of strong/weak retention.

## CNN rescoring

CNN rescoring is on by default on every surface and uses GNINA's distilled
`fast` model, which costs about 1.1x empirical scoring. Set `cnn="default"` /
`--cnn default` for Rustina's bundled GNINA 1.3 three-model ensemble
(marginally better, 1.3-1.9x slower), `skip_cnn=True` / `--skip-cnn` to rank by
the empirical score alone, or pass a native weight file path.

On the controlled PoseBusters protocol the CNN is the largest single measured
lever: **+13.6 points, p = 0.00032**.

`cnn_rotations` (CLI `--cnn-rotations`, 1-24, default 1) averages the CNN pose
score over that many orientations of the voxel grid, the equivalent of GNINA's
`--cnn_rotations`. Only the rescore is affected: the search, the poses, and
every empirical number are identical at any setting. The default of 1 is
deliberate. Measured over 361 redocking targets by rescoring stored poses at
all 24 orientations, averaging changed the chosen pose on 4 of 361 targets and
never improved one, while costing about 7 s per target on six cores. See
`CLAUDE.md` for the measurement.

## Installation

Install the Python package:

```bash
python -m pip install rustina
```

Ligand preparation is an optional extra because it installs RDKit, molscrub,
Meeko, SciPy, and Gemmi:

```bash
python -m pip install "rustina[prep]"
```

For a source checkout, build a production extension with:

```bash
maturin develop --release
```

## Python API

Inputs to docking and scoring are prepared PDBQT paths or raw PDBQT strings.

```python
import rustina

poses = rustina.dock(
    receptor="data/receptor.pdbqt",
    ligand="data/ligand.pdbqt",
    cx=-14.0,
    cy=18.0,
    cz=-15.0,
    sx=14.0,
    sy=18.0,
    sz=15.0,
    runs=32,
    threads=8,
    seed=42,
)
print(poses[0]["affinity"], poses[0]["cnn_score"])
```

`num_modes` caps the number of ranked poses returned, exactly as `--num-modes`
caps the models the CLI writes. It is applied after ranking, so it never
changes the search or which pose comes first.

The low-level CLI and Python API require an explicit pocket center and size.
The product-style benchmark driver defaults to ligand-sized dynamic boxes,
using the crystallographic heavy-atom extent plus adaptive padding
`max(4 Å, 8 Å - 0.2 Å × torsions)`. The fixed-25 release qualification remains
explicitly pinned to its published protocol.

Screen multiple ligands while reusing the receptor grid:

```python
results = rustina.screen(
    receptor="data/receptor.pdbqt",
    ligands=["ligand-1.pdbqt", "ligand-2.pdbqt"],
    cx=-14.0,
    cy=18.0,
    cz=-15.0,
    sx=14.0,
    sy=18.0,
    sz=15.0,
    threads=8,
)
```

Experimental negative-image screening compiles typed pocket hotspots once,
generates rigid pharmacophore alignments, and selects/refines them with the
atom-specific Vina grids. It is a retrieval prototype, not a replacement for
final docking:

```python
pocket = rustina.prepare_negative_image(
    "data/receptor.pdbqt",
    cx=-14.0, cy=18.0, cz=-15.0,
    sx=14.0, sy=18.0, sz=15.0,
)
matches = pocket.screen(["ligand-1.pdbqt", "ligand-2.pdbqt"], threads=6)

# Pose-generation diagnostics: refine and serialize up to 32 distinct basins.
candidates = pocket.screen_candidates(
    "ligand-1.pdbqt", limit=16, refine_candidates=32
)
```

Each successful record reports physical grid energy, pharmacophore coverage,
clash diagnostics, pre/post-refinement transforms, selection provenance, and a
ready-to-score PDBQT pose. Malformed ligands return an isolated
`{"error": ...}` record. Use `scripts/benchmark_negative_image.py` for
versioned teacher-recall experiments; do not interpret the prototype score as
experimental binding affinity.

Before promoting RIFT as a docking prescreen, run
`scripts/run_rift_pose_qualification.py`. It evaluates randomized crystal,
single-ETKDG, and eight-ETKDG inputs on the fixed-box Astex diagnostic panel,
reporting raw proposal RMSD separately from grid-refined RMSD and physical
PoseBusters validity. Qualification and retrieval benchmark drivers print
configuration/target progress as they run and keep resumable raw results under
`scratch/`.

Optional ligand preparation accepts SMILES, structure files/blocks, or an
RDKit molecule:

```python
states = rustina.prepare_ligand("CC(=O)Nc1ccc(O)cc1", random_seed=42)
print(states[0]["pdbqt"])
```

### Template docking

Template docking is a supported, opt-in v1 feature for congeneric series when a
ligand with known coordinates is already in the receptor coordinate frame. It requires the
`prep` extra because maximum-common-substructure matching is performed with
RDKit:

```python
poses = rustina.dock_reference(
    receptor="receptor.pdbqt",
    ligand="query.sdf",
    reference="co-crystal-ligand.sdf",
    cx=-14.0,
    cy=18.0,
    cz=-15.0,
    sx=14.0,
    sy=18.0,
    sz=15.0,
    runs=8,
    seed=42,
)
print(
    poses[0]["relaxed_core_rmsd"],
    poses[0]["physical_score"],
    poses[0]["reference_satisfied"],
)
```

The mapped heavy-atom core is restrained by a soft flat-bottom potential during
search, followed by a short unrestrained relaxation. Results report constrained
and relaxed core RMSD, physical and guided scores, the restraint penalty,
mapping identity, both MCS coverage fractions, and whether the relaxed core
remains within the default 1.0 Angstrom satisfaction threshold.

Automatic guidance requires at least six mapped heavy atoms and 50% query
coverage. An explicit `atom_map` can define a smaller anchor. Rustina fails
rather than silently switching to free docking when automatic guidance does not
meet these gates. `screen_reference()` applies the same workflow to a series
while reusing one receptor grid.

Up to eight symmetry-distinct MCS mappings are evaluated deterministically.
The exact requested run budget is distributed across them globally; mappings do
not multiply the run count. Candidates are merged and deduplicated before one
CNN rescore using the requested `cnn_pool_size`.

#### Timing smoke test

A release-mode `5SAK_ZRY` methyl-analog smoke test used four total runs, four
mapping hypotheses, empirical scoring, matched pose-pool budgets, and three
seeds. Median wall times were:

| CPU threads | Normal docking | Template docking | Difference |
|---:|---:|---:|---:|
| 2 | 0.711 s | 0.732 s | +3% |
| 8 | 0.620 s | 0.682 s | +10% |

The two-thread end-to-end template call, including MCS generation and final
relaxation, took 0.758 s median. These numbers characterize one small smoke
case, not expected performance across ligand series. The solved-congeneric
qualification requirements are documented in
[`docs/benchmarks/REFERENCE_DOCKING_PROTOCOL.md`](docs/benchmarks/REFERENCE_DOCKING_PROTOCOL.md).

Template docking assumes the reference and target receptor use the same
coordinate frame. It does not align receptor structures or provide shape-only
or pharmacophore guidance. `5SAK_ZRY` demonstrates restraint behavior and
physical validity only; its methyl analog has no experimental pose and is not
an accuracy benchmark.

`rustina.build_profile()` returns `debug` or `release`. Performance results are
valid only when the actually imported extension reports `release`.

## Command line

```bash
rustina dock \
  --receptor data/receptor.pdbqt \
  --ligand data/ligand.pdbqt \
  --output docked.pdbqt \
  --cx -14 --cy 18 --cz -15 \
  --sx 14 --sy 18 --sz 15 \
  --search-mode qvina2-race \
  --scoring vina \
  --runs 32 --threads 8 --seed 42 \
  --num-modes 9

rustina score \
  --receptor data/receptor.pdbqt \
  --ligand data/ligand.pdbqt \
  --scoring ad4 \
  --cnn
```

Batch docking prepares the receptor once and docks every ligand in sorted order
into an output directory with a `summary.csv`:

```bash
rustina dock --batch ligands/ --output results/ \
  --receptor data/receptor.pdbqt \
  --cx -14 --cy 18 --cz -15 --sx 14 --sy 18 --sz 15
```

GPU mode streams the batch through the shared pipelined submission; CPU modes
dock serially with runs parallelized across threads. Each ligand stays
seeded-deterministic.

Use `rustina dock --help` and `rustina score --help` for the complete supported
options.

## Reproducibility and benchmarks

Docking is deterministic when `seed` is provided. Release benchmarks must
record the Rustina version, build profile, input dataset revision, complete
arguments, CPU model, and raw per-target results. A benchmark arm must name its
search mode and scoring function explicitly rather than relying on a default: a
default is not a protocol, and when it moves the arm either changes meaning
silently or stops running. Historical pre-v1 research is preserved by the
`pre-v1-research` tag; v1 benchmark qualification lives under
`docs/benchmarks/`.

### PoseBusters + Astex Benchmark

All numbers below come from a single binary (`97d0dd73`), seed 42, budget 8,
393 targets (308 PoseBusters + 85 Astex), pinned to six physical cores on a
Ryzen 5 5600X. Wilson 95% intervals are on the all-target denominator; at ~390
targets the minimum detectable difference is roughly ±4 points, so do not read
a smaller gap as real without a paired test.

Superseded numbers: this section previously reported 79.1% over "485 evaluated
targets" and 68.6% over "433" — both larger than the 393 targets that exist.
Those rows pooled records from several different binaries into one figure. See
`docs/benchmarks/MERGE_QUALIFICATION_RESULTS.md`.

#### Controlled protocol — the headline

Generated ETKDGv3+UFF start conformers and the pinned fixed 25 Å box;
crystallographic coordinates are evaluation-only. This is the protocol
comparable to published PoseBusters results.

| Search mode | CNN | RMSD<=2A (95% CI)  | PB-valid | Mean s |
|-------------|------|--------------------|----------|--------|
| `rustina` (research policy) | fast | 53.7% (48.7-58.6)  | 95.7%    | 2.18   |
| `qvina2`    | fast | 63.9% (59.0-68.5)  | 96.9%    | 4.83   |

On this protocol the `qvina2` compatibility mode outperforms the research
`rustina` policy by 10.2 points — but it also takes 2.2x as long at the same
`--runs 8`, so an unknown part of that gap is simply more search. A
budget-matched comparison has not been run.

#### Experimental protocol — a labelled upper bound, not a headline

Crystal ligand conformers fed as docking input and a box sized from the crystal
ligand extent. Both leak the answer, so this runs about 28 points high. It is
published because it is the corpus used for offline re-ranking research, not
because it is a fair number.

| Search mode | CNN | RMSD<=2A | Top-5 <=2A | PB-valid | Any pose <=2A | Mean s |
|---------|------|----------|------------|----------|---------------|--------|
| `rustina` | fast | 81.4%    | 92.6%      | 98.0%    | 97.5%         | 1.84   |

The gap between 81.4% top-1 and the 97.5% sampling ceiling is the engine's real
limitation: near-native poses are usually generated and then ranked below #1.

Rustina is research software. Docking scores and predicted poses are not a
substitute for experimental evidence or clinical decision-making.

## Rust development checks

```bash
cargo fmt --all --check
cargo clippy --all-targets --locked -- -D warnings
cargo test --all-targets --locked
```

The test profile uses basic optimization (`opt-level = 1`) for the numerical
kernels and limited debug information. Debug assertions and integer overflow
checks remain enabled. The first test build compiles a separate set of optimized
artifacts; subsequent runs reuse them. For full variable/type information in a
debugger, use `CARGO_PROFILE_TEST_DEBUG=2 cargo test ...`. For unoptimized stepping,
also set `CARGO_PROFILE_TEST_OPT_LEVEL=0`.

Two GPU tests are budgeted rather than exhaustive, because both were gating the
whole suite's wall time on a single test:

```bash
# Full 1248-grid GPU/CPU CNN parity sweep (default is 4 chunks of 48 poses).
RUSTINA_GPU_PARITY_CHUNKS=26 cargo test --locked gpu_cnn_matches_cpu

# Engine throughput probe -- diagnostics only, asserts nothing.
cargo test --locked gpu_cnn_engine_throughput_bench -- --ignored --nocapture
```

Run the full parity sweep before releasing a change to the GPU CNN kernels.

The `scoring_bench` example measures fixed-work grid scoring independently of
the number of optimizer evaluations needed by a search:

```bash
cargo build --release --locked --example scoring_bench
RAYON_NUM_THREADS=1 taskset -c 0 target/release/examples/scoring_bench \
  --receptor receptor.pdbqt --ligand ligand.pdbqt \
  --cx 0 --cy 0 --cz 0 --size 25 --spacing 0.375 --scoring ad4
```

Supply the target's pocket center. The JSON records per-evaluation timings and
a fingerprint of maps, scores and gradients for exact before/after comparisons.
Use release binaries and an idle, pinned physical core for timing.

## Attribution and license

Rustina is MIT licensed. Its algorithms and bundled CNN parameters build on
AutoDock Vina, Smina, QuickVina2, GNINA, Vina-GPU, and related published work.
See [THIRD_PARTY_NOTICES.md](THIRD_PARTY_NOTICES.md) and
[CITATION.cff](CITATION.cff) for provenance and citations. Model conversion is
documented in [docs/models.md](docs/models.md).

