diff --git a/.binder/environment.yml b/.binder/environment.yml index 810c274..47c765c 100644 --- a/.binder/environment.yml +++ b/.binder/environment.yml @@ -2,6 +2,7 @@ name: openfe-notebooks channels: - jaimergp/label/unsupported-cudatoolkit-shim - conda-forge + - openeye dependencies: - MDAnalysis - click @@ -29,3 +30,11 @@ dependencies: - typing_extensions - gufe==0.6.* - openfe==0.6.* + ## needed for perses + - openmoltools + - cloudpathlib + - dask + - distributed + - openeye-toolkits + + diff --git a/networks/Preparing AlchemicalNetworks.ipynb b/networks/Preparing AlchemicalNetworks.ipynb new file mode 100644 index 0000000..ea1de86 --- /dev/null +++ b/networks/Preparing AlchemicalNetworks.ipynb @@ -0,0 +1,844 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4ed54201-aa1f-4ba9-bdb7-1c7f5a2cbd05", + "metadata": {}, + "source": [ + "# Preparing `AlchemicalNetwork`s for use with `fah-alchemy`" + ] + }, + { + "cell_type": "markdown", + "id": "6da80e66-6e17-407f-9ba4-a25ded4777bc", + "metadata": {}, + "source": [ + "`fah-alchemy` is a platform for evluating the free energy differences between chemical systems in an alchemical network.\n", + "This notebook will illustrate how to build alchemical networks suitable for submission to a deployed `fah-alchemy` instance." + ] + }, + { + "cell_type": "markdown", + "id": "9d6eaa9e-daa0-4fee-8e72-d1116a1d4954", + "metadata": {}, + "source": [ + "`fah-alchemy` works in terms of `gufe` objects; the `gufe` module defines the data model for `AlchemicalNetwork`s and all objects they are composed of. We'll import the classes of objects we'll use in this tutorial here." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "27b4aa61-ca82-44fa-be84-c802f1083a29", + "metadata": {}, + "outputs": [], + "source": [ + "# suppress `numba` warnings, if present\n", + "from numba.core.errors import NumbaWarning\n", + "import warnings\n", + "\n", + "warnings.simplefilter('ignore', category=NumbaWarning)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "6bb3f4a4-2135-494a-8365-9ef229dd124d", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "LICENSE: Could not open license file \"oe_license.txt\" in local directory\n", + "LICENSE: N.B. OE_LICENSE environment variable is not set\n", + "LICENSE: N.B. OE_DIR environment variable is not set\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n" + ] + } + ], + "source": [ + "from gufe import AlchemicalNetwork, Transformation, ChemicalSystem\n", + "from gufe.components import ProteinComponent, SmallMoleculeComponent, SolventComponent\n", + "\n", + "from openff.units import unit" + ] + }, + { + "cell_type": "markdown", + "id": "99190957-890d-4729-ac72-bef35fbb212f", + "metadata": {}, + "source": [ + "## Sample network from `openfe-benchmark`" + ] + }, + { + "cell_type": "markdown", + "id": "c7362410-7045-42a6-8249-93a1e832d8c9", + "metadata": {}, + "source": [ + "We'll use a sample network in `openfe-benchmark` for demonstration purposes. The sources can be found here: https://github.com/OpenFreeEnergy/openfe-benchmarks\n", + "\n", + "In particular, we'll use the `tyk2` network. We'll extract ligands manually from the ligand SDF, and the protein target from its PDB." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "08f7ddfd-4fd5-4bef-8bf9-b2a286b094c6", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING:root:Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.\n" + ] + } + ], + "source": [ + "from importlib import resources\n", + "from rdkit import Chem\n", + "\n", + "from openfe_benchmarks import tyk2" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "d3d7946f-6b31-4dfd-b337-7a41d717e391", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tyk2_system = tyk2.get_system()\n", + "tyk2_system" + ] + }, + { + "cell_type": "markdown", + "id": "7aff6317-abe0-40df-99db-9d1eb121b20e", + "metadata": {}, + "source": [ + "The connections for the network are defined here; we'll use these for building up our own `AlchemicalNetwork`." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "ec682c29-12cf-46dc-aea7-8b7e45324529", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[('lig_ejm_31', 'lig_ejm_50'),\n", + " ('lig_ejm_46', 'lig_jmc_23'),\n", + " ('lig_ejm_31', 'lig_ejm_55'),\n", + " ('lig_ejm_31', 'lig_ejm_48'),\n", + " ('lig_ejm_31', 'lig_ejm_54'),\n", + " ('lig_ejm_31', 'lig_ejm_47'),\n", + " ('lig_ejm_31', 'lig_ejm_46'),\n", + " ('lig_ejm_46', 'lig_jmc_27'),\n", + " ('lig_ejm_46', 'lig_jmc_28'),\n", + " ('lig_ejm_42', 'lig_ejm_43'),\n", + " ('lig_ejm_31', 'lig_ejm_42'),\n", + " ('lig_ejm_45', 'lig_ejm_55')]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tyk2_system.connections" + ] + }, + { + "cell_type": "markdown", + "id": "5d3285c5-0b7e-461f-a761-dabd01182775", + "metadata": {}, + "source": [ + "## Define `ChemicalSystem`s for network nodes" + ] + }, + { + "cell_type": "markdown", + "id": "90ddf9d5-71d0-44ed-8eea-5ea4d512f615", + "metadata": {}, + "source": [ + "An `AlchemicalNetwork` features `ChemicalSystem`s as nodes and `Transformation`s as directed edges between nodes. We'll start by defining the nodes for our network.\n", + "\n", + "A `ChemicalSystem` is made of one or more `Component`s. These can be one of `ProteinComponent`, `SmallMoleculeComponent`, or `SolventComponent`, and potentially others as needed. This design allows for memory efficient representation of large networks with perhaps hundreds or thousands of nodes, but perhaps far fewer variants in proteins, ligands, etc." + ] + }, + { + "cell_type": "markdown", + "id": "0b64f1c8-0fa0-4b10-baba-7014f9de9440", + "metadata": {}, + "source": [ + "### Define `Component`s for a given `ChemicalSystem`" + ] + }, + { + "cell_type": "markdown", + "id": "90fd672d-ce78-48de-811a-7490ada1ea3b", + "metadata": {}, + "source": [ + "Let's start by assembling the ligands. These are defined as `SmallMoleculeComponent`s, and can be initialized with RDKit molecules. \n", + "\n", + "We'll read a multimolecule SDF from `openfe-benchmarks` and create a `SmallMoleculeComponent` for each ligand in the file:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "7880e75b-d28d-411c-b5c1-3aafbed31099", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[SmallMoleculeComponent(name=lig_ejm_31),\n", + " SmallMoleculeComponent(name=lig_ejm_42),\n", + " SmallMoleculeComponent(name=lig_ejm_43),\n", + " SmallMoleculeComponent(name=lig_ejm_45),\n", + " SmallMoleculeComponent(name=lig_ejm_46),\n", + " SmallMoleculeComponent(name=lig_ejm_47),\n", + " SmallMoleculeComponent(name=lig_ejm_48),\n", + " SmallMoleculeComponent(name=lig_ejm_50),\n", + " SmallMoleculeComponent(name=lig_ejm_54),\n", + " SmallMoleculeComponent(name=lig_ejm_55),\n", + " SmallMoleculeComponent(name=lig_jmc_23),\n", + " SmallMoleculeComponent(name=lig_jmc_27),\n", + " SmallMoleculeComponent(name=lig_jmc_28)]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "with resources.path('openfe_benchmarks.data',\n", + " 'tyk2_ligands.sdf') as fn:\n", + " ligands_sdf = Chem.SDMolSupplier(str(fn), removeHs=False)\n", + " ligands = [SmallMoleculeComponent(rdkit_ligand) for rdkit_ligand in ligands_sdf]\n", + "\n", + "ligands" + ] + }, + { + "cell_type": "markdown", + "id": "d0f8e9a5-80a1-4eef-86e9-b47a9ba1f9c5", + "metadata": {}, + "source": [ + "We'll also load our protein into a `ProteinComponent`:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "8ffcd2d6-47ec-46ee-bfe8-8182915f01bd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "ProteinComponent(name=tyk2)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "with resources.path('openfe_benchmarks.data',\n", + " 'tyk2_protein.pdb') as fn:\n", + " protein = ProteinComponent.from_pdb_file(str(fn), name='tyk2')\n", + "\n", + "protein" + ] + }, + { + "cell_type": "markdown", + "id": "e653c208-1f69-4f08-a38e-55b7123292d9", + "metadata": {}, + "source": [ + "We'll also need at least one `SolventComponent` to encode our choice of solvent and counterions, with concentration:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "3c47166d-143b-41f0-9dce-7805e56d7191", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "SolventComponent(name=O, Na+, Cl-)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solvent = SolventComponent(positive_ion='Na', \n", + " negative_ion='Cl',\n", + " neutralize=True, \n", + " ion_concentration=0.15*unit.molar)\n", + "solvent" + ] + }, + { + "cell_type": "markdown", + "id": "9306e96d-f257-4b71-83e1-9814701570cd", + "metadata": {}, + "source": [ + "The `SolventComponent` doesn't actually perform any actual solvation (packing water molecules, ions); that is performed just before simulation time during `Protocol` execution." + ] + }, + { + "cell_type": "markdown", + "id": "ea4dfd8d-4ab6-44c8-a268-abb3d1505803", + "metadata": {}, + "source": [ + "Each of the ligands have been pre-docked into the protein and aligned to their common scaffold. It is important to recognize that any processing required to prepare ligand and protein structures for alchemical free energy calculations should be done *before* the steps we are taking here." + ] + }, + { + "cell_type": "markdown", + "id": "b0b11357-2179-496b-9bff-89303e1c9c33", + "metadata": {}, + "source": [ + "### Build the `ChemicalSystem`s" + ] + }, + { + "cell_type": "markdown", + "id": "c314a1a2-afcc-4390-9ed0-16fae3d48c6e", + "metadata": {}, + "source": [ + "We can now construct the `ChemicalSystem`s we want represented in our network. Since we are planning to perform relative binding free energy (RBFE) calculations, we'll define both *complex* and *solvent* variants for each ligand." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "10edc8b3-5f45-4e37-b1d5-bd508601e352", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'lig_ejm_31': ChemicalSystem(name=lig_ejm_31_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_31), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_42': ChemicalSystem(name=lig_ejm_42_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_42), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_43': ChemicalSystem(name=lig_ejm_43_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_43), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_45': ChemicalSystem(name=lig_ejm_45_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_45), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_46': ChemicalSystem(name=lig_ejm_46_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_46), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_47': ChemicalSystem(name=lig_ejm_47_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_47), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_48': ChemicalSystem(name=lig_ejm_48_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_48), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_50': ChemicalSystem(name=lig_ejm_50_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_50), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_54': ChemicalSystem(name=lig_ejm_54_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_54), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_ejm_55': ChemicalSystem(name=lig_ejm_55_complex, components={'ligand': SmallMoleculeComponent(name=lig_ejm_55), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_jmc_23': ChemicalSystem(name=lig_jmc_23_complex, components={'ligand': SmallMoleculeComponent(name=lig_jmc_23), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_jmc_27': ChemicalSystem(name=lig_jmc_27_complex, components={'ligand': SmallMoleculeComponent(name=lig_jmc_27), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)}),\n", + " 'lig_jmc_28': ChemicalSystem(name=lig_jmc_28_complex, components={'ligand': SmallMoleculeComponent(name=lig_jmc_28), 'solvent': SolventComponent(name=O, Na+, Cl-), 'protein': ProteinComponent(name=tyk2)})}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "complexed = {l.name: ChemicalSystem(components={'ligand': l, \n", + " 'solvent': solvent, \n", + " 'protein': protein}, \n", + " name=f\"{l.name}_complex\") \n", + " for l in ligands}\n", + "complexed" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "c8a0b60d-bfae-4f99-9f4d-656c982433b9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'lig_ejm_31': ChemicalSystem(name=lig_ejm_31_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_31), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_42': ChemicalSystem(name=lig_ejm_42_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_42), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_43': ChemicalSystem(name=lig_ejm_43_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_43), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_45': ChemicalSystem(name=lig_ejm_45_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_45), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_46': ChemicalSystem(name=lig_ejm_46_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_46), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_47': ChemicalSystem(name=lig_ejm_47_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_47), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_48': ChemicalSystem(name=lig_ejm_48_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_48), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_50': ChemicalSystem(name=lig_ejm_50_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_50), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_54': ChemicalSystem(name=lig_ejm_54_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_54), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_ejm_55': ChemicalSystem(name=lig_ejm_55_water, components={'ligand': SmallMoleculeComponent(name=lig_ejm_55), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_jmc_23': ChemicalSystem(name=lig_jmc_23_water, components={'ligand': SmallMoleculeComponent(name=lig_jmc_23), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_jmc_27': ChemicalSystem(name=lig_jmc_27_water, components={'ligand': SmallMoleculeComponent(name=lig_jmc_27), 'solvent': SolventComponent(name=O, Na+, Cl-)}),\n", + " 'lig_jmc_28': ChemicalSystem(name=lig_jmc_28_water, components={'ligand': SmallMoleculeComponent(name=lig_jmc_28), 'solvent': SolventComponent(name=O, Na+, Cl-)})}" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solvated = {l.name: ChemicalSystem(components={'ligand': l, \n", + " 'solvent': solvent}, \n", + " name=f\"{l.name}_water\") \n", + " for l in ligands}\n", + "solvated" + ] + }, + { + "cell_type": "markdown", + "id": "7566fda4-6c8d-4003-8601-45f05face65a", + "metadata": {}, + "source": [ + "We now have all our network nodes defined. Next, we need to define the `Transformation`s that we wish to perform between them." + ] + }, + { + "cell_type": "markdown", + "id": "7bc3a65a-979c-4114-b948-f7315796709d", + "metadata": {}, + "source": [ + "## Define `Transformation`s between `ChemicalSystem`s as network edges" + ] + }, + { + "cell_type": "markdown", + "id": "0a92b4c1-fffb-410c-9c76-177fe3ea0a87", + "metadata": {}, + "source": [ + "A `Transformation` is a directed edge between two `ChemicalSystem`s. It includes a `Protocol` parameterized with `Settings`, and if optionally a `ComponentMapping`. \n", + "\n", + "The `Protocol` defines the actual computational method used to evaluate the `Transformation` to yield estimates for the free energy difference between the `ChemicalSystem`s.\n", + "\n", + "The `ComponentMapping` defines the atom mapping(s) between corresponding `Component`s in the two `ChemicalSystem`s. This is often critical for relative binding free energy calculations, since the choice of mapping can heavily influence convergence of the resulting estimates." + ] + }, + { + "cell_type": "markdown", + "id": "ace9de53-2527-47f5-8aae-98d9adef2c1d", + "metadata": {}, + "source": [ + "### Define the `Protocol` used for `Transformation` evaluation" + ] + }, + { + "cell_type": "markdown", + "id": "0338be2a-b861-4534-bc78-49c2e00ca2ac", + "metadata": {}, + "source": [ + "For this example, we'll use the same `Protocol` for all our `Transformation`s, with identical `Settings` for each." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "183b23a6-1839-4586-bbbe-60b2defd2f01", + "metadata": {}, + "outputs": [], + "source": [ + "from perses.protocols.nonequilibrium_cycling import NonEquilibriumCyclingProtocol" + ] + }, + { + "cell_type": "markdown", + "id": "9e686c7c-9ef7-4e97-b384-63334eeb894e", + "metadata": {}, + "source": [ + "Any given `Protocol` features a `default_settings` method, which can be used to get the default settings that are specific to that `Protocol`:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "d7bac455-ddc1-4735-96d7-5c6681ab3cbc", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "{'num_replicates': 1,\n", + " 'lambda_functions': {'lambda_sterics_core': 'lambda',\n", + " 'lambda_electrostatics_core': 'lambda',\n", + " 'lambda_sterics_insert': 'select(step(lambda - 0.5), 1.0, 2.0 * lambda)',\n", + " 'lambda_sterics_delete': 'select(step(lambda - 0.5), 2.0 * (lambda - 0.5), 0.0)',\n", + " 'lambda_electrostatics_insert': 'select(step(lambda - 0.5), 2.0 * (lambda - 0.5), 0.0)',\n", + " 'lambda_electrostatics_delete': 'select(step(lambda - 0.5), 1.0, 2.0 * lambda)',\n", + " 'lambda_bonds': 'lambda',\n", + " 'lambda_angles': 'lambda',\n", + " 'lambda_torsions': 'lambda'},\n", + " 'softcore_LJ_v2': True,\n", + " 'interpolate_old_and_new_14s': False,\n", + " 'timestep': 4.0 ,\n", + " 'neq_splitting': 'V R H O R V',\n", + " 'eq_steps': 1000,\n", + " 'neq_steps': 100,\n", + " 'platform': 'CUDA',\n", + " 'save_frequency': 100}" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "protocol_settings = NonEquilibriumCyclingProtocol.default_settings()\n", + "protocol_settings.dict()" + ] + }, + { + "cell_type": "markdown", + "id": "d7898feb-419a-45cc-8837-53cd583d130b", + "metadata": {}, + "source": [ + "These can be edited, e.g. with:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "998ef869-7086-4c24-9c53-2d4e9553b6b9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "protocol_settings.save_frequency = 200" + ] + }, + { + "cell_type": "markdown", + "id": "8af105c1-08fe-4f49-a422-a6d6307695d5", + "metadata": {}, + "source": [ + "We'll construct our full `Settings` for our chosen `NonEquilibriumCyclingProtocol`, which will include the more general `ThermoSettings` and `ForcefieldSettings` as well:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "19abd504-1e44-413e-8740-42413163a25a", + "metadata": {}, + "outputs": [], + "source": [ + "from openff.units import unit\n", + "from gufe.settings.models import (\n", + " Settings, \n", + " ThermoSettings, \n", + " OpenMMSystemGeneratorFFSettings,\n", + ")\n", + "from perses.protocols.settings import NonEqCyclingSettings\n", + "\n", + "settings = Settings(\n", + " settings_version=0,\n", + " forcefield_settings=OpenMMSystemGeneratorFFSettings(),\n", + " thermo_settings=ThermoSettings(temperature=300*unit.kelvin),\n", + " protocol_settings=protocol_settings,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "70e558b5-93e4-4c91-8ff0-bc5b99646a71", + "metadata": {}, + "source": [ + "We can now produce a parameterized `NonEquilibriumCyclingProtocol` instance:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "bd1d0a0d-00c9-44a1-820c-186b4ee05334", + "metadata": {}, + "outputs": [], + "source": [ + "protocol = NonEquilibriumCyclingProtocol(settings)" + ] + }, + { + "cell_type": "markdown", + "id": "57239762-43b9-4501-b3ad-9eebbcd64a20", + "metadata": {}, + "source": [ + "### Build the `Transformation`s" + ] + }, + { + "cell_type": "markdown", + "id": "8ccdec56-519e-4293-86c4-60b0255cfcad", + "metadata": {}, + "source": [ + "We can now construct the `Transformation`s we want represented in our network. We'll use the predefined connections from the `tyk2` system from above as the basis for our choices here, but you could use any network planner of your choice to generate connections and use those instead." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "abdcc8cb-6460-42f5-9a3b-f360041fd1e4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[('lig_ejm_31', 'lig_ejm_50'),\n", + " ('lig_ejm_46', 'lig_jmc_23'),\n", + " ('lig_ejm_31', 'lig_ejm_55'),\n", + " ('lig_ejm_31', 'lig_ejm_48'),\n", + " ('lig_ejm_31', 'lig_ejm_54'),\n", + " ('lig_ejm_31', 'lig_ejm_47'),\n", + " ('lig_ejm_31', 'lig_ejm_46'),\n", + " ('lig_ejm_46', 'lig_jmc_27'),\n", + " ('lig_ejm_46', 'lig_jmc_28'),\n", + " ('lig_ejm_42', 'lig_ejm_43'),\n", + " ('lig_ejm_31', 'lig_ejm_42'),\n", + " ('lig_ejm_45', 'lig_ejm_55')]" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tyk2_system.connections" + ] + }, + { + "cell_type": "markdown", + "id": "b1b65f19-15bb-49df-88e8-ddf8adccf8e1", + "metadata": {}, + "source": [ + "**TODO: need to add mappings for each edge; these would be included in the `Transformation` creations below.**" + ] + }, + { + "cell_type": "markdown", + "id": "cfd3fe16-1666-4b78-8263-becff81dfef1", + "metadata": {}, + "source": [ + "Since we are planning to perform relative binding free energy (RBFE) calculations, we'll define both *complex* and *solvent* variants for each `Transformation`:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "9c1b4bf5-8fcb-4ce8-8eba-12ad3c4c93f8", + "metadata": {}, + "outputs": [], + "source": [ + "complexed_transformations = [Transformation(stateA=complexed[edge[0]], \n", + " stateB=complexed[edge[1]], \n", + " protocol=protocol) \n", + " for edge in tyk2_system.connections]" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "2e278a01-432f-4fb5-a004-a76eff100058", + "metadata": {}, + "outputs": [], + "source": [ + "solvated_transformations = [Transformation(stateA=solvated[edge[0]], \n", + " stateB=solvated[edge[1]], \n", + " protocol=protocol) \n", + " for edge in tyk2_system.connections]" + ] + }, + { + "cell_type": "markdown", + "id": "36c79a08-d4aa-4388-a0af-be51906f7672", + "metadata": {}, + "source": [ + "## Create the `AlchemicalNetwork`" + ] + }, + { + "cell_type": "markdown", + "id": "6ebc1174-b348-4e44-97f7-444c44c8d6d0", + "metadata": {}, + "source": [ + "An `AlchemicalNetwork` is simply the combination of `ChemicalSystem`s (nodes) and `Transformation`s (directed edges) that we want to evaluate $\\Delta G$s for. This data structure functions as a declaration of what you want to compute, and is the central object on which systems like `fah-alchemy` operate. \n", + "\n", + "We'll finish here by creating an `AlchemicalNetwork` from the collection of objects we've built so far." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "b45baada-efa4-46f7-a2cb-6e24a9309f1e", + "metadata": {}, + "outputs": [ + { + "ename": "TypeError", + "evalue": "unhashable type: 'Settings'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[24], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m network \u001b[38;5;241m=\u001b[39m \u001b[43mAlchemicalNetwork\u001b[49m\u001b[43m(\u001b[49m\u001b[43medges\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43msolvated_transformations\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mcomplexed_transformations\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 2\u001b[0m \u001b[43m \u001b[49m\u001b[43mnodes\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43msolvated\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvalues\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mcomplexed\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvalues\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 3\u001b[0m \u001b[43m \u001b[49m\u001b[43mname\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mtyk2_relative_benchmark\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 4\u001b[0m network\n", + "File \u001b[0;32m~/.conda/envs/openfe-notebooks/lib/python3.9/site-packages/gufe/tokenization.py:55\u001b[0m, in \u001b[0;36m_GufeTokenizableMeta.__call__\u001b[0;34m(cls, *args, **kwargs)\u001b[0m\n\u001b[1;32m 54\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__call__\u001b[39m(\u001b[38;5;28mcls\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m---> 55\u001b[0m instance \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__call__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 56\u001b[0m \u001b[38;5;66;03m# add to registry if not already present\u001b[39;00m\n\u001b[1;32m 57\u001b[0m TOKENIZABLE_REGISTRY\u001b[38;5;241m.\u001b[39msetdefault(instance\u001b[38;5;241m.\u001b[39mkey, instance)\n", + "File \u001b[0;32m~/.conda/envs/openfe-notebooks/lib/python3.9/site-packages/gufe/network.py:31\u001b[0m, in \u001b[0;36mAlchemicalNetwork.__init__\u001b[0;34m(self, edges, nodes, name)\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\n\u001b[1;32m 26\u001b[0m \u001b[38;5;28mself\u001b[39m,\n\u001b[1;32m 27\u001b[0m edges: Optional[Iterable[Transformation]] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m 28\u001b[0m nodes: Optional[Iterable[ChemicalSystem]] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m 29\u001b[0m name: Optional[\u001b[38;5;28mstr\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m 30\u001b[0m ):\n\u001b[0;32m---> 31\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_edges: FrozenSet[Transformation] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mfrozenset\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43medges\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mif\u001b[39;00m edges \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;28mfrozenset\u001b[39m()\n\u001b[1;32m 32\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_nodes: FrozenSet[ChemicalSystem]\n\u001b[1;32m 34\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_name \u001b[38;5;241m=\u001b[39m name\n", + "File \u001b[0;32m~/.conda/envs/openfe-notebooks/lib/python3.9/site-packages/gufe/transformations/transformation.py:111\u001b[0m, in \u001b[0;36mTransformation.__hash__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 110\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__hash__\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m--> 111\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mhash\u001b[39;49m\u001b[43m(\u001b[49m\n\u001b[1;32m 112\u001b[0m \u001b[43m \u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 113\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_stateA\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 114\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_stateB\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 115\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_mapping\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 116\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_name\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 117\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_protocol\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 118\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 119\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/.conda/envs/openfe-notebooks/lib/python3.9/site-packages/gufe/protocols/protocol.py:122\u001b[0m, in \u001b[0;36mProtocol.__hash__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 121\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__hash__\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m--> 122\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mhash\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;18;43m__class__\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;18;43m__name__\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_settings\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[0;31mTypeError\u001b[0m: unhashable type: 'Settings'" + ] + } + ], + "source": [ + "network = AlchemicalNetwork(edges=(solvated_transformations + complexed_transformations), \n", + " nodes=(list(solvated.values()) + list(complexed.values())),\n", + " name=\"tyk2_relative_benchmark\")\n", + "network" + ] + }, + { + "cell_type": "markdown", + "id": "262c5908-640d-45dc-9ca4-6d912428b8ac", + "metadata": {}, + "source": [ + "That's it! We simply toss in all `Transformation`s (edges) and `ChemicalSystem`s (nodes) we want included in this `AlchemicalNetwork`, and optionally give it a name that means something to us (it need not be unique, but can be used to query for network(s) from `fah-alchemy` later)." + ] + }, + { + "cell_type": "markdown", + "id": "43849b09-5031-4eb0-989c-ceca2110f736", + "metadata": {}, + "source": [ + "We could have chosen here to leave the `nodes` argument off, since every `ChemicalSystem` we included was already represented among the `edges`, but we show it here for completeness. In this way, it's possible to include `ChemicalSystem`s in the network that aren't connected via any `Transformation`s to others, though in practice there isn't much utility in this." + ] + }, + { + "cell_type": "markdown", + "id": "4e4bfff8-0072-4ad7-9b5e-3ac3f6bdc551", + "metadata": {}, + "source": [ + "### Optional: Run a `Protocol` locally" + ] + }, + { + "cell_type": "markdown", + "id": "954ce90a-2c7f-4ffc-a218-e5c8298d9ce5", + "metadata": {}, + "source": [ + "We can run our parameterized `NonEqulibriumCyclingProtocol` locally as a way to check if things are working as we expect. We'll pick one of our `Transformation`s out from our `AlchemicalNetwork`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "751b25e3-4d1e-4d5e-a870-12cc12021f73", + "metadata": {}, + "outputs": [], + "source": [ + "transformation = list(network.edges)[0]" + ] + }, + { + "cell_type": "markdown", + "id": "7b068982-1e5d-478d-9adc-6f86059ebbbd", + "metadata": {}, + "source": [ + "We'll generate a `ProtocolDAG` that encodes the actual operations to perform in order to execute the `Protocol`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdcc8b9f-65aa-4f8b-8661-1fe4eeb997b3", + "metadata": {}, + "outputs": [], + "source": [ + "protocoldag = transformation.create()" + ] + }, + { + "cell_type": "markdown", + "id": "1fea6b7d-0ba8-4804-be0c-6e3bf935bc78", + "metadata": {}, + "source": [ + "And we'll run it locally, in-process. This will run each `ProtocolUnit` in the `ProtocolDAG` in series, in dependency order:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "271b0164-9e0e-4578-b637-81124903e455", + "metadata": {}, + "outputs": [], + "source": [ + "from gufe.protocols.protocoldag import execute_DAG" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58c4e234-a112-4ba7-9eec-e55b64e39ba1", + "metadata": {}, + "outputs": [], + "source": [ + "protocoldagresult = execute_DAG(protocoldag)" + ] + }, + { + "cell_type": "markdown", + "id": "f6e222af-ae1c-4437-aaac-0fe0f5d9fb16", + "metadata": {}, + "source": [ + "The above will raise an exception if at any point execution failed." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.16" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}