Bio 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

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

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

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

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

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

Common Errors

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