fep/sampling: biologist-actionable MBAR-failure messages

When MBAR raises internally (singular matrix, column-sum-zero,
NaN, convergence stall) the exception text is meaningful to
stat-mech experts and opaque to wet-lab biologists. Translate
the common failure modes to physical causes ("adjacent λ-windows
don't overlap") and concrete parameter bumps ("try --n-windows
21 --n-production-steps 25000; GPU recommended"), so someone who
just burned 6 hours of GPU time knows what knob to turn.

- sampling.py: _biologist_reason_for_mbar_error() catches column-
  sum / w_nk / overlap / singular / convergence / NaN / bounds
  phrases and emits bumped-parameter hints, capped at 31 windows
  and 100 k prod steps so the advice is realistic, not scary.
- sample_alchemical_windows now wraps estimate_free_energy and
  populates result.reason via the translator (previously the raw
  pymbar exception bubbled up as "binding sampling failed: ...").
- Test: 10/10 regression over column-sum, singular, NaN, bounds,
  fallback, and the suggestion-cap boundaries.
- CI: wired into smoke.yml next to the existing sampling smoke.

Happy path: existing sampling smoke still PASSes at smoke params.

Author: Rockman6

.github/workflows/smoke.yml

@@ -160,6 +160,9 @@ jobs:
       - name: fep sampling smoke (Layer 1.3, methane vacuum Langevin+MBAR)
         run: python -u tests/fep/test_sampling_smoke.py
 
+      - name: fep MBAR-error translator (biologist-actionable failure messages)
+        run: python -u tests/fep/test_mbar_error_translation_smoke.py
+
       - name: fep end-to-end smoke (Layer 1.3, methane hydration ΔG pipeline)
         run: python -u tests/fep/test_hydration_dg_smoke.py
 

src/fep/sampling.py

@@ -285,7 +285,15 @@ def _apply_alchemical_params(ctx, comp_state):
             result.ghmc_acceptance.append(0.0)
     del _eval_ctx, _eval_int
 
-    dG_kT, ddG_kT = estimate_free_energy(u_kn, N_k)
+    try:
+        dG_kT, ddG_kT = estimate_free_energy(u_kn, N_k)
+    except Exception as e:
+        result.reason = _biologist_reason_for_mbar_error(
+            e, n_windows=n_windows,
+            n_production_steps=n_production_steps)
+        result.n_samples_per_window = n_samples_per
+        result.wall_seconds = time.time() - t0
+        return result
     kT = 0.0019872041 * temperature_K
     result.ok = True
     result.n_samples_per_window = n_samples_per
@@ -295,8 +303,65 @@ def _apply_alchemical_params(ctx, comp_state):
     return result
 
 
+def _biologist_reason_for_mbar_error(
+    exc: BaseException,
+    *,
+    n_windows: int,
+    n_production_steps: int,
+) -> str:
+    """Translate pymbar/MBAR internal exceptions into biologist-
+    actionable guidance. The raw pymbar text says things like
+    "Column sum W_nk = 0" which is not meaningful to anyone who
+    hasn't read the MBAR paper. The translation names the physical
+    cause (adjacent λ-windows don't overlap, sampling too short)
+    and suggests exact parameter bumps.
+    """
+    msg = str(exc).lower()
+    # Suggested bumps: 2× windows, 3× prod steps, capped so we
+    # don't print absurdly large numbers that scare off users on
+    # laptops. These are the Milestone-A-tier numbers that
+    # reliably converge streptavidin-class binders on GPU.
+    suggest_windows = min(max(n_windows * 2, 21), 31)
+    suggest_prod = min(max(n_production_steps * 3, 25000), 100000)
+    hint = (
+        f" → try --n-windows {suggest_windows} "
+        f"--n-production-steps {suggest_prod} "
+        "(Milestone-A-tier sampling; GPU recommended)")
+
+    # Pattern-match the distinctive pymbar failure strings.
+    if ("column sum" in msg) or ("w_nk" in msg) or (
+            "overlap" in msg):
+        return (
+            "adjacent λ-windows don't overlap — sampling was "
+            "insufficient for MBAR to reweight between states. "
+            "This is common for tight binders (|ΔG| > 10 kcal/mol) "
+            "at smoke-tier parameters." + hint)
+    if ("singular" in msg) or ("convergence" in msg):
+        return (
+            "MBAR self-consistent iteration did not converge — "
+            "typically caused by states with near-identical "
+            "reduced potentials (redundant λ schedule) or NaN/Inf "
+            "samples." + hint)
+    if ("bounds" in msg) or ("negative uncertainty" in msg):
+        return (
+            "MBAR returned non-physical uncertainty — too few "
+            "independent samples per window." + hint)
+    if ("nan" in msg) or ("inf" in msg):
+        return (
+            "reduced-potential matrix contains NaN/Inf — the "
+            "integrator blew up on at least one window. Try a "
+            "smaller timestep (--timestep-fs 0.5) or verify the "
+            "input structure has no steric clashes." + hint)
+    # Fallback: preserve the original text but frame it clearly.
+    return (
+        f"MBAR failed with: {str(exc)[:200]}." + hint)
+
+
 def estimate_free_energy(u_kn, N_k) -> tuple[float, float]:
-    """Run MBAR and return (ΔG_λ=0→λ=1, uncertainty) in kT."""
+    """Run MBAR and return (ΔG_λ=0→λ=1, uncertainty) in kT.
+
+    Raises on convergence / overlap failure. Callers should wrap
+    with _biologist_reason_for_mbar_error to translate."""
     from pymbar import MBAR
 
     mbar = MBAR(u_kn, N_k)

tests/fep/test_mbar_error_translation_smoke.py

@@ -0,0 +1,129 @@
+"""Regression tests for the biologist-readable MBAR-failure
+translator in src/fep/sampling.py.
+
+The underlying problem: pymbar's built-in error strings ("Column
+sum W_nk = 0 for state 3 and 8 other columns") are meaningful to
+statistical mechanics experts and meaningless to a wet-lab biologist
+running `cellsim fep-binding` for the first time. After a 6-hour
+GPU burn, the biologist needs to know (a) *why* it failed and (b)
+*what parameter to change* — not what a reweighting matrix is.
+
+These tests pin the translator behaviour so we never regress to
+the raw pymbar string."""
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+REPO = Path(__file__).resolve().parents[2]
+sys.path.insert(0, str(REPO))
+
+from src.fep.sampling import _biologist_reason_for_mbar_error
+
+
+def _msg(exc_text: str, **kwargs) -> str:
+    defaults = dict(n_windows=11, n_production_steps=1500)
+    defaults.update(kwargs)
+    return _biologist_reason_for_mbar_error(
+        Exception(exc_text), **defaults)
+
+
+def test_column_sum_triggers_overlap_message():
+    out = _msg(
+        "Column sum W_nk = 0 for state 0 and 4 other columns; "
+        "check that the states have sufficient overlap")
+    assert "λ-windows don't overlap" in out, out
+    assert "tight binders" in out, out
+    assert "--n-windows" in out, out
+    assert "--n-production-steps" in out, out
+
+
+def test_w_nk_phrase_triggers_overlap_message():
+    # Some pymbar versions lowercase the matrix name.
+    out = _msg("w_nk matrix has columns that do not overlap")
+    assert "λ-windows don't overlap" in out
+
+
+def test_singular_matrix_triggers_convergence_message():
+    out = _msg("Singular matrix encountered in MBAR iteration")
+    assert "did not converge" in out or "convergence" in out.lower(), out
+    assert "--n-windows" in out
+
+
+def test_nan_reduced_potentials_mentions_integrator():
+    out = _msg("Reduced potential contains NaN values")
+    assert "NaN" in out or "integrator" in out.lower(), out
+    assert "timestep" in out.lower(), out
+
+
+def test_bounds_error_mentions_samples():
+    out = _msg("BoundsError: squared uncertainty is negative")
+    assert "uncertainty" in out.lower()
+    assert "samples per window" in out.lower() or (
+        "samples" in out.lower())
+
+
+def test_fallback_preserves_original_text():
+    out = _msg("Mystery error the translator hasn't seen")
+    # Fallback must include a truncated copy of the original so
+    # debug info isn't lost, but also still suggest bumped params.
+    assert "Mystery error" in out
+    assert "--n-windows" in out
+
+
+def test_suggested_windows_bumped_at_least_21():
+    # Tight-binder-tier floor: 21 windows (Milestone-A practice).
+    out = _msg("column sum 0", n_windows=11, n_production_steps=1500)
+    assert "--n-windows 22" in out or "--n-windows 21" in out, out
+    assert "--n-production-steps" in out
+
+
+def test_suggested_prod_capped_at_100k():
+    # Don't print scare numbers — 100k is a reasonable ceiling
+    # for the suggestion text regardless of user input.
+    out = _msg("column sum 0", n_windows=11,
+               n_production_steps=50000)
+    # Bump would be 3× = 150000, but capped at 100000.
+    assert "--n-production-steps 100000" in out, out
+
+
+def test_suggested_windows_capped_at_31():
+    out = _msg("column sum 0", n_windows=21,
+               n_production_steps=1500)
+    # 2× = 42, capped at 31.
+    assert "--n-windows 31" in out, out
+
+
+def test_gpu_recommended_for_tight_binders():
+    out = _msg("Column sum W_nk = 0")
+    assert "GPU" in out
+
+
+if __name__ == "__main__":
+    funcs = [
+        test_column_sum_triggers_overlap_message,
+        test_w_nk_phrase_triggers_overlap_message,
+        test_singular_matrix_triggers_convergence_message,
+        test_nan_reduced_potentials_mentions_integrator,
+        test_bounds_error_mentions_samples,
+        test_fallback_preserves_original_text,
+        test_suggested_windows_bumped_at_least_21,
+        test_suggested_prod_capped_at_100k,
+        test_suggested_windows_capped_at_31,
+        test_gpu_recommended_for_tight_binders,
+    ]
+    fails = []
+    for f in funcs:
+        try:
+            f()
+            print(f"[PASS] {f.__name__}")
+        except AssertionError as e:
+            print(f"[FAIL] {f.__name__}: {e}")
+            fails.append(f.__name__)
+        except Exception as e:
+            import traceback
+            traceback.print_exc()
+            print(f"[ERROR] {f.__name__}: {e}")
+            fails.append(f.__name__)
+    print(f"{len(funcs) - len(fails)}/{len(funcs)} PASS")
+    sys.exit(0 if not fails else 1)