molforge.wrappers.freeenergy

freeenergy

Concrete endpoint free-energy engines and result parsers.

Wrappers for the external MM/PB(GB)SA tools — Amber's MMPBSA.py and gmx_MMPBSA — that implement :class:molforge.freeenergy.MMGBSAEngine. Currently exposes the Amber output parser; the engines that build inputs and invoke the tools are layered on top.

AmberMMGBSA

AmberMMGBSA(
    *,
    mmpbsa_executable: str = "MMPBSA.py",
    antemmpbsa_executable: str = "ante-MMPBSA.py",
    igb: int = 5,
    strip_mask: str | None = ":WAT,HOH,Na+,Cl-,K+",
    verbose: bool = False,
)

Bases: MMGBSAEngine

Endpoint free energy via Amber's MMPBSA.py.

Post-processes an MD trajectory: splits the complex topology into complex/receptor/ligand with ante-MMPBSA.py (using masks derived from the receptor / ligand selections), runs MMPBSA.py, and parses its FINAL_RESULTS_MMPBSA.dat.

This engine orchestrates the tools; it does not parameterize a system. It needs an Amber topology (prmtop) and a matching trajectory file on disk — either passed explicitly or carried by a trajectory produced by :class:molforge.wrappers.md.AMBER (whose metadata records the run directory holding system.prmtop and prod.nc). Anything else raises a clear error rather than trying to build a topology.

Parameters:

Name	Type	Description	Default
mmpbsa_executable	str	MMPBSA.py binary name or path.	'MMPBSA.py'
antemmpbsa_executable	str	ante-MMPBSA.py binary name or path.	'ante-MMPBSA.py'
igb	int	Generalized Born model index passed to MM/GBSA runs.	5
strip_mask	str | None	Amber mask of atoms stripped from the complex before splitting (solvent/ions), or None to strip nothing.	':WAT,HOH,Na+,Cl-,K+'
verbose	bool	If true, stream tool stdout/stderr instead of capturing.	False

run

run(
    trajectory: Trajectory,
    *,
    receptor: Selection,
    ligand: Selection,
    solvent_model: str = "gb",
    prmtop: str | PathLike[str] | None = None,
    trajectory_file: str | PathLike[str] | None = None,
    start_frame: int = 1,
    end_frame: int | None = None,
    interval: int = 1,
    salt_conc: float = 0.0,
    idecomp: int = 0,
    print_res: str = "within 6",
    **_kwargs: object,
) -> FreeEnergyResult

Estimate ΔG_bind from trajectory with MM/GBSA or MM/PBSA.

Results are cached on the run's :class:~molforge.core.Provenance; an identical repeat returns the cached result without invoking the tools.

Parameters:

Name	Type	Description	Default
trajectory	Trajectory	The ensemble to average over; its topology defines the complex and its metadata may locate the Amber inputs.	required
receptor	Selection	Selection identifying the receptor atoms.	required
ligand	Selection	Selection identifying the ligand atoms.	required
solvent_model	str	"gb" (MM/GBSA, default) or "pb".	'gb'
prmtop	str | PathLike[str] | None	Explicit Amber complex topology; overrides metadata.	None
trajectory_file	str | PathLike[str] | None	Explicit trajectory file; overrides metadata.	None
start_frame	int	First frame to analyze (1-based).	1
end_frame	int | None	Last frame to analyze; defaults to the trajectory length.	None
interval	int	Stride between analyzed frames.	1
salt_conc	float	Salt concentration (mol/L).	0.0
idecomp	int	Per-residue decomposition scheme (0 = off, default; 1/2 per-residue). When set, the result carries a :attr:FreeEnergyResult.decomposition.	0
print_res	str	Which residues to decompose (&decomp print_res, default residues within 6 Å of the interface); only used when idecomp is set.	'within 6'

Returns:

Name	Type	Description
A	FreeEnergyResult	class:FreeEnergyResult from the tool's final block, with
	FreeEnergyResult	provenance attached (and a decomposition if idecomp set).

