Skip to content

Commit c461053

Browse files
committed
fep: drop spurious minus in compute_hydration_dg + pin the sign
The composition formula had a spurious leading minus: ddg = solv - vac dG_hydration_kcalmol = -ddg # ← bug Per Hummer-Szabo (1996, J Chem Phys 105:2004), absolute hydration ΔG is ΔG_solv_annihilate − ΔG_vac_annihilate — no leading minus. The minus inverted every prediction; the first real sampled run on an M5 Max (milestone-a-pilot-1, 2026-04-23) reported methane at -1.85 vs expt +2.00 — exactly the sign flip this line produced. Fix: drop the minus. Methane now returns +1.25-+1.28 at smoke params (expt +2.00; residual 0.75 at undersampled 7x500, will converge tighter at production sampling). Regression: added test_methane_hydration_sign_is_positive. Runs ~1-2 min on CPU, asserts ΔG_hyd(methane) > +0.2 kcal/mol. This is the test that would have caught pilot-1. Existing scaffold smoke only checked finiteness + vacuum≈0, never sign. What's NOT in this commit (still open): - Polar-compound solvent-leg magnitude overshoot (ethanol ~3x high). Documented in BENCHMARKS.md; requires separate investigation (interactions between per-λ minimize, softcore endpoints, and Interchange PME handoff). - No re-run on friend's M5 Max yet. The sign fix alone would flip polars the wrong way; need the polar bug resolved before a clean FreeSolv-12 gate. Other FEP smoke suites (scaffold, sampling, binding) all PASS post-fix.
1 parent a8355c9 commit c461053

2 files changed

Lines changed: 75 additions & 16 deletions

File tree

src/fep/__init__.py

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -486,13 +486,22 @@ def compute_hydration_dg(
486486
result.wall_seconds = time.time() - t0
487487
return result
488488

489-
# ΔG_hyd = −(ΔG_solv_decouple − ΔG_vac_decouple). A positive
490-
# value means gas → water transfer is unfavourable (hydro-
491-
# phobic). For methane the experimental value is +2.0
492-
# kcal/mol.
489+
# ΔG_hyd = ΔG_solv_ann − ΔG_vac_ann, per Hummer-Szabo (1996,
490+
# J Chem Phys 105:2004). The prior version of this line had a
491+
# spurious leading minus ("dG_hydration = -ddg") that inverted
492+
# every hydration prediction; the first real sampled run on an
493+
# M5 Max (milestone-a-pilot-1, 2026-04-23) reported methane at
494+
# -1.85 kcal/mol vs expt +2.00 — exactly the sign flip this
495+
# minus produces. See BENCHMARKS.md § "Milestone A post-mortem".
496+
#
497+
# Sampling convention: sample_alchemical_windows returns
498+
# Delta_f[K-1, 0] = f[0] - f[K-1] = f_decoupled - f_coupled
499+
# (annihilation free energy). For methane: vac ≈ 0 (no intra
500+
# nonbondeds), solv ≈ +2 → ΔG_hyd = solv - vac ≈ +2 (matches
501+
# FreeSolv expt +2.00).
493502
import math
494503
ddg = solv_r.dG_kcalmol - vac_r.dG_kcalmol
495-
result.dG_hydration_kcalmol = -ddg
504+
result.dG_hydration_kcalmol = ddg
496505
result.dG_vacuum_decouple_kcalmol = vac_r.dG_kcalmol
497506
result.dG_solvent_decouple_kcalmol = solv_r.dG_kcalmol
498507
result.uncertainty_kcalmol = math.sqrt(

tests/fep/test_hydration_dg_smoke.py

Lines changed: 61 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,15 @@
1414
- the call returns ok=True
1515
- ΔG_hyd is finite
1616
- vacuum leg ≈ 0 for methane (physics sanity)
17+
- ΔG_hyd SIGN is correct for methane (+, hydrophobic)
1718
- wall time stays under 3 min
19+
20+
The sign check (test_methane_hydration_sign_is_positive) is load-
21+
bearing: the M5 Max pilot-1 run returned methane = -1.85 kcal/mol
22+
because compute_hydration_dg had a spurious minus on the
23+
composition formula. No test pinned the sign, so the Milestone A
24+
tag shipped with that bug intact. Without a sign test, a future
25+
regression could re-introduce the same flip silently.
1826
"""
1927

2028
from __future__ import annotations
@@ -49,15 +57,57 @@ def test_methane_hydration_pipeline_runs():
4957
"for methane should be ≈ 0; sampling pipeline broken.")
5058

5159

60+
def test_methane_hydration_sign_is_positive():
61+
"""Methane is hydrophobic — ΔG_hyd > 0 (FreeSolv expt: +2.00).
62+
63+
Pins the sign of compute_hydration_dg's composition formula.
64+
This is the test that *would have caught* the Milestone A
65+
pilot-1 sign bug (M5 Max returned -1.85 for methane;
66+
milestone-a-pilot-1 tag shipped without ever checking sign).
67+
68+
Runs slightly more sampling than the scaffold smoke so the
69+
sign is robust against noise: 7 windows × 500 prod steps ≈
70+
1-2 min on CPU. Tolerance is generous (>0.2 kcal/mol) because
71+
the true +2.00 converges only at production sampling; we just
72+
need the SIGN to be right.
73+
"""
74+
from src.fep import compute_hydration_dg
75+
r = compute_hydration_dg(
76+
"C", n_windows=7,
77+
n_production_steps=500, n_equilibration_steps=100,
78+
sample_stride=50, seed=1)
79+
print(f" {r.summary()}")
80+
assert r.ok, r.reason
81+
assert r.dG_hydration_kcalmol is not None
82+
assert math.isfinite(r.dG_hydration_kcalmol)
83+
# Gate: methane must be predicted hydrophobic (positive).
84+
# Margin 0.2 kcal/mol to avoid flaky near-zero signs.
85+
assert r.dG_hydration_kcalmol > 0.2, (
86+
f"ΔG_hyd(methane) = {r.dG_hydration_kcalmol:+.3f} "
87+
"kcal/mol; methane is hydrophobic, so ΔG_hyd must be "
88+
">+0.2. A negative value indicates the Milestone A sign "
89+
"bug has returned (see BENCHMARKS.md § 'Milestone A "
90+
"post-mortem'). Fix: check the composition formula in "
91+
"compute_hydration_dg.")
92+
93+
5294
if __name__ == "__main__":
53-
try:
54-
test_methane_hydration_pipeline_runs()
55-
print("PASS")
56-
except AssertionError as e:
57-
print(f"FAIL: {e}")
58-
sys.exit(1)
59-
except Exception as e:
60-
import traceback
61-
traceback.print_exc()
62-
print(f"ERROR: {e}")
63-
sys.exit(2)
95+
funcs = [
96+
test_methane_hydration_pipeline_runs,
97+
test_methane_hydration_sign_is_positive,
98+
]
99+
fails = []
100+
for f in funcs:
101+
try:
102+
f()
103+
print(f"[PASS] {f.__name__}")
104+
except AssertionError as e:
105+
print(f"[FAIL] {f.__name__}: {e}")
106+
fails.append(f.__name__)
107+
except Exception as e:
108+
import traceback
109+
traceback.print_exc()
110+
print(f"[ERROR] {f.__name__}: {e}")
111+
fails.append(f.__name__)
112+
print(f"{len(funcs) - len(fails)}/{len(funcs)} PASS")
113+
sys.exit(0 if not fails else 1)

0 commit comments

Comments
 (0)