neq

Run non-equilibrium forward and reverse sampling.

Functions:

• run_neq –

  Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly

run_neq

run_neq(
    simulation: Simulation,
    coords_0: State,
    coords_1: State,
    protocol: SwitchingProtocol,
) -> tuple[float, float]

Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly pulled along an alchemical pathway as described by Ballard and Jarzynski [1] (Figure 3) and Gapsys et al [2].

Both the forward and reverse directions will be simulated.

References

[1] Ballard, Andrew J., and Christopher Jarzynski. "Replica exchange with nonequilibrium switches: Enhancing equilibrium sampling by increasing replica overlap." The Journal of chemical physics 136.19 (2012): 194101.

[2] Gapsys, Vytautas, et al. "Large scale relative protein ligand binding affinities using non-equilibrium alchemy." Chemical Science 11.4 (2020): 1140-1152.

Returns:

• tuple[float, float] –

  The forward and reverse work values.

Source code in absolv/neq.py

def run_neq(
    simulation: openmm.app.Simulation,
    coords_0: openmm.State,
    coords_1: openmm.State,
    protocol: absolv.config.SwitchingProtocol,
) -> tuple[float, float]:
    """Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly
    pulled along an alchemical pathway as described by Ballard and Jarzynski [1]
    (Figure 3) and Gapsys et al [2].

    Both the forward and reverse directions will be simulated.

    References:
        [1] Ballard, Andrew J., and Christopher Jarzynski. "Replica exchange with
        nonequilibrium switches: Enhancing equilibrium sampling by increasing replica
        overlap." The Journal of chemical physics 136.19 (2012): 194101.

        [2] Gapsys, Vytautas, et al. "Large scale relative protein ligand binding
        affinities using non-equilibrium alchemy." Chemical Science 11.4 (2020):
        1140-1152.

    Returns:
        The forward and reverse work values.
    """

    forward_potentials = _simulate(
        simulation, coords_0, protocol, reverse_direction=False
    )
    reverse_potentials = _simulate(
        simulation, coords_1, protocol, reverse_direction=True
    )

    forward_work = (forward_potentials[:, 1] - forward_potentials[:, 0]).sum()
    reverse_work = (reverse_potentials[:, 1] - reverse_potentials[:, 0]).sum()

    return forward_work, reverse_work