Raises:

Type	Description
ValueError	If a selection is empty/splits a residue, or the Amber topology / trajectory file can't be located.
MMGBSAEngineNotInstalledError	If the tools aren't installed.

GromacsMMGBSA

GromacsMMGBSA(
    *,
    executable: str = "gmx_MMPBSA",
    igb: int = 5,
    verbose: bool = False,
)

Bases: MMGBSAEngine

Endpoint free energy via gmx_MMPBSA (the GROMACS sibling).

Post-processes a GROMACS MD trajectory: resolves the receptor / ligand selections to .ndx index groups, runs gmx_MMPBSA, and parses its FINAL_RESULTS_MMPBSA.dat. Like the Amber engine it orchestrates the tool rather than parameterizing a system: it needs a GROMACS structure (.tpr) and trajectory (.xtc/.trr) on disk — passed explicitly or carried by a trajectory produced by :class:molforge.wrappers.md.GROMACS (whose metadata records the run directory holding md.tpr, md.xtc and topol.top). A topology (.top) is used with -cp when available. Results are cached on the run's :class:~molforge.core.Provenance.

Parameters:

Name	Type	Description	Default
executable	str	gmx_MMPBSA binary name or path.	'gmx_MMPBSA'
igb	int	Generalized Born model index passed to MM/GBSA runs.	5
verbose	bool	If true, stream tool stdout/stderr instead of capturing.	False

run

run(
    trajectory: Trajectory,
    *,
    receptor: Selection,
    ligand: Selection,
    solvent_model: str = "gb",
    structure: str | PathLike[str] | None = None,
    trajectory_file: str | PathLike[str] | None = None,
    topology: str | PathLike[str] | None = None,
    start_frame: int = 1,
    end_frame: int | None = None,
    interval: int = 1,
    salt_conc: float = 0.0,
    idecomp: int = 0,
    print_res: str = "within 6",
    **_kwargs: object,
) -> FreeEnergyResult

Estimate ΔG_bind from trajectory with gmx_MMPBSA.

An identical repeat returns the cached result without invoking the tool.

Parameters:

Name	Type	Description	Default
trajectory	Trajectory	The ensemble to average over; its topology defines the atom numbering and its metadata may locate the GROMACS inputs.	required
receptor	Selection	Selection identifying the receptor atoms.	required
ligand	Selection	Selection identifying the ligand atoms.	required
solvent_model	str	"gb" (MM/GBSA, default) or "pb".	'gb'
structure	str | PathLike[str] | None	Explicit GROMACS structure (.tpr); overrides metadata.	None
trajectory_file	str | PathLike[str] | None	Explicit trajectory (.xtc/.trr).	None
topology	str | PathLike[str] | None	Explicit GROMACS topology (.top); passed with -cp when present.	None
start_frame	int	First frame to analyze (1-based).	1
end_frame	int | None	Last frame to analyze; defaults to the length.	None
interval	int	Stride between analyzed frames.	1
salt_conc	float	Salt concentration (mol/L).	0.0
idecomp	int	Per-residue decomposition scheme (0 = off, default; 1/2 per-residue). When set, the result carries a :attr:FreeEnergyResult.decomposition.	0
print_res	str	Which residues to decompose (&decomp print_res, default residues within 6 Å of the interface); only used when idecomp is set.	'within 6'

Returns:

Name	Type	Description
A	FreeEnergyResult	class:FreeEnergyResult from the tool's final block, with
	FreeEnergyResult	provenance attached (and a decomposition if idecomp set).

Raises:

Type	Description
ValueError	If a selection is empty, or the structure / trajectory can't be located.
MMGBSAEngineNotInstalledError	If gmx_MMPBSA isn't installed.

absolute_binding_free_energy

absolute_binding_free_energy(
    complex_leg: FreeEnergyResult,
    solvent_leg: FreeEnergyResult,
    *,
    restraint_correction: FreeEnergyResult | float = 0.0,
    method: str = "ABFE",
    metadata: Mapping[str, object] | None = None,
) -> FreeEnergyResult

