Recipe — FEP Point-Mutation Algorithm
P1 recipe — verbatim pseudocode from 2007-chipot-free-energy-calculations-book §2.8.6 (Chipot), p. ~60. Use this when computing the relative free-energy effect of a single residue mutation in a protein (e.g., E1659A in STRC, A1078C disulfide design) using the dual-topology paradigm. The dual-topology approach pairs with Recipe — Soft-Core Potential for Alchemical End Points to avoid end-point catastrophes.
Pseudocode (verbatim §2.8.6)
(a) Build the topologies representative of state 0 and state 1, and establish an exclusion list to prevent atoms that are not common to 0 and 1 from interacting. (b) Generate an ensemble of configurations that are representative of the reference state, λ. (c) For each configuration, evaluate the potential energy using the reference-state Hamiltonian, U(x; λ = 0). (d) Repeat the same calculations using the Hamiltonian of the target state. (e) For each configuration, evaluate the potential energy difference using (2.48). (f) Compute the ensemble average ⟨exp{−β[U(x; λ + Δλ) − U(x; λ)]}⟩_λ, from which the free energy difference ΔA = A(λ + δλ) − A(λ) can be derived. (g) Increment λ and go to stage (b).
— [Chipot & Pohorille 2007, p. 60, §2.8.6 algorithm steps a–g]
Equation (2.48), the dual-topology Hamiltonian:
U(x; λ) = λ U₁(x) + (1 − λ) U₀(x)
where U₀ contains the reference-state side-chain interactions and U₁ contains the target-state ones. Both topologies coexist on the same backbone atoms, never seeing each other through an exclusion list. (See Single-Topology vs Dual-Topology Alchemical Paradigms.)
Default parameter choices for STRC alchemical work
Lifted from book §2.8 + Fig. 2.7 caption (argon hydration calibration, p. ~46):
| Parameter | Default | Source |
|---|---|---|
| Number of λ-windows | 21 (uneven width near 0 and 1) | Fig. 2.7 caption [Chipot 2007] |
| Equilibration per window | 40 ps | Fig. 2.7 caption |
| Production per window | 400 ps | Fig. 2.7 caption |
| MD time step | 2 fs | Fig. 2.7 caption |
| Water model | TIP3P | Fig. 2.7 caption |
| Long-range electrostatics | particle-mesh Ewald (PME) | Fig. 2.7 caption |
| Van der Waals cutoff | 10 Å | Fig. 2.7 caption |
| Decouple electrostatic vs vdW transformations | yes | §2.8.4 (“strongly recommended”) |
| Window-narrowing near λ→0, 1 | yes | §2.8.5 (“end-point catastrophe” mitigation) |
The book validates this protocol on argon hydration: forward (creation) ΔA = +2.11 kcal/mol; reverse (annihilation) ΔA = −2.08 kcal/mol; experiment 2.002 kcal/mol at 298 K. Hysteresis 0.03 kcal/mol — the protocol is converged at this scale of perturbation.
Per-step efficiency upgrades (apply unless reason not to)
- Decouple electrostatics first, then vdW (§2.8.4) — annihilate point charges in window 1, then vdW in window 2. Rationale: simultaneous decoupling causes vanishing atoms to clash with nonzero residual charge. Two narrower P(ΔU) distributions converge faster than one broad joint one.
- Use BAR, not exponential averaging, in step (f) — replace ⟨exp(−βΔU)⟩ with the Bennett iterative solution. See Recipe — Bennett Acceptance Ratio Estimator. Reduces variance ~2–10×.
- Run forward + backward at every window (“simple overlap sampling,” §2.9.1) — pairs naturally with BAR; gives free-energy estimate plus immediate hysteresis check.
- Hamiltonian hopping (§2.9.2) — if windows show ergodicity issues, run all λ states in parallel and attempt swap moves with Metropolis acceptance every N_sample steps.
STRC application
- h01 phase5 —
pharmacochaperone_phase5_mmpbsa.py: this script currently runs MM-PBSA at WT and E1659A endpoints. To strengthen Phase 5 from “MM-PBSA scoring” to “rigorous alchemical scoring,” wrap step (a)–(g) with NAMD/GROMACS using the topology of the docked top-3 ligand series. Use this recipe + Recipe — Soft-Core Potential for Alchemical End Points + Recipe — Bennett Acceptance Ratio Estimator together. - h26 phase 1d — A1078C / S1080C / S1579C cysteine designs: dual-topology FEP between native and Cys side chains. Note: cysteine side chain has more atoms than serine — the dual-topology approach makes this trivial; the single-topology approach would require disappearing-atom handling.
Limitations to flag in any script docstring
- Step (b) “ensemble of configurations” assumes the reference-state MD is converged — not free; check ergodicity (STRC AF3 Static Pocket Blindness to Loop Dynamics is the failure mode for h01 if loop dynamics aren’t sampled).
- The exponential average in step (f) is biased and high-variance; always upgrade to BAR.
- Bonded-term contributions (bond lengths and angles changing during mutation) typically cancel in bound vs. free states (§2.8.4) — flag in the script if the mutation alters a sterically-constrained bond near a clash.
Connections
[source]2007-chipot-free-energy-calculations-book[part-of]free-energy-methods[see-also]Recipe — Soft-Core Potential for Alchemical End Points[see-also]Recipe — Bennett Acceptance Ratio Estimator[see-also]Single-Topology vs Dual-Topology Alchemical Paradigms[see-also]Phase Space Overlap and FEP Sampling[see-also]MFEP Two-Stage Strategies Table[applies]index[applies]index