Rustina · scoring programme

Beyond the Argmin


Six work items that extract ranking signal from the energy landscape Rustina already computes and currently discards. Most need no training at all; the ones that do need between one and a few hundred parameters, not a new network.

Headroom 97.5% pool ceiling vs ~62% top-1 Baselines F1x 61.8% · qvina2+CNN 61.3% Corpus SPINDR 35,349 / 68 / 225

The thesis

Every scoring module in src/interaction_ot, negative_image, pharmacophore_match, ranker, cnn — scores a pose. So does Vina, Vinardo, and GNINA's CNN. Competing there means out-training people with a thousand times the compute.

But the exhaustive FFT computes the score at every lattice translation for every enumerated rotation, and then keeps the argmin and throws the rest away. That discarded landscape is a null model, an entropy estimate, and a specificity test — none of which any pose-scorer can see even in principle, and all of which are already paid for. The programme below is the systematic recovery of that signal.

Correction — GBDT was mis-scoped

I argued you'd over-reached to deep learning on the grounds that ~400 targets can't support a GNN. Wrong: SPINDR packed v4 holds 35,349 train complexes. At that scale ranker v2 is entirely defensible and the over-parameterisation argument does not apply.

GBDT stays in the plan (W6) with a different job: a diagnostic baseline that measures how much of ranker v2's lift comes from the learned representation versus features a tree could have read, and the fusion layer where the new landscape features cash out. The real data constraint is elsewhere — the 68-complex val split is far too thin for hyperparameter selection, which shapes the protocol below.

Data discipline

SplitCountRole in this programme
train35,349All fitting, via grouped k-fold CV. Model selection happens here, not on val.
val68Smoke check only. Too small to select on — a 3-point swing is one complex.
test225Touched once per milestone. Never for tuning.

Join key is spindr_id — it maps manifest-test.json onto the packed test split 225/225, whereas the packed id field (1ia9__1__1.B__1.G, a PLINDER system id) overlaps zero. Use the wrong key and you will silently train on test.

Two standing rules. Every fit reports its CV fold spread, not just its mean — with only 68 val complexes the temptation to chase noise is severe. And before any number from this programme is quoted, fix the hardcoded pb_valid = True at scripts/evaluate_pure_equifield.py:102 and route validity through the real posebusters_top1 from run_spindr_matrix.py. Several of these items trade enthalpy for geometry; unmeasured validity would hide exactly that.

Shared prerequisites

S1 Pose-pool corpus over train ~4 CPU-h

Every ranking item needs (target, pose, symmetry-corrected RMSD) triples. Run the F3 arm (512 orientations, 0.375 Å, 1.42 s median) over train, retaining the top 32 poses per target rather than 8 — the extra poses cost nothing and give the ranker negatives near the decision boundary.

First
Check whether scratch/ranker-v2-guarded-seed42.json and the re_export_corpus_depth.py pipeline already carry pools with RMSDs. If so, reuse and skip the regeneration entirely.
Store
Poses plus the landscape summary from W1 — recomputing the field later would cost as much as the docking run itself.
Budget
35,349 × 1.42 s ÷ 6 cores ≈ 2.3 h, plus RMSD evaluation.
S2 One feature-extraction entry point ~1 day

Every item below emits features. Give them one home so the Rust scorer and the Python trainer cannot drift apart — the train/inference mismatch that killed EquiField (relu(2.5−d)²·5.0 in the trainer against 4.0·(3.4−d)² in equifield.rs:577) came from exactly this having no single owner.

pub struct PoseFeatures { /* ~60 named f32 fields */ }

pub fn extract(
    grids:     &ReceptorGrids,
    landscape: &LandscapeSummary,   // W1
    hydration: &HydrationField,     // W2
    pose:      &Pose,
) -> PoseFeatures;

Serialise as a stable named schema, never a bare float array. Emit the same struct from the Rust path and the Python export, and unit-test that they agree bit-for-bit on a fixed complex.

Work items

W1 Landscape functionals — basin free energy, normalisation, confidence 1 param~3 days

The flagship. Vina's entire entropy model is −0.0585 × N_rot: a torsion count. Replace it with the actual configurational integral over the field you already have.

Z(R) = Σ_T exp(−β·E(T,R)) · h³
F    = −kT · ln Σ_R w_R · Z(R)

