Free energy calculations

Version Compatibility

Reference examples tested with: OpenFE 1.7+, OpenMM 8.1+, GROMACS 2024+, AMBER pmemd 22+, alchemlyb 2.1+, pymbar 4.0+, RDKit 2024.09+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: openfe --version; gmx --version; pmemd.cuda --version

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Free Energy Calculations

Predict binding affinity differences (RBFE) or absolute binding affinities (ABFE) using alchemical free-energy methods. FEP+ (Schrödinger) is the commercial industry standard; OpenFE (Open Free Energy) is the open-source reference. Modern best practice achieves 1-2 kcal/mol RMSE vs experimental for well-set-up RBFE on rigid receptors. Boltz-2 affinity module (Wohlwend 2025) approaches FEP accuracy at 1000x speed on benchmarks, but FEP remains gold standard for production lead optimization.

For docking input poses, see chemoinformatics/virtual-screening. For pose validation before FEP, see chemoinformatics/pose-validation. For ML alternatives, see chemoinformatics/ml-docking-rescoring.

FEP Method Taxonomy

Method	Cost / pair	Accuracy	Use case	Fails when
FEP+ (Schrödinger)	hours-days GPU	1-2 kcal/mol RMSE	Commercial lead opt	License cost
OpenFE RBFE	hours-days GPU	comparable to FEP+	Open-source RBFE	Setup automation less mature
OpenFE ABFE	days GPU	2-3 kcal/mol RMSE	Absolute affinity	Slower; more setup care
GROMACS RBFE	hours-days GPU	1-2 kcal/mol	Power users, custom setup	Manual setup is error-prone
AMBER pmemd RBFE	hours-days GPU	1-2 kcal/mol	Tradition; force-field maturity	Manual setup
FEP-SPell-ABFE	days GPU	2-3 kcal/mol	Automated ABFE	Limited adoption
QligFEP v2.1	minutes-hours	1.5-3 kcal/mol	Q-based ligand FEP	Less standard
MM/PBSA	minutes	3-5 kcal/mol RMSE	Endpoint, fast	Limited accuracy; entropy missing
MM/GBSA	minutes	3-5 kcal/mol RMSE	Endpoint, faster than PBSA	Same caveats
Boltz-2 affinity	seconds GPU	0.66 Pearson on FEP subset	ML alternative; 1000x faster	Novel chemotypes
ALEPB / EE-AMBER	days	1-2 kcal/mol	Specialized	Limited tools

Decision: For lead-optimization SAR validation, OpenFE RBFE (open) or FEP+ (commercial) is the standard. For prospective discovery, MM/GBSA is a fast first-pass (3-5 kcal/mol RMSE); use FEP for top 10-50 candidates.

Decision Tree by Scenario

Scenario	Recommended workflow
Rank close analogs (R-group SAR)	RBFE via OpenFE (cycle: lig1↔lig2↔lig3)
Cross-scaffold ranking	ABFE per ligand; or coordinated RBFE with star network
Lead optimization 10-50 compounds	RBFE; perturbation-graph design
Single ligand affinity	ABFE (no reference needed)
Quick first-pass on top 1k	MM/GBSA after docking
Novel scaffold prospective	Boltz-2 affinity + FEP confirmation on top
Selectivity (target vs off-target)	RBFE on both proteins; report delta-delta-G
Allosteric vs orthosteric	ABFE comparable; check pose stability with MD
Ions / metal centers	Specialized force field (ZAFF, MCPB.py); not standard FEP

Relative Binding Free Energy (RBFE) Setup

Goal: Calculate delta-delta-G between two ligands (lig1 → lig2) in pocket.

Approach: Alchemical transformation lig1 → lig2 in both bound state (pocket + ligand + water) and unbound state (ligand + water alone). Thermodynamic cycle:

delta(delta-G_binding) = (delta-G_lig1->lig2 in pocket) - (delta-G_lig1->lig2 in solvent)

# OpenFE simplified setup (real usage requires complete protocol setup)
from openfe import SmallMoleculeComponent, ProteinComponent, SolventComponent
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol

protein = ProteinComponent.from_pdb_file('receptor.pdb')
ligA = SmallMoleculeComponent.from_sdf_file('ligand_A.sdf')
ligB = SmallMoleculeComponent.from_sdf_file('ligand_B.sdf')
solvent = SolventComponent()

protocol = RelativeHybridTopologyProtocol(
    RelativeHybridTopologyProtocol.default_settings()
)

OpenFE's setup automates: mapping atoms between ligands (LOMAP), building hybrid topology, generating lambda windows, equilibration, production MD.

Lambda Window Scheduling

Stage	Lambda values	Purpose
Decoupling (vdW)	0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0	Turn off ligand vdW
Charging (Coulomb)	0.0, 0.25, 0.5, 0.75, 1.0	Turn off ligand partial charges
Restraint (ABFE only)	0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0	Boresch-style restraints

