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 free-energy differences (RBFE) or standard binding free energies (ABFE) using alchemical methods. FEP+ is a commercial workflow and OpenFE is an open-source framework. Accuracy and cost vary substantially with system, perturbation, force field, setup, sampling, and evaluation design; report the protocol and benchmark relevant to the intended decision. The Boltz-2 report includes benchmark-specific comparisons with FEP methods but does not replace prospective validation on the project chemistry.

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 congeneric lead-optimization questions, evaluate a validated RBFE protocol and perturbation network. Use endpoint methods only for decisions supported by a project-specific benchmark. Candidate counts and escalation gates should follow compute budget, uncertainty, and prospective validation rather than a universal top-N rule.

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()
)

The protocol object does not itself choose an atom mapping or define a simulation. Create or inspect a mapping (Kartograf is the OpenFE 1.7 CLI default; LOMAP is also supported), construct bound and solvent Transformation objects, create their protocol DAGs, and run them through openfe quickrun or the documented Python execution interface. Always inspect the selected mapping before running.

Lambda Window Scheduling

The 12-20 windows and 5-20 ns per-window ranges are repository starting ranges, not universal prescriptions. Select and extend them from overlap, exchange, and replicate-convergence diagnostics for the system; total cost therefore varies substantially.

Enhanced Sampling

REST2 (Replica Exchange with Solute Tempering) is one enhanced-sampling approach used in some FEP workflows. It scales selected interactions to improve barrier crossing, but suitability and implementation are engine- and protocol-specific.

In FEP+, REST2 region typically includes:

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

FEP+ can use a configured REST2 region. OpenFE's RelativeHybridTopologyProtocol uses Hamiltonian replica exchange across its lambda states by default; that is not the same as REST2, and OpenFE does not automatically apply REST2. Use only enhanced-sampling modes supported and documented by the selected protocol and version.

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 import concat
from alchemlyb.parsing import gmx
from alchemlyb.estimators import MBAR
from alchemlyb.postprocessors.units import to_kcalmol

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

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

MBAR.delta_f_ and d_delta_f_ are dimensionless (in kT) until explicitly converted. The parser shown above reads GROMACS XVG files. For other engines, use the engine-specific parser supported by the installed alchemlyb version.

Cycle Closure Analysis

For a single directed thermodynamic cycle, the signed closure residual is the sum of its edges and should be consistent with zero within uncertainty. An RMS closure statistic requires residuals from multiple cycles and a stated aggregation convention.

def cycle_closure_residual(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

Interpret each closure residual relative to propagated edge uncertainties, replicate behavior, shared-edge correlations, network topology, and the decision supported. If reporting RMS across cycles, state which cycles were included and avoid treating correlated cycles as independent observations.

Absolute Binding Free Energy (ABFE)

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

Goal: Estimate the standard binding free energy of a single ligand prospectively. Conversion to an equilibrium dissociation constant requires an explicit standard-state convention; ABFE does not generically predict an assay Ki.

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

Use OpenFE's documented AbsoluteBindingProtocol workflow: construct the ligand and complex chemical systems, select and inspect the restraint setup, create the corresponding Transformation objects, serialize them with Transformation.to_json(), and execute each transformation with openfe quickrun. Do not substitute an ad hoc absolute-free-energy CLI; OpenFE does not provide that command.

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

ABFE cost relative to RBFE depends on the protocols, number of legs/windows/repeats, and convergence requirements; estimate it from the explicit campaign plan.

MM/PBSA, MM/GBSA Endpoint Methods

Lower-cost endpoint alternatives whose usefulness must be established on a matched benchmark:

# 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: Exploratory ranking when a matched retrospective benchmark shows the endpoint method supports the intended decision. Do not transfer generic correlation ranges across targets or protocols.

Force Field Selection

For OpenFE 1.7, the versioned documentation shows OpenFF 2.1.1 for the ligand and Amber-family protein/water XMLs including ff14SB and TIP3P. Inspect and serialize the actual protocol settings because defaults change between releases.