Close a double-decoupling cycle into an absolute binding ΔG.

Absolute binding free energy (ABFE) by double decoupling annihilates the ligand's interactions with its environment in two phases — restrained in the complex, and free in solvent — and adds a standard-state restraint correction. The thermodynamic cycle gives

ΔG_bind = ΔG_solvent − ΔG_complex + restraint_correction

where the two legs are decoupling free energies (coupled → non-interacting): a strong binder is hard to decouple from the complex (large positive ΔG_complex), so ΔG_solvent − ΔG_complex comes out negative — favorable — as it should.

Unlike the relative cycle, this yields an absolute ΔG_bind, so the result is a :class:FreeEnergyResult that ranks directly.

Convention notes:

• Pass the legs as decoupling free energies. If your λ schedule runs the other way (coupling), negate them (or reverse the λ order before :func:from_alchemlyb).
• restraint_correction is a signed contribution added as-is — supply it with the sign your protocol uses (e.g. a Boresch standard-state correction). It may be a float or a :class:FreeEnergyResult (whose uncertainty then propagates).

Parameters:

Name	Type	Description	Default
complex_leg	FreeEnergyResult	ΔG of decoupling the (restrained) ligand in the complex.	required
solvent_leg	FreeEnergyResult	ΔG of decoupling the ligand in solvent.	required
restraint_correction	FreeEnergyResult | float	Signed standard-state / restraint term.	0.0
method	str	Value for :attr:FreeEnergyResult.method.	'ABFE'
metadata	Mapping[str, object] | None	Extra items merged into the result's metadata.	None

Returns:

Name	Type	Description
A	FreeEnergyResult	class:FreeEnergyResult with the absolute ΔG_bind in kcal/mol
	FreeEnergyResult	(components is None) and the three terms' errors propagated
	FreeEnergyResult	in quadrature.

from_alchemlyb

from_alchemlyb(
    estimator: Any,
    *,
    temperature: float | None = None,
    energy_unit: str | None = None,
    method: str | None = None,
    metadata: Mapping[str, object] | None = None,
) -> FreeEnergyResult

Ingest a fitted alchemlyb estimator into a :class:FreeEnergyResult.

Reads estimator.delta_f_ / estimator.d_delta_f_ and, unless overridden, labels the result with the estimator's class name ("MBAR", "BAR", "TI", …).

Parameters:

Name	Type	Description	Default
estimator	Any	A fitted alchemlyb estimator (anything exposing delta_f_ and d_delta_f_).	required
temperature	float | None	See :func:from_delta_f.	None
energy_unit	str | None	See :func:from_delta_f.	None
method	str | None	Overrides the estimator-class-name default.	None
metadata	Mapping[str, object] | None	Extra items merged into the result's metadata; the estimator name is always recorded under "estimator".	None

Returns:

Name	Type	Description
A	FreeEnergyResult	class:FreeEnergyResult in kcal/mol.

Raises:

Type	Description
AttributeError	If estimator lacks delta_f_ / d_delta_f_.
ValueError	Propagated from :func:from_delta_f.

from_delta_f

from_delta_f(
    delta_f: Any,
    d_delta_f: Any,
    *,
    temperature: float | None = None,
    energy_unit: str | None = None,
    method: str = "FEP",
    metadata: Mapping[str, object] | None = None,
) -> FreeEnergyResult

Ingest alchemlyb delta_f_ / d_delta_f_ matrices.

Takes the full first-state → last-state transformation (the [0, -1] corner) and its error, converting to kcal/mol.

Parameters:

