fep

Modify OpenMM systems to depend on FEP lambdas.

Notes

Currently the implementation here is very simple and only assumes morphing full ligands to ideal states. This is sufficient for the current use case of implementing separated topology calculations, but may be expanded in future if full relative FEP calculations are to be supported.

Modules:

• femto –

  A comprehensive toolkit for predicting free energies

Functions:

• apply_fep –

  Modifies an OpenMM system so that different interactions can be scaled by

Attributes:

• LAMBDA_VDW_LIGAND_1 –

  The global parameter used to scale the vdW interactions of ligand 1.

• LAMBDA_VDW_LIGAND_2 –

  The global parameter used to scale the vdW interactions of ligand 2.

• LAMBDA_CHARGES_LIGAND_1 –

  The global parameter used to scale the electrostatic interactions of

• LAMBDA_CHARGES_LIGAND_2 –

  The global parameter used to scale the electrostatic interactions of

LAMBDA_VDW_LIGAND_1 module-attribute

LAMBDA_VDW_LIGAND_1 = 'lambda_vdw_lig_1'

The global parameter used to scale the vdW interactions of ligand 1.

LAMBDA_VDW_LIGAND_2 module-attribute

LAMBDA_VDW_LIGAND_2 = 'lambda_vdw_lig_2'

The global parameter used to scale the vdW interactions of ligand 2.

LAMBDA_CHARGES_LIGAND_1 module-attribute

LAMBDA_CHARGES_LIGAND_1 = 'lambda_charges_lig_1'

The global parameter used to scale the electrostatic interactions of ligand 1.

LAMBDA_CHARGES_LIGAND_2 module-attribute

LAMBDA_CHARGES_LIGAND_2 = 'lambda_charges_lig_2'

The global parameter used to scale the electrostatic interactions of ligand 2.

apply_fep

apply_fep(
    system: System,
    ligand_1_idxs: set[int],
    ligand_2_idxs: set[int] | None,
    config: FEP,
)

Modifies an OpenMM system so that different interactions can be scaled by corresponding lambda parameters.

Notes

• All intra-molecular interactions of alchemical ligands are currently replaced with exceptions and are not scaled by lambda parameters. This will lead to slightly different energies from the original system as cutoffs are not applied to them by OpenMM.
• LJ vdW interactions between ligands and the rest of the system will be replaced with a Beutler-style soft-core LJ potential with a-b-c of 1-1-6 and alpha=0.5.
• Ligand indices must correspond to all atoms in the ligand, alchemically modifying part of a molecule is not yet supported.

Parameters:

• system (System) –

  The system to modify in-place.

• ligand_1_idxs (set[int]) –

  The indices of the ligand atoms whose interactions should be scaled by lambda.

• ligand_2_idxs (set[int] | None) –

  The indices of the ligand atoms whose interactions should be scaled by 1 - lambda.

• config (FEP) –

  Configuration options.

Source code in femto/fe/fep.py

def apply_fep(
    system: openmm.System,
    ligand_1_idxs: set[int],
    ligand_2_idxs: set[int] | None,
    config: femto.fe.config.FEP,
):
    """Modifies an OpenMM system so that different interactions can be scaled by
    corresponding lambda parameters.

    Notes:
        * All intra-molecular interactions of alchemical ligands are currently replaced
          with exceptions and are not scaled by lambda parameters. This will lead to
          slightly different energies from the original system as cutoffs are not
          applied to them by OpenMM.
        * LJ vdW interactions between ligands and the rest of the system will be
          replaced with a Beutler-style soft-core LJ potential with a-b-c of 1-1-6 and
          alpha=0.5.
        * Ligand indices **must** correspond to **all** atoms in the ligand,
          alchemically modifying part of a molecule is not yet supported.

    Args:
        system: The system to modify in-place.
        ligand_1_idxs: The indices of the ligand atoms whose interactions should be
            scaled by lambda.
        ligand_2_idxs: The indices of the ligand atoms whose interactions should be
            scaled by 1 - lambda.
        config: Configuration options.
    """

    forces_by_type = collections.defaultdict(list)

    for force in system.getForces():
        if type(force) not in _SUPPORTED_FORCES:
            raise ValueError(
                f"Force type {type(force)} is not supported when alchemical modifying "
                f"a system."
            )
        forces_by_type[type(force)].append(force)

    updated_forces = []

    for force_type, forces in forces_by_type.items():
        if len(forces) != 1:
            raise NotImplementedError("only one force of each type is supported.")

        force = forces[0]

        if force_type == openmm.NonbondedForce:
            nonbonded_forces = _apply_nonbonded_lambdas(
                force, ligand_1_idxs, ligand_2_idxs, config
            )
            updated_forces.extend(nonbonded_forces)
        else:
            updated_forces.append(copy.deepcopy(force))

    for i in reversed(range(system.getNumForces())):
        system.removeForce(i)

    for force in updated_forces:
        if force is not None:
            system.addForce(force)