One exp-and-sum pass over arrays correlate.rs already produced — negligible beside the FFT itself. Then watershed the field by steepest descent into basins, integrate Z per basin, and rank basins by F rather than by E_min. A broad basin beats a sharp one of equal depth, which is the right physics: crystal poses sit in redundant, robust basins while empirical false minima are knife-edge artefacts of one lucky contact.

New file
src/landscape.rs, consuming correlate::TranslationField
Emits
F_basin, basin volume, occupancy p_i = Z_i/ΣZ_j; pocket-normalised z (E_min − μ)/σ against the landscape's own distribution; landscape entropy −Σp ln p; rank gap between basin 1 and 2.
Params
β alone. Grid-search on train CV against top-1 accuracy.
Basins across rotations
Watershed is per-rotation. Merge basins whose minima fall within 1 Å RMSD — the retained-pose clustering already does this, so reuse it rather than inventing a second convention.

The pocket-normalised z may be the highest-value line here and is nearly free: Vina scores scale with ligand size, which is why absolute scores are near-useless for cross-ligand ranking. Normalising against this ligand's own landscape in this pocket cancels both ligand-size and pocket-size bias. Expect the effect in screening enrichment (vs_enrichment.py) more than in single-target top-1.

W2 Hydration field — enclosure × water affinity ~4 params~3 days

Where Vina is weakest. Its hydrophobic term is a distance-based contact count with no notion of enclosure or water structure, yet the hydrophobic effect is driven by the entropy of displaced water, not by carbon proximity.

Check first: ReceptorGrids already holds a per-atom-type map, so the O/OA map may already be W(x), the water-oxygen affinity at every voxel. If so this item needs no new correlation at all — only the enclosure term and the combination rule.

Enclosure B(x) comes from erosion.rs: the sliding-window minima machinery already computes how deeply a voxel sits inside the free-space mask. Then classify:

  • high B, low W → trapped, poorly-coordinated water → displacement is strongly favourable. This is the hydrophobic hotspot, properly derived.
  • high B, high W → conserved water → penalise displacement, unless the displacing ligand atom is a donor/acceptor that can replace its H-bonds.
New file
src/hydration.rs
Params
Two thresholds and two weights. Fit on train CV.
Cost
Receptor-side and built once — amortised across an entire screen, which is where it matters most.
Validation
Before touching docking: do predicted conserved-water sites coincide with crystallographic waters in the train receptors? That is a direct, honest check of the physics with no scoring involved.
W3 Specificity by counterfactual pocket 0 params~4 days

Empirical functions reward nonspecific bulk — a greasy ligand that sticks anywhere scores well everywhere. Ask instead how well this ligand scores elsewhere on the same protein, and score the difference:

specificity = E_site − E_best_elsewhere

Uniquely affordable here because the FFT already sweeps translations exhaustively; a sampling engine would pay for a second full docking run. Decoy-free and training-free.

Cost control
Run the counterfactual sweep on a coarse 1.0 Å grid with a reduced rotation set. You need a crude "does it stick elsewhere" signal, not a second accurate dock — this is what keeps the item from costing 10–50× a normal run.
Emits
Absolute gap, percentile rank of the true site among all surface patches, and the number of competing sites within 2 kcal/mol.
Risk
Large proteins with several genuine pockets will show a small gap for correct poses. Report the percentile alongside the gap so the ranker can learn the distinction rather than being misled by it.
W4 Torsion strain prior as a unary DP potential 0 params~2 days

Nothing in src/ looks like a strain library, and this one compounds with the pose-quantisation work directly.

Histogram observed torsion angles by bond type across the 35,349 train natives — no external dataset needed, you already hold it. Define strain(τ) = −ln p(τ | bond type), smoothed as a von Mises mixture so the potential and its derivative are continuous.

The payoff is structural: a per-torsion prior is a unary potential, and unary potentials are exactly what the tropical DP already handles exactly. It costs nothing in torsion_field.rs, it biases enumeration toward conformers that are chemically real, and PoseBusters explicitly checks internal energy — so top-1 and PB validity improve from the same table lookup.