Name	Type	Description	Default
delta_f	Any	The estimator's delta_f_ matrix (a pandas DataFrame, or any 2-D array-like). λ states must be ordered, so the corner is the whole transformation.	required
d_delta_f	Any	The matching d_delta_f_ error matrix.	required
temperature	float | None	Temperature in K, needed only when the unit is kT and the matrix doesn't carry a temperature in .attrs. An explicit value overrides .attrs.	None
energy_unit	str | None	"kT" (default when unknown), "kcal/mol", or "kJ/mol". An explicit value overrides .attrs; otherwise the matrix's .attrs['energy_unit'] is used.	None
method	str	Value for :attr:FreeEnergyResult.method (e.g. "FEP", "TI"); :func:from_alchemlyb fills the estimator name.	'FEP'
metadata	Mapping[str, object] | None	Extra items merged into the result's metadata.	None

Returns:

Name	Type	Description
A	FreeEnergyResult	class:FreeEnergyResult with delta_g / uncertainty in
	FreeEnergyResult	kcal/mol and components set to None — an alchemical ΔG is
	FreeEnergyResult	a single number, not an MM/GBSA-style term breakdown.

Raises:

Type	Description
ValueError	If a matrix isn't 2-D, the unit is unknown, or a kT matrix has no temperature (neither argument nor .attrs).

relative_binding_free_energy

relative_binding_free_energy(
    complex_leg: FreeEnergyResult,
    solvent_leg: FreeEnergyResult,
    *,
    reference: str,
    other: str,
) -> DeltaDeltaG

Close a relative-FEP thermodynamic cycle into a binding ΔΔG.

A relative FEP perturbation transforms one ligand into another (reference → other) along two legs: bound to the receptor (the complex leg) and free in solution (the solvent leg). The thermodynamic cycle gives the relative binding free energy

ΔΔG_bind = ΔG_complex − ΔG_solvent
         = ΔG_bind(other) − ΔG_bind(reference)

so a single perturbation's two ingested legs become the binding ΔΔG between the two ligands — the quantity relative FEP actually reports. A single leg's ΔG is not a binding affinity; only the cycle is.

Both legs must be the same perturbation direction (reference → other) and in the same units (kcal/mol, as returned by :func:from_alchemlyb / :func:from_delta_f).

Parameters:

Name	Type	Description	Default
complex_leg	FreeEnergyResult	ΔG of the reference → other transformation in the complex.	required
solvent_leg	FreeEnergyResult	ΔG of the same transformation in solvent.	required
reference	str	Label of the reference ligand.	required
other	str	Label of the ligand it is perturbed into.	required

Returns:

Name	Type	Description
A	DeltaDeltaG	class:~molforge.freeenergy.DeltaDeltaG: the signed ΔΔG_bind
	DeltaDeltaG	(negative means other binds more tightly) with the two legs'
	DeltaDeltaG	errors propagated in quadrature. Treat a difference within its
	DeltaDeltaG	uncertainty as a tie.

build_mmpbsa_input

build_mmpbsa_input(
    *,
    solvent_model: str = "gb",
    start_frame: int = 1,
    end_frame: int,
    interval: int = 1,
    salt_conc: float = 0.0,
    igb: int = 5,
    verbose: int = 1,
    idecomp: int = 0,
    dec_verbose: int = 0,
    print_res: str = "within 6",
) -> str

Build the mmpbsa.in namelist text for an MM/PB(GB)SA run.

Shared by the Amber and gmx_MMPBSA engines — both consume the same &general / &gb / &pb (and optional &decomp) namelists.

Parameters:

Name	Type	Description	Default
solvent_model	str	"gb" (writes a &gb namelist, MM/GBSA) or "pb" (writes a &pb namelist, MM/PBSA).	'gb'
start_frame	int	First trajectory frame to analyze (1-based).	1
end_frame	int	Last trajectory frame to analyze (inclusive).	required
interval	int	Stride between analyzed frames.	1
salt_conc	float	Salt concentration in mol/L (saltcon for GB, istrng for PB).	0.0
igb	int	Generalized Born model index (GB only; 5 = OBC-II).	5
verbose	int	MMPBSA verbose level.	1
idecomp	int	Per-residue decomposition scheme. 0 (default) writes no &decomp namelist; 1 or 2 request a per-residue decomposition (1 folds 1-4 terms into internal, 2 into electrostatic/vdw). 3/4 are the pairwise schemes.	0
dec_verbose	int	&decomp output level (0 = total contributions).	0
print_res	str	Which residues to decompose (&decomp print_res); the default selects residues within 6 Å of the interface.	'within 6'

