In-Silico Pipeline (EBFC Gen 2.0 Zero-Lab Proof)

SSOT Documents — Read These First

Document	What it covers
docs/01_03_EBFC_Enzymatic_Bio_Fuel_Cell.md §3.4	Pipeline spec, TRL gate, L1-L4 definitions, artifact table
docs/protocols/ebfc/in_silico/PIPELINE_STATUS.md	Live status: running/queued/completed scripts, decision matrix
docs/protocols/ebfc/in_silico/SUMMARY.md	All L1-L4 results in one page
docs/protocols/ebfc/in_silico/L3_quantum_chemistry.md	DFT details, cascade methods comparison, ΔSCF, L3b cathode
docs/protocols/ebfc/in_silico/L1_protein_architecture.md	AlphaFold 3 results, d_FAD distance
docs/00_02_Academic_Integration_and_IP.md §1	Мінаєв (DFT) + ЧМА Бушуєва (enzymes/EIS) validation; xylem sap (bio hub) — реєстр партнерів
docs/00_02_Academic_Integration_and_IP.md §2	Publication plan — Стаття 1 (honest reframe 2026-06-05)
docs/00_07_Action_Plan_Tracker.md	HW.5.IS section — operational task status
tools/in_silico/README.md	Setup, quickstart, GPU notes, GAFF explanation

After completing ANY task in this pipeline — update ALL the SSOT docs above.

Script Dependency Graph

Parameterization (CPU, ~minutes):
  02 FAD → 03 GEN → 04 CSO → 05 CLB → 06 PPy → 07 PVI → 08 SBMA
           ↓
L2 MD (GPU):            ↓ SMILES change → rerun ALL downstream
  10 (baseline) ← 02,03
  11 (full matrix) ← 02-05    → 11* (10ns extended)
  12 (temp sweep) ← 02-05     → 4 temperatures
  13 (PSBMA diffusion) ← 08   → D_eff
  14 (xylem sap) ← 02-05      → 6 species

L3 DFT anode (CPU):
  20 (FAD) ──┐
  21f (Os B3LYP dimethyl, os_complex.json) ─┼── 22 (cascade verdict) ← rerun after 21f
  21d (Os ωB97X) ──┘
  28 (tunneling pathway) ← PDB only, no DFT deps
  29 (Nelsen λ) ← standalone (FADH₂•⁺ pathological → metal hops = ③) · 29b ← rescues 29 (FADH⁻/FADH• couple → anode inner-sphere λ_i 0.39 eV)
  21e (Os mediator Hammett series ①) ← 21b geometry
  32 (PCET E°) → 33 (PCET cascade) ← lumiflavin
  34 (micro-solvation ② cluster-continuum) ← 21b geom + hexaaqua + aqua/bis-Im speciation → 34b (ωB97X ΔSCF cross-check)

L3b DFT cathode (CPU):
  23 (ZIF clusters) → 24 (hopping t_ij, 3 pairs) → 25 (k_ET vs λ ③) ← 35 (metal λ, Nelsen 4-pt)
  24b (FO-DFT two-state coupling) = rigor upgrade of 24's crude t_ij → re-ran 25 [CHEM.14 ✅: t_ij 0.00546 eV + 0.18 eV site-gap → borderline robust to coupling, ×10⁵ excluded]

L4 Kinetics (CPU, seconds):
  30 (delta_t) → 30b (Monte Carlo) → 31 (EIS) → 40 (validation)

Bridge:
  27 (MD→DFT ensemble, FAD HOMO) ← 11 (DCD trajectory) + DFT
  28b (CHEM.16 tunneling ensemble ✅) ← 11 (DCD) + 28 (Beratan-Onuchic over frames) → β·d 2.02±0.13, thermally robust

Critical Rules

0. 🔴 conda env list in a non-interactive shell answers about ITSELF, not about the machine — ask the FILESYSTEM. The conda shell function is only defined by an interactive init, so a tool-driven Bash gets a stub that reports just base; reading that as «the env is not installed here» is a measurement substitution, and it cost a wrong verdict once (2026-09-08: silken_md was present all along). Roster: ls ~/miniforge3/envs/. Run a script without activation: ~/miniforge3/envs/silken_md/bin/python scripts/NN_*.py (same for silken_lint, silken_ml, pymol_tmp). ⚠️ Regenerating a cache is safe ONLY if you prove it: back the JSON up, re-run, then diff FIELD-BY-FIELD against the backup — a run that changes a number you did not intend is otherwise indistinguishable from one that changes only the text you edited.

