free-energy-md

Free-Energy MD

Description

This skill covers quantitative free-energy calculations from molecular dynamics: alchemical methods (free energy perturbation, thermodynamic integration, BAR, MBAR), potential of mean force from umbrella sampling, absolute and relative binding/hydration/solvation free energies, thermodynamic cycles, lambda schedules and soft-core potentials, overlap diagnostics, standard-state corrections, and statistical uncertainty analysis using pymbar and alchemlyb. GROMACS provides the simulation engine with the lambda free-energy infrastructure; PLUMED provides CV-based PMF calculations; pymbar/alchemlyb provide statistically optimal analysis. Invoke this skill when a quantitative ΔG, ΔΔG, or PMF is the target observable and cannot be obtained from plain MD averages.

Domain Context

Free-energy calculations in MD operate on two distinct paradigms:

Alchemical methods compute ΔG by coupling a Hamiltonian H(λ) that interpolates between two physical end-states (λ = 0 and λ = 1) through unphysical intermediate states. The Zwanzig relation gives ΔG = −k_B T ln⟨exp(−ΔU/k_B T)⟩₀, which converges only when the configuration spaces of the two states overlap sufficiently. BAR and MBAR exploit bidirectional sampling across multiple lambda windows to give statistically optimal estimates. No physical trajectory connects the two end-states; the path is a computational construction, but the endpoint free-energy difference is exact in the limit of convergence.

PMF methods compute F(ξ) = −k_B T ln P(ξ) directly along a physical collective variable (CV) by enhanced sampling. The PMF is a projection of the true free energy landscape onto one or two coordinates; its accuracy depends on the quality of the CV as a reaction coordinate.

Several physical constraints govern applicability:

• Phase-space overlap requirement: For FEP (Zwanzig), the configurations sampled at λ = 0 must include configurations representative of λ = 1. A metric: if exp(−ΔU/k_B T) is dominated by a handful of frames (effective sample size N_eff ≪ N), the estimate is unreliable. BAR/MBAR relax this by pooling data from all intermediate windows.
• Soft-core potentials at lambda endpoints: At λ → 0 or λ → 1, a fully decoupled particle has a core that overlaps with solvent, creating numerical singularities in dU/dλ and instability in the integrator. Soft-core (Beutler) potentials replace the Lennard-Jones core with a regularized form that remains finite at λ = 0. All production alchemical FEP requires soft-core LJ.
• Standard-state correction for binding: Absolute binding free energies ΔG_bind are referenced to a standard state of 1 M (1 mol/L). Restraints that confine the ligand during alchemical annihilation impose a loss of translational/rotational entropy that must be corrected analytically (Boresch correction). Omitting this correction can introduce errors of 10–30 kJ/mol. [EXPERT REVIEW NEEDED]
• Thermodynamic cycle closure: Free energies are state functions. The sum of ΔG around a closed thermodynamic cycle must equal zero. A non-zero cycle closure error ≥ 1 kJ/mol indicates either a numerical problem (insufficient sampling in a window) or a force-field inconsistency. Use cycle closure as a convergence diagnostic.
• Force-field accuracy is the dominant error source. A perfectly converged FEP calculation with a poorly parameterized force field gives a precisely wrong answer. Experimental validation of the method on a known binding affinity series is required before applying it to novel systems.