Bond typing
Key on the two central atoms' AutoDock types plus heavy-atom degree and aromaticity. Coarse but sufficient; fall back to a pooled distribution when a type has < 50 observations.
Artifact
release/torsion_prior.bin, hashed into release/model-provenance.toml like every other committed binary.
Validation
Held-out natives must score lower strain than their own docked decoys. If they don't, the typing is too coarse — check before wiring it into the DP.
W5 Buried unsatisfied polar atoms 2 params~1 day

The cheapest item on the list and one of the strongest known discriminators of real binding modes. Empirical functions reward formed H-bonds and stay silent on the converse: desolvating a polar atom and giving it no partner costs several kcal/mol.

hbond_repair.rs already fixes donor/acceptor assignment — which matters here, since a receptor prepared without polar hydrogens presents every backbone nitrogen as an acceptor and would corrupt this count exactly as it corrupts the hbond term. Run the repair first, then count.

Definition
A ligand polar atom is unsatisfied when B(x) exceeds the burial threshold and no complementary partner sits within standard H-bond distance and angle.
Emits
Separate counts for ligand-side and receptor-side unsatisfied atoms — burying a receptor donor is a distinct and equally real cost.
Reuse
B(x) from W2. Sequence W2 first or share the enclosure computation.
W6 Fusion layer, diagnostic baseline, and calibrated confidence ~4 days

Three jobs, one model. LambdaMART over ~60 features, groups keyed by target.

  1. Fusion. Where W1–W5 cash out. Combine them with the existing Vina, Vinardo, CNN and ranker-v2 scores. The new features are orthogonal to what a pose scorer can perceive, and orthogonal signals fuse well — which the existing z-mean fusion already demonstrates in a cruder form.
  2. Diagnostic. Train once on physical features alone, without ranker-v2 output. The gap against ranker v2 measures how much its lift comes from the learned representation versus features a tree could read directly. Either answer is useful: a small gap says simplify, a large gap says the GNN is earning its keep.
  3. Pocket conditioning. Don't build a separate mixture-of-experts — feed pocket descriptors (buriedness distribution, polar/apolar ratio, metal presence, volume) as features and let the trees learn the conditioning. Same effect, one model, no extra machinery.
Protocol
Grouped 5-fold CV on train for both hyperparameters and early stopping. Val (68) as a smoke check. Test (225) once.
Determinism
Fixed seed, sorted feature order, no hash-map iteration in the funnel — the bit-identical invariant in CLAUDE.md applies to this path too.
Confidence
Emit P(top-1 ≤ 2 Å), calibrated by isotonic regression on CV out-of-folds. Report a reliability diagram and ECE, not just accuracy.

That last deliverable is the one with no competitor: no mainstream docking engine tells you when to distrust its answer, and W1's landscape entropy is a genuinely informative input to it.

Sequencing

Week 1S1 corpus (check for reuse first) · W5 · W4 histogram build. All independent of the gate work.
Week 2S2 feature schema · W2 hydration, starting with the "is the O map already W(x)?" check.
Week 3W1 landscape functionals — the flagship, and the one to protect time for.
Week 4W3 specificity sweep · W4 wired into the DP once the pose-quantisation gate reports.
Week 5W6 fusion, diagnostic, calibration. First and only test-set read.

W1, W2, W4 and W5 each produce features that stand alone, so ablate them individually before fusing. A feature that doesn't move CV top-1 on its own is unlikely to earn its complexity inside the ensemble, and finding that out early is cheaper than debugging it later.

Expectations, honestly

I can't promise numbers. What I can say is where the headroom demonstrably is: the pool holds a near-native pose for 97.5% of targets while top-1 selects one about 62% of the time, so roughly 35 points are lost purely to ranking. Recovering 5–10 of those across the whole programme would be a strong result and would put Rustina clearly ahead of qvina2+CNN on its own benchmark.

The items most likely to disappoint are W3, where multi-pocket proteins may wash out the signal, and the fusion step, where correlated features often deliver far less together than their individual ablations suggest. The items most likely to surprise on the upside are the pocket-normalised z in W1 — for screening rather than pose prediction — and W5, which is a day's work for a discriminator the empirical functions simply omit.

One structural risk worth stating plainly: several of these items reward burial and contact, and a ranker that learns to love buried poses will quietly degrade PoseBusters validity while top-1 improves. That is precisely why the hardcoded pb_valid has to be fixed before the first measurement, not after the last.