mstk.sfe.SFEManager

class mstk.sfe.SFEManager(system, mol_id, n_lambda=16)

Alchemical system for solvation free energy calculation.

Takes an mstk System and builds an OpenMM system with alchemical forces for coupling a single molecule into its environment. Two-phase protocol: LJ with soft-core first, then Coulomb.

Parameters:

  • system (mstk.simsys.System) – The molecular system (topology + force field).

  • mol_id (int) – Index of the molecule to decouple (0-based).

  • n_lambda (int, optional) – Number of lambda windows (default 16).

omm_system

The alchemical OpenMM system, ready to create a Simulation.

Type:

openmm.System

schedule

Lambda schedule. Each element is (lambda_coul, lambda_vdw).

Type:

list of (float, float)

mol_atoms

Global atom indices of the decoupled molecule.

Type:

list of int

topology

Reference to the original topology.

Type:

mstk.topology.Topology

Examples

>>> # Create alchemical manager
>>> sfe = SFEManager(system, mol_id=0, n_lambda=16)
>>>
>>> # Create simulation with the alchemical OpenMM system
>>> integrator = mm.LangevinMiddleIntegrator(300*unit.kelvin, 1/unit.ps, 0.002*unit.ps)
>>> sim = app.Simulation(top.to_omm_topology(), sfe.omm_system, integrator)
>>> sim.context.setPositions(top.positions)
>>>
>>> # Set lambda state for a specific window
>>> sfe.set_lambda_state(sim.context, 5)
>>>
>>> # Evaluate dU to all windows at current configuration
>>> U = []
>>> for k in range(16):
...     sfe.set_lambda_state(sim.context, k)
...     U.append(sim.context.getState(getEnergy=True).getPotentialEnergy())
>>>
>>> # After collecting dU from all windows, run MBAR
>>> dG, error, dG_adj, dG_adj_err = SFEManager.compute_mbar(u_kn, N_k, 300.0)

Methods

__init__(system, mol_id[, n_lambda])

compute_mbar(u_kn, N_k, temperature)

Run MBAR analysis on a reduced energy matrix.

default_lambda_schedule(n_lambda)

Generate a lambda schedule with n_lambda windows.

set_lambda_state(context, window)

Set the alchemical state of the context by window index.
set_lambda_state(context, window)

Set the alchemical state of the context by window index.

Parameters:

  • context (openmm.Context) –

  • window (int) – Index into self.schedule.

static default_lambda_schedule(n_lambda)

Generate a lambda schedule with n_lambda windows.

First ~3/4 of windows couple LJ (lambda_vdw 0->1, lambda_coul=0). Remaining ~1/4 couple Coulomb (lambda_coul 0->1, lambda_vdw=1).

Parameters:

n_lambda (int) – Total number of lambda windows.

Returns:

schedule – Each element is (lambda_coul, lambda_vdw).

Return type:

list of (float, float)

static compute_mbar(u_kn, N_k, temperature)

Run MBAR analysis on a reduced energy matrix.

Parameters:

  • u_kn (np.ndarray, shape (n_lambda, n_total)) – Reduced potential energy matrix. u_kn[k, n] = beta * dU to state k for sample n.

  • N_k (np.ndarray, shape (n_lambda,)) – Number of samples from each state.

  • temperature (float) – Temperature in K.

Returns:

  • dG (float) – Total free energy difference between first and last window in kJ/mol.

  • error (float) – Statistical uncertainty in kJ/mol.

  • dG_adj (np.ndarray, shape (n_lambda - 1,)) – Free energy difference between adjacent windows in kJ/mol.

  • dG_adj_err (np.ndarray, shape (n_lambda - 1,)) – Uncertainty of dG_adj in kJ/mol.