Modern best practice uses 12-20 lambda windows per leg. Sampling at each window: 5-20 ns. Total simulation time per pair: 1-5 GPU-days.

REST2 Enhanced Sampling

REST2 (Replica Exchange with Solute Tempering) is the de facto standard for FEP enhanced sampling. Scales solute-solute and solute-solvent interactions; allows ligand to overcome local minima.

In FEP+, REST2 region typically includes:

• The entire ligand
• Flexible binding-site loops
• Catalytic / ionic residues with high pKa shift potential

In OpenFE, REST2 is automatically applied to the alchemical region.

MBAR/BAR Analysis

After production simulation, extract delta-G via MBAR (Multistate Bennett Acceptance Ratio) or BAR (Bennett Acceptance Ratio). MBAR uses data from all windows simultaneously; BAR uses adjacent windows.

from alchemlyb.parsing import gmx
from alchemlyb.estimators import MBAR
import pandas as pd

u_nks = []
for window in range(12):
    df = gmx.extract_u_nk(f'window_{window}.xvg', T=300)
    u_nks.append(df)

u_nk = pd.concat(u_nks)
mbar = MBAR().fit(u_nk)
print(f'delta-G: {mbar.delta_f_.iloc[0, -1]:.2f} +/- {mbar.d_delta_f_.iloc[0, -1]:.2f} kcal/mol')

alchemlyb is the standard analysis package for FEP results from GROMACS, AMBER, OpenMM.

Cycle Closure Analysis

Thermodynamic cycles must close (sum of edges = 0). Cycle closure error = root-mean-square error across closed cycles.

def cycle_closure(rbfe_results, cycle):
    # cycle is list of edges, each (lig_i, lig_j, delta_g, sd)
    total = sum(d_g for _, _, d_g, _ in cycle)
    total_var = sum(sd**2 for _, _, _, sd in cycle)
    return total, total_var ** 0.5

Acceptable cycle closure: < 0.5 kcal/mol RMS. Higher indicates insufficient sampling or force-field issues.

Absolute Binding Free Energy (ABFE)

ABFE computes delta-G of binding for a single ligand (no reference compound).

Goal: Compute Kd or Ki for a single ligand prospectively.

Approach: Decouple ligand from solvated state and from pocket-bound state separately; correction terms for analytical end states.

openfe absolute-free-energy run \
  --protein receptor.pdb \
  --ligand ligand.sdf \
  --output abfe_results/ \
  --n-lambda-charge 5 \
  --n-lambda-vdw 12 \
  --n-lambda-restraint 7

ABFE is harder than RBFE: requires Boresch-style restraints to keep ligand near pocket during decoupling. Restraint contribution must be analytically corrected.

ABFE cost: ~3x RBFE cost per ligand.

MM/PBSA, MM/GBSA Endpoint Methods

Fast (<1 hour) alternative; lower accuracy (3-5 kcal/mol RMSE):

# MM/GBSA via AMBER MMPBSA.py
MMPBSA.py -i input.in -cp complex.parm7 -rp receptor.parm7 \
          -lp ligand.parm7 -y trajectory.nc

Sample input:

&general
  startframe = 100, endframe = 1000, interval = 10
&gb
  igb = 5
&pb
  istrng = 0.150

Use case: Rank order top 100 docking poses; MM/GBSA correlates ~0.5-0.7 with experimental binding; better than docking score (0.3-0.5) but worse than FEP (0.7-0.9).

Force Field Selection

Force field	Use for	Notes
OPLS4 (Schrödinger)	FEP+ default	Commercial; well-tested
OpenFF SAGE 2.x	OpenFE default	Open-source modern
GAFF2	AMBER FEP	Use for ligand only; protein FF14SB
GAFF	Legacy	Replaced by GAFF2
CGenFF	CHARMM-style FEP	CHARMM force-field family
ANI-2x	Mixed QM/MM	Experimental for FEP
MACE-OFF	Modern ML force field	Promising for FEP, limited tooling

For OpenFE: defaults are SAGE 2.1.0 for ligand, FF14SB for protein, TIP3P for water. Override only if benchmarking.

Per-Tool Failure Modes

Insufficient sampling

Trigger: Lambda windows simulated < 5 ns each; ligand has slow rotamer change.

Mechanism: REST2 helps but isn't a panacea; some conformational changes take 100s of ns.

Symptom: Replicates disagree by > 1 kcal/mol; cycle closure > 1 kcal/mol RMS.

Fix: Increase per-window sampling to 10-20 ns; add REST2 to additional residues; check if pose is genuinely stable.

Force-field artifacts

Trigger: Charged ligand or charged pocket residue.