Returns:

Type	Description
str	The input-file text, ready to write to mmpbsa.in.

Raises:

Type	Description
ValueError	On an unknown solvent_model, an invalid frame range / interval, or an out-of-range idecomp.

parse_mmpbsa_dat

parse_mmpbsa_dat(
    text: str, *, solvent_model: str = "gb"
) -> FreeEnergyResult

Parse MMPBSA.py output text into a :class:FreeEnergyResult.

Parameters:

Name	Type	Description	Default
text	str	Full contents of a FINAL_RESULTS_MMPBSA.dat file.	required
solvent_model	str	"gb" to read the Generalized Born section (MM/GBSA, default) or "pb" for Poisson–Boltzmann (MM/PBSA).	'gb'

Returns:

Type	Description
FreeEnergyResult	A result whose delta_g and uncertainty come from the
FreeEnergyResult	DELTA TOTAL row (Average and Std. Err. of Mean) and whose
FreeEnergyResult	components hold the per-term breakdown. entropy is
FreeEnergyResult	None — the entropy section is not parsed here — so
FreeEnergyResult	delta_g is the enthalpic binding total. Frame count and the
FreeEnergyResult	DELTA TOTAL standard deviation are recorded in metadata.

Raises:

Type	Description
ValueError	If solvent_model is unknown, the requested section is absent, or a required term row is missing.

parse_mmpbsa_decomp

parse_mmpbsa_decomp(
    text: str, *, section: str = "delta"
) -> Decomposition

Parse a MMPBSA.py per-residue decomposition.

Reads FINAL_DECOMP_MMPBSA.dat (written when idecomp is set) and returns the per-residue contributions from the requested species' "Total Energy Decomposition" block. The default "delta" section is the binding contribution (complex − receptor − ligand) — the one that answers which residues drive the affinity.

Parameters:

Name	Type	Description	Default
text	str	Contents of FINAL_DECOMP_MMPBSA.dat.	required
section	str	Which species block to read — "delta" (default), "complex", "receptor", or "ligand".	'delta'

Returns:

Name	Type	Description
A	Decomposition	class:~molforge.freeenergy.Decomposition over the residues in
	Decomposition	that block, in report order.

Raises:

Type	Description
ValueError	If section is unknown or the block is absent.

selection_to_amber_mask

selection_to_amber_mask(
    topology: Protein, selection: Selection
) -> str

Convert a molforge selection to an Amber residue mask.

The selection (a field-filter mapping like {"entity_type": "ligand"}, forwarded to :meth:~molforge.core.AtomArray.where, or a boolean atom mask) is resolved against topology and expressed as a residue-number mask over the topology's sequential 1-based residue numbering — e.g. ":1-120" or ":121". This is the numbering Amber's ambmask uses, and the form MM/PB(GB)SA expects for splitting a complex into receptor and ligand.

Parameters:

Name	Type	Description	Default
topology	Protein	Structure the selection is resolved against (the complex).	required
selection	Selection	Field filters or a boolean atom mask.	required

Returns:

Type	Description
str	An Amber residue mask string beginning with ":".

Raises:

Type	Description
ValueError	If the selection matches no atoms, or splits a residue (endpoint masks must cover whole residues).

from_cinnabar

from_cinnabar(
    source: Any,
    *,
    computational_only: bool = True,
    method: str = "FEP (network)",
    metadata: Mapping[str, object] | None = None,
) -> dict[str, FreeEnergyResult]

Ingest cinnabar's per-ligand absolute estimates.

Parameters:

