BespokeFit: Using bespoke force field parameters with OpenFE protocols#

This tutorial gives a step-by-step guide on the use of OpenFF-BespokeFit generated force field parameters with OpenFE protocols. Here we will focus on relative binding free energy (RBFE) calculations, but this strategy can be used with any OpenMM-based OpenFE protocol.

Here we will assume you have succesfully planned your RBFE campaign using OpenFE by following the showcase or RBFE tutorial and will be using the TYK2 test system in this example with the planned network provided for you at assets/ligand_network.graphml. You should have also already generated a set of bespoke parameters for your ligand series following the BespokeFit production guide and combined the parameters into a single SMIRNOFF style offxml file following the gathering results guide an example bespoke force field is provied for you at assets/bespoke-1.3.0-default.offxml.

Reading input files#

Before creating the OpenFE simulation inputs we first we need to load up our planned ligand network of transformations and our bespoke force field file.

[ ]:

from gufe import LigandNetwork
from openff.toolkit import ForceField

# load our network file
ligand_network = LigandNetwork.from_graphml(
    open("assets/ligand_network.graphml", "r").read()
)
# load the bespoke force field
bespoke_force_field = ForceField("assets/bespoke-1.3.0-default.offxml")

Running the simulations#

We are now ready to build a new set of transformations which can be executed locally using the OpenFE quickrun CLI. We will be following the python API tutorial to create the transformations, you should check that documentation for more details on the next steps.

[2]:

# create the OpenFE RBFE protocol using our bespoke force field
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol
import openfe


# create the default protocol settings
settings = RelativeHybridTopologyProtocol.default_settings()
# add our new force field as a string
# this avoids the need to move the file around when executing the transformations
settings.forcefield_settings.small_molecule_forcefield = bespoke_force_field.to_string()

# create the protocol
protocol = RelativeHybridTopologyProtocol(settings)

# create the solvent and protein components
solvent = openfe.SolventComponent()
protein = openfe.ProteinComponent.from_pdb_file("assets/tyk2_protein.pdb")

# follow the tutorial to create the AlchemicalNetwork
transformations = []
for mapping in ligand_network.edges:
    for leg in ['solvent', 'complex']:
        # use the solvent and protein created above
        sysA_dict = {'ligand': mapping.componentA,
                     'solvent': solvent}
        sysB_dict = {'ligand': mapping.componentB,
                     'solvent': solvent}

        if leg == 'complex':
            sysA_dict['protein'] = protein
            sysB_dict['protein'] = protein

        # we don't have to name objects, but it can make things (like filenames) more convenient
        sysA = openfe.ChemicalSystem(sysA_dict, name=f"{mapping.componentA.name}_{leg}")
        sysB = openfe.ChemicalSystem(sysB_dict, name=f"{mapping.componentB.name}_{leg}")

        prefix = "rbfe_"  # prefix is only to exactly reproduce CLI

        transformation = openfe.Transformation(
            stateA=sysA,
            stateB=sysB,
            mapping={'ligand': mapping},
            protocol=protocol,  # use protocol created above
            name=f"{prefix}{sysA.name}_{sysB.name}"
        )
        transformations.append(transformation)

network = openfe.AlchemicalNetwork(transformations)

We can now write out each of the transformations to disk for independent execution:

[3]:

import pathlib
# first we create the directory
transformation_dir = pathlib.Path("transformations")
transformation_dir.mkdir(exist_ok=True)

# then we write out each transformation
for transformation in network.edges:
    transformation.dump(transformation_dir / f"{transformation.name}.json")

Recap#

So to recap the workflow can be reduced to the following steps:

Plan the RBFE network
Create a single SMIRNOFF style force field with all of the bespoke parameters for the network using the BespokeFit combine CLI
Store the force field as a string in the OpenFE protocol under the settings.forcefield_settings.small_molecule_forcefield field
Use these settings to create the protcol and create the AlchemicalNetwork following the normal steps

Hopefully its clear that this strategy can be applied to any bespoke parameters you wish to add to the force field not just those from BespokeFit, simply edit your base SMIRNOFF style force field using the OpenFF-Toolkit, set it in the protocol and simulate!