1. Shared lib is SSOT — all constants, paths, banner() from lib/. Never redefine locally (19 local banners eliminated in refactoring).

2. SMILES fixes cascade — fixing a SMILES in 02-08 → rerun ALL downstream MD using that ligand. Check GAFF cache too.

3. DFT: NEVER two heavy jobs on same CPU — cache thrashing makes both 2× slower (observed: a 92-min solo job took 190+ min alongside another). Small molecules (H₂O, lumiflavin) OK to parallelize; MD-on-GPU alongside a CPU DFT job is fine (different hardware) — the hazard is two CPU-bound DFT jobs sharing cache.

4. No density_fit() for Os/Ce — auto-generated aux basis for heavy metals is 3× SLOWER than standard integrals.

5. level_shift=0.3 for open-shell transition metals — UKS on Os(III), Co-Ce, Ce without it → SCF oscillates forever (60h+ observed).

6. MD NaN at NVT ramp — 10K pre-relaxation (1000 steps at 10K) + start NVT from 50K + ramp step 10K (not 5K). Especially at high T (313K) and with many ligands. Use maxIterations=10000 for minimization. ⚠️ Scope: the production MD scripts (10/11/12/14/15), NOT the CI smoke gate — rule 8 two lines below forbids raising the cap there, and for a measured reason. The two read as contradictory only if you take this line as repo-wide.

7. MD trajectories gitignored — only JSON/PNG summaries in cache/dft/ and cache/kinetics/ committed.

8. CI smoke test — in_silico_smoke.yml: CPU-only, padding 0.5nm, 500 min iterations, 40 min timeout. Deterministic anti-NaN gate (b8b506f): seed=42 on BOTH the velocity draw and the Langevin thermostat + a 1000-step 10K pre-relax before the 298K run (the script-12-14 house pattern). 🔴 That gate is NOT sufficient, and the claim that used to stand here — «the flaky NaN was unseeded RNG … re-running the same code passes» — is REFUTED by measurement (run 34261017635, 2026-09-08): both stochastic sources were seeded and it NaN'd anyway. Seeding makes the run reproducible GIVEN the positions; it cannot fix positions that are still clashing. The real root is the truncated minimisation: the AF3 model solvates to 311 918 atoms even at padding 0.5 nm (extended tails, not a compact globule), minimizeEnergy(maxIterations=500) on that box burned 1617 s of a 40-min budget and left PE at −5.34M kJ/mol, and the 10 K pre-relax then explodes. ⛔ Do NOT "just re-run" and do NOT raise the iteration cap without measuring — at 500 iterations the job already spends 2/3 of its timeout. The axis is BOX SIZE (crop the model / smaller padding / drop disordered tails), not RNG. ✅ Fixed 2026-09-09 — 01_smoke_test_water_box.py now drops the disordered N-terminal tail (residues 1-26, AF3 pLDDT < 50) before solvation (~69K atoms vs ~310K; see the script's own "N-terminal truncation" docstring section for the measurement). in_silico_smoke.yml itself was NOT touched. Verify a touch here locally with the exact CI env (SILKEN_FORCE_PLATFORM=CPU SILKEN_WATER_PADDING=0.5 SILKEN_MIN_ITERATIONS=500 SILKEN_MD_STEPS=100).