Name	Type	Description	Default
source	Any	Either a cinnabar FEMap (with absolute values already generated via generate_absolute_values()) or the DataFrame its get_absolute_dataframe() returns. A FEMap is detected by that method and read through it.	required
computational_only	bool	Keep only calculated ligands, dropping any experimental reference rows (the computational column). If that column is absent, all rows are kept.	True
method	str	Value for :attr:FreeEnergyResult.method.	'FEP (network)'
metadata	Mapping[str, object] | None	Extra items merged into every result's metadata.	None

Returns:

Type	Description
dict[str, FreeEnergyResult]	{label: FreeEnergyResult} in kcal/mol (components is
dict[str, FreeEnergyResult]	None). Wrap in :class:~molforge.freeenergy.FreeEnergyRanking
dict[str, FreeEnergyResult]	to rank. Because the absolute values share one network offset,
dict[str, FreeEnergyResult]	ranking and every pairwise ΔΔG are unaffected by it.

Raises:

Type	Description
ValueError	If the label / ΔG / uncertainty columns can't be found.

parse_gmx_mmpbsa_dat

parse_gmx_mmpbsa_dat(
    text: str, *, solvent_model: str = "gb"
) -> FreeEnergyResult

Parse gmx_MMPBSA output text into a :class:FreeEnergyResult.

Parameters:

Name	Type	Description	Default
text	str	Full contents of a gmx_MMPBSA FINAL_RESULTS_MMPBSA.dat.	required
solvent_model	str	"gb" (MM/GBSA, default) or "pb" (MM/PBSA).	'gb'

Returns:

Type	Description
FreeEnergyResult	A result whose delta_g / uncertainty come from the
FreeEnergyResult	ΔTOTAL row (Average and the final SEM column) and whose
FreeEnergyResult	components hold the per-term breakdown; entropy is
FreeEnergyResult	None. Frame count and the ΔTOTAL sample standard deviation
FreeEnergyResult	are recorded in metadata.

Raises:

Type	Description
ValueError	If solvent_model is unknown, the requested section is absent, or a required term row is missing.

parse_gmx_mmpbsa_decomp

parse_gmx_mmpbsa_decomp(
    text: str, *, section: str = "delta"
) -> Decomposition

Parse a gmx_MMPBSA per-residue decomposition.

gmx_MMPBSA writes the same FINAL_DECOMP_MMPBSA.dat structure as MMPBSA.py (it reuses that writer), so this reads the requested species' "Total Energy Decomposition" block just like :func:~molforge.wrappers.freeenergy.parse_mmpbsa_decomp. The default "delta" section is the per-residue binding contribution.

gmx's delta rows carry an extra Location column (the residue's spot in the receptor/ligand topology); it is dropped, and the residue keeps its complex-numbering resname resnum label.

Parameters:

Name	Type	Description	Default
text	str	Contents of FINAL_DECOMP_MMPBSA.dat.	required
section	str	Which species block to read — "delta" (default), "complex", "receptor", or "ligand".	'delta'

Returns:

Name	Type	Description
A	Decomposition	class:~molforge.freeenergy.Decomposition over the residues in
	Decomposition	that block, in report order.

Raises:

Type	Description
ValueError	If section is unknown or the block is absent.

selection_to_ndx_group

selection_to_ndx_group(
    topology: Protein,
    selection: Selection,
    name: str,
    *,
    per_line: int = 15,
) -> str

Render a molforge selection as a GROMACS .ndx index group.

gmx_MMPBSA identifies the receptor and ligand by index-group number (-cg <receptor> <ligand>), so a selection is resolved against the topology and written as a named group of 1-based atom numbers — the numbering GROMACS uses, which matches the topology's atom order (as produced by the GROMACS MD wrapper).

Parameters:

Name	Type	Description	Default
topology	Protein	Structure the selection is resolved against.	required
selection	Selection	Field filters or a boolean atom mask.	required
name	str	Group name, written as [ name ].	required
per_line	int	Atom numbers per line (GROMACS wraps long groups).	15

Returns:

Type	Description
str	The .ndx group block (header line plus wrapped atom numbers).

Raises:

Type	Description
ValueError	If the selection matches no atoms.