View on GitHub

Download Notebook

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. Similar 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.

See the User Guide: Alchemical Networks for more background information.

This notebook creates an Hybrid Topology RBFE AlchemicalNetwork as an example. For other protocol types, see our step-by-step tutorials on running SepTop RBFE and ABFE calculations.

Start with a LigandNetwork

This cookbook assumes you’ve already loaded a LigandNetwork.

To learn how to build a LigandNetwork, see Planning a Ligand Network:

import openfe
import rdkit.Chem

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.KartografAtomMapper()],
    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:

from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol

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:

from openff.units import unit

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

Create the Alchemical Network

Manually

To create an AlchemicalNetwork 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 = {
    # Specify the components common to all systems in this leg
    "solvent": {'solvent': solvent},
    # Specify the components common to all systems in this leg
    "complex": {'solvent': solvent, 'protein': protein}
}

# create a list of transformations

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,
            name=f"rbfe_{system_a.name}_{system_b.name}"
        )

        transformations.append(transformation)

Finally, combine the transformations into an Alchemical Network:

alchemical_network = openfe.AlchemicalNetwork(transformations)
alchemical_network

<AlchemicalNetwork-130d122e6027d856bebe9a480b35f363>

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.to_json(transformations_dir / f"{n}_{transformation_name}.json")