Mechanism: GAFF2/SAGE may misparameterize unusual functional groups (perfluoro, charged sulfonate near Asp/Glu).

Symptom: Large RBFE error (>2 kcal/mol) for a specific transformation.

Fix: Visual inspection; check ligand topology with rdkit; consider non-bonded fix or fragment-specific parameters.

Mapping ambiguity

Trigger: Two ligands differ in scaffold (not just R-groups).

Mechanism: LOMAP atom mapping may not find good correspondence; results from ambiguous mappings unreliable.

Symptom: Mapping score low; large dummy-atom count; cycle closure errors.

Fix: Manual mapping using OpenFE's editor; or use ABFE per ligand instead of RBFE.

Restraint contribution wrong (ABFE)

Trigger: Boresch restraint applied to flexible region of ligand.

Mechanism: Analytical restraint correction assumes harmonic potential at well-defined minimum.

Symptom: ABFE differs by >3 kcal/mol from experiment systematically.

Fix: Choose Boresch restraint atoms from rigid ligand core; not flexible side chains.

MM/GBSA -- bias from entropy missing

Trigger: Comparing ligands of very different size.

Mechanism: MM/GBSA misses entropy contribution; larger ligands appear more favorable.

Symptom: Larger ligands always rank higher.

Fix: Use MM/GBSA only for within-series ranking; supplement with FEP for cross-size.

Boltz-2 affinity -- chemotype OOD

Trigger: Novel chemotype outside training distribution.

Mechanism: Boltz-2 is trained on PDBbind + ChEMBL; novel scaffolds extrapolate.

Symptom: Boltz-2 affinity and FEP affinity disagree.

Fix: Use Boltz-2 as a screen; validate top candidates with FEP. Treat Boltz-2 confidence band carefully.

Reconciliation: FEP+ vs OpenFE

Aspect	FEP+	OpenFE
Force field	OPLS4 (proprietary)	SAGE 2.1.0 (open)
Workflow	Schrödinger GUI	Python CLI/API
Atom mapping	Automated LOMAP-style	LOMAP via openmm-tools
Reported accuracy	1-2 kcal/mol RMSE	Comparable; emerging benchmarks
Cost	Schrödinger license	Free + compute time
Decision	Commercial team default	Open-source / academic / cost-sensitive

For new groups, OpenFE 1.7+ is the recommended starting point; FEP+ is the gold standard for established pharma pipelines.

Common Errors

Symptom	Cause	Fix
Lambda window simulation diverges	Bad initial pose	Re-relax pose with MM minimization first
Cycle closure > 1 kcal/mol	Insufficient sampling	Increase per-window time; check replicate convergence
MBAR returns NaN	Insufficient overlap between windows	Add intermediate lambda windows
Restraint contribution wrong	Boresch atoms on flexible region	Choose 3 atoms on rigid ligand core
GROMACS REST2 setup wrong	Hot region not specified	-rest2-hot 'protein and resi 100-110' style selection
ABFE under-estimates by ~3 kcal/mol	Forgetting analytical correction	Apply delta-G_restraint_correction term
MM/GBSA rmsd doesn't match docking	Different trajectory frames	Compute MM/GBSA on MD-relaxed pose

References

• Mey et al., Living J. Comput. Mol. Sci. 2:18378 (2020) -- alchemical free energy best practices.
• Wang et al., J. Am. Chem. Soc. 137:2695 (2015) -- FEP+ method.
• Open Free Energy (OpenFE) consortium 2023+ -- OpenFE framework. Cite the current release via the OpenFE Zenodo DOI (https://github.com/OpenFreeEnergy/openfe); the earlier "Henderson 2023 Comput Phys Commun" attribution could not be verified.
• Cournia et al., J. Chem. Inf. Model. 60:4153 (2020) -- RBFE for lead optimization.
• Aldeghi M et al -- ABFE protocols (consult current literature; the earlier "Aldeghi 2018 J Cheminform 10:43" citation could not be verified — Aldeghi's 2018 ABFE work appeared as a Methods in Molecular Biology book chapter).
• Wohlwend et al. (2025) -- Boltz-2 affinity prediction.
• Shirts & Chodera, J. Chem. Phys. 129:124105 (2008) -- MBAR.
• Bennett, J. Comput. Phys. 22:245 (1976) -- BAR.

Related Skills

• chemoinformatics/virtual-screening - Source poses for FEP input
• chemoinformatics/pose-validation - PoseBusters-validate before FEP
• chemoinformatics/conformer-generation - Generate ligand 3D for FEP setup
• chemoinformatics/molecular-standardization - Standardize ligand before FEP
• chemoinformatics/ml-docking-rescoring - Boltz-2 affinity as alternative
• chemoinformatics/qsar-modeling - Surrogate models for high-throughput