View on GitHub

Download Notebook

Run in Colab

Create an Alchemical Network

The final setup step is to compile all the different bits of information into a description of a single simulation campaign. This description takes the form of an Alchemical Network. Similarly to LigandNetwork, the AlchemicalNetwork class is a graph of all the transformations in the campaign; however, in an AlchemicalNetwork, these transformations include all the information needed to perform the transformation. By contrast, a LigandNetwork includes only the ligands themselves.

Setup

import openfe, rdkit.Chem
from openff.units import unit
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol

This cookbook assumes you’ve already loaded a LigandNetwork. For more information, see Generate a Ligand Network Automatically:

ligand_network = openfe.ligand_network_planning.generate_minimal_spanning_network(
    ligands=[
        openfe.SmallMoleculeComponent(mol)
        for mol in rdkit.Chem.SDMolSupplier(
            "assets/somebenzenes.sdf",
            removeHs=False,
        )
    ],
    mappers=[openfe.setup.LomapAtomMapper()],
    scorer=openfe.lomap_scorers.default_lomap_score,
)

This cookbook assumes you’ve already loaded a Protocol. For more information, see Choose and Configure a Protocol:

protocol = RelativeHybridTopologyProtocol(RelativeHybridTopologyProtocol.default_settings())

This cookbook assumes you’ve already loaded a solvent and all the other chemical components of the system, including any proteins or cofactors. For more information, see Loading your data into Components:

solvent = openfe.SolventComponent(
    ion_concentration=0.15 * unit.molar
)
protein = openfe.ProteinComponent.from_pdb_file(
    "assets/t4_lysozyme.pdb"
)

Create the Alchemical Network

Automatically

The LigandNetwork.to_rbfe_alchemical_network() method makes constructing alchemical networks for relative binding free energy calculations very simple:

alchemical_network_auto = ligand_network.to_rbfe_alchemical_network(
    solvent=solvent,
    protein=protein,
    protocol=protocol,
)

Manually

If your needs are not catered to by the above method, you can instead loop over the LigandNetwork edges and manually create the Transformation objects for each of them. This gives you full control over the entire network. For more information, see Under the Hood:

# In an RBFE, each edge includes two "legs":
# one for the ligand complexed to the protein,
# and the other for the ligand free in solution
legs = {
    "solvent": {
        # Specify the components common to all systems in this leg
        'solvent': solvent,
    },
    "complex": {
        # Specify the components common to all systems in this leg
        'solvent': solvent,
        'protein': protein,
    }
}
transformations = []

for mapping in ligand_network.edges:
    for leg, common_components in legs.items():
        system_a = openfe.ChemicalSystem(
            {
                'ligand': mapping.componentA,
                **common_components,
            },
            name=f"{mapping.componentA.name}_{leg}"
        )

        system_b = openfe.ChemicalSystem(
            {
                'ligand': mapping.componentB,
                **common_components,
            },
            name=f"{mapping.componentB.name}_{leg}"
        )

        transformation = openfe.Transformation(
            stateA=system_a,
            stateB=system_b,
            mapping=mapping,
            protocol=protocol,
            # Using the same name as to_rbfe_alchemical_network()
            name=f"easy_rbfe_{system_a.name}_{system_b.name}"
        )

        transformations.append(transformation)

Finally, combine the transformations into an Alchemical Network:

alchemical_network = openfe.AlchemicalNetwork(transformations)

We can confirm that this produces identical results to the previous strategy:

assert alchemical_network == ligand_network.to_rbfe_alchemical_network(
    solvent=solvent,
    protein=protein,
    protocol=protocol,
)

Write the Alchemical Network to Disk

While the AlchemicalNetwork class itself has no on-disk representation, the Transformation instances that compose it do. Write an alchemical network to disk by iterating over its edges. For more information, see Dumping a Transformation to JSON:

from pathlib import Path

transformations_dir = Path("transformations")
transformations_dir.mkdir(exist_ok=True)

for n, transformation in enumerate(alchemical_network.edges):
    transformation_name = transformation.name or transformation.key
    transformation.dump(transformations_dir / f"{n}_{transformation_name}.json")