9. 🔴 Declaring an anchor dimension: name the REFERENT, not the number — and if canon gives a RANGE, say which end you took and why. Nothing binds these scripts to canon geometry and, measured 2026-09-09 (00_07 HW.45), a value gate was refused with its price named, so this line is the only carrier. The reason a gate cannot help: both live divergences carry numbers that ARE true canon values, just of a different referent — 54's D_BUS = 1.3 is the cathode channel where 01_01 §1.4 freezes the rod at Ø1.0 (rod 1.0 + 2×0.15 liner = 1.3), and its L_ANODE_GYROID = 30 is the lower end of the canon 30–50 while the CEM mirror and sibling 58 both carry 40. A pin comparing a constant to canon is green on both. ⊕ One of the two was CLOSED 2026-09-10 and you need that before you grep: the D_BUS case is fixed — the rod diameter now has a single home (lib/constants.py D_BUS_ROD_MM) that 54/55/58 import, so the divergence is not findable in the tree any more, and the remedy was not a gate but ONE HOME. L_ANODE_GYROID = 30 vs the CEM's 40 is still live and is the instance to read. ⚠️ The class is also older than the scripts: the 1D thermal_penetration.json entered as a cache file with no script at all, its only geometry field inside the canon range, and its defect was a dimension the model did not have (no PEEK break) — so «the numbers check out» has never been the right question here. 🔴 Third live instance, and it is the one that shows the rule has a BLIND HALF (2026-09-11, 00_07 HW.33): the constant took NEITHER end of the range — it took a number from the home's PROSE. 52_z_stack_tolerance.RF_ANT_TI_CLEARANCE_MIN = 12.0 cites 02_01 §5.3; that section's normative table says ≥ 8 (10–15 desirable) with λ/40 = 8.6 as the ground, and its only 12 is the OUTCOME of a proposed two-deck board stack. So «which end did you take» has a third answer — «neither, I took a worked example» — and it is the hardest to spot, because the number IS in the cited section. ⊕ Two properties worth carrying. (a) It crossed the machine halves and doubled: tools/cad's Cem.RfClearanceMinMm = 12f carries the same mirror, and the picogk skill had already listed that field among its «silent defaults» with the verdict «all correct today» — i.e. BOTH halves checked the ADDRESS and neither read the CLAUSE. (b) git log -S prices the consensus: canon row 2026-05-16, the 12 a month later in TWO commits of one day into nine homes — one date, ONE witness (00_05 §5). 🔑 So the rule grows a step: name the referent, name which end of the range — and if the number is not an END, say what it IS (a worked example · a design point · a vendor figure) and whose. 🔴 FIFTH face, 2026-09-12, and it is a new axis rather than a new instance: the referent can be wrong by AXIS while every end is quoted correctly. 50's docstring and its printed line both named wood at 9–16 GPa as the isoelasticity target — the correct ends of the correct range of the LONGITUDINAL modulus, for a part canon says is loaded TRANSVERSELY at 0.5–1.5 GPa (01_01 §5.1). A ~10× bar, stated on every run, and it had travelled: the same figure stood in the vendor-facing RFQ acceptance row. ⚠️ Carrier gap this exposes and does not close: the fix lives in a docstring and two print() calls — 50 writes no note into its cache and test_doc_cache_sync has no row for it, so a console line is the one carrier in this skill's list with no machine reader at all. 🔴 And the MIRROR of the same rule, measured 2026-09-11 on 55_bus_mechanical: our OWN threshold quoted as if it were an external norm. INFINITE_LIFE_SF = 2.0 is declared in that script with the words «the SF line this script calls infinite life» — an honest self-declaration — and the tracker was citing it as the bar an alloy must clear, which reads as an industry standard and borrows authority the number does not have. Same axis, opposite direction: #9 catches a canon range flattened into a floor; this catches a house constant inflated into a standard. Practice: when a verdict turns on a threshold, say WHOSE it is in the same sentence — ours, the vendor’s, or the standard’s (and then name the standard). 🔴 Third face of the same coin, measured 2026-09-11 on that same script: its PROSE names a quantity it never COMPUTES, and the prose is what travels. 55_bus_mechanical carries seven wear/fretting mentions — docstring, comment, and the verdict strings it PRINTS and writes into its cache — because the 2026-09-11 ratification made wear the main axis. All seven are prose; one of them says so outright («The discriminating costs are NOT computed here»). So a reader who greps the script for the subject of the verdict finds it richly discussed and concludes the model answers it. ⛔ It did not: wear had NO model anywhere in the tree, so «rated for 20 years» was a claim with no instrument (00_07 HW.34). ✅ Since 2026-09-12 55 §wear_budget BOUNDS it — the rate the pair may have and still keep the wall — while the specific wear rate itself is still measured nowhere. 🔑 Reflex for these scripts specifically, because their verdict strings ride into the cache and from there into canon and the tracker: a verdict/caveats string may state what the script MEASURED and what it DID NOT — and when it names the second, it must name it as absent, not merely discuss it. The cheapest check before quoting any script's conclusion onward: grep it for the quantity and ask whether the hits are assignments or sentences. 🔴 SIXTH face, 2026-09-12, and it is the only one where the defect is a FLAG rather than a number: a constant that keeps calling itself an ASSUMPTION after the assumption was ratified. 55_bus_mechanical.LINER_LENGTH_ASSUMED_MM entered with a comment saying canon freezes the liner WALL and «says NOTHING about the axial extent», and wrote liner_length_is_assumed: true into the cache. The axial ⚖️ was ratified the SAME DAY (01_01 §1.4, 2026-09-12 — the tube covers the channel end to end and its lower end protrudes into the PEEK gap), so within hours the flag was declaring a FROZEN DIM soft, and the number it guarded excluded the ratified protrusion (17 instead of 18 mm, every ΔT row understated ~6 %). ⚠️ Why this one evades the whole rule above: the referent was NAMED correctly and the value WAS a true canon figure — the channel run really is 17 — so «name the referent / which end of the range» is satisfied and says nothing. What was wrong is the MODALITY: the comment asserted an open question that had closed. 🔑 Practice, and it is cheap because the flag is the tell: a constant that declares its own subject OPEN must be re-read whenever that subject is ratified — and the ratification lands in CANON, never in this file, so nothing here will red. Grep your own constants for ASSUMED/is_assumed/«open ⚖️» after any verdict in their domain; a flag that calls a spec soft is worse than no flag, because downstream it reads as «this may still move». ⊕ Same pass, same file, and worth carrying because the GREP could not find it: the protrusion fix earlier that day moved the first-contact range 4.0–6.5 → 5.03–8.45 mm, and the FMEA register kept the old pair. ⊕ That new pair was itself a ROW-MIX — the bare-column minimum (5.03) beside the bonded-column maximum (8.45); the bonded range is 5.21–8.45 (corrected 2026-09-13), so the same class struck one step later, inside its own correction — and the whole pair was retired 2026-09-14 as free-shape crossings, not stations (the «contact station» bullet under When Modifying). That fix's own tracker leg prescribed taking the perimeter by a run on the token 36, and this number shares no token with the span it derives from — a derived figure has no lexical link to its cause, so the perimeter of a span fix is «what else is derived from this span», never a grep.

Practice: every geometric constant cites either a tools/cad/cem/*.json field (canon-gated by cem_canon_sync) or a canon row; prefer importing from tools/in_silico/lib/constants.py over a local literal (51 held D_SHAFT = 11.0 and its own copy of the ISO 286 deviations in a file that already imported the same Ø11 as R_INTERFACE_M — and the copy repeated the r6-under-s6 table read, fixed 2026-09-14); and a deliberately-off value carries its caveat in every carrier — docstring · the comment over the constant · the caveats/verdict string that rides into the cache · its PIPELINE_STATUS row · its SUMMARY section. The 2026-09-08 sweep reached 55 and missed 54 because it swept by FILE; the unit is the OCCURRENCE. 🔴 Fourth live instance, 2026-09-12, and it is the rule's own practice line violated by the rule's own example file: 55_bus_mechanical.L_FREE_UNSUP = 36.0 decomposed itself in its comment as «gap 6 + cathode bore ~14 + flange/pad standoff ~16», and that third term cited NOTHING — the CEM stack gives 23 mm and the pad is the rod's own end face, so the span was over-stated 1.57× and every unsupported SF under-stated by the same factor. A constant that decomposes itself in prose reads as derived and is not. ✅ FIXED 2026-09-12 — the three spans are read from tools/cad/cem/zone2_sleeve.json + tools/cad/cem/cathode_flange.json at runtime now, so they move with the manifests; the one term with no JSON field (Zone1InsertionMm, an HW.8 placeholder living only in Cem.cs) is quoted AS a by-value crossing rather than promoted into a CEM field, because promoting it would canonise a placeholder. 🔑 What the fix is worth carrying for is not the number but WHAT MOVED WITH IT: three published verdicts reversed — every alloy now clears SF 2 bare on the welded branch (was 4 of 6), the seam's break-even k went 0.426 → 0.554 so our own AS_PRINTED_DERATE marker NO LONGER covers the joint, and the conformal films stopped bearing on the wall over the WHOLE µ sweep — they miss it at the lowest µ and bear on the rest (clearance_regime.partly_gap_limited_branches), i.e. the shorter span made the rod stiff enough to change which configuration the lining argument is about. ⚠️ For two days the cache said «never reaches the wall» for them and listed them in free_cantilever_sf_describes_these: the regime label's last branch caught «not on every µ» instead of «on no µ» (fixed 2026-09-14, pinned label-vs-contacts in test_bus_mechanical_weld_seam). A categorical label derived from a sweep must say which quantifier it uses — ALL, SOME or NONE — or its else-branch will pick one for you. (⊕ All three movements were re-derived on the contact equilibrium 2026-09-14: the seam k is gone rather than moved, and the regime keys are touches_down_on_mus / touches_down_on_every_mu now — the quantifier lesson survived the re-derivation, the numbers did not.) A span constant is not a digit: it is the configuration every verdict in the file is about. ⊕ And a FOURTH face of the same coin, cheaper to catch than the other three: a WORD with no measurer, which travels exactly like a number. «blind bore» is this script's own prose; it reached five doc homes while CathodeFlange.cs cuts the channel THROUGH, and a blind bore cannot pass a conductor, so it was never physically possible. L/D 12.6 through is a different machining class from L/D 12.6 blind, and that ratio is what carries «the vendor picks the operation». Practice: a descriptive adjective in a script comment is a claim about the part, and it needs a source exactly like a number does. 🔑 The positive worked example of face 2 now lives in the same file, so use it rather than re-deriving: WELD_KNOCKDOWN_MEASURED = None keeps the missing input VISIBLY missing, and the only comparison the block draws is against AS_PRINTED_DERATE — OUR constant, named as ours — instead of a borrowed weld figure. That is what «say whose threshold it is» looks like when applied BEFORE the fact rather than after.

10. 🔴 The L4 kinetics anchor is not what its comment says — do not cite it, and do not «fix» it one constant at a time (00_07 HW.5.IS, 2026-09-13). J_MAX_25C = 494e-6 is labelled «dgrGcGDH + Os-polymer (Zafar 2012, PMC3275720)», but that paper gives 494 ± 17 for NATIVE GcGDH at 20 mM glucose (deglycosylated: 520 ± 20; phosphate pH 7.4, graphite, flow), i.e. an OPERATING POINT — script 30 uses it as the Michaelis-Menten ASYMPTOTE, so under the source's own conditions the model returns 247. KM_GLUCOSE = 20 has no primary (read values for GcGDH 10.1–19.0 mM; ≈ 10 at pH 5.5), and test_cache_integrity pins both constants, so 40's «compare with measured Km» checks a guess. The runtime m(delta_t) thresholds are 600/7200 s placeholders, so no money rests on these numbers — the cascade is in-silico (30 → delta_t_lookup.json → 30b · 31 · 40 · 57 → SUMMARY L4 · PIPELINE_STATUS · 01_03 §3.4) and it moves as ONE re-run with a field diff of every cache.

DFT Gotchas (Hard-Won Lessons)

• PySCF no SDD — use lanl2dz for Cu/Co, stuttgart_rsc for Ce (Ce not in lanl2dz)
• wb97x-d not supported — use wb97x (range separation is the main fix, dispersion ~0.05 eV)
• ωB97X Koopmans orbital energies ≠ redox potentials — RSH gives accurate IPs but LUMO systematically too high for inter-molecular comparisons. Use ΔSCF (total energies) instead. B3LYP Koopmans works better due to error cancellation.
• Adiabatic ΔSCF — composite approach: geom opt at B3LYP/def2-SVP, SP at ωB97X/def2-TZVP. Saves orders of magnitude vs full ωB97X opt.
• Cl on flat PES — geometry optimization never converges GAU displacement criterion for Cl in Os complex. Programmatic octahedral geometry sufficient (LUMO diff < 0.002 eV after 30 cycles).
• Spin parity — odd electrons → odd spin (2S). Auto-detect: spin = mol.nelectron % 2 as fallback.
• PCET with H₃O⁺/PCM — PCM oversolvates small ions (H₃O⁺ by ~7 eV). Don't use for proton transfer corrections. Need explicit water for meaningful PCET.
• FAD in MD topology — GAFF renames FAD to "UNK", all atoms have "x" suffix. 86 atoms total, 53 heavy. Full FAD has odd electron count — set charge=1 for even.
• Os mediator speciation matters (② / script 34) — chloro vs aqua vs bis-imidazole shifts E°(Os III/II) by ~0.5 eV. On the dimethyl device mediator (+309 mV, Zafar 2012; OS-RECOMPUTE 2026-06-17) the ② frame is a chloro(+1/+2)↔{aqua,bis-Im}(+2/+3) differential-solvation bracket: chloro +0.21 (3 Cl⁻-waters, lower) ↔ bis-Im +0.55 / aqua +0.49 (upper), on the [Os(H₂O)₆] n6→n18 +0.98 eV benchmark; the +2/+3 couples carry the larger group-8 PCM bias. Decompose the cascade gap into speciation + solvation + the 4,4'-dimethyl substituent ① (+0.142 Koopmans / +0.149 adiabatic — method-dependent, NOT a stale +0.146); don't lump it. chloro↔+2/+3 bracket is functional-ROBUST; the internal aqua↔bis-Im order is functional-SENSITIVE (34b: ωB97X aqua>bis-Im, B3LYP-dimethyl bis-Im>aqua, <0.15 eV). No density_fit for Os, level_shift=0.3 for the Os(III) UKS doublet. (Pre-recompute plain-bpy values −0.91/−0.61/−0.40 + the "+200 mV / aqua>bis-Im>chloro" framing are superseded.)
• lo.PM / lo.Boys crash (PySCF lib.einsum version bug) — ValueError: not enough values to unpack (expected 4, got 3) in pipek.py. For the 2-orbital FO-DFT localisation (24b) skip PySCF lo entirely: diagonalise the metal-projected 2×2 Mulliken population matrix in the {i,j} MO basis → rotation R → H_ab = off-diagonal of Rᵀ·diag(εᵢ,εⱼ)·R (Mulliken-Hush diabatisation, pure numpy; F is diagonal = ε in the orthonormal MO basis).

MD Gotchas (Hard-Won Lessons)

• 10K pre-relaxation mandatory — 500-1000 steps at 10K before NVT ramp. Without it → NaN on ~50% of runs with multiple ligands.
• Fibonacci sphere placement — deterministic (seed=42) but can create bad contacts at specific positions. 313K (40°C) particularly vulnerable — skip if NaN persists after 3 attempts (3/4 temps sufficient).
• GAFF matching — GAFFTemplateGenerator matches by graph structure. One Molecule per unique chemical species is enough.
• L2 10ns RMSD ~4 Å is normal — AF3 structures relax 3-5 Å under AMBER ff14SB for large enzymes. Check Rg (radius of gyration) — if stable → protein folded, RMSD is just conformational relaxation. Full equilibration needs 20-50 ns.
• 25GB DCD files — use stride=10 or stride=100 when loading with mdtraj. Full load kills memory.
• PBC unwrap for ensemble graph analysis (CHEM.16 / 28b) — a PBC-wrapped protein/cofactor splits a contact graph (artificial >cutoff gaps) → Dijkstra returns NaN. make_molecules_whole() makes each molecule whole but leaves a SEPARATE non-covalent cofactor (FAD) in a different periodic image → still disconnected (verified 1/15 frames). Use traj.image_molecules(inplace=True) (default anchor = largest molecule = protein) to co-locate everything into the protein's image (15/15 frames). Apply on the FULL topology, before atom_slice.

Cascade Verdict Summary