STRC Research

Recipe — Bennett Acceptance Ratio Estimator

4 min read

Recipe — Bennett Acceptance Ratio Estimator

P1 recipe — Eq. (5.50) from 2007-chipot-free-energy-calculations-book §5.7.4 (Hummer) and Eq. (2.51)–(2.52) from §2.9.1 (Chipot). Use BAR whenever forward and reverse work (or potential-energy) samples are available between two states; it is the maximum-likelihood / minimum-variance estimator and is strictly better than naive exponential averaging in essentially every case relevant to STRC work. Pairs with Recipe — FEP Point-Mutation Algorithm step (f) and with any nonequilibrium pulling protocol (Jarzynski/Crooks).

The implicit equation (Eq. 5.50, verbatim)

∑_{i=1}^{N_f}  1 / [ 1 + (N_f / N_b) · exp(β(W_i − ΔA)) ]
=
∑_{i=1}^{N_b}  1 / [ 1 + (N_b / N_f) · exp(β(W_i + ΔA)) ]

where:

  • N_f, N_b = number of forward / backward trajectories (or FEP samples).
  • W_i on the left = work along forward trajectory i (state 0 → state 1).
  • W_i on the right = work along reversed trajectory i; sign convention as in Eq. (5.33), so reversed work has the opposite sign.
  • β = 1/(k_B T).
  • ΔA = the unknown free-energy difference, solved iteratively (Newton–Raphson).

Equivalent FEP form (Eq. 2.51, §2.9.1):

exp(−β ΔA_{i, i+1})
= ⟨ {1 + exp[β(ΔU_{i, i+1} − C)]}^{−1} ⟩_i
/ ⟨ {1 + exp[−β(ΔU_{i, i+1} − C)]}^{−1} ⟩_{i+1}
· exp(−β C)

with the constant C set by:

C = ΔA_{i, i+1} + (1/β) · ln(n_i / n_{i+1})

— [Chipot & Pohorille 2007, Eq. (2.51)–(2.52) p. 65; Eq. (5.50) p. ~225]

Why use BAR over exponential averaging

Per book §5.8 (1D Brownian-dynamics calibration, Fig. 5.3):

  • Exponential average (Eq. 5.44) at N=100: bimodal distribution of estimated ΔA, strong bias, large variance. Only at N≥1,000 does the histogram center.
  • BAR (Eq. 5.50) at N=100: distribution is centered at the true ΔA (here = 0). Wider than the cumulant estimator only because the test system is highly symmetric — for asymmetric perturbations (Fig. 5.4) BAR strongly outperforms cumulant estimators which acquire systematic error.
  • BAR is maximum-likelihood-optimal under the assumption of statistically independent samples — a duplicate result was obtained by Shirts et al. via maximum likelihood (cited in book §5.7.4).

Bottom line: replace ⟨exp(−β ΔU)⟩ with BAR in every FEP/NEW analysis unless there is a specific reason not to.

Protocol

  1. Run the simulation at λ=0 and λ=1 (or for NEW: forward and reverse pulling protocols) to collect:
    • N_f forward work / energy-difference samples {W_i^f} or {ΔU_{0→1}^i}
    • N_b backward samples {W_i^b} or {ΔU_{1→0}^i}
  2. Make an initial guess ΔA_0 (e.g., from the simple average of forward and backward exponential estimates).
  3. Solve Eq. (5.50) for ΔA via Newton–Raphson. The pymbar library implements this directly (pymbar.bar); GROMACS gmx bar and NAMD’s ParseFEP also embed it.
  4. Estimate the BAR variance from the curvature of the log-likelihood at the converged ΔA — pymbar returns this as σ_BAR by default.
  5. Sanity-check using the Crooks fluctuation theorem (book §5.4): the forward and backward work histograms should cross at W = ΔA. If they don’t, sampling is incomplete.

Cautions (book §5.7.4 + §5.8)

  • Sample independence is required — if forward samples are strongly correlated and outnumber the backward, the BAR estimate becomes distorted (book gives the particle-insertion-vs-removal example). Decorrelate by subsampling at the autocorrelation time, or use MBAR for multi-state extension.
  • Particle insertion vs removal (or analogously: side-chain creation vs annihilation) is the canonical case where BAR can fail; flag in any STRC alchemical script if N_f >> N_b or vice versa.
  • For highly asymmetric perturbations BAR still works but the variance grows; staging into intermediates (see MFEP Two-Stage Strategies Table) is the better route.

How to use this in STRC

  • h01 phase5 — pharmacochaperone_phase5_mmpbsa.py: when alchemically scoring a pharmacochaperone candidate, pipe NAMD/GROMACS .fepout files through pymbar.bar rather than NAMD’s built-in ParseFEP exponential. This step alone reduces ΔΔG variance by ~2–10× (book §5.8 calibration).
  • h26 phase 1d — disulfide-cysteine designs: each Ser→Cys or Ala→Cys mutation produces forward and reverse FEP samples; BAR gives the per-design ΔΔA_dimer with proper error bars.
  • h09 hydrogel (any future SMD pulling work): when running steered MD across the assembly coordinate, collect forward and reverse pulls and analyse with BAR (the NEW form of Eq. 5.50). The Crooks crossing test is also free.

Connections

Graph View

Backlinks