bio-free-energy-calculations

What this skill does

## 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)
```

```python
# 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.

```python
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.

```python
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.

```bash
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):

```bash
# 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 fiel