STRC Research

Recipe — Soft-Core Potential for Alchemical End Points

4 min read

Recipe — Soft-Core Potential for Alchemical End Points

P1 recipe — equation lifted verbatim from 2007-chipot-free-energy-calculations-book §2.8.5, p. ~59 (Chipot). Use this whenever the alchemical FEP transformation creates or annihilates atoms (which is essentially always, except for charge-only mutations). Without soft-core scaling, dual-topology runs blow up at λ → 0 or λ → 1 because residual partial charges interact with vanishing van der Waals cores at zero separation (“end-point catastrophe”). This recipe pairs with Recipe — FEP Point-Mutation Algorithm.

The soft-core vdW form (Eq. 2.49, verbatim)

For an appearing or disappearing atom i interacting with an unaltered atom j:

U_ij^softvdw(r_ij; λ) = 4 ε_ij λⁿ {
    1 / [α_vdW (1 − λ)² + (r_ij / σ_ij)⁶]²
  − 1 / [α_vdW (1 − λ)² + (r_ij / σ_ij)⁶]
}

where:

  • α_vdW is a positive constant (the soft-core parameter),
  • σ_ij and ε_ij are the standard Lennard-Jones parameters from the macromolecular force field,
  • the role of α_vdW (1 − λ)² in each denominator is to eliminate the singularity of the vdW interaction. As λ → 1 the term vanishes and the standard Lennard-Jones form is recovered; as λ → 0 the soft-core regularises the well, and the derivative dU/dλ stays bounded.
  • The exponent n is typically 1 (linear) or 2; book treats n as an exponent governing the scaling of the amplitude.

Introduction of this soft-core potential results in bounded derivatives of the potential energy function when λ tends towards 0.

— [Chipot & Pohorille 2007, p. 59, §2.8.5]

Decoupled electrostatic + vdW protocol (book §2.8.4 recommendation)

When using soft-core, run the alchemical transformation in two stages, not one:

  1. Stage E — turn off ligand/side-chain electrostatic charges (q_i: q_i^(0) → 0) at fixed full vdW. λ = 0 → 1.
  2. Stage V — turn off ligand/side-chain vdW interactions using the soft-core form above. λ = 0 → 1.

Total ΔA = ΔA_E + ΔA_V. The book justifies this:

A number of authors have opted for decoupling the mutation of the electrostatic and the van der Waals contributions, taking advantage of the fact that free energy is a state function… P(ΔU) for individual stages are narrower, and, consequently, easier to evaluate if electrostatic and van der Waals energies are modified individually rather than concomitantly.

— [Chipot & Pohorille 2007, p. 58, §2.8.4]

Default α_vdW values

Book describes α_vdW qualitatively as “a positive constant” and does not tabulate it. Values in published implementations (cite the engine, not the book):

  • NAMD / NAMD3 — α_vdW = 0.5 (linear); user-tunable via alchVdwShiftCoeff.
  • GROMACS — sc-alpha = 0.5, sc-power = 1, sc-sigma = 0.3 nm are the de facto defaults.
  • AMBER — scalpha = 0.5, scbeta = 12.0 (vdW vs Coulomb shift).

For h01 phase5, default to NAMD’s 0.5 unless evidence of remaining instability — flag this in the script as an engine choice, not derived from Chipot 2007.

Failure modes flagged in the book

  • Even with very narrow windows (δλ ≃ 10⁻⁵) the catastrophe occurs without a soft-core (§2.8.5). Narrowing windows alone is not a substitute for soft-core scaling.
  • Flexibility of the mutated group multiplies the problem if more than one conformation contributes — combine with replica-exchange in λ if seen.

How to use this in STRC

  • h01 phase5 alchemical scoring — turn this on by default whenever any pharmacochaperone candidate is annihilated/created in the binding pocket. Without it, the bound-state ΔA leg of the cycle in 2007-chipot-free-energy-calculations-book Fig. 2.8 becomes unreliable.
  • h26 disulfide cysteines — Ser → Cys gains an SH; Ala → Cys gains an SH and a methylene. Always uses soft-core in the dual-topology paradigm.
  • Documentation in scripts — when adding any FEP step, add a comment: # soft-core λ-scaling: NAMD alchVdwShiftCoeff=0.5 (engine default); rationale Chipot & Pohorille 2007 §2.8.5 eq 2.49.

Connections

Graph View

Backlinks