Golden path: input formats -> docking -> xtb -> DFT/QM-MM (2026-09-18)

Status of each claim is marked MEASURED (with source) or ESTIMATED (with reasoning). Two companion docs: golden-path-iteration-1.md (the survey that scoped this) and golden-path-qmmm-pipeline.md (the full QM/MM tier design).

STATUS AS OF 2026-09-19 — what has landed on main since this was written

C3 (transition-state search) is MERGED (#106) and now DEMONSTRATED under QM/MM embedding, not merely reachable. The TS cost model was also corrected from 2*6N + n_steps to 2*(6N+1) + (n_steps+1) (#110) -- a constant +3, so no conclusion moves, but the tables below are each 3 low.

This document was written against a main where several of its blockers were live. Ten PRs have merged since. VERIFIED against origin/main just now:

STILL OPEN, and these are the real remaining gaps:

• No saddle search. CLOSED AND MERGED 2026-09-19 (#106): ferric_scf::saddle::find_saddle, P-RFO with a Bofill update. C0-C5 is complete, and it is now WIRED TO AND DEMONSTRATED ON the QM/MM evaluator -- NH3 umbrella inversion under point-charge embedding converges in 5 steps to exactly one imaginary mode (qmmm_saddle_converges.rs), wired to the QM/MM evaluator. IRC LANDED 2026-09-19: irc::follow_irc walks mass-weighted steepest descent off the saddle in both directions, so "which two minima does this connect?" is now answerable. NH3 inversion lands at +0.8044/-0.8044 Bohr pyramidalisation -- opposite sides, degenerate, barrier 0.016078 Ha both ways. Runs under QM/MM embedding through the same closure find_saddle takes, so the path and the search cannot disagree about the field. Section 4 has the full scope.
• No analytic dispersion gradients. CLOSED 2026-09-19 (#99). ferric_d3::d3bj_gradient is implemented and validated against finite difference to 1e-12..1e-13 on four systems, with the step-size scan showing textbook O(h^2) convergence. task="optimize" and with_gradient=True now WORK; MEASURED end to end on Ar2/PBE/STO-3G from a short start, 6.7917 Bohr uncorrected vs 6.6418 with d3bj -- the attraction shortens the bond, which is the direction that shows the correction reached the GRADIENT and not only the energy. task="frequencies" and run_frequencies(xc=..., dispersion=...) also WORK on a closed-shell KS reference: the Hessian is the finite difference of the KS + dispersion gradient, validated against second differences of the D3(BJ) energy and of the full SCF + MBD@rsSCS energy (see Validation). Open-shell references are refused.
• QM/MM dispersion. D3/D4/XDM/VV10 are all QM-atom-pairwise; MM point charges carry none. Needs LJ terms in ferric-mm.
• Pose noise. MEASURED per-pose sd 29.07 kcal/mol against 1-2 kcal/mol substituent effects, and funnel.py keys one row per MOLECULE so it cannot express an ensemble. This gates the QM tier, not the cheap tiers.
• Site dependence. MEASURED within/between ratio 0.94-0.95 on two independent constructions: WHERE a group goes matters as much as WHICH group. The pipeline's unit must be the (substituent, SITE) pair.

0a. THE ANSWER TABLE -- start here

Every row was EXECUTED on origin/main on 2026-09-19 to produce the cost in it; a row nobody could run is not in this table. Costs are one call, warm process, single-threaded, on this box -- read them as orders of magnitude, and see "the first call costs 24x" below before planning a campaign from them.

Tier 1 is now MEASURED here, not only cited. Every previous pass recorded that docking could not be re-measured on this box because mk_prepare_receptor.py was "missing". It was not -- it ships in the venv's bin/ and was simply not on PATH, so the check that found it absent was a PATH artifact rather than a missing dependency. Run 2026-09-20 against the 7LCJ pocket at exhaustiveness 4:

prepare_receptor itself is 3.1 s, once per receptor. The 26.4 s/ligand figure from RESULTS.md M11 sits inside this curve at drug scale and is CONFIRMED independently. Cost scales ~N^1.5 in atom count, so quoting one number hides a 16x spread across the sizes a campaign actually sees -- the row now gives three points.

pocket_polarization is the plot for a SINGLE pose (added 2026-09-20). site_substituent_heatmap ranks substituents across sites and says nothing about one calculation. compute_binding_energy also returns charges_vacuum and charges_field, and nothing plotted them -- so the one output that distinguishes an embedded result from a number went unlooked at. dq = q_field - q_vacuum is the pocket pushing electrons around the ligand: MEASURED on danuglipron/7LCJ, max |dq| = 0.037 e, summing to zero to 1e-13. The plot annotates that sum, because a NONZERO total means the two SCFs were not the same molecule and the interaction energy alone cannot show it.

binding_energy is TWO SCFs, not one tier. RUN end to end 2026-09-20: danuglipron's cryo-EM pose (71 atoms) in the 7LCJ pocket (6458 charges), STO-3G/RHF, 137 s wall, giving delta_e_kcal_mol = -17.41. The row used to say "tier 3/4 above", which is a pointer rather than a cost and does not add up from the tier rows: the call derives the pocket charges (2.66 s) and then runs the ligand SCF TWICE, embedded and in vacuum. At def2-SVP, scale by the tier-4 basis factor.

The number is for ONE POSE and does not rank analogues -- compute_binding_energy's own docstring carries the five closed protocols and the 4.07 kcal/mol ddE noise floor against 1-2 kcal/mol substituent effects. Pass that floor to site_substituent_heatmap(noise_floor=).

The DFT row is quoted at the DEFAULT basis now. tier4_dft's default is def2-svp, not STO-3G, and at 9 atoms that is 2.64 s vs 0.75 -- 3.5x (re-measured 2026-09-20, n=3, tight). Every STO-3G figure in tiers.py is labelled as such and none is wrong, but the leading number was the one a reader budgets with, and it was not the one a default call costs. Worth recording how this nearly went the other way: a bare run_dft(mol, sto3g, functional="pbe") takes 0.69 s while tier4_dft on the same molecule takes 2.64, which reads as ~2 s of wrapper overhead. Profiling put 2.667 of 2.667 s inside run_dft -- there is no wrapper cost at all, the two calls simply used different bases.

assess_smiles reaches the NETWORK by default. include_web=True adds the admetlab3 and protox3 providers, measured at 1.6 s/analogue against 54 ms for the offline rdkit-alerts path alone. Which one you pay decides whether a 1000-analogue sweep takes a minute or half an hour, and whether it runs offline at all.

The substitution row was RE-MEASURED 2026-09-20 (benzoic acid, {F, Cl}, the configuration that gives exactly 7 proposals): 248 ms on the FIRST call in a process, 7.6 ms on every call after it. The old flat "216 ms" was a cold-start number quoted as a per-call cost, which over-states a sweep of 100 analogues by ~30x. Both are given because both are real and they answer different questions -- budget the first call once, the rest at 7.6 ms.

The two tier rows were RE-MEASURED 2026-09-20 through the tier functions themselves (aspirin, 21 atoms, n=9 after a warm-up call): FF median 9.0 ms (min 8.8, max 10.3), GFN2 median 39.0 ms (min 35.4, max 44.0). They had read 32 ms and 53 ms with no provenance, and the FF figure was 3.5x high -- tiers.py measures 8.2 ms at 19 atoms and 21.6 at 34, which cannot bracket 32 ms at 21. Both rows now carry the size breakdown, because a single number hides that tier 2 and tier 3 scale differently with atom count. test_the_answer_table_tier_costs_agree_with_tiers_py pins them to the source.

Two rows carry a caveat that outweighs their cost.

Ranking by binding energy is not licensed: all four pose protocols are closed and the best available ddE noise is ~4.07 kcal/mol against effects of 1-2 (RESULTS.md M4-M14). site_substituent_heatmap(noise_floor=...) greys out every cell inside that limit precisely so a figure cannot imply otherwise. A paired estimator (M17) measures 0.221-0.615 in gas-phase MMFF and is PROVISIONAL -- untested in a pocket.

Toxicology discriminates between MOTIFS, not between substituents that leave the motif alone. MEASURED liability scores: benzene 0.000, aspirin 0.083 (phenol ester), nitroaromatic 0.167, catechol 0.208 (PAINS + NIH + BMS), Michael acceptor 0.250 -- a chemically sensible ordering. But a halogen scan around an unchanged scaffold gives every analogue the SAME score as the parent (13 aspirin proposals, all 0.0833), because the alerts are substructure matches and F/Cl/Me on a ring do not hit one. That is correct behaviour and it bounds the use: the gate REMOVES a liability-bearing motif from the set; it does not order the survivors. Every endpoint says so in its own note -- "a rank-only liability density, NOT a probability of toxicity".

xtb is a gate, not a ranker: Spearman 0.011 against DFT over a 3 kcal/mol span (M16, n=20, 95% CI [-0.434, +0.451]). Use it to separate the anion from the neutral, not to order two conformers.

The last two rows are the ones a catalyst user reaches for, and both exist because a number alone was not enough. run_optimize reports converged, steps and a final energy, which cannot tell "ran out of steps near a minimum" from "walked uphill and oscillated" -- optimization_trace draws the climb with its magnitude. And min_link_to_charge_distance() says a frontier is bad without saying WHERE: qmmm_partition shows the cut, the link atoms and the offending charge. Pass it expect_min_angstrom= -- it cross-checks your coordinates against that accessor and refuses a mismatch, because link_atom_positions() is ANGSTROM while point_charges() is BOHR and converting both errs by 1.89x in the SAFE direction.

"Has a plot" means the plot answers THAT question, not that a figure exists. The transition-state row is the example: a TS search produces an imaginary MODE (a 3N vector), and imaginary_mode shows whether it displaces the reacting atoms -- which is the second, non-optional half of confirming a saddle.

0. The headline defect: the docked pose is thrown away — FIXED 2026-09-19

RESOLVED. tools/pipeline/funnel.py::_harvest_geometry writes the key and is called from the driver (funnel.py:217). Chain re-verified on origin/main: tier1_dock returns symbols + coords_angstrom in its payload, _harvest_geometry writes context["geometry"][id], tiers._embedded reads it. The docked pose now reaches tiers 3 and 4.

Two things from the diagnosis are worth keeping, which is why the original write-up is left below rather than deleted.

FIRST, WHERE the write had to go. _run_stage dispatches through a ProcessPoolExecutor, so a tier that mutates context mutates a PER-WORKER COPY: the write is lost under the parallel path while appearing to work serially. Harvesting from the returned results in the driver is the only placement that holds for both. That generalises to any state a tier tries to pass forward.

SECOND, THE GREP BELOW CANNOT SEE ITS OWN FIX. Re-run today it still reports one read and no writes, because the write is context.setdefault("geometry", {})[...] and setdefault matches none of the three alternatives in that pattern. A grep that would have to be rewritten to notice the bug being fixed is not a verification you can re-run -- and this one was labelled "VERIFIED by grep, both directions".

The original 2026-09-18 write-up follows.

VERIFIED by grep, both directions, on 2026-09-18:

$ grep -rnE '\["geometry"\]|\.get\("geometry"|"geometry":' tools/ experiments/ --include=*.py | grep -v test
tools/pipeline/tiers.py:84:    cached = context.get("geometry", {}).get(iso.canonical)

One read, zero writes -- AS OF THE SURVEY. Fixed since; see the heading. funnel.py:184 now writes context.setdefault("geometry", {})[...] from _harvest_geometry. The paragraphs below describe the DEFECT as found, because the reasoning is what makes the fix legible; they are not the current state.

The consequence is concrete. tier1_dock DOES return the pose (tiers.py:181-185, payload carries symbols and coords_angstrom), and _embedded DOES look for a cached geometry (tiers.py:84) -- but no code connects the two. So _embedded falls through to tier2_forcefield, which re-embeds from SMILES with ETKDG in free solution.

Tier 1 costs 26.4 s/ligand (RESULTS.md M11, exhaustiveness 4). That entire spend was discarded, and tiers 3 and 4 then scored a gas-phase conformer that had never seen the pocket.

(Cost figures here cite RESULTS.md, not a file.py:NN line: a citation that resolves to a module doc comment is not a measurement, and the cost table below says so at length.) This violates the funnel's own stated premise -- _embedded's docstring says the cache exists so "tiers 3 and 4 would not be scoring DIFFERENT geometries of the same candidate".

This is independent of QM/MM. It is the single highest-value fix in the pipeline and should land before any new tier is added, because every downstream number is currently computed on the wrong geometry.

VERIFIED end to end, with and without the fix (2026-09-19)

Ran the REAL run_funnel with a geometry-producing tier 1 and a tier-4 stand-in that calls tiers._embedded -- the actual accessor tiers 3 and 4 use, not a proxy for it. The docked pose was set to an unmistakable [(9.9,9.9,9.9), (8.8,8.8,8.8)].

The pre-fix arm is the important half: it is not that the pose arrives slightly perturbed, it is that a DIFFERENT MOLECULE GEOMETRY arrives, re-embedded from SMILES with a different atom count. Anything computed downstream of that describes a conformer the pocket never saw.

This is the A/B that makes the defect concrete rather than inferred from a grep, and it also confirms the fix is not vacuous -- the no-fix arm genuinely fails.

Where the fix goes (not inside a tier)

_run_stage (funnel.py:111-115) dispatches through a ProcessPoolExecutor, so a tier mutating context in a worker mutates a COPY -- the change never returns. Any fix that writes context["geometry"] from inside tier1_dock will silently do nothing under the parallel path while appearing to work serially.

The results DO return to the driver (funnel.py:158), so the harvest belongs in run_funnel's stage loop, between stages:

results = _run_stage(stage, population, context)
# harvest any geometry a tier produced, so later tiers reuse it
for r in results:
    if r.ok and r.payload and "coords" in r.payload:
        context.setdefault("geometry", {})[r.candidate_id] = {
            "symbols": r.payload["symbols"], "coords": r.payload["coords"],
        }

Anchor test before changing anything: assert that with tier 1 in the stack, tier 3 receives the DOCKED coordinates, not ETKDG's. The trivial-limit anchor is a funnel with no docking stage, where behaviour must be byte-identical to today.

0b. QUICKSTART -- the pipeline in code you can paste

Everything below this section describes the pipeline and costs it. This one RUNS it. Every snippet was executed on 2026-09-19 and its real output is shown; none is illustrative.

# A. any input format -> a ferric Molecule (Angstrom, seeded ETKDG for SMILES)
from tools.structure import from_smiles, read
mol = from_smiles("CC(=O)Oc1ccccc1C(=O)O", seed=0xF00D)   # aspirin: 21 atoms
mol = read("ligand_with_hydrogens.pdb")   # | .sdf | .mol2 | .xyz | .pqr | .gro
#   ^ a PLACEHOLDER path -- substitute your own file. Everything below runs
#     as written; this line is the only one that needs editing.
#   `read` RETURNS A MOLECULE. `read_structure` does NOT -- it stops at a
#   `Structure` (symbols/coords/charge/multiplicity/source) and never imports
#   ferric, which is what you want when inspecting a file without pulling in
#   the extension. Call `.to_molecule()` on it if you need one.
#
#   This doc named `read_structure` here and called it "the same entry point"
#   as `from_smiles`. It is not: pasting that line gives a TypeError on the
#   next `.symbols()`, because a Structure exposes symbols as a FIELD and a
#   Molecule as a METHOD. Verified 2026-09-19 by running it.
#
#   THE PDB MUST ALREADY HAVE EXPLICIT HYDROGENS. An ordinary
#   crystallographic PDB does not, and `read` refuses it with a
#   StructureError -- it does not protonate. That is correct (a species
#   missing its hydrogens is not the molecule), but it means "x.pdb" is not
#   a generic example, hence the filename above.
#
#   A RECEPTOR is the separate case: it goes through derive_pocket_charges,
#   which runs pdb2pqr and DOES protonate. That does not change `read`.

# B. enumerate analogues at every matching site
from tools.pipeline.substitution import propose_substitutions, relative_descriptors
parent = "c1ccccc1C(=O)O"
props = propose_substitutions(parent, {"F": "F", "Cl": "Cl", "Me": "C"})
#   -> 10 proposals, and props[0] is the PARENT. That is the anchor, not a bug:
#      the parent rides through every stage so scores can be reported as ddE.

# C. gate RELATIVE to the parent, never on absolutes
relative_descriptors(props[0].smiles, "c1ccccc1C(=O)O")
#   -> (-1.42e-14, 0.0, 0.0)   the parent against itself -- NOT exactly zero.
#      The distinction is real and worth knowing:
#        relative_descriptors(s, s)                     -> exactly (0.0,0.0,0.0)
#        relative_descriptors(props[0].smiles, s)       -> -1.42e-14 in dMW
#      because props[0].smiles is the CANONICAL form ("O=C(O)c1ccccc1") of the
#      input spelling ("c1ccccc1C(=O)O"). Same molecule, different atom order,
#      so the MW float sum lands one ulp apart. Compare ddE against a
#      tolerance, never `== 0.0`, whenever either side has been round-tripped
#      through a canonicaliser.

# D. liability flags (published alert sets -- NOT a probability of harm)
from tools.tox.alerts import RdkitAlertsProvider
provider = RdkitAlertsProvider()          # build ONCE: 47 ms vs 9.4 ms/molecule
eps = provider.fetch("CC(=O)Oc1ccccc1C(=O)O")
#   -> 13 endpoints, e.g. alert_brenk = 0.333

# E. the pocket half, composed -- all four tiers share ONE signature
#    (iso, context) -> TierResult, so run_funnel chains them.
from tools.pipeline import run_funnel, Stage
from tools.pipeline.tiers import tier1_dock, tier2_forcefield, tier3_gfn2, tier4_dft
from tools.campaign.hierarchy import Tier

stages = [
    Stage(Tier.FORCE_FIELD,   tier2_forcefield, keep=2, name="ff"),
    Stage(Tier.SEMIEMPIRICAL, tier3_gfn2,       keep=2, name="xtb"),
    Stage(Tier.QUANTUM,       tier4_dft,        keep=1, name="dft"),
]
# the funnel takes Isomers, and section B produced SubstitutionProposals --
# this conversion is the one line between them.
from tools.isomers.model import Isomer

candidates = [
    Isomer(
        smiles=p.smiles,
        kind="substitution",
        transform=p.label,
        parent_smiles=parent,
    )
    for p in props
]
rep = run_funnel(candidates, stages, {"seed": 0xF00D, "basis": "sto-3g"})
#   tier 1 is omitted above only because it needs the `docking` extra and a
#   receptor; add Stage(Tier.EMPIRICAL, tier1_dock, ...) with
#   context["receptor_pdbqt"] and ["box_center"] to run it.

MEASURED 2026-09-19, two small candidates, STO-3G, one process:

funnel wall: 1.58 s
  FORCE_FIELD    in=2 out=2 failed=0   0.03 s
  SEMIEMPIRICAL  in=2 out=2 failed=0   0.04 s
  QUANTUM        in=2 out=1 failed=0   1.51 s     <- 96% of the wall
survivors: ['CC(=O)O']   dft = -225.76133078 Ha

96% of the wall in the last tier on TWO candidates is the funnel's whole argument in one line, and it gets worse with candidate count: the cheap tiers scale with the population, tier 4 scales with what reaches it. That is why keep= matters more than any per-call cost in this note.

Note the DFT total (-225.76133078) is 4.18 mHa BELOW the SCF energy the ladder logs (-225.7571497). That difference is the D3(BJ) correction, -2.62 kcal/mol -- tier4_dft passes dispersion="d3bj" by default. It is a real contribution at chemical-accuracy scale, and it is why run_end.energy and a method's result.total are different numbers.

For the pocket half (docking, xtb, DFT) the entry points are tools.docking.vina_dock.dock_ligand, tools.campaign.fit.pose_fit and tools.active_site.binding_energy.compute_binding_energy; the composition is tools.pipeline.run_funnel, and tools/pipeline/tests/test_golden_path_smoke.py runs two real tiers through it end to end in about 7 seconds.

Before you rank anything with the output, read "WHICH METHOD FOR WHICH QUESTION" below. The pipeline will happily produce a ddE ordering that the measurements do not support -- five pose protocols have been tried and the best available noise is 4.07 kcal/mol against substituent effects of 1-2.

1. Input formats: what can enter

Rust Molecule has only load_xyz / load_xyz_with_charge / parse_xyz (crates/ferric-core/src/mol.rs). Everything else is Python-side.

tools/structure (added 2026-09-18, PR #91) funnels every format through Molecule.from_xyz_string, so a PDB and an XYZ of the same geometry give bit-identical Bohr coordinates rather than drifting through two parsers.

What each entry path COSTS (MEASURED 2026-09-19)

The format table above says what can enter; these say what it costs. Min of 1-3 reps, single-threaded, on this box:

The remaining five paths, measured 2026-09-19 so the table covers what section 1 says can enter rather than only the three that had been timed. All one molecule (aspirin, 21 atoms) so the number is the PARSER and not the size, min of 5, single-threaded:

Every file-parsing path is 0.1-0.3 ms and none of them matters. The only entry cost worth planning around is SMILES (~25-80x a file read, because ETKDG generates a conformer rather than reading one) and the one-off PDB->PQR at 2.66 s. A campaign that reads 1000 SDF files spends 0.3 s total on parsing.

THE FIRST from_smiles CALL COSTS 24x THE REST. MEASURED on aspirin:

first call   214.66 ms
min of 7       9.38 ms
median         9.75 ms

RDKit import plus ETKDG warm-up, paid once per PROCESS. The 8.5 ms in the table above is steady state and is right (re-measured: ethanol 2.38 ms, aspirin 8.52 ms) -- but a ONE-MOLECULE run pays 215 ms, not 8.5, and a per-molecule subprocess pays it every time.

This is why the per-item costs in this note are the wrong unit for planning a campaign: 1000 molecules in one process is 9.6 s of embedding, and 1000 subprocesses is 215 s. The difference is entirely warm-up, and it does not appear in any per-call figure.

END-TO-END through the real tier functions, aspirin (21 atoms), one process:

SMILES -> 3D (first call)   205.76 ms   <- warm-up dominates
xyz -> Molecule               0.21 ms
tier 2 MMFF                  24.67 ms
tier 3 GFN2-xTB              50.26 ms

Tier 2 at 24.67 ms sits between the table's 21.6 ms @ 34 atoms and the ~73 ms projected @ 71, so the projection holds at this size. Tier 3 at 50.26 ms is within the 0.05-0.152 s band.

Note the pdb row: a PDB that already has hydrogens reads at file speed. It is the crystal PDB with no hydrogens that is refused -- see the refusal note below, which is about protonation, not about the format being slow.

Four orders of magnitude separate them, and the ordering is the point: the PDB path is ~300x the SMILES path and ~9000x an xyz read. It is also a ONE-OFF per target -- PocketCharges is derived once and reused across the whole ensemble (that is the reason the type exists), so 2.66 s amortises to nothing over a 1000-analogue campaign and is a real cost for a one-molecule run.

Note from_smiles at 8.5 ms here vs 214 ms/proposal for embed_proposals in the cheap-stage table: same ETKDG machinery, ~25x apart, because a drug-sized analogue is far harder to embed than the small test molecule timed here. Quote the 214 ms for campaign planning.

A refusal worth knowing about before you hit it. read_structure REJECTS a crystal PDB with StructureError: ... no hydrogens. That is correct -- a PDB from the PDB has no hydrogens, and silently treating it as a QM molecule would hand the solver a species that does not exist. A receptor goes through derive_pocket_charges (which runs pdb2pqr and protonates), not through read_structure. The error names the problem, but the two paths are easy to confuse on first use.

PQR and the chain-ID column: a non-bug, CHECKED (2026-09-19). pqr_parser hard-requires exactly 10 whitespace fields and reads coordinates at fields[5:9]. PQR files that carry a chain ID have ELEVEN fields and shift the coordinates to fields[6:10], so such a file is rejected with Unexpected PQR field count (11, expected 10) -- and both repo fixtures are 10-field, so the 11-field layout is untested.

That looks like a gap and is not one. MEASURED: pdb2pqr30 drops the chain ID, emitting 10 fields even from a PDB whose ATOM records carry chain A. Fed a chain-bearing PDB through run_pdb2pqr, all 16 output records came back 10-field. The parser matches its only producer in this pipeline, and the 11-field layout is not reachable through derive_pocket_charges.

It IS reachable if someone hands you a PQR from another tool (APBS's own writers, some Amber paths). The failure is then a clean error naming the field count, not a silent misparse -- coordinates read from the wrong columns would be far worse. Pinned by test_an_eleven_field_pqr_is_refused_not_misparsed.

Unit hazard, worth stating once: Python geometry entry is Angstrom; point_charges and QmmmSystem.point_charges() are BOHR; PocketCharges holds Bohr. Mixing them is a silent 1.89x error, not a crash.

Charge and multiplicity are never inferred. No common structure format records multiplicity (2S+1) at all. ferric's parity check catches an odd-electron species left at the default singlet, but CANNOT catch an even-electron triplet (O2, many carbenes) -- that converges cleanly to a state that does not exist.

1b. Dependency status, CHECKED rather than assumed (2026-09-19)

Every stage below names a tool. Whether those tools are actually present is a separate question from whether the code is written, and it is the one that decides if a stage runs today.

ONE FALSE ALARM WORTH RECORDING, because the same check will mislead the next person: import pdb2pqr FAILS (ModuleNotFoundError) while pdb2pqr30 is on PATH. That is not a broken install -- pdb2pqr_runner.py is a deliberate SUBPROCESS wrapper around the CLI ("Subprocess wrapper around the PDB2PQR CLI", line 1), so the module never needs to import. Checking importability would have reported a working stage as broken.

So the docking branch's tooling is complete. The blockers listed elsewhere in this document are about PHYSICS and DATA FLOW (missing dispersion, pose noise), not about missing software.

2. Cost estimates

WHICH METHOD FOR WHICH QUESTION (the table this note was missing)

Everything below this heading costs methods. This one says which to REACH FOR, and -- more usefully -- what resolution each can actually deliver on this campaign. A method that is cheap and cannot answer your question is not a bargain.

The row that matters most is the one with no method. Four pose protocols have been measured (RESULTS.md M4-M13) and the best available ddE noise is 4.07 kcal/mol against substituent effects of 1-2. Every OTHER row in this table is a green light; that one is not, and no amount of tier-4 DFT fixes it, because the error is in the pose ensemble and not the electronic structure.

Two rows are worth reading together. "Where does this ligand sit?" is a green light and "which pose is best?" is a red one, from the SAME tool. Docking generates the right answer among its candidates and cannot pick it out. That is not a defect to fix -- it is the empirical reason tiers 2-4 exist.

Embedding is free only while the charge set is small (CORRECTED 2026-09-19)

The row above read "~1.0x, embedding is essentially free". That is true for a handful of charges and FALSE for a real pocket. MEASURED, benzene/STO-3G against the 7LCJ pocket (derive_pocket_charges, 6458 charges, net -1.000 e):

(The sub-1.0 ratios at 10-100 are SCF iteration-count noise on a 0.2 s baseline, not a speed-up from adding charges. A separate warm run of the same pair gave 5.61x at 6458, so read the large-N cost as 3-6x rather than a single figure.)

So both numbers are right about different things, and the old row generalized the small one. A whole pocket is a real cost. The PDB->charges step itself is 2.46 s for 7LCJ and is paid once, not per candidate.

AND THE ANALOGUE MUST BE IN THE POCKET, which does not happen by itself. embed_proposals returns ETKDG conformers centred on the ORIGIN; the pocket sits at its crystal coordinates. MEASURED on 7LCJ: the embedded analogue's centroid is (0.00, 0.00, 0.00) A and the pocket's is (124.3, 148.3, 116.9) -- 226 A apart. Feeding those coordinates straight to an embedded SCF is not an error, it is a confident dE of -0.001 to +0.005 kcal/mol, i.e. a gas-phase answer wearing a QM/MM label.

This is the concrete shape of "harvest the docked pose" (section 0). The chain SMILES -> propose_substitutions -> embed_proposals -> run_rhf runs end to end and is WRONG without a placement step between the embed and the score. A dE of essentially zero against a charged pocket is the tell -- see the single-charge control under G3.

THERE ARE TWO PATHS TO A GEOMETRY AND ONLY ONE IS PLACED. Verified in the source, because the two look interchangeable from a call site:

A DOCKED POSE IS UNITED-ATOM, so it is not a QM geometry as it stands. Vina merges nonpolar hydrogens into their carbons: MEASURED, aspirin docks as 14 atoms where 21 went in. tier1_dock now re-hydrogenates the pose via united_atom.restore_hydrogens before handing it on, and fails the candidate if it cannot -- because the failure downstream is silent in the worst place:

AUDITED afterwards, because the obvious question is whether any other hop does this: it does not. read_structure (.sdf/.pdb/.xyz), from_smiles, tier2_forcefield and embed_proposals all preserve every atom, pinned by test_no_geometry_hop_silently_changes_the_MOLECULE. Docking was the only one, and only because PDBQT is a united-atom format.

funnel._harvest_geometry exists to carry the first into context["geometry"] so tiers 3 and 4 score the DOCKED pose instead of re-embedding. Use the funnel, or take coords_angstrom off the pose yourself. The embed path is for enumeration and gas-phase work; it is not a substitute for docking. embed_proposals' docstring SAYS so, with the 226 A figure.

But DO NOT truncate with a naive distance cut. That was the obvious next move and it is measured here because it does not work. Keeping charges within r of the probe, 7LCJ, water/STO-3G, error against the full 6458-charge answer:

The error is NOT monotone in r -- 1.82 -> -0.46 -> 0.04 -> 0.06 -> 0.92 -> 0.76 -> 0.88 -- so "use a bigger radius" does not buy accuracy, and a 15 A cut is worse than a 10 A one. A sphere through a protein also cuts residues in half: the net charge kept wanders from +0.02 to +3.55 e against the full pocket's clean -1.000, and a spurious monopole is exactly the kind of error that does not decay with distance.

What is NOT established: that the error is CAUSED by the net charge. Spearman over these 7 points gives rho = -0.07 for |err| vs |net q| and +0.04 vs radius, and at n=7 the smallest detectable |rho| is ~0.75 -- neither comes close. The non-monotonicity and the charge wander are both MEASURED; the link between them is a hypothesis this sweep cannot test. A truncation scheme that cuts on whole RESIDUES (keeping each one neutral) is the standard fix and is untested here.

Until then: use the whole pocket and pay the 3-6x, or validate your own cut against it. err vs full is cheap to compute -- one extra SCF.

WHICH FUNCTIONAL AND BASIS (the other half of "which method")

The table above says which METHOD KIND to reach for. It never said which FUNCTIONAL or BASIS, which is the choice a tier-4 user actually makes -- tiers.py defaults to PBE/def2-svp and nothing explained why or when to depart from it.

Grades below are from wiki/VALIDATION.md, which is the authority; the worst measured error against PySCF is quoted rather than a tolerance, because a guard band says what a test permits and not what the code does.

All five are validated on cc-pVDZ and def2-SVP only, four small molecules. Larger systems and other bases are unverified, which is a scope limit and not a prediction of failure.

The gradient story is narrower than the energy story, and it is the gradient that a geometry optimization or a saddle search depends on:

• meta-GGA gradients are s/p-shell only -- ferric's AO Hessians do not go past 6-31G, so a SCAN optimization at cc-pVDZ is out of reach. d shells carry a ~6e-5 residual that the GGA path SHARES, so it is an AO-Hessian limit rather than a meta-GGA one.
• there is no meta-GGA f_xc Newton kernel, so those SCFs fall back to DIIS.
• there is no analytic Hessian for anything. Every Hessian in this note is finite-differenced at 6N+1 gradients, which is what makes the TS and IRC budgets what they are.

Basis costs ~10x, and almost every timing in this note is STO-3G. MEASURED 2026-09-19, RHF wall time on one process:

The EXPONENT is nearly unchanged -- tail-fitted 4.69 (STO-3G) against 4.10 (def2-SVP), global 3.48 against 4.04 -- so the basis is a near-constant MULTIPLIER over this range, not a steeper curve. That makes the correction easy and reusable: multiply any STO-3G figure in this note by ~10 to get what def2-svp costs, and keep the scaling exponent.

It matters because def2-svp is the tier-4 DEFAULT while STO-3G is what the end-to-end timings use. A campaign budgeted from the STO-3G rows is budgeted an order of magnitude light.

Basis. def2-svp is the tier-4 default and the larger of the two validated sets. STO-3G appears throughout this note because it is what the end-to-end timings use -- it is a demonstration basis, NOT a production one, and no number measured at STO-3G should be read as an accuracy claim.

ECP. def2-ECP is Proven (narrow) for RHF -- Xe 2e-12 Ha, I2 1.2e-6 -- so heavy elements are reachable, on three systems and one basis.

MEASURED (quoted with source)

Budget from the ex=4 row, not the ex=32 one. M11 measured that across an 8x range of exhaustiveness the mean redock RMSD moved 0.097 A -- SMALLER than the 0.131 A between-seed SEM -- and ex=32 had the WORST mean of the four levels tried. So ex=32 costs 6.8x for no accuracy, and the ~2 min row is kept only because tiers.py:11 still documents it.

Note the parallel efficiency while you are here: giving up 11 of 12 cores costs 4.1x, so Vina's internal parallelism runs at 34%. Fan-out across ligands only wins above ~4 workers, which is a narrower claim than "the machine sits idle at low exhaustiveness". | 2 | MMFF94 (tier2_forcefield = embed + optimize) | 2.2 ms @ 9 atoms, 8.2 @ 19, 21.6 @ 34; ~73 ms projected @ 71 | MEASURED 2026-09-19 | | 2 | MMFF94 ~1 ms/pose | superseded: that figure cited tiers.py:12, which is the same doc comment -- a circular citation, never a measurement, and it described a single point rather than embed+optimize | | | 3 | GFN2-xTB via tier3_gfn2 | 0.05-0.15 s @ 9-19 atoms | MEASURED 2026-09-19 | | 4 | ferric DFT via tier4_dft | 0.66 s @ 9, 8.7 s @ 19 (STO-3G); 96.1 s @ 32 (def2-SVP, ~450 bf) | MEASURED 2026-09-19 / RESULTS.md | | 4 | ferric DFT | 612.4 s @ 71 atoms, STO-3G/PBE (~234 bf), 18 iters, converged | RESULTS.md |

DO NOT DERIVE A SCALING LAW FROM THOSE TWO ROWS. They differ in BASIS as well as size, and in the unhelpful direction: the BIGGER system used the SMALLER basis. Fitting them gives p = 2.32, which UNDERSTATES pure N-scaling because part of the size increase was paid for by a cheaper basis.

Measured directly, and the confound turns out to be small. A DFT N-sweep at FIXED basis (PBE/STO-3G, alkanes C2-C4, 2026-09-19) gives a tail exponent of 2.32 -- exactly what the confounded pair gives. Two independent routes to the same number:

The last row is a DIFFERENT KIND of point and is here to stop a future session misreading it. Measured through tier4_dft: 0.25 s @ 3 atoms, 1.02 @ 6, 2.68 @ 9 -- so the table's "1.0 s @ 6 atoms" is exact. The OVERALL 3->9 exponent is 2.16, which looks like it undercuts the 2.3-2.6 band. It does not: these sizes are far below the 9-71 range the band is fitted on, and fitting the TAIL (6->9) gives 2.38, inside it. Averaging in the flat small-N start pulls the exponent DOWN, the same artefact this document flags for the QM-radius sweep below ("gives 2.51 by averaging in the flat small-R start").

So p ~ 2.3-2.6 is the honest band, and the basis confound moves the answer less than the RHF-vs-DFT difference does. Still quote the fixed-basis numbers, because the agreement is what makes the confounded pair usable rather than the other way round -- and note all three are 2-3 point tail fits, which is indicative, not decisive.

One measurement artifact worth naming: C1 came in at 5.30 s against C2's 0.37 s -- a first-call warm-up (basis parse, grid construction), not chemistry. Tail-fitting excludes it automatically, which is the reason the repo's protocol says to fit the tail rather than the whole series.

ESTIMATED (reasoning stated, do not quote as measured)

No QM/MM wall-time measurement exists anywhere in the repo -- the 9 Rust and 29 Python qmmm tests are PySCF CORRECTNESS tests on H2/ethane/water with no timings. So:

• Tier 3.5 (in a pocket field): ~2x a gas-phase single point. RETRACTED 2026-09-18 -- MEASURED at ~1.0x, so the estimate was wrong by about 2x. See the measurement immediately below. I had reasoned "MM charges add a one-body term" and then guessed 2x anyway; the same reasoning actually PREDICTS a small overhead, because a one-body term is cheap. The guess did not follow from the argument I gave for it. QUALIFIED 2026-09-19: that ~1.0x holds up to ~1000 charges. A WHOLE pocket (7LCJ, 6458 charges) costs 3.4-5.6x -- the one-body term is cheap PER CHARGE and there are thousands of them. See "Embedding is free only while the charge set is small". The original retraction stands at the scale it was measured; it just does not generalize to an untruncated pocket.

• Tier 5 (QM/MM DFT optimization): ~17-35 h in a real pocket (a ~5.8 h gas-phase floor x the 3-6x embedding cost); was "1.5-17 h". The 10x width came entirely from an unmeasured BFGS step count. MEASURED 2026-09-18: ethane with a link atom across the C-C cut (the catalyst shape) converges in 34 steps -- and in 34 at BOTH STO-3G and 6-31G, so the step count is set by the GEOMETRY, not the basis, as BFGS theory predicts. 34 x the 612 s single point = ~5.8 h. Still ESTIMATED, because the MULTIPLICAND is a 71-atom DFT single point measured on a different system; only the multiplier is now measured. CAVEAT, and it is the load-bearing one: 34 steps is a near-rigid CH3 relaxing a few bond lengths. A floppy ligand in a pocket has far more soft degrees of freedom and will take more. Treat 34 as a FLOOR for the step count, not a typical value -- one geometry is not a distribution. SECOND CAVEAT, added 2026-09-19: the 612 s multiplicand is a GAS-PHASE single point. A QM/MM optimization pays the embedding cost at EVERY step, and with a whole pocket that is not free. MEASURED, benzene/STO-3G against the 7LCJ pocket (6458 charges), same optimizer, same convergence:

  vacuum       4 steps, 0.278 s/step
  full pocket  4 steps, 1.574 s/step     -> 5.7x
  

  The STEP COUNT is unchanged, so the two factors MULTIPLY rather than trade off: ~5.8 h becomes ~17-35 h at 3-6x. Both caveats push the same way, so read 5.8 h as a hard floor built from a gas-phase multiplicand and a near-rigid multiplier -- not as an estimate of a real catalyst job.

  THIRD CAVEAT, and it cuts the other way (2026-09-20): the 3-6x embedding factor DOES NOT TRANSFER to a large QM region. It was measured on benzene, 12 atoms. MEASURED on danuglipron, 71 atoms, same 6458-charge pocket, RHF/STO-3G:

  system	atoms	vacuum	embedded	ratio	ABSOLUTE overhead
  benzene	12	0.278 s	1.574 s	5.7x	1.30 s
  danuglipron	71	43.1 s	49.1 s	1.14x	6.0 s

  The pocket contributes roughly FIXED work -- 6458 charges folded into hcore once -- so the absolute overhead grows slowly (1.3 -> 6.0 s) while the QM SCF grows fast (0.28 -> 43 s), and the RATIO therefore FALLS with QM size. Applying benzene's 5.7x at 71 atoms overstates the embedding cost by ~5x. Use the absolute overhead, not the ratio, and never carry a ratio measured at one QM size to another.

  A worked consequence: an RHF/STO-3G QM/MM optimization of this 71-atom ligand in the full pocket MEASURED 77.8 s/step (6 steps, 466.9 s), so the 34-step floor is ~0.73 h, not 5.8. The 5.8 h figure is a KS-DFT number -- tier 4's 612 s multiplicand is tier4_dft, i.e. DFT with a grid, not the 49.1 s RHF single point measured here. Quote the one that matches the method you are actually running.

Every tiers.py:NN citation in this table pointed at a DOC COMMENT

Three of the cost rows cited tiers.py:11/13/14 or vina_dock.py:7. Every one of those lines is the module header's own cost table -- the citations pointed at prose, not at a measurement, and the golden path and the docstring were quoting each other.

A file.py:NN reference survives the check "is this sourced?", which is what let three unmeasured figures sit in a MEASURED table. They also rot: correcting the tier-2 line shifted the numbering, so tiers.py:13 and :14 came to point at the wrong rows entirely.

All three are now measured THROUGH THE TIER FUNCTIONS (not through the underlying library, which is a different cost):

Tier 3's old ~0.5 s was the right order. Tier 2's ~1 ms was 20x low. Both were unmeasured; the difference is luck, not diligence.

When a cost in this table cites a file and line, open it. If the line is a doc comment, the number is unsourced.

Tier 2 is 20x costlier than claimed, and it changes nothing (MEASURED 2026-09-19)

The ~1 ms/pose figure cited tiers.py:12 -- which is the same doc comment. A circular citation, never a measurement, and it described a single MMFF point while tier2_forcefield does embed + optimize.

Tail exponent 1.66, projecting to ~73 ms at danuglipron's 71 atoms -- about 20x the claimed figure.

And the conclusion is unchanged. Tier 2 runs on the ~10% that survive docking, so at N=1000 it is 0.4% of the campaign against docking's 78%. Recorded because the next person to notice the discrepancy should not have to re-measure it to find out it does not matter: a wrong number that changes no decision is still worth fixing once, and worth marking as decision-neutral so it is not fixed twice.

The xtb -> DFT ratio is ~200x at drug scale, not 1000x (MEASURED 2026-09-19)

"DFT costs ~1000x" appeared twice as a STOP-HERE decision criterion. Measured on the SAME molecules (the only comparison that means anything -- the table's 0.5 s and 96-612 s rows came from different systems), GFN2-xTB relax vs PBE/STO-3G single point:

The ratio GROWS with size, as N^1.17 on the last two points, because DFT scales ~N^2.3 while xtb is much flatter. Projecting:

So the ~1000x figure is right only above ~270 atoms -- roughly 4x larger than anything this pipeline runs. At danuglipron's 71 atoms it overstates by ~5x.

The decision does not change: 200x is still a decisive reason to stop at tier 3 for a coarse sort. What changes is that the number is now measured, and anyone budgeting a campaign from it is out by 5x rather than in the right place by luck.

The CHEAP stages, measured at last (2026-09-19)

The cost table covered the quantum tiers and said nothing about the four stages that run before them -- enumeration, descriptors, embedding, toxicology -- even though those are what a substitution campaign spends its first hour on. All MEASURED on danuglipron (73 proposals from 8 substituent groups), min of 3-5 reps, single-threaded:

What this changes about where to worry. For 1000 analogues the whole cheap half is 1000 x (2.8 + 1.7 + 214 + 9.4) ms = ~3.8 minutes, of which embedding is 94%. Everything upstream of docking is free at campaign scale; the only cheap-tier stage worth optimizing is the ETKDG embed, and only if the campaign is much larger than 1000.

A docstring correction. RdkitAlertsProvider's own comment says per-molecule construction "dominates the runtime of a batch". MEASURED the ratio is 5x (47 ms build vs 9.4 ms/molecule), so constructing per molecule would cost 6x a batch, not orders of magnitude. Building it once is still right; the stated reason overstates the effect.

End to end: DOCKING dominates a substitution campaign, not DFT (2026-09-19)

Composing the measured per-stage numbers over a 10x-per-tier funnel (all -> 10% docked -> 10% xtb -> 1% DFT), single-threaded:

Shares are scale-invariant at this funnel ratio. Recomputed 2026-09-19 with M11's MEASURED 26.4 s/ligand for tier 1 rather than the ~20 s estimate used first:

The correction STRENGTHENS the conclusion rather than softening it: docking is 79% of the campaign, DFT 18%. Quote the measured row of the two above -- and read the basis caveat below before quoting either, because both rows are STO-3G and the split inverts at the default basis.

The funnel RUN end to end, at both bases (2026-09-20)

Not a model of the shares -- the actual pipeline, 10 substitution candidates of benzoic acid through dock -> FF -> xtb -> DFT against the 7LCJ pocket, keeping 6/4/2/1. Zero failures at either basis, same survivor (O=C(O)c1cccc(F)c1):

Two things this settles. At STO-3G the measured split is 73/27, which reproduces the modelled share above independently -- different system, different candidate count, same answer. And the basis caveat is REAL, not defensive: at the default the ranking inverts, DFT becomes the majority of the run, and "docking dominates, not DFT" stops being true.

So the headline holds only with the basis attached. Docking is the budget at STO-3G; DFT is the budget at def2-svp, which is what tier4_dft runs unless you say otherwise. The absolute docking cost is unchanged between the rows (32.7 s both times) -- it is DFT that moves.

That inverts the intuition this pipeline was designed around. DFT is the most expensive thing PER CALL by five orders of magnitude (6e+2 s vs 1e-5 s), and it is still only 18% of the campaign, because the funnel has already cut the population 100x by the time it runs. Docking is 79% -- it is cheap per pose and runs on EVERYTHING, 20 poses each.

AND IT INVERTS AGAIN AT THE PRODUCTION BASIS (2026-09-19). Both rows above use the 612 s DFT point, which is STO-3G. def2-svp -- the tier-4 DEFAULT -- costs ~10x that (measured under "Basis costs ~10x"), and nothing else in the table moves:

So "docking dominates, not DFT" is true of a DEMONSTRATION basis and false of the default. The M11 conclusion below stands for what it measured -- tier 1 at production exhaustiveness -- but do not carry the 79/18 split into a production budget.

This is the same conclusion M11 reached from the other direction ("the funnel spent 2.6x more than it needed to", and the fix was tier-1 effort and fan-out, not tier 4). Two independent routes to "tier 1 is the budget" is worth more than either alone.

Practical consequence. The lever is exhaustiveness and worker fan-out at tier 1, both already measured (M11: ex=4 matches ex=32's accuracy at a quarter the cost; 10 workers x cpu=1 gives 6.2x). Optimizing the DFT tier -- the instinctive target -- can win at most 18%.

ESTIMATED, with the inputs labelled: the per-stage costs are MEASURED (above, and the hierarchy table), the funnel RATIOS are a design choice, and the DFT 600 s is a 71-atom single point from a different system. Change the ratios and the shares move; the ordering is robust to anything reasonable.

Transition-state search costs 2 Hessians + n_steps (MEASURED, 2026-09-19)

New entry: until ferric_scf::saddle landed there was no saddle search to cost. crates/ferric-scf/tests/saddle_cost.rs counts the actual calls rather than timing them, because a call count is a property of the algorithm while a wall time is a property of this box.

total = n_hessian * (6N + 1)  +  (n_steps + 1)   (gradient evaluations)

CORRECTED 2026-09-19 (was 2*6N + n_steps). harmonic_frequencies takes one energy-and-gradient at the UNDISPLACED geometry before its displacement loop, and n_gradient_evaluations is zeroed AFTER that call -- so the reported count is 6N while the true cost is 6N+1. find_saddle likewise takes one gradient before its first step. A constant +3 for a default search, so every conclusion below is unchanged and the tables are each 3 low.

Worth recording WHY it survived review, because the failure is reusable: the cost test supplies its own ANALYTIC Hessian closure and therefore never calls harmonic_frequencies at all. The documented model described the finite-difference path while every assertion measured a synthetic one. The test looked like it pinned the cost model and pinned something else. Fixed by a_finite_difference_hessian_costs_6n_plus_one_gradients, which runs the real FD Hessian.

MEASURED on an analytic surface whose saddle is known in closed form (hessian_recalc_every = 0, the default):

So a search on N = 20 atoms costs 242 + n_steps + 1 gradient evaluations, and the two Hessians dominate until n_steps exceeds ~240. That is the whole reason hessian_recalc_every defaults to 0; MEASURED, setting it to 1 doubles the Hessian work (2 -> 4 on this surface, i.e. 240 -> 480 gradient-equivalents at N=20).

The multiplicand, measured at three sizes (2026-09-19). The TS cost model 2*(6N+1) + (n_steps+1) had a measured MULTIPLIER and an unmeasured MULTIPLICAND -- one gradient, taken from a 71-atom DFT single point on a different system. Measured directly on linear alkanes, RHF/STO-3G, single-threaded:

Last-two exponent p = 2.58. THREE POINTS IS NOT A SCALING MEASUREMENT -- quote it as indicative, not as an exponent, and note the tail was fitted rather than the whole series (a global fit averages in the flat N=5->8 start).

Projecting t(N) = t(11)*(N/11)^2.58:

The step count barely matters and the BASIS dominates. Going 30 -> 100 steps at N=20 moves the total 26%; going N=20 -> 40 moves it 11x. And every row above is STO-3G, the cheapest basis there is -- a real catalyst at def2-SVP or better is orders above these. Treat the table as the N-SCALING SHAPE, not as wall times.

HOW to size the QM region

C0 says the QM region "sets the cost" and the section below says to size it first. This is HOW, and it is the first decision a catalyst user makes.

QmSelection offers three ways, and the choice matters:

What a radius actually buys, MEASURED 2026-09-20 (danuglipron in the 7LCJ pocket, seeded on one ligand atom, from Python):

That refusal is the thing to know before you sweep a radius. A pocket point charge enters the structure as symbol "X" with z = 0; it has no basis functions, so it can never be quantum, and the sphere reaches one as soon as it leaves the ligand. The error now says so and names the knob:

QM atom 156 has symbol "X", which is not an element, so it cannot be in
the QM region -- a bare-charge site (z = 0) has no basis functions. ...
If you selected by radius, REDUCE qm_radius_angstrom until the sphere
holds only real atoms, or list the QM atoms explicitly with qm_indices.

It used to stop at "which is not an element", which is true, names the index, and still leaves you guessing whether the seed or the radius was wrong.

Two different enums share the name WithinRadius, and they answer opposite questions. QmSelection::WithinRadius picks which atoms are QUANTUM; MoveMm::WithinRadius(f64) picks which MM atoms are allowed to MOVE during an optimization. Confusing them gives a QM region of the wrong size or a frozen pocket, and neither fails loudly. Note also that MoveMm's radius is measured ONCE at the starting geometry -- a set re-evaluated as atoms move would change the coordinate vector's length mid-optimization.

What a radius costs. MEASURED on a uniform shell model (3 ligand atoms plus a 2-Bohr lattice), counting only -- no SCF:

 radius(Bohr)  QM atoms   TS budget 2*(6N+1)+31 grads   rel. DFT cost N^2.3
     3.0           29                  381                     1.0x
     4.0           50                  633                     3.5x
     5.0           95                 1173                    15.3x
     6.0          145                 1773                    40.5x
     8.0          333                 4029                   274x
    10.0          551                 6645                   873x

TWO exponents compound here, which is why this is the expensive knob. Atom count grows as R^2.63 (tail-fitted over the last three points; the global fit gives 2.51 by averaging in the flat small-R start, and a uniform shell would give exactly 3 by volume). Then DFT cost grows as N^2.3 on top. So 3 -> 10 Bohr is ~870x the cost PER GRADIENT while the gradient COUNT also grows 17x.

A uniform shell is the pessimistic case and this is a COUNTING model, not a chemistry one. A real pocket is not uniform -- solvent is stripped, the protein is not a lattice, and whole-residue completion changes the boundary -- so read the EXPONENT and the shape, not the absolute atom counts. What transfers is that radius is the dominant lever and that a Bohr is not a small unit here.

That is the practical guidance the cost model was missing: size the QM region first (golden path C0 already says it "sets the cost" -- this is by how much), and do not spend effort shaving P-RFO steps.

Putting a number on a catalyst TS. Using the same 612 s DFT single point the tier-5 estimate uses, and treating a gradient as ~1 single point:

TS search, N = 20, n_steps = 30    273 gradients   ~46 h
TS search, N = 20, n_steps = 100   343 gradients   ~58 h
IRC, both branches                 142 gradients   ~24 h
-------------------------------------------------------
TS + IRC, n_steps = 30             415 gradients   ~71 h

WALL CLOCK, measured against main 2026-09-19. The budget above is in gradient evaluations, which is the right unit for projecting -- but the whole chain had never been TIMED. NH3 umbrella inversion, STO-3G, one process, OPENBLAS_NUM_THREADS=1:

C3 saddle search     1.42 s    9 steps, is_transition_state() = True
C5 IRC (both)        2.87 s    63 + 63 steps, both converged
C6 barrier                     11.141 / 11.141 kcal/mol
------------------------------------------------------------
TOTAL C3 -> C6       4.29 s

The IRC is 67% of it -- 2x the search, not the "more than half" the gradient count predicts, because its steps are plain gradients while the search pays for two finite-difference Hessians AND the steps. Both ratios say the same thing for planning: budget the pair.

This is a 4-atom molecule at the cheapest basis, so treat it as proof the chain RUNS end to end from Python, not as a catalyst estimate. The projection below is the estimate.

The IRC is not a rounding item. At 142 gradients (MEASURED: 71 per branch on NH3 inversion, two branches) it adds more than half the TS search again, and a catalyst study needs it -- without it "exactly one imaginary mode" says the geometry is A saddle, not that it is the one connecting your reactant and product. Budget the pair, not the search alone.

ESTIMATED, and the multiplicand is the load-bearing weakness -- it is a 71-atom DFT single point measured on a different system. What is MEASURED is the multiplier (2*(6N+1) + n_steps + 1 for the search, ~71 gradients per IRC branch) and the 6N+1 Hessian cost. Note the step count matters much less than it does for a minimization: going from 30 to 100 steps moves the TS total by 26%, because the fixed 242-gradient Hessian cost swamps it.

What is NOT in this budget, so it is not mistaken for a full study: the reactant and product optimizations that precede the search, a frequency run at each endpoint for ZPE, and any conformational search over the QM region. Each is its own multiple of the same single-point cost.

The floor caveat, same as the BFGS one below. The step count above comes from an analytic two-atom surface. A real catalyst TS has soft degrees of freedom it does not. Treat any step count from this test as a FLOOR.

An EMBEDDED saddle search, demonstrated end to end (2026-09-19)

examples/qmmm_saddle.rs originally showed only a REFUSAL -- H2 in an MM field has no saddle, so find_saddle declines. That is worth showing, but it does not demonstrate the workflow WORKS: code that rejected everything would print the same thing. crates/ferric-scf/tests/qmmm_saddle_converges.rs now pins the positive half.

NH3 umbrella inversion under point-charge embedding: converges in 5 steps, exactly 1 imaginary mode, is_transition_state() = true, z spread 0.0075 Bohr (i.e. planar, which is the physically right answer).

Two findings came out of building it, both of which cost real time:

1. An MM field that BREAKS THE SYMMETRY DEFINING THE SADDLE does not make the search harder -- it removes the target.

An antisymmetric pair puts a CONSTANT force along z, so planar NH3 stops being a stationary point at all. The search is right not to converge -- there is nothing there -- but it reads exactly like a solver bug. The tell that it is not step starvation: raising max_steps 60 -> 200 moved the energy by 2e-8.

Before debugging an embedded saddle search that will not converge, evaluate the gradient AT the symmetric geometry under the field. If it does not vanish, the field removed the saddle.

2. Relax every coordinate EXCEPT the one under study first. A textbook 1.01 A N-H left a bond-stretch gradient of 5.2e-3 that the convergence test (g_max 3e-4) rightly refuses; the STO-3G planar optimum is 1.006 A (g_max 6.1e-4). The search looked broken in a coordinate with nothing to do with the umbrella.

TS candidates REJECTED, recorded so they are not retried: linear H3+ and linear H2O are both SECOND-order saddles (2 imaginary, -1068 and -2328.7 cm^-1, each doubly degenerate). The bend of a linear molecule comes in a perpendicular pair, so "the linear form of a bent molecule" is almost never a transition state. NH3 is the smallest unambiguous closed-shell TS that is not already the starting geometry.

The negative half of that test is MUTATION-VERIFIED and is what makes the pair meaningful: setting external_potential: None passes the positive case -- the gas-phase saddle is right there -- and FAILS the antisymmetric one. Without it, a find_saddle that ignored the embedding entirely would look correct.

Point-charge embedding is essentially FREE (MEASURED, 2026-09-18)

Vacuum vs point-charge-embedded RHF on the SAME molecule (water), 8 MM charges of +/-0.5 e on a 6-Bohr shell, OPENBLAS_NUM_THREADS=1, min of 9 reps on an idle box (load ~1.0):

0.94-1.06x across three bases -- indistinguishable from free. The iteration count moved once (8 -> 9 at STO-3G) and not at all in the other two, so the "weak assumption" the retracted estimate worried about does not bite here.

The physics is unsurprising in hindsight: hcore_with_external folds the point-charge term into hcore ONCE before the SCF loop (see CLAUDE.md's external_potential notes), so embedding adds a fixed one-body cost and nothing per-iteration. Sub-1.0 ratios are timing noise, not speedups -- at 3 reps the STO-3G row read 0.74x, which is what prompted going to 9.

Consequence for planning: the MM environment is not what costs you. The QM region is, and specifically its GRADIENT peak -- see the next section. Do not shrink a QM region to "afford the embedding"; there is nothing to afford.

CAVEAT on scope, now RESOLVED by the sweep below: 8 charges on water measures the MECHANISM, not a protein-scale run.

...and it stays linear out to pocket scale (MEASURED, 2026-09-18)

Same QM region (water/cc-pVDZ), sweeping the MM charge count on an 8-Bohr Fibonacci shell. Min of 5 reps, OPENBLAS_NUM_THREADS=1, idle box.

Per-charge cost converges to 0.075-0.078 ms and stays there across three decades. A log-log fit over the last three points gives slope 0.993 -- linear to within 0.7%, which is what a one-body hcore term must be.

The SCF iteration count is 11 at EVERY charge count, vacuum included. The embedding field does not make the SCF harder to converge; that was the "weak assumption" behind the retracted 2x estimate, and it simply does not occur.

Extrapolating at 0.078 ms/charge (the 1000-charge row predicts 174 ms against a measured 176 ms, so the fit holds):

So embedding is free at a few hundred charges and becomes a real but modest cost at whole-protein scale -- and it is LINEAR, so it never overtakes the QM region's own N^3-N^4 scaling. Choose the QM region on its gradient peak (next section); choose the MM cutoff on physics, not on cost.

The DFT GRADIENT peak, not the SCF, is what sizes a QM/MM job (MEASURED)

New measurement, 2026-09-18, serial on an idle box (gradient_pool_peak_measurement.rs, 4 passed):

Two things follow for anyone sizing a QM/MM job.

1. Budget for the GRADIENT, not the SCF. The gradient peak is 27-92x the SCF peak that precedes it. A QM region whose SCF fits comfortably can still fail at the gradient step, which is the step every optimization and every one of the 6N+1 frequency displacements needs. Sizing a QM region from a single-point SCF is the wrong measurement.

2. alkane_16 is already at the edge. Its grid is TRUNCATED -- 405,038 of 412,500 points -- i.e. the unbatched path does not fit at 212 basis functions on this box. That is well below drug scale (danuglipron is ~73 atoms), so this is not a corner case.

Batching (PR #92) caps the gradient peak at the budget instead: benzene/cc-pVDZ 1.860 -> 0.250 GB (7.4x), benzene/6-31G 1.099 -> 0.250 GB (4.4x). The two water rows are unchanged at 1.0x by design -- batching engages only when the working set does not fit.

MEASUREMENT CAVEAT, learned the hard way: these tests carry #[ignore = "production-scale measurement; run serially on a quiet box"]. Run concurrently they FAIL, and the failure looks exactly like a code defect -- it is not. Satisfy the precondition (--test-threads=1, idle box) before believing either the numbers or a failure.

Frequencies / TS verification cost 6N+1 gradients (MEASURED)

harmonic_frequencies central-differences the ANALYTIC gradient, so a full Hessian is 6N+1 gradient evaluations (6N displaced plus one undisplaced, which n_gradient_evaluations does not count), each requiring its own converged SCF. There is no analytic Hessian to fall back on -- see section 4.

VERIFIED TWICE, and the second pass CORRECTED the first. Reading the loop gives for b in 0..n_coord over 3N coordinates with TWO energy_and_gradient calls inside (+delta, -delta), plus one at the undisplaced geometry before the loop -- from which I concluded 6N+1. That was wrong by one. The live counter FrequencyResult.n_gradient_evaluations reports exactly 6N:

The undisplaced call supplies result.energy and is not counted as a gradient evaluation; frequencies.rs:194 documents the field as 6N.

RE-CONFIRMED 2026-09-20 on four sizes, against merged main, and recorded here because the counter alone invites exactly the wrong correction:

The counter is 6N at every size, and frequencies.rs:260 makes one energy_and_gradient call before the 0..n_coord loop's two per coordinate. So BOTH numbers are right for different questions: budget 6N+1 SCFs, expect the counter to say 6N. Someone reading only the counter will "fix" the 6N+1 figures in this document and understate every Hessian budget by one gradient.

Lesson worth keeping: reading a loop is better than trusting a docstring, but an EXPOSED COUNTER beats both -- it cannot drift from what the code did. energy_and_gradient takes (mol, basis_name, op, scf_config, reference) -- no density argument, so every one of those calls runs a complete SCF from the default guess.

Scaling the one MEASURED DFT point (612.4 s at 71 atoms, STO-3G/PBE):

• a 20-atom QM region -> 120 gradient calls
• a 40-atom QM region -> 240 gradient calls

ESTIMATED, and the multiplier is the honest part: each call is a fresh SCF, so the cost is 6N x (one converged SCF + one gradient), not 6N x (one gradient).

VERIFIED: frequencies.rs contains no guess-reuse or restart machinery -- grep for guess / initial_density / restart / previous finds only unrelated test comments. Every displaced SCF therefore starts from the DEFAULT GUESS, not from the converged undisplaced density, even though a delta-Bohr displacement perturbs the density negligibly. Seeding each displacement from the undisplaced converged density is a self-contained optimization that should cut the iteration count substantially on every one of the 6N calls. Unmeasured, but the mechanism is not in doubt: it is the same reason geometry optimizers carry the density between steps.

Two practical consequences for a catalyst workflow. First, TS verification is not a cheap afterthought appended to an optimization; on a realistic QM region it can dominate the whole job. Second, this is the strongest argument for wanting libint2 deriv_order=2: an analytic Hessian replaces 6N SCFs with one CPKS solve. That is a real speedup of something that WORKS, not an unblock -- unlike the missing saddle search, which blocks the workflow outright.

One structural result that IS solid

QM/MM memory is set by the QM region alone. MM point charges carry no grid points, so embedding a ligand in a whole protein costs the same AO cache as the ligand in vacuum. This follows directly from cost.py's model (grid scales with ATOM COUNT of the QM region) and is why QM/MM is affordable where a larger QM region is not.

But read it with the gradient measurement below: the quantity that must fit is the GRADIENT peak of the QM region, which is 27-92x its SCF peak, and a 212 basis-function region already truncates the grid unbatched. "QM/MM is affordable" means the MM environment is nearly free -- it does NOT mean a generous QM region is. The two statements are often conflated, and only the first is true.

Where the budget should NOT go (MEASURED, and counter-intuitive)

tools/campaign/hierarchy.py rule 6, from RESULTS.md M11: raising Vina's exhaustiveness from 4 to 32 cost 6.8x and improved the top score by 0.005 kcal/mol -- against a scoring function whose published RMSE is ~2.5 kcal/mol, a gain ~500x smaller than its own error bar. Redock RMSD across an 8x effort range moved 0.097 A against a 0.131 A between-seed SEM, so it was not resolvable either.

The sensitive variable was the STARTING CONFORMER (0.75-1.24 A across ETKDG seeds), so three seeds at the cheap setting beat one seed at the expensive one -- cheaper AND better sampled.

Two consequences for the cost table above. First, the "~2 min/ligand" tier-1 row is exhaustiveness 32; at exhaustiveness 4 with three seeds you get a better answer for less. Second, and more general: the measured tier costs tell you what a run WILL cost, not where to spend. Find the sensitive variable by measuring, then spend there.

Why the funnel exists at all (MEASURED)

Also from hierarchy.py: Vina reproduced the crystal pose at 0.95 A in ~2 minutes; the best any tier-3 method managed in 62 minutes was 2.41 A. But r(vina_score, pose RMSD) = +0.461, and only 4 of 20 poses were under 2.0 A -- so the cheap score FINDS the right pose and cannot reliably RANK it. That asymmetry is the empirical justification for tiers 2-4 existing: if tier 1 could pick its own best pose, no rescoring would be needed.

This is also the cautionary tale for the geometry defect in section 0. The danuglipron campaign ran tier 3 alone on free-solution conformers -- no tier 1, so it never searched -- and spent four rounds of increasingly careful statistics on a metric fed geometries 2.2-4.1 A from the binding mode. The unwritten context["geometry"] reproduces exactly that failure inside a pipeline that LOOKS like it has a docking tier.

Which dispersion model, once one exists (researched 2026-09-19)

D3(BJ) is IMPLEMENTED (#99, merged), energy and analytic gradient both. The comparison behind that choice, because "add dispersion" has four plausible answers and they are not equivalent:

Accuracy is not the discriminator. Published S66 MAE (kcal/mol): BLYP-D3 0.19 vs BLYP-XDM 0.19; B3LYP-D3 0.28 vs B3LYP-XDM 0.22. The 2026 GMTKN55 paper states D3(BJ), XDM(BJ) and XDM(Z) "perform similarly". Anyone claiming a clear winner on dimer binding energies is overreading the data.

Name collision worth knowing. D3(BJ) uses Becke-Johnson DAMPING but Grimme's C6 table. XDM is Becke & Johnson's own MODEL, deriving C6/C8/C10 from the exchange-hole dipole moment. They are different things that share names.

ferric is unusually close to XDM. Each ingredient VERIFIED present:

So XDM Step 3 is already implemented and tested here. The gap is assembling grad^2-rho from Hessians that exist, plus a scalar Newton solve per grid point. The reference implementation, postg (Otero-de-la-Roza & Johnson), is GPL-3.0 and a PROGRAM not a library -- no permissively-licensed XDM library exists in any language, and it is in neither PySCF nor Psi4. That argues for a native implementation rather than against one.

A dispersion correction does NOT fix QM/MM. D3/D4/XDM/VV10 are all QM-atom-pairwise. ferric feeds MM atoms in as point charges carrying no dispersion at all, so QM-MM dispersion is a SEPARATE gap -- probably LJ terms across the boundary in ferric-mm. Easy to assume away, so stated explicitly.

The label "DFT + dispersion" -- FIXED and MERGED 2026-09-19

crates/ferric-d3 implements two-body D3(BJ) natively, Z=1..103, with ZERO new dependencies. Tier 4 now defaults to dispersion="d3bj", so the label three places carry is finally true.

INDEPENDENTLY VERIFIED by me, not taken from the implementer's report: same water geometry through ferric and through live simple-dftd3 1.6.0.

Machine precision. 21 Rust tests pass. The reference tables are GENERATED by generate_tables.py from upstream s-dftd3 Fortran rather than hand-copied, so they are reproducible and auditable; CONTRIBUTING.md gained a separate "Derived data (not linked)" section, because the crate links nothing but its tables still derive from LGPL-3.0 source -- an obligation the existing linking table could not express.

WHAT IS STILL MISSING, stated because a dispersion correction invites the assumption that everything dispersive is now handled:

• No analytic gradients. CLOSED (#99), and this line was wrong even before that -- task="optimize" was REFUSED, never silently uncorrected. Recorded because a stale "silently wrong" claim is worse than a stale "missing" one: it invites a reader to distrust results that were never produced. The gradient's load-bearing subtlety, since it is easy to reimplement wrongly: C6_AB is interpolated by both atoms' COORDINATION NUMBERS, so moving atom X changes C6 for pairs that do not contain X. Omitting that chain rule costs ~1e7x in FD error on real molecules (water 1.3e-6, CF2ONH 3.1e-4) and EXACTLY NOTHING on a dimer -- Ar2 is unchanged to 1.46e-13. A dimer-only test cannot see it.
• No ATM three-body term. Absent rather than approximated -- there is no knob that does nothing. MEASURED contribution rises with system size: water 0.0001% -> benzene 0.1003% of the two-body energy. Extrapolating that trend to drug-scale ligands is explicitly NOT done.
• QM/MM dispersion is still unfixed. D3 is QM-atom-pairwise; the QM-to-MM-point-charge gap needs LJ terms in ferric-mm.

The label "DFT + dispersion" was wrong before that

hierarchy.py:15 and vina_dock.py:8 both label tier 4 "DFT + dispersion". Grep for D3 / D3(BJ) / Grimme / dftd3 across ferric-dft, ferric-scf and ferric-python finds NO dispersion correction. VV10 exists inside specific functionals, not as an add-on for PBE. ferric_d3.py wraps the external dftd3 package but is a loose helper, not wired into the tier.

Dispersion is the DOMINANT attractive term in ligand binding, so this is not a cosmetic mislabel: tier 4 as run is missing the physics its own label claims.

3. The golden path, as an ordered list

Decision points are marked. Steps 1-6 are available today; step 7 is blocked (see section 4).

1. Ligand in. SMILES -> tools.structure.from_smiles; a file -> tools.structure.read. State charge and multiplicity explicitly.
2. Receptor in. PDB -> pdb2pqr (active_site/pdb2pqr_runner) to add hydrogens and assign MM charges. Do NOT feed a crystallographic PDB directly; tools.structure refuses one with no hydrogens for this reason.
3. Tier 1, dock. ~2 min/ligand. The pose is a HYPOTHESIS -- run redock_rmsd against a known bound pose on your target first, or you do not know the search works.
4. Harvest the pose into context["geometry"] -- see section 0. Without this, everything below scores a gas-phase conformer.
5. Tier 3, xtb. 0.05-0.152 s (MEASURED, 9-19 atoms). DECISION: stop here? If you are rank-ordering many ligands and only need a coarse sort, GFN2 is often enough. Going to DFT costs ~200x per candidate at this scale (MEASURED; ~1000x only above ~270 atoms).
6. Tier 3.5/4, DFT. DECISION: gas phase or embedded? Gas-phase DFT on a docked pose ignores the pocket electrostatics entirely. Use the pocket field (context["point_charges"], already consumed at tiers.py:196 and tiers.py:249) when the pocket is charged or polar. DECISION: do you trust the number? See the pose-noise caveat below.
7. Catalyst / barrier work. BLOCKED AVAILABLE since 2026-09-19 -- ferric.run_saddle -> run_frequencies -> run_irc, all three taking point_charges=/external_field= so the whole chain runs on ONE surface. Section 4 has the cost model and section 3b(b) the decision procedure. DECISION: is your QM region right? That is the choice that sets the cost -- see "how to size the QM region", and note the frontier trap: keeping the host MM charge across a covalent cut puts a point charge 0.443 A from the link atom and the optimization DIVERGES. Use delete-host at minimum.

The pose-noise caveat, which decides whether step 6 is worth running

funnel.py:162 keys results on iso.canonical -- ONE row per MOLECULE. But QM/MM binding energy is a property of a POSE. RESULTS.md M5/M6 MEASURED a per-pose sd of 29.07 kcal/mol and a within-molecule/between-molecule spread ratio of 228 vs 41 kcal/mol; M6 computes that resolving a 0.25 kcal/mol gap needs ~108,000 poses per candidate.

A one-pose-per-ligand QM/MM tier therefore reports a number whose noise is orders of magnitude above the distinction being asked of it. Decide the pose treatment (ensemble? Boltzmann weight? best-N?) BEFORE writing the adapter, because the data structure the funnel uses cannot currently express it.

3b. Decision procedure: what a QM/MM workflow should DO

The two workflows are NOT the same problem and must not share a recipe.

(a) Docking / binding affinity

The question is a RELATIVE energy between ligands, so error cancellation does most of the work and the QM region can be modest.

G0. POSE-QUALITY GATE -- before scoring anything.
    tools.campaign.align.pose_quality_gate(aligned, threshold=2.0 A)
    Redock a KNOWN complex; if no pose clears the bar, STOP. Nothing
    downstream can recover, and this is precisely the check whose absence
    cost the danuglipron campaign four measurement rounds.
G1. Dock (tier 1), then HARVEST the pose into context["geometry"] (section 0).
G2. Rank with GFN2-xTB in the pocket field (tier 3 + point_charges).
    DECISION: if you only need a coarse sort, STOP HERE. DFT costs ~200x at
    danuglipron scale (MEASURED, see below).
G3. QM region = the ligand. Pocket = MM point charges. No link atoms needed
    when the cut does not cross a covalent bond -- which for a non-covalent
    ligand it does not. This is the case ferric handles cleanly today.
    HOW MUCH does the field matter? MEASURED, water/STO-3G vs vacuum:

        one -0.5 charge at 3.2 A      -5.97 kcal/mol
        one -0.5 charge at 2.1 A     -12.89
        one -1.0 (Asp-like) at 2.6 A -17.39

    TENS of kcal/mol for a charged residue in contact range -- far larger
    than the substituent effects a campaign tries to resolve. Embedding is
    not a refinement here; omitting it changes the answer.

    BUT CHECK YOUR POCKET MODEL IS EXERTING A FIELD AT ALL. A symmetric or
    antisymmetric charge arrangement can cancel almost exactly at the
    ligand: MEASURED -0.002 kcal/mol for a +-0.4 pair at +-4.2 A, against
    -5.97 for a single -0.5 at 3.2 A. A near-zero embedding shift usually
    indicts the MODEL rather than showing the pocket does not matter, so
    compare against a single-charge control before concluding it is free.
G4. Report a DIFFERENCE (dG_bind between ligands, or vs a reference ligand),
    never an absolute. The absolute carries the full method error; the
    difference is what error cancellation protects.

What one G0-G4 pass COSTS, per ligand. The catalyst branch has a closed form (2*(6N+1) + (n_steps+1) gradients); this branch had per-stage numbers scattered across the note and no way to add them up. MEASURED 2026-09-19, one process, single-threaded:

**Docking is the same column as the formula below, so the table and the budget have one scope. 26.4 s/ligand at exhaustiveness 4 on 12 cores (RESULTS.md M11); it is quoted per ligand rather than per atom because Vina's cost is driven by the search, not by the atom count over this range.

*The ethanol FF number is LARGER than aspirin's on a SMALLER molecule because it is the first call in the process -- the RDKit/ETKDG warm-up documented under "the first from_smiles call costs 24x". Read 12-23 ms as the steady state.

So the budget for N ligands through G0-G3 is

N * (26.4 s docking + ~0.02 s FF + ~0.05 s xtb) + N_survivors * DFT

and DFT is the only term whose exponent hurts: 2.6 -> 62.6 -> 265.3 s across 9 -> 21 -> 34 atoms. Fitted on ATOM COUNT the exponent is 3.0 (21->34), 3.75 (9->21), 3.48 globally -- steeper than the ~N^2.3 measured elsewhere in this note, and the difference is the axis, not a contradiction: that figure is in BASIS FUNCTIONS at fixed basis, and these three molecules differ in composition as well as size, so atom count is the cruder axis. Quote whichever you fit, and say which. Everything before DFT is flat by comparison.

The practical consequence depends on the BASIS, and that is easy to get backwards. At the STO-3G numbers in the table, docking 100 ligands costs 44.0 min and DFT on 10 survivors at 34 atoms costs 44.2 -- comparable, which says "choose how many reach tier 4 before optimizing what tier 4 does".

At the tier-4 DEFAULT basis that conclusion inverts. def2-svp costs ~10x STO-3G (measured above), so the same 10 survivors cost 442 min against docking's 44:

So for a production run, tier 4 IS the budget and making it cheaper is where the work is. The STO-3G reading is right only for a demonstration basis. This is the same trap as reading any STO-3G row here as a production cost.

No pose treatment works (RESULTS.md M4-M14). Ensemble, Boltzmann weight and best-N were each measured:

So funnel.py:162's one-row-per-MOLECULE keying is not the blocker it was written up as -- no pose treatment the data structure could express resolves a 1-2 kcal/mol substituent effect. G2 is the last step whose output is trustworthy. Its own "STOP HERE if you only need a coarse sort" is now the recommendation rather than an option, and G4's dG_bind difference is reportable only when the gap is large (>~5 kcal/mol, i.e. outside the measured noise), not for lead optimisation.

This does NOT weaken G0-G3: the pose is found reliably (M9, 0.95 A redock) and the coarse sort works. It bounds what G4 may claim.

G1-G4 VERIFIED END TO END on merged main (2026-09-18)

Not a plan -- run against origin/main after PR #91 landed the readers:

PDB -> Molecule : 3 atoms, 10 electrons, ['O', 'H', 'H']
vacuum RHF      : -74.72406147 Ha, converged=True, 14 it
embedded RHF    : -74.74153533 Ha, converged=True, 12 it
G4 difference   : dE = -0.017474 Ha (-10.97 kcal/mol)

Every hop the docking branch claims is real today: tools.structure.read takes a PDB (explicit charge/multiplicity, no guessing), hands a Molecule to run_rhf, the same QM region runs again with point_charges=, and the answer reported is a DIFFERENCE rather than an absolute. Both SCFs converged.

The script is /tmp scratch, not committed -- it is four calls and is reproduced above in full effect. What matters is that it was RUN, so the (a) branch below is a description of working code rather than an intention.

The (b) branch is too, as of 2026-09-20. run_saddle landed on 2026-09-19, and every C-step was then run against merged main as two smoke tests on two systems, not one continuous C1-C5 chain -- see "Every C-step RUN against merged main" under (b).

The same thing from the CLI, no Python (2026-09-19)

G3 no longer requires writing a script. A [qmmm] TOML section drives the same path -- the QM region becomes the molecule that is solved, the MM region becomes the external potential it is solved in:

[qmmm]
pqr = "testdata/molecules/water_na.pqr"
qm_indices = [0, 1, 2]            # or: qm_seeds = [0], qm_radius_angstrom = 1.5
# link_bonds = [[0, 3]]           # required when the cut crosses a covalent bond
# boundary_scheme = "delete-host" # default; also "keep", "rc", "rcd"

RUN, not described (examples/water-qmmm.toml, water + one Na+ at 4 A, STO-3G):

[ferric] QM/MM: 3 QM atoms, 1 MM charges from testdata/molecules/water_na.pqr
embedded : -74.9653197421 Ha
vacuum   : -74.9629466809 Ha
G4 dE    : -0.002373 Ha = -1.489 kcal/mol

A geometry comes from the PQR, not from [molecule] xyz. The xyz key is still accepted and still ignored when [qmmm] is present -- a PQR carries BOTH coordinates and MM charges, and an xyz carries no charges, so the PQR has to win. This bites when computing the vacuum reference for G4: deleting the [qmmm] section makes the run fall back to the xyz, and if that file holds a different geometry you get a different molecule. Here water.xyz is an optimized HF/cc-pVDZ structure and gives -74.9631468000, which is NOT the vacuum energy of the embedded geometry and would put a 0.13 kcal/mol error straight into the difference.

For a G4 difference, take the vacuum number at the SAME geometry -- write the PQR's QM atoms out as an xyz, or call run_rhf twice with and without point_charges=. charge and multiplicity under [molecule] DO still apply, to the QM region.

(b) Catalyst optimization

The question is a BARRIER, i.e. a SADDLE POINT. Error cancellation does not save you, and the QM region must contain the reacting bonds.

C0. QM region MUST contain every bond that breaks or forms, plus any residue
    donating/accepting a proton or coordinating the metal. This is bigger
    than a ligand-only region and sets the cost.
    AND: check the MM field does not break a symmetry that DEFINES your
    saddle. MEASURED -- an antisymmetric charge pair makes planar NH3
    non-stationary (max|g_z| 3.8e-3 vs ~1e-16 in gas phase), so the search
    fails for want of a target and looks like a solver bug.
C1. The cut WILL cross covalent bonds, so link atoms are mandatory:
    .with_link_atoms(bonds, DEFAULT_LINK_SCALE) and a boundary-charge scheme
    (Z1/RC/RCD). ferric has all of these, PySCF-validated.
C2. Optimize the reactant and product complexes -- ferric CAN do this
    (optimize_qmmm is a minimizer).
C3. FIND THE TRANSITION STATE. AVAILABLE AND MERGED (#106, 2026-09-19) from
    BOTH languages: ferric_scf::saddle::find_saddle (Rust) and
    ferric.run_saddle(mol, basis, ...) (Python). The Python binding
    matters because the whole tools/ pipeline is driven from Python --
    without it C3 existed in a language the pipeline does not speak.
    See section 4.
C4. Verify the TS: n_imaginary == 1, AND the imaginary mode must point along
    the reaction coordinate (one imaginary frequency is necessary, not
    sufficient -- a methyl rotor gives one too).
    WHICH OBJECT: `n_imaginary` and `is_transition_state()` are on
    **SaddleResult** (from `run_saddle`). **FrequencyResult** (from
    `run_frequencies`) has no `n_imaginary` -- it exposes `frequencies` as a
    PROPERTY, not a method, and you count the negatives yourself. Writing
    `harmonic_frequencies -> n_imaginary()` is neither object's API and raises
    AttributeError.
    COMPLETE since #97: `PyFrequencyResult.normal_modes` is a real
    #[pyo3(get)] accessor on main (VERIFIED against origin/main
    2026-09-19), so both halves are reachable from Python.
C5. Barrier = E(TS) - E(reactant), with ZPE from the same frequency run.

C0-C5 VERIFIED reachable from PYTHON, end to end (2026-09-19)

The steps landed one at a time across several PRs, and the failure mode is a procedure that READS as complete while one step lives only in Rust. That happened twice: C3 until run_saddle was bound, and C4's mode vectors until #97 -- after which this document carried a stale "MODE VECTORS are Rust-only" caveat for a day.

So it is now asserted by execution rather than by reading:

C0/C1 QM region + link atoms     ferric.QmmmSystem
C2    optimize reactant/product  ferric.run_optimize_qmmm
C3    FIND the transition state  ferric.run_saddle
C4    verify                     ferric.run_frequencies
C5    WHICH minima does it join? ferric.run_irc          <- added 2026-09-19
C6    barrier                    IrcResult.forward_barrier() / reverse_barrier()

C3 -> C5 EXECUTED from Python (2026-09-19), NH3 umbrella inversion at
STO-3G. SCOPE: this runs `run_saddle` and `run_irc` only -- it does NOT
build a QmmmSystem (C0/C1), call run_optimize_qmmm (C2), or call
run_frequencies (C4). Those have their own coverage, so this is not a
C0-C5 chain test:

  saddle   converged, n_imaginary = 1, is_transition_state() = True
  IRC      -0.4257 / +0.4257 A pyramidalisation, both branches converged
  barrier  11.142 kcal/mol, symmetric to 3 decimals

The SIGN is the load-bearing check, not the energy. NH3's two pyramidal
minima are mirror images, so a walk that went the same way twice -- the
most likely direction bug -- gives the same energy with the SAME SIGN.
Only the sign test catches it, and it is MUTATION-VERIFIED through the
Python layer.

Pinned by crates/ferric-python/tests/test_saddle.py, and mutation-tested: renaming a checked attribute fails the test, so it is not a tautology over hasattr. It asserts REACHABILITY only -- each step has its own correctness tests; what this catches is a step quietly leaving the language tools/ is written in.

...and C3 now RUNS under embedding, not just reachable (2026-09-19)

Reachability was the weaker claim, and the QM/MM example made it weaker still: it demonstrated only a REFUSAL (H2 in an MM field has no saddle, so find_saddle declines). A refusal alone does not show the procedure works -- code that rejected everything would print the same thing.

crates/ferric-scf/tests/qmmm_saddle_converges.rs closes that. NH3 umbrella inversion under point-charge embedding converges in 5 steps to exactly one imaginary mode, is_transition_state() = true, z spread 0.0075 Bohr -- i.e. planar, the physically right answer. C3 is now demonstrated on the embedded surface, which is the surface a catalyst question is actually asked on.

The trap it surfaced, which belongs in C0. An MM field that BREAKS THE SYMMETRY DEFINING THE SADDLE removes the target rather than making the search harder: with an antisymmetric charge pair, max|g_z| at the planar geometry is 3.8e-3 against ~1e-16 in gas phase, so planar NH3 is not a stationary point at all. find_saddle correctly fails and it reads like a solver bug. See the transition-state cost section for the table and the diagnostic.

Every C-step RUN against merged main (2026-09-20)

Not one end-to-end run on one system -- two smoke tests on two systems, and the distinction matters because the sizes and the surfaces differ. Labelling this "C1-C5 end to end" would claim a continuity these runs do not have.

C1/C2/C4 -- ethane, QM = one CH3, covalent cut with a link atom, STO-3G. Embedded throughout, one MM field:

C1  QM 5 atoms ['C','H','H','H','H']; 3 MM charges; min link-charge 1.304 A
C2  optimize   converged=True  steps=4   E=-39.72650708            [0.1 s]
C4  freqs@min  9 modes, n_imag=0   counter=30 gradients, 31 SCF calls  [0.8 s]
C4  normal_modes reachable from Python: True

C4's counter=30 is n_gradient_evaluations; the BUDGET is 31 SCF calls (6N+1, the extra one undisplaced and uncounted). Both numbers appear here deliberately -- see the accounting section above for why quoting only the counter understates a Hessian by one.

C3/C5 -- planar NH3 inversion, STO-3G. A separate system, because the ethane methyl has no saddle to find:

C3  vacuum saddle   converged=True   n_imag=1  is_TS=True   E=-55.43766531
C3  in the MM field converged=False  n_imag=1  is_TS=False  (no stationary point)
C5  IRC             from the VACUUM saddle geometry and its imaginary mode,
                    run with point_charges= -- returns an IrcResult

The C3 pair is the point, and it is not a solver failure. The same search converges in vacuum and does not in the field: a symmetry-breaking MM field makes planar NH3 non-stationary, so there is no saddle left to find. Running the vacuum case is the discriminator -- without it, converged=False reads as a broken optimizer. See C0's symmetry warning.

The field IS being applied, and the honest way to show that is at ONE geometry with two Hamiltonians, since the field search has no stationary point to quote an energy from:

(The shift is taken at ONE geometry. Subtracting the two saddle energies instead gives +0.635 kcal/mol, but that compares two DIFFERENT geometries, one of them not stationary, so it is not a field shift.)

The basis argument is a BasisSet for energies and a NAME for geometry changes. Checked across the entry points 2026-09-20:

Not arbitrary -- a call that MOVES the nuclei has to rebuild the basis at each new geometry, so it needs the name rather than a prepared set. But nothing in either signature says which it wants, and passing the wrong one is a TypeError about PyString conversion that reads like a bug in your code rather than a convention. It cost a run here.

Three more API details cost a run each here, all now fixed in the C-steps above: OptimizeResult.mol is a METHOD (o.mol()), FrequencyResult.frequencies is a PROPERTY (no parentheses), and n_imaginary is on SaddleResult, not on FrequencyResult. run_saddle is closed-shell only -- a doublet guess is refused with a clear message rather than silently solved.

These are smoke tests, not the regression net. The chain is pinned by test_the_whole_embedded_chain_runs_and_the_barrier_moves (crates/ferric-python/tests/test_saddle.py), which asserts the same split -- the vacuum converges, the field does not -- and uses irc.saddle_energy as a field-detector.

C1 and C4 VERIFIED to work (2026-09-18)

Run against the merged extension, so these are not claims:

C1 bare cut       : QM has 4 atoms            (ethane, QM = one CH3)
C1 with link atom : QM has 5 atoms  -> ADDED  symbols ['C','H','H','H','H']
C4 H2 at minimum  : 1 mode, n_imaginary = 0   frequencies [5018.8] cm-1

C1 (QmmmSystem(...).with_link_atoms([(0, 4)])) caps a cut C-C bond with an H, exactly as the catalyst path requires -- the cut across a covalent bond is the step that distinguishes a catalyst QM region from a ligand one. C4's n_imaginary returns 0 at a MINIMUM, which is the verifier's negative control: a verifier that cannot report "this is not a saddle" cannot report "this is" either.

So C0-C2 and C4-C5 are working code. C3 was the ONLY gap, established by grep (no dimer / NEB / P-RFO / eigenvector-following anywhere) rather than by this script -- a negative cannot be demonstrated by running something.

C3 CLOSED AND MERGED 2026-09-19 (ferric_scf::saddle, #106), and demonstrated under QM/MM embedding rather than only in the gas phase. The chain C0-C5 is complete. What that does and does not mean is in section 4. The IRC gap is closed (irc::follow_irc); what remains is the analytic Hessian.

API INCONSISTENCY worth knowing before writing a workflow: run_rhf takes a BasisSet OBJECT (ferric.BasisSet.bundled("sto-3g")), while run_frequencies takes a basis-name STRING. Passing the wrong one is a clean TypeError, not a silent failure, but it costs a round trip. The frequency entry point is run_frequencies, NOT harmonic_frequencies -- the latter is the Rust name and is not what pyo3 exports.

C4 now COMPLETABLE from Python (2026-09-19, PR #97)

FrequencyResult.normal_modes is bound, so a Python workflow can finally do what C4 requires -- count imaginary modes AND inspect what one displaces. Demonstrated on water/STO-3G:

3 modes, n_imaginary = 0  -> MINIMUM
softest mode (2049.4 cm-1) displacements:
    O0  0.0016
    H1  0.0125
    H2  0.0125

The bend moves both hydrogens symmetrically while the oxygen barely moves -- physically right, and the kind of check that distinguishes "one imaginary mode" from "one imaginary mode ALONG THE REACTION COORDINATE". A methyl rotor also gives exactly one imaginary frequency; only the vector tells them apart.

So the catalyst path's verification half is complete: C0-C2 and C4-C5 all work from Python. C3 -- FINDING the saddle -- remains the single blocker, and it is a missing capability (no dimer/NEB/P-RFO anywhere), not a missing binding.

API caveat, VERIFIED 2026-09-18: FrequencyResult.normal_modes is NOT exposed in the pyo3 bindings... SUPERSEDED 2026-09-19. #97 exposed it. PyFrequencyResult.normal_modes is a real #[pyo3(get)] accessor (crates/ferric-python/src/lib.rs:1768 on origin/main, re-verified 2026-09-19), returning Vec<Vec<f64>> in the documented (mode, 3N) layout. A Python workflow can now both COUNT imaginary modes and INSPECT them, so C4 is complete from Python.

If a local check disagrees, check the loaded extension before the source: the .so symlinked into .venv points at the MAIN checkout's target/release, so a worktree can be testing a stale build. That is what made this caveat look current when I first re-read it today.

C3 does not exist. SUPERSEDED 2026-09-19 by ferric_scf::saddle::find_saddle (P-RFO), and wired to QM/MM in examples/qmmm_saddle.rs. C0-C5 is complete in principle; section 4 states what that does and does not mean. A catalyst workflow no longer has to import its transition state from another code -- though doing so and using ferric to verify (C4) and compute the barrier at a better level (C5) remains a perfectly good option, and is the cheaper one when a TS is already in hand.

4. Catalyst optimization: the search gap is CLOSED, the cost one is not

2026-09-19. The blocker was C3 -- no saddle search anywhere in the tree. That is now implemented.

What landed

ferric_scf::saddle::find_saddle (MERGED, #106): partitioned rational function optimization. It partitions the Hessian eigenspace and solves a separate RFO step in each -- maximize along one followed mode, minimize in the orthogonal complement (Banerjee/Adams/Simons/ Shepard, JPC 89, 52 (1985)). Between steps the Hessian is carried by a Bofill update, chosen because it does NOT preserve positive definiteness; BFGS would drive out the negative eigenvalue the whole search depends on.

Why it could not be a flag on optimize.rs: that is a MINIMIZER, and its quasi-Newton update is kept positive definite on purpose. No step size turns a minimizer into a saddle finder.

Where the Hessian comes from

P-RFO does NOT use hessian.rs: hessian.rs::rhf_hessian is a documented stub that ALWAYS returns Err (see four paragraphs below). The working Hessian is frequencies.rs's central-differenced analytic gradient, and that is what saddle.rs calls.

The cost, which is now the binding constraint

The Hessian is 6N+1 gradient evaluations -- the counter reports 6N and omits the undisplaced call (MEASURED via the gradient counter: H2 = 12, water = 18 -- exactly 6N). For a 20-atom QM region that is 120 gradients for ONE Hessian. Rebuilding it every step is not a search, it is a Hessian benchmark, which is why hessian_recalc_every defaults to 0 (build once, then Bofill). A realistic catalyst run is therefore:

1 Hessian (6N+1 gradients) + ~20-60 P-RFO steps (1 gradient each)

so the Hessian dominates at small N and the steps dominate past roughly N = 10. That ratio, not the algorithm, is what sizes a catalyst job now.

What is still NOT available

• No analytic Hessian. libint2 deriv_order=2 is absent from this build (VERIFIED from ~/.local/include/libint2/config.h: INCLUDE_ERI 1, INCLUDE_ONEBODY 1). Raising it means re-running libint's generation stage, a once-off out-of-band build, not a cmake flag. Every Hessian here is finite difference.

• No IRC. CLOSED 2026-09-19. irc::follow_irc answers it. Note what it does NOT do: it stops on a gradient threshold and does not CONFIRM the endpoint is a minimum -- that needs a Hessian there, another 6N+1 gradients per side, and IrcBranch::converged reports which stopping condition fired so the caller can decide. P-RFO finds a first-order saddle; it does not prove which reaction it belongs to.

• No reaction-coordinate constraint / relaxed scan. MoveMm freezes whole MM atoms and QM atoms are always free (free_atom_indices, qmmm.rs:1815), so there is still no constrained-scan route to a starting guess.

• Not wired to QM/MM. DONE 2026-09-19, and it needed less than expected. crates/ferric-scf/examples/qmmm_saddle.rs runs the whole chain on an embedded system: SCF -> analytic embedded gradient -> finite-difference embedded Hessian -> projection -> saddle step.

  The part I expected to block it did not: a QM/MM Hessian needs no new machinery. frequencies::harmonic_frequencies already threads config.external_potential into the same rhf_gradient(.., ext) / ks_gradient_closed(.., ext) calls (frequencies.rs:649), so an embedded Hessian is one ordinary call.

  Three limitations remain, and they are why this ships as an EXAMPLE rather than a library entry point:

  • MM charges are FIXED (to_external_potential() evaluated once). A barrier computed this way omits MM relaxation along the reaction coordinate. optimize_qmmm rebuilds the field per step (qmmm.rs:2066) precisely because that matters when MM atoms are free.
  • Link atoms do not track the frontier as the QM region distorts. optimize_qmmm shares this.
  • Making it a library function means extracting optimize_qmmm's 172-line inline evaluator closure -- a refactor with its own risk, and its own PR.

  Reachable from Python since 2026-09-19: ferric.run_saddle(mol, basis, xc=, multiplicity=, max_steps=, trust_radius=, follow_mode=, delta=) -> SaddleResult, with is_transition_state() as a method so converged alone cannot be read as a TS. The refusal crosses the FFI boundary with its reason intact -- VERIFIED on H2 at 0.74 A, which returns "the projected Hessian at the starting geometry has NO negative eigenvalue (lowest = 9.612869e-1)" rather than an opaque failure.

• Cartesian only. Internal-coordinate P-RFO converges in fewer steps on floppy systems.

• The mode is not checked for being the RIGHT one. Exactly one imaginary frequency means first-order saddle, not "saddle for the reaction you meant" -- a methyl rotor gives one too. SaddleResult::imaginary_mode returns the vector so a caller can check; the module does not pretend to.

What it refuses to do, on purpose

• Starting with no negative projected eigenvalue is a HARD ERROR naming the lowest eigenvalue. P-RFO from a minimum's basin has nothing to climb and would otherwise return a minimum labelled as a transition state.
• is_transition_state() requires gradient convergence AND exactly one imaginary mode. Convergence alone is satisfied by every stationary point.

Honest status line

ferric can SEARCH for a transition state, CONFIRM one, and now FOLLOW the reaction path off it in both directions. The QM/MM wiring for the search exists and is demonstrated (NH3 inversion under point-charge embedding converges in 5 steps to exactly one imaginary mode). What remains missing is the ANALYTIC Hessian -- every Hessian here is finite-differenced from analytic gradients at 6N+1 evaluations. So a catalyst workflow knows what to do, and the remaining work is cost, not capability.

5. Ordered next actions

0. Wire saddle::find_saddle to the QM/MM evaluator. DONE 2026-09-19 -- examples/qmmm_saddle.rs, verified running end to end on an embedded system. The catalyst workflow can now run. What remains is promoting it from an example to a library entry point (needs optimize_qmmm's evaluator extracted) and lifting the fixed-MM-field approximation; see section 4.

1. Harvest the docked pose into context["geometry"] in run_funnel's stage loop (section 0). DONE (#93) -- funnel._harvest_geometry writes the key.

2. Fix the "DFT + dispersion" label. DONE (#99) -- tier4_dft passes dispersion="d3bj" (native ferric-d3) by default, so the label matches what the tier computes. QM/MM dispersion is still absent (section on dispersion above).

3. Decide the pose treatment (section 3) before any QM/MM adapter.

4. Add the tier 3.5 producer: derive_pocket_charges -> context. Both quantum tiers ALREADY consume context["point_charges"], so this is a producer, not a new capability.

5. Fix relax_pose_in_pocket's movable pocket: it builds its QmmmSystem from (q,x,y,z) with symbol "X" (pose_relaxation.py:188), but move_mm != "none" needs an MmTopology in the same atom order, and PocketCharges carries no element symbols. move_mm="within"/"all" type-checks and is unreachable in practice.

6. Fix the drifted docstring: qmmm.rs:26-34 says "no Lennard-Jones QM-MM term", but qmmm_mm_terms (qmmm.rs:1578) computes one at qmmm.rs:1674-1704. The code is right; the comment is stale.

7. Expose FrequencyResult.normal_modes to Python. [DONE, #97: FrequencyResult.normal_modes is in the bindings, which is what lets a Python caller check find_saddle's imaginary mode points along the reaction coordinate, step C4's second half.] Without it a Python workflow can count imaginary modes but not check one points along the reaction coordinate, so it cannot complete TS verification (step C4). Small, self-contained, and a prerequisite for any Python-driven catalyst work.

8. Seed each finite-difference displacement from the undisplaced converged density. VERIFIED that frequencies.rs has no restart machinery, so all 6N displaced SCFs start cold from the default guess despite being a delta-Bohr perturbation apart. Self-contained, and it pays off on every frequency run -- which is every TS verification.

9. Surface tools/ in site/src/SUMMARY.md. [DONE: the pipeline is published under "End-to-end applications", "Toxicity screening" and the "Project notebooks" section of SUMMARY.md.]