Relative Binding Free Energy is the validation step after docking and rescoring. Where MM-GBSA Rescoring gives an end-point estimate with no sampling, this tool samples alchemical intermediate states in explicit water and reports free-energy differences with their uncertainty. It runs OpenFE on OpenMM, on GPU hardware only.
Choosing a mode
| Mode | Question it answers | Inputs |
|---|---|---|
rbfe (default) | How much tighter does ligand B bind than ligand A? | Receptor PDB, docked poses, ligand pairs |
ahfe | How favourable is moving each ligand from vacuum into water? | 3D ligands only |
rest2 | Does a docked pose stay put, and where do binding-site waters sit? | Receptor PDB, one docked pose |
How it works
rbfe. For each pair a → b you nominate, a hybrid topology morphs a into b along a series of lambda windows, once bound to the receptor (complex leg) and once free in water (solvent leg). Windows swap configurations by Hamiltonian replica exchange, and MBAR turns the samples into a free energy for each leg:
ΔΔG(a → b) = ΔG(complex) − ΔG(solvent) = ΔG_bind(b) − ΔG_bind(a)
in kcal/mol; negative means b binds more tightly. The atoms of a and b are matched with kartograf, so pairs should share a scaffold. All pairs run in one job, so closed cycles in your network give a cycle-closure residual, and a maximum-likelihood fit of the network gives a ΔG for each ligand: relative (centred on zero), or anchored to an experimental value for one reference ligand.
ahfe. Each ligand's charges are switched off and then its van der Waals contacts are removed, in water and in vacuum, giving its hydration free energy.
rest2. Replica exchange with solute scaling: the ligand and every residue within 5 Å of it are simulated as if heated, replica by replica, while the rest of the system stays at the physical temperature. Only snapshots of the unscaled replica are reported.
Ligands use the OpenFF Sage force field (openff-2.2.0), the protein ff14SB, and water TIP3P. Each ligand is charged once, before anything runs, so both ends of every pair share the same charges. To reuse charges from the Small-Molecule Force-Field Parameterizer, bind its systems to ligand_systems in a workflow; otherwise AM1-BCC runs once per ligand.
Inputs
protein: a prepared receptor PDB with hydrogens (for example from the Protein Preparer);rbfeandrest2.poses: docked poses as SDF, for example from GNINA or AutoDock Vina. Each ligand is named by its SDF title, and the first pose of each name is used.edges(rbfe): the pairs, as[ligand_a, ligand_b]names. Choosing pairs that differ by a small change is your call; the tool does not design the network.ligands(ahfe): 3D ligands as SDF.
Sampling settings are the lambda windows per leg, equilibration and production time per window (or per replica), independent repeats, the temperature, and for rest2 the number of replicas and the hottest effective solute temperature. Their defaults and ranges are listed below.
Outputs
ddg.csv: one row per pair with ΔΔG, its uncertainty, both leg ΔGs, the diagnostics below, the cycle-closure residual, andconverged/status.dg.csv: one ΔG per ligand (network fit or hydration).rest2.csvandframes/: protein–ligand snapshots with the ligand's RMSD from the docked pose and the number of waters within 3.5 Å of it.summary.json(settings, versions, warnings) and the full OpenFE results.
In a workflow, pairs carry rbfe_ddg_kcal_mol, ligands rbfe_dg_kcal_mol or hydration_dg_kcal_mol, and both fe_uncertainty_kcal_mol, so filter steps can rank on them.
Convergence
Each result is checked three ways: its uncertainty (MBAR, or the spread across repeats), the agreement of forward and reverse estimates on half the data, and the overlap between neighbouring lambda windows. A result that fails any check keeps its value but is marked converged=false, with a warning in summary.json; it is never averaged in silently. Treat those values as unconverged and rerun with longer sampling or more windows.
Limits
- One job has a fixed budget of simulated time across all its pairs, windows, repeats and nanoseconds; split a larger network across jobs. Pairs a job cannot start within its time limit are reported as not run, never estimated.
- This is hit validation for tens of pairs, not screening.
- REST2 ligand RMSD is not corrected for symmetry: a flipped ring reads as a large RMSD. Look at the frames before concluding a pose is lost.
Not covered
Absolute protein–ligand binding free energies, membrane or cofactor-dependent targets, QM/MM, and protein–protein binding. For a quick physics re-rank without sampling, use MM-GBSA Rescoring; to run plain MD of a complex, use OpenMM MD.