Fix all 13 BUG_AUDIT findings in the FEP binding pathway (Milestone B)#10
Open
Rockman6 wants to merge 10 commits into
Open
Fix all 13 BUG_AUDIT findings in the FEP binding pathway (Milestone B)#10Rockman6 wants to merge 10 commits into
Rockman6 wants to merge 10 commits into
Conversation
#3) The DDM cycle must ADD the confinement work +kT·ln(V_std/V_harm) (the cost of localising a decoupled ligand from the 1 M standard-state volume into the restraint volume) to ΔG_bind. The helper returned the *release* free energy −kT·ln(V_std/V_harm) and the call site added it, flipping the correction. At the default k this biased every absolute ΔG_bind by ≈ +2·5.27 ≈ 10.5 kcal/mol toward over-binding — enough to turn a non-binder into an apparent binder and to corrupt every absolute-Kd consumer downstream (bridge.binding_to_hill → Campaign-2 rate constants). ΔΔG within a congeneric series was insulated (the constant cancels in B−A). Fix flips the sign inside _harmonic_restraint_free_energy_kcalmol (BUG_AUDIT Option B) and leaves the call-site "+" intact. Verified against openmmtools' HarmonicRestraintForce.compute_standard_state_ correction: our additive term == −(openmmtools SSC), pinned by a new regression test. Updated the two existing closed-form tests to the corrected positive convention. Docstring now states the r0 = 0 limitation (BUG_AUDIT #9/#10) instead of the false "independent of r0" claim; r0-aware form lands with the restraint-geometry fix. The restraint-on-real-bound leg (#4) is still assumed ≈ 0 — flagged inline; absolute ΔG_bind is not yet trustworthy until that leg lands.
…AUDIT #2, #8) #2 seed plumbing: the `seed` param was passed through sample_alchemical_windows but setRandomNumberSeed() was never called on any GHMC integrator, so OpenMM used system-clock entropy and no two runs — even at seed=1 — were reproducible. Reproducibility claims in BENCHMARKS.md were false and the tightened methane gate was exposed to CI flakiness. Add _seed_ghmc_move(), which wraps GHMCMove._get_integrator to setRandomNumberSeed on the integrator it builds (equil=seed, prod=seed+1; replica-exchange move=seed). Seed 0 → nonzero so seed=0 is still deterministic. #8 n_windows guard: _split_lambda_schedule silently clamped n_windows<3 to 3 while the caller kept the original value, desyncing the u_kn matrix and crashing with an out-of-bounds write. Now raises ValueError at the schedule builder and at sample_alchemical_windows entry; the public compute_hydration_dg / compute_absolute_binding_dg return a clean non-ok envelope (no wasted OpenFF build) instead of an opaque traceback. Regression tests (no MD): seed lands on the integrator, seed 0 maps nonzero, helper is idempotent; n_windows<3 rejected at the schedule builder, the sampler, and both public entry points. Full fep suite 122 passed.
…(BUG_AUDIT #5) DDM subtracts the solvent leg from the complex leg, so the bulk-water decoupling contribution only cancels if BOTH legs decouple the ligand from an identical water model at an identical nonbonded cutoff. The amber14 complex leg used amber14/tip3pfb at a 1.0 nm PME cutoff, but the solvent leg was routed through the pure-SMIRNOFF _build_alchemical_ legs (tip3p.offxml at Sage's 0.9 nm switched cutoff) — mismatched on both water model and cutoff, leaving a systematic 1-3 kcal/mol offset on every absolute ΔG_bind. Add _build_solvent_alchemical_system_amber14 (mirrors the amber14 complex builder minus the protein) and route the solvent leg by force_field_path: amber14 → the new matched builder, smirnoff → the existing SMIRNOFF builder (which also stays the validated hydration solvent leg, untouched). The new builder floors solvent-leg padding at 1.2 nm: the 1.0 nm cutoff needs box > 2.0 nm, and a small ligand at tight padding (methane at 0.8 nm → ~1.9 nm box) otherwise trips OpenMM's "cutoff > half box" at sampler minimisation. Box size doesn't affect a neutral ligand's decoupling ΔG, so DDM cancellation is unaffected. Caught by the opt-in end-to-end sampled smoke, which now runs to completion: ΔG_bind(methane/1ubq) = +6.91 kcal/mol (ΔG_R = +5.27, sign-fix live). Tests: pin solvent-leg cutoff == complex-leg cutoff (1.0 nm) and box span > 2×cutoff. Full fep scaffold suite 123 passed.
The binding complex builders pass allow_undefined_stereo=True, but ligand_hydration_fep (__init__.py:219) and _build_alchemical_legs (:374) used the default (False). A SMILES with an ambiguous stereo centre (e.g. an exocyclic C=N, or any undefined tetrahedral centre) therefore built the complex leg fine and then raised UndefinedStereochemistryError on the solvent leg — an asymmetric mid-run failure on any such input. Set allow_undefined_stereo=True on both. RDKit assigns stereo deterministically for a fixed SMILES, so the two independent legs build the same molecule; and for a truly undefined centre the solvent leg's choice is irrelevant to ΔG (enantiomers solvate identically in achiral water). Regression: alanine (CC(N)C(=O)O) now builds through both _build_alchemical_legs and the hydration scaffold without raising.
…ion (BUG_AUDIT #7/#9/#10/#11) #7 restraint geometry: the ligand is PLACED at the binding-site centre but the CoM restraint pulled it toward the Cα-shell anchor centroid with r0=0 — a different point. The placed pose therefore started far up the harmonic wall (methane/1ubq: offset 0.42 nm ≈ 88 kcal/mol at the default k) and was yanked out of the docked pose during minimisation, so the "restrained decoupled" state was restrained at the wrong location. Fix: set r0 = ‖placement − anchor centroid‖ so the restraint minimum coincides with the pose (both amber14 and smirnoff builders). An explicit positive restraint_r0_nm still overrides. #9/#10/#11 r0-aware correction: the standard-state correction used the r0=0 isotropic-Gaussian volume V_harm = (2π·kT/k)^{3/2}, valid only at r0=0. Add _harmonic_restraint_volume_nm3 with the exact closed form for the radial harmonic ∫₀^∞ 4πr²·exp(−β·½·k·(r−r0)²)dr (via erf), verified against numerical integration to machine precision; it reduces exactly to the Gaussian at r0=0 (so the r0=0 tests are unchanged) and grows with r0. compute_absolute_binding_dg now passes the builder's actual r0 into the correction. For methane/1ubq (r0=0.42) the correction drops from +5.27 to +1.48 kcal/mol. Docstring no longer claims "independent of r0" (#11). Tests: correction is r0-aware + monotonically shrinking in r0; geometry helper returns the placement→anchor distance; the built pose sits at ~0 (not ~88) kcal/mol restraint energy. Full fep scaffold suite 128 passed; opt-in sampled smoke runs to completion (ΔG_R=+1.48 live). NOTE: the restraint-on-real-bound leg (BUG_AUDIT #4) is still assumed ≈ 0 — its magnitude is much reduced now that the restraint is aligned with the pose (#7), but a residual remains; addressed next.
…ign (BUG_AUDIT #12/#13)
#13: _sample_via_replica_exchange created a tempfile.mkdtemp scratch
dir with a .nc storage file and never removed it — one leak per call,
which fills $TMPDIR over a multi-compound bench. Wrap the sampler +
analysis in try/finally that closes both reporters and rmtree's the
scratch dir on every exit path (success, run failure, MBAR failure).
#12: estimate_free_energy's docstring said it returns ΔG_{λ=0→λ=1},
but the code returns Delta_f[K-1,0] = ΔG_{λ=1→λ=0} = +ΔG_decoupling.
Callers rely on the decoupling convention; the inverted docstring is
the same artefact that produced the reverted c461053 hydration sign
bug. Docstring now states the actual direction and warns against
flipping the sign.
Regression: the RE end-to-end smoke now also asserts no cellsim_rex_
scratch dirs leak across the call. 5/5 replica-exchange tests pass.
…cit (BUG_AUDIT #4) The DDM cycle keeps the CoM restraint on at every λ of the complex leg, so the "coupled restrained" endpoint is not the physical (unrestrained) bound state. The work between them — ΔG_restraint_on_real = G(coupled,restrained) − G(coupled,unrestrained) — was silently set to zero, so any absolute-Kd consumer (bridge.binding_to_hill → Campaign-2) trusted a value that omitted it. This commit: - Adds sample_restraint_coupling_dg(), which samples the leg by scaling lambda_restraint 1→0 at full coupling (state 0 = ON, state K-1 = OFF, so MBAR Delta_f[K-1,0] = G(on)−G(off) = ΔG_restraint_on_real). Uses a GHMCIntegrator with the seed plumbed. - Composes ΔG_bind = −(ΔG_dec_cx − ΔG_dec_sv) + ΔG_R+std − ΔG_restraint_on_real (derivation-checked sign: the restrained state is less stable, so subtracting the positive leg makes ΔG_bind correctly more favourable). - Adds dG_restraint_on_real_kcalmol + restraint_on_real_included to BindingDGResult. The leg is opt-in (include_restraint_on_real_leg / --restraint-on-real-leg on dg+bench); when off, the result FLAGS itself PROVISIONAL rather than silently assuming zero. ΔΔG cancels the term entirely. Its magnitude is much smaller now that the restraint is aligned with the pose (#7) — methane/1ubq smoke gives +0.25 kcal/mol. Validated to run end-to-end (3-leg cycle completes, field/flag populated); the PRODUCTION magnitude and the k-invariance check (BUG_AUDIT #4's verification) must be confirmed on the GPU Milestone-B run before absolute ΔG_bind is signed off. Full fep suite 137 passed.
… fixes Review of the BUG_AUDIT fix branch surfaced two real issues in the fixes themselves — fixing them rather than shipping the overclaim: 1. Seed reproducibility (#2) was overclaimed. _seed_ghmc_move seeds the integrator returned by GHMCMove._get_integrator, but apply() steps the integrator returned by the ContextCache's get_context — a separate instance. Verified empirically: compute with seed=1 twice still gives different ΔG on BOTH the Reference (deterministic) and OpenCL platforms, so it is the cache indirection, not platform threading. The seed plumbing is a correct prerequisite but NOT sufficient. Corrected the docstrings, call-site comments, and the test-module claim to say so; real reproducibility needs a direct-Context sampler refactor (the pattern sample_restraint_ coupling_dg already uses). Guidance: report FEP as mean±SD over seeds until then. (#8 n_windows guard in the same commit is solid.) 2. The restraint-on-real "PROVISIONAL" note (#4) was written to result.reason on an ok=True result. reason is the failure channel: BindingDGResult.summary() only prints it when not ok, so the caveat was invisible to anyone reading the summary, while still landing in the bench CSV reason column where it could be read as a failure. Now the provisional status is carried solely by the machine-readable restraint_on_real_included flag and surfaced in summary() (visible) — reason stays empty on success. Tests updated to assert the honest behavior (seed lands on the built integrator; provisional shows in summary() not reason). Fast + opt-in sampled smokes pass.
The Milestone-B ABSOLUTE gate (|ΔG_bind(biotin) − (−18.3)| < 3) needs the closed cycle, so the streptavidin runner now passes --restraint-on-real-leg alongside --replica-exchange. This makes the GPU run push-button for validating the audit fixes end-to-end on real data — the one step that must confirm the absolute number before Milestone-B sign-off. EGFR (ranking / ΔΔG) intentionally omits the leg since it cancels in B−A and would just cost ~50% more wall.
…AUDIT #2) Follow-up to the review finding that seeding alone does not make runs reproducible. Root-caused it: the non-determinism is the PLATFORM, not the seed. On the fast Metal/OpenCL/CUDA/CPU platforms runs scatter; on the deterministic Reference platform the openmmtools GHMC integrator is seed-INDEPENDENT (verified: seed=1 and seed=2 give bit-identical trajectories) — and the sampler's per-window minimise + cross-λ eval contexts were silently created with no platform arg, so they ran on the fast platform even when GHMC was pinned elsewhere. Fix: `sample_alchemical_windows(..., deterministic=True)` pins the Reference platform on EVERY context (initial + per-window minimise, GHMC propagation via a fresh un-shared ContextCache passed to apply(), and the eval context) so no state leaks between runs and nothing runs on a non-deterministic backend. Validated bitwise-reproducible: two runs give identical ΔG (0.29943760 == 0.29943760); the default (fast, non-deterministic) path is byte-for-byte unchanged — all 137 prior tests still pass. Threaded `deterministic` through sample_restraint_coupling_dg, compute_hydration_dg, and compute_absolute_binding_dg. Corrected the #2 docstrings/test-module notes (reproducibility is now delivered via deterministic mode, not "future work"). New witness test runs on the tiny gas-phase ethanol leg so the single-threaded Reference platform stays fast (~7 s). Use deterministic mode for CI / a reproducibility witness; production stays fast + reports mean ± SD over runs.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes every finding in
BUG_AUDIT.md(2026-05-31, 13 confirmed: 7 HIGH / 3 MEDIUM / 3 LOW). None had been fixed before this branch — the post-audit commits only added replica exchange, which is orthogonal to the sign/cycle/FF bugs.Verification
tests/fepsuite: 137 passed (baseline was 110; every fix ships a regression test).tests/bridge tests/chem tests/uq: 13 passed (no cross-module breakage).Findings → commits
db1e79aHarmonicRestraintForceSSC.ed33bd0cb521711812032allow_undefined_stereo=Trueon the solvent/hydration legs (was an asymmetric mid-run crash).1047e15737f1b7estimate_free_energy; replica-exchange temp-dir leak.eded459--restraint-on-real-leg); when off, the result is explicitly flagged provisional (no longer silent). ΔΔG cancels it.Self-review (max effort) found and fixed 2 defects in the fixes
ecbb196+ streptavidin-runner commit:_get_integratoris a prerequisite but not sufficient — openmmtools'ContextCachesteps a separately-managed integrator. Confirmed empirically:seed=1twice differs on both the deterministic Reference and OpenCL platforms. Claim, docstrings, and test corrected to be honest; true reproducibility needs a direct-Context sampler refactor (tracked follow-up). Interim guidance: report FEP as mean ± SD over seeds.reasonon anok=Trueresult, making it invisible insummary()and misleading to CSV consumers. Moved to therestraint_on_real_includedflag +summary().Status / merge intent
|ΔG_bind − (−18.3)| < 3,--replica-exchange --restraint-on-real-leg, both wired intoscripts/run_binding_streptavidin_gpu.sh). The buggy code landed near −29.The prior untracked
run/fep/streptavidin_*result predates these fixes and should be re-run on this branch before any absolute number is reported.