fep

Prepare OpenMM systems for FEP calculations.

Functions:

• apply_fep –

  Generate a system whereby a number of the molecules can be alchemically

• set_fep_lambdas –

  Set the values of the alchemical lambdas on an OpenMM context.

apply_fep

apply_fep(
    system: System,
    alchemical_indices: list[set[int]],
    persistent_indices: list[set[int]],
    custom_alchemical_potential: str | None = None,
) -> System

Generate a system whereby a number of the molecules can be alchemically transformed from a base chemical system.

Notes

• Currently only OpenMM systems that have:

• vdW + electrostatics in a single built-in non-bonded force

• electrostatics in a built-in non-bonded force, vdW in a custom non-bonded force and vdW 1-4 interactions in a custom bond force
• all of the above sans any electrostatics

are supported. * By default a soft-core version of the LJ potential with a-b-c of 1-1-6 and alpha=0.5 that can be scaled by a global lambda_sterics parameter will be used for alchemical-chemical vdW interactions embedded in an OpenMM NonbondedForce while the energy expression set on a CustomNonbondedForce will be be modified to have the form "lambda_sterics*({original_expression})".

Parameters:

• system (System) –

  The chemical system to generate the alchemical system from

• alchemical_indices (list[set[int]]) –

  The atom indices corresponding to each molecule that should be alchemically transformable. The atom indices must correspond to all atoms / v-sites in each molecule as alchemically transforming part of a molecule is not supported.

• persistent_indices (list[set[int]]) –

  The atom indices corresponding to each molecule that should not be alchemically transformable.

• custom_alchemical_potential (str | None, default: None ) –

  A custom expression to use for the potential energy function that describes the chemical-alchemical intermolecular interactions. See the Notes for information about the default value.

  The expression must include "lambda_sterics".

Source code in absolv/fep.py

def apply_fep(
    system: openmm.System,
    alchemical_indices: list[set[int]],
    persistent_indices: list[set[int]],
    custom_alchemical_potential: str | None = None,
) -> openmm.System:
    """Generate a system whereby a number of the molecules can be alchemically
    transformed from a base chemical system.

    Notes:
        * Currently only OpenMM systems that have:

          - vdW + electrostatics in a single built-in non-bonded force
          - electrostatics in a built-in non-bonded force, vdW in a custom
            **non-bonded** force and vdW 1-4 interactions in a custom **bond** force
          - all of the above sans any electrostatics

          are supported.
        * By default a soft-core version of the LJ potential with a-b-c of 1-1-6
          and alpha=0.5 that can be scaled by a global `lambda_sterics` parameter
          will be used for alchemical-chemical vdW interactions embedded in an
          OpenMM ``NonbondedForce`` while the energy expression set on a
          ``CustomNonbondedForce`` will be be modified to have the form
          ``"lambda_sterics*({original_expression})"``.

    Args:
        system: The chemical system to generate the alchemical system from
        alchemical_indices: The atom indices corresponding to each molecule that
            should be alchemically transformable. The atom indices **must**
            correspond to  **all** atoms / v-sites in each molecule as alchemically
            transforming part of a molecule is not supported.
        persistent_indices: The atom indices corresponding to each molecule that
            should **not** be alchemically transformable.
        custom_alchemical_potential: A custom expression to use for the potential
            energy function that describes the chemical-alchemical intermolecular
            interactions. See the Notes for information about the default value.

            The expression **must** include ``"lambda_sterics"``.
    """

    system = copy.deepcopy(system)

    (
        nonbonded_force,
        custom_nonbonded_force,
        custom_bond_force,
    ) = _find_nonbonded_forces(system)

    if nonbonded_force is not None:
        _add_electrostatics_lambda(nonbonded_force, alchemical_indices)

    if custom_nonbonded_force is not None:
        for alchemical_force in _add_custom_vdw_lambda(
            custom_nonbonded_force,
            alchemical_indices,
            persistent_indices,
            custom_alchemical_potential,
        ):
            system.addForce(alchemical_force)

    elif nonbonded_force is not None:
        for alchemical_force in _add_lj_vdw_lambda(
            nonbonded_force,
            alchemical_indices,
            persistent_indices,
            custom_alchemical_potential,
        ):
            system.addForce(alchemical_force)

    else:
        raise NotImplementedError

    return system

set_fep_lambdas

set_fep_lambdas(
    context: Context,
    lambda_sterics: float | None = None,
    lambda_electrostatics: float | None = None,
)

Set the values of the alchemical lambdas on an OpenMM context.

Parameters:

• context (Context) –

  The context to update.

• lambda_sterics (float | None, default: None ) –

  The (optional) value of the steric lambda.

• lambda_electrostatics (float | None, default: None ) –

  The (optional) value of the electrostatics lambda.

Source code in absolv/fep.py

def set_fep_lambdas(
    context: openmm.Context,
    lambda_sterics: float | None = None,
    lambda_electrostatics: float | None = None,
):
    """Set the values of the alchemical lambdas on an OpenMM context.

    Args:
        context: The context to update.
        lambda_sterics: The (optional) value of the steric lambda.
        lambda_electrostatics: The (optional) value of the electrostatics lambda.
    """

    if lambda_sterics is not None:
        assert (
            0.0 <= lambda_sterics <= 1.0
        ), f"`{LAMBDA_STERICS}` must be between 0 and 1"
        context.setParameter(LAMBDA_STERICS, lambda_sterics)

    if lambda_electrostatics is not None:
        assert (
            0.0 <= lambda_electrostatics <= 1.0
        ), f"`{LAMBDA_ELECTROSTATICS}` must be between 0 and 1"

        context.setParameter(LAMBDA_ELECTROSTATICS, lambda_electrostatics)