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.

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.

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

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:

python -m pip install rustina

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

python -m pip install "rustina[prep]"

For a source checkout, build a production extension with:

maturin develop --release

Python API

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

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:

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:

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:

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:

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:

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.

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

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:

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.

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.

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

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:

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

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 and CITATION.cff for provenance and citations. Model conversion is documented in docs/models.md.