Demo - Using alchemiscale to evaluate a relative binding free energy network¶
This notebook details the process of running a relative binding free energy calculation using the openfe toolkit and the execution platform alchemiscale.
The notebook is broken down into the following section:
Creating an
AlchemicalNetworkSubmitting an
AlchemicalnetworktoalchemiscaleMonitoring the progress of an
alchemiscalesubmissionGathering results and analyzing them using
cinnabar
Credits:
Irfan Alibay (@ialibay) for initial structure and content
David Dotson (@dotsdl) for expansion and editing for use as a tutorial
Creating an AlchemicalNetwork for the tyk2 benchmark system¶
We first illustrate how to build AlchemicalNetwork objects,
from a starting point of chemical models stored in sdf and pdb files.
An AlchemicalNetwork is used to represent an entire network of calculations,
and is composed of many smaller objects:
An
AlchemicalNetworkcomposed ofeach node a
ChemicalSystemeach containing many components, such as
SmallMoleculeComponent,ProteinComponentinternally each Component usually wraps an RDKit representation
each directed edge a
Transformation, containingtwo
ChemicalSystems, the ‘A’ and ‘B’ sidezero or more
Mappingobjects relating these two sidesa
Protocoldefining the computational method to be applied to other items
import openfe
import gufe
from openff.units import unit
from rdkit import Chem
print("gufe version: ", gufe.__version__)
print("openfe version: ", openfe.__version__)
gufe version: 1.0.0
openfe version: 1.0.1+0.g48dcbb26.dirty
import logging
logger = logging.getLogger()
logger.setLevel(logging.WARN)
Define ChemicalSystems for network nodes¶
We’ll start by defining the nodes for our network.
A ChemicalSystem is made of one or more Components. 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.
Reading Ligands¶
The ligands are concatenated in a single sdf file, we’ll read these using RDKit.
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.
ligands = [
openfe.SmallMoleculeComponent(m) for m in Chem.SDMolSupplier('ligands.sdf', removeHs=False)
]
ligands
[SmallMoleculeComponent(name=lig_ejm_54),
SmallMoleculeComponent(name=lig_jmc_23),
SmallMoleculeComponent(name=lig_ejm_47),
SmallMoleculeComponent(name=lig_jmc_27),
SmallMoleculeComponent(name=lig_ejm_46),
SmallMoleculeComponent(name=lig_ejm_31),
SmallMoleculeComponent(name=lig_ejm_42),
SmallMoleculeComponent(name=lig_ejm_50),
SmallMoleculeComponent(name=lig_ejm_45),
SmallMoleculeComponent(name=lig_jmc_28),
SmallMoleculeComponent(name=lig_ejm_55),
SmallMoleculeComponent(name=lig_ejm_43),
SmallMoleculeComponent(name=lig_ejm_48)]
# The OpenFE RBFE protocol does not yet support charge changes (to be changed very soon!)
# Here we raise an error if a change in formal charge is identified
base_formalcharge = Chem.rdmolops.GetFormalCharge(ligands[0].to_rdkit())
for mol in ligands[1:]:
if Chem.rdmolops.GetFormalCharge(mol.to_rdkit()) != base_formalcharge:
errmsg = f"The molecule {mol} has a different formal charge than the first ligand molecule"
raise ValueError(errmsg)
Reading the protein¶
The protein is supplied as a PDB file, readable via the ProteinComponent.from_pdb_file class method.
protein = openfe.ProteinComponent.from_pdb_file('protein.pdb', name='tyk2')
protein
ProteinComponent(name=tyk2)
Defining the solvent¶
We’ll also need at least one SolventComponent to encode our choice of solvent and counterions, with concentration.
The concentration is defined as having units supplied by openff.units, this package is used to avoid confusion.
The SolventComponent doesn’t actually perform any actual solvation (packing water molecules, ions); that is performed just before simulation time during Protocol execution.
solvent = openfe.SolventComponent(positive_ion='Na',
negative_ion='Cl',
neutralize=True,
ion_concentration=0.15*unit.molar)
solvent
SolventComponent(name=O, Na+, Cl-)
Build the ChemicalSystems¶
We can now construct the ChemicalSystems 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.
This produces a dictionary mapping the ligand name to the ChemicalSystem that contains that ligand.
There are two dictionaries, for complexed and solvated ligands respectively.
complexed = {l.name: openfe.ChemicalSystem(components={'ligand': l,
'solvent': solvent,
'protein': protein},
name=f"{l.name}_complex")
for l in ligands}
complexed
{'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)}),
'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)})}
solvated = {l.name: openfe.ChemicalSystem(components={'ligand': l,
'solvent': solvent},
name=f"{l.name}_solvent")
for l in ligands}
solvated
{'lig_ejm_54': ChemicalSystem(name=lig_ejm_54_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_54), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_jmc_23': ChemicalSystem(name=lig_jmc_23_solvent, components={'ligand': SmallMoleculeComponent(name=lig_jmc_23), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_47': ChemicalSystem(name=lig_ejm_47_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_47), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_jmc_27': ChemicalSystem(name=lig_jmc_27_solvent, components={'ligand': SmallMoleculeComponent(name=lig_jmc_27), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_46': ChemicalSystem(name=lig_ejm_46_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_46), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_31': ChemicalSystem(name=lig_ejm_31_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_31), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_42': ChemicalSystem(name=lig_ejm_42_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_42), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_50': ChemicalSystem(name=lig_ejm_50_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_50), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_45': ChemicalSystem(name=lig_ejm_45_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_45), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_jmc_28': ChemicalSystem(name=lig_jmc_28_solvent, components={'ligand': SmallMoleculeComponent(name=lig_jmc_28), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_55': ChemicalSystem(name=lig_ejm_55_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_55), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_43': ChemicalSystem(name=lig_ejm_43_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_43), 'solvent': SolventComponent(name=O, Na+, Cl-)}),
'lig_ejm_48': ChemicalSystem(name=lig_ejm_48_solvent, components={'ligand': SmallMoleculeComponent(name=lig_ejm_48), 'solvent': SolventComponent(name=O, Na+, Cl-)})}
Define Transformations between ChemicalSystems as network edges¶
A Transformation is a directed edge between two ChemicalSystems. It includes a Protocol parameterized with Settings, and optionally a ComponentMapping.
The Protocol defines the actual computational method used to evaluate the Transformation to yield estimates for the free energy difference between the ChemicalSystems.
The ComponentMapping defines the atom mapping(s) between corresponding Components in the two ChemicalSystems. This is often critical for relative binding free energy calculations, since the choice of mapping can heavily influence convergence of the resulting estimates.
Define the Protocol used for Transformation evaluation¶
For this example, we’ll use the same Protocol for all our Transformations, with identical Settings for each.
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol
Any given Protocol has a default_settings() method, which can be used to get the default settings that are specific to that Protocol:
protocol_settings = RelativeHybridTopologyProtocol.default_settings()
protocol_settings.dict()
{'forcefield_settings': {'constraints': 'hbonds',
'rigid_water': True,
'hydrogen_mass': 3.0,
'forcefields': ['amber/ff14SB.xml',
'amber/tip3p_standard.xml',
'amber/tip3p_HFE_multivalent.xml',
'amber/phosaa10.xml'],
'small_molecule_forcefield': 'openff-2.1.1',
'nonbonded_cutoff': 1.0 <Unit('nanometer')>,
'nonbonded_method': 'PME'},
'thermo_settings': {'temperature': 298.15 <Unit('kelvin')>,
'pressure': 0.9869232667160129 <Unit('standard_atmosphere')>,
'ph': None,
'redox_potential': None},
'protocol_repeats': 3,
'solvation_settings': {'solvent_model': 'tip3p',
'solvent_padding': 1.2 <Unit('nanometer')>},
'partial_charge_settings': {'partial_charge_method': 'am1bcc',
'off_toolkit_backend': 'ambertools',
'number_of_conformers': None,
'nagl_model': None},
'lambda_settings': {'lambda_functions': 'default', 'lambda_windows': 11},
'alchemical_settings': {'softcore_LJ': 'gapsys',
'explicit_charge_correction_cutoff': 0.8 <Unit('nanometer')>,
'endstate_dispersion_correction': False,
'use_dispersion_correction': False,
'softcore_alpha': 0.85,
'turn_off_core_unique_exceptions': False,
'explicit_charge_correction': False},
'simulation_settings': {'equilibration_length': 1.0 <Unit('nanosecond')>,
'production_length': 5.0 <Unit('nanosecond')>,
'minimization_steps': 5000,
'time_per_iteration': 1 <Unit('picosecond')>,
'real_time_analysis_interval': 250 <Unit('picosecond')>,
'early_termination_target_error': 0.0 <Unit('kilocalorie_per_mole')>,
'real_time_analysis_minimum_time': 500 <Unit('picosecond')>,
'sampler_method': 'repex',
'sams_flatness_criteria': 'logZ-flatness',
'sams_gamma0': 1.0,
'n_replicas': 11},
'engine_settings': {'compute_platform': None},
'integrator_settings': {'timestep': 4 <Unit('femtosecond')>,
'langevin_collision_rate': 1.0 <Unit('1 / picosecond')>,
'barostat_frequency': 25 <Unit('timestep')>,
'remove_com': False,
'reassign_velocities': False,
'n_restart_attempts': 20,
'constraint_tolerance': 1e-06},
'output_settings': {'checkpoint_interval': 250 <Unit('picosecond')>,
'forcefield_cache': 'db.json',
'output_indices': 'not water',
'checkpoint_storage_filename': 'checkpoint.chk',
'output_filename': 'simulation.nc',
'output_structure': 'hybrid_system.pdb'}}
protocol_settings
{'alchemical_settings': {'endstate_dispersion_correction': False,
'explicit_charge_correction': False,
'explicit_charge_correction_cutoff': <Quantity(0.8, 'nanometer')>,
'softcore_LJ': 'gapsys',
'softcore_alpha': 0.85,
'turn_off_core_unique_exceptions': False,
'use_dispersion_correction': False},
'engine_settings': {'compute_platform': None},
'forcefield_settings': {'constraints': 'hbonds',
'forcefields': ['amber/ff14SB.xml',
'amber/tip3p_standard.xml',
'amber/tip3p_HFE_multivalent.xml',
'amber/phosaa10.xml'],
'hydrogen_mass': 3.0,
'nonbonded_cutoff': <Quantity(1.0, 'nanometer')>,
'nonbonded_method': 'PME',
'rigid_water': True,
'small_molecule_forcefield': 'openff-2.1.1'},
'integrator_settings': {'barostat_frequency': <Quantity(25, 'timestep')>,
'constraint_tolerance': 1e-06,
'langevin_collision_rate': <Quantity(1.0, '1 / picosecond')>,
'n_restart_attempts': 20,
'reassign_velocities': False,
'remove_com': False,
'timestep': <Quantity(4, 'femtosecond')>},
'lambda_settings': {'lambda_functions': 'default', 'lambda_windows': 11},
'output_settings': {'checkpoint_interval': <Quantity(250, 'picosecond')>,
'checkpoint_storage_filename': 'checkpoint.chk',
'forcefield_cache': 'db.json',
'output_filename': 'simulation.nc',
'output_indices': 'not water',
'output_structure': 'hybrid_system.pdb'},
'partial_charge_settings': {'nagl_model': None,
'number_of_conformers': None,
'off_toolkit_backend': 'ambertools',
'partial_charge_method': 'am1bcc'},
'protocol_repeats': 3,
'simulation_settings': {'early_termination_target_error': <Quantity(0.0, 'kilocalorie_per_mole')>,
'equilibration_length': <Quantity(1.0, 'nanosecond')>,
'minimization_steps': 5000,
'n_replicas': 11,
'production_length': <Quantity(5.0, 'nanosecond')>,
'real_time_analysis_interval': <Quantity(250, 'picosecond')>,
'real_time_analysis_minimum_time': <Quantity(500, 'picosecond')>,
'sampler_method': 'repex',
'sams_flatness_criteria': 'logZ-flatness',
'sams_gamma0': 1.0,
'time_per_iteration': <Quantity(1, 'picosecond')>},
'solvation_settings': {'solvent_model': 'tip3p',
'solvent_padding': <Quantity(1.2, 'nanometer')>},
'thermo_settings': {'ph': None,
'pressure': <Quantity(0.986923267, 'standard_atmosphere')>,
'redox_potential': None,
'temperature': <Quantity(298.15, 'kelvin')>}}
# Here we make it so that each simulation only features a single repeat
# We then do multiple repeats by running each simulation multiple times
protocol_settings.protocol_repeats = 1
We can now produce a parameterized RelativeHybridTopologyProtocol instance:
protocol = RelativeHybridTopologyProtocol(protocol_settings)
Build the Transformations¶
We can now construct the Transformations we want represented in our network.
We’ll use openfe’s implementation of the LomapAtomMapper to produce mappings between ligands, the default LOMAP scorer for assessing the quality of the mappings produced, and the minimal spanning tree network planner. However, you could use any network planner of your choice to generate the network connections and mappings and use those instead.
mapper = openfe.LomapAtomMapper(max3d=1.0, element_change=True)
scorer = openfe.lomap_scorers.default_lomap_score
network_planner = openfe.ligand_network_planning.generate_minimal_spanning_network
# Note creating a minimum spanning tree requires an all-to-all mapping generation
# This can be quite time consuming...
ligand_network = network_planner(
ligands=ligands,
mappers=[mapper],
scorer=scorer
)
from openfe.utils.atommapping_network_plotting import plot_atommapping_network
plot_atommapping_network(ligand_network)
Since we are planning to perform relative binding free energy (RBFE) calculations, we’ll define both complex and solvent variants for each Transformation:
for mapping in ligand_network.edges:
print(mapping.componentA.name, mapping.componentB.name)
lig_ejm_46 lig_ejm_31
lig_ejm_55 lig_ejm_43
lig_ejm_31 lig_ejm_50
lig_ejm_31 lig_ejm_48
lig_ejm_47 lig_ejm_31
lig_ejm_42 lig_ejm_50
lig_jmc_23 lig_ejm_46
lig_ejm_42 lig_ejm_43
lig_jmc_23 lig_jmc_28
lig_jmc_27 lig_jmc_28
lig_ejm_54 lig_ejm_55
lig_ejm_31 lig_ejm_45
complexed_transformations = []
solvated_transformations = []
for mapping in ligand_network.edges:
ligA_name = mapping.componentA.name
ligB_name = mapping.componentB.name
ligA = complexed[ligA_name]['ligand']
ligB = complexed[ligB_name]['ligand']
complexed_transformations.append(
openfe.Transformation(stateA=complexed[ligA_name],
stateB=complexed[ligB_name],
mapping={'ligand': mapping},
protocol=protocol,
name=f"{complexed[ligA_name].name}->{complexed[ligB_name].name}")
)
solvated_transformations.append(
openfe.Transformation(stateA=solvated[ligA_name],
stateB=solvated[ligB_name],
mapping={'ligand': mapping},
protocol=protocol,
name=f"{solvated[ligA_name].name}->{solvated[ligB_name].name}")
)
Create the AlchemicalNetwork¶
An AlchemicalNetwork is simply the combination of ChemicalSystems (nodes) and Transformations (directed edges) that we want to evaluate \(\Delta G\)s for. This data structure functions as a declaration of what you want to compute.
We’ll finish here by creating an AlchemicalNetwork from the collection of objects we’ve built so far.
network = openfe.AlchemicalNetwork(edges=(solvated_transformations + complexed_transformations),
nodes=(list(solvated.values()) + list(complexed.values())),
name="tyk2_relative_benchmark")
network
<AlchemicalNetwork-901e26f212324c47939436bc9b64a964>
That’s it! We simply toss in all Transformations (edges) and ChemicalSystems (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) later).
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 ChemicalSystems in the network that aren’t connected via any Transformations to others, though in practice there isn’t much utility in this.
Submitting an AlchemicalNetwork to alchemiscale¶
import alchemiscale
print("alchemiscale version: ", alchemiscale.__version__)
alchemiscale version: 0.5.0
from alchemiscale import AlchemiscaleClient, Scope, ScopedKey
from getpass import getpass
identity_key = getpass()
asc = AlchemiscaleClient('https://api.alchemiscale.org', 'ddotson', identity_key)
asc
AlchemiscaleClient('https://api.alchemiscale.org')
asc.list_scopes()
[<Scope('*-*-*')>]
asc.query_networks(scope=Scope('ddotson'))
[]
AlchemicalNetworks and the objects associated with them live in a specific Scope. Because my user identity is able to use any Scope from the above ('*-*-*'), I can choose whatever I want. AlchemicalNetworks submitted to the same Scope will share Transformations they have in common, along with results for those Transformations; AlchemicalNetworks in different Scopes will share nothing.
How you choose to use the Scopes you have access to is entirely up to you. I’ll choose this is my Scope for this network:
scope = Scope('ddotson', 'tyk2', 'demo')
And then we’ll submit the AlchemicalNetwork to that Scope with:
an_sk = asc.create_network(network, scope)
We get back a ScopedKey for the `AlchemicalNetwork:
an_sk
<ScopedKey('AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo')>
This uniquely identifies the AlchemicalNetwork in the database, and allows us to perform operations later on it. We can get the ScopedKeys for all AlchemicalNetworks that exist in Scopes accessible to us with:
asc.query_networks(scope=Scope('ddotson'))
[<ScopedKey('AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo')>]
We can likewise retrieve an AlchemicalNetwork that we previously submitted with:
network_again = asc.get_network(an_sk)
Executing Transformations by creating and actioning Tasks¶
Submitting an AlchemicalNetwork with create_network only defines it on the server; it does not initiate any calculations. For this, we need to create and action Tasks.
To obtain a free energy estimate \(\Delta G\) for each Transformation in the AlchemicalNetwork, we’ll create a Task for each one:
tf_sks = asc.get_network_transformations(an_sk)
tasks = asc.create_transformations_tasks(tf_sks)
tasks
[<ScopedKey('Task-ef9ccf8712d6458582e075a5be87d06e-ddotson-tyk2-demo')>,
<ScopedKey('Task-f6a91ccc6f9e4bcbace3331a0109c9da-ddotson-tyk2-demo')>,
<ScopedKey('Task-cccc389fc316494994cb6c8f834aa006-ddotson-tyk2-demo')>,
<ScopedKey('Task-6205c49da54647a491ae11bce9105288-ddotson-tyk2-demo')>,
<ScopedKey('Task-7003ead7ec9a4ff18965067add5e80fb-ddotson-tyk2-demo')>,
<ScopedKey('Task-de763d3ed6f045f596e96373ae231cfd-ddotson-tyk2-demo')>,
<ScopedKey('Task-a2505069149d47738f2b30a47f24d860-ddotson-tyk2-demo')>,
<ScopedKey('Task-eee82adb878f4701bd61ce765abc5da6-ddotson-tyk2-demo')>,
<ScopedKey('Task-4e9d22efa6364a96854fd2e45dfe1939-ddotson-tyk2-demo')>,
<ScopedKey('Task-3fade5ba936a4d1fb967c695137aa210-ddotson-tyk2-demo')>,
<ScopedKey('Task-320fc3cb408a4e5d9cea7304cf73d7f1-ddotson-tyk2-demo')>,
<ScopedKey('Task-191cb8e4a03c4595a249cf1b7b8771db-ddotson-tyk2-demo')>,
<ScopedKey('Task-112dcba14c0342fba611c9079c0c6004-ddotson-tyk2-demo')>,
<ScopedKey('Task-1200e42c2a8a4d43b1ae39528cb48e8d-ddotson-tyk2-demo')>,
<ScopedKey('Task-35577a6ccebc44f58a129f0ad061d7b7-ddotson-tyk2-demo')>,
<ScopedKey('Task-859860913fa54db68d28ed2e2ec7e494-ddotson-tyk2-demo')>,
<ScopedKey('Task-0d67a73e932f47a094e9c1111a0833bc-ddotson-tyk2-demo')>,
<ScopedKey('Task-5e85cac6182b48258d572d037b52bc71-ddotson-tyk2-demo')>,
<ScopedKey('Task-6ec4facb5ba54ebaa9905bacb0e39a5c-ddotson-tyk2-demo')>,
<ScopedKey('Task-91e70ed59c084d3388df9009ab1553e2-ddotson-tyk2-demo')>,
<ScopedKey('Task-21b0cba98a594da4a15f02c5ab20ed62-ddotson-tyk2-demo')>,
<ScopedKey('Task-ce83f523df324508a6e6cfbce1c34393-ddotson-tyk2-demo')>,
<ScopedKey('Task-acbad8588c5b4f0c962f347785277c94-ddotson-tyk2-demo')>,
<ScopedKey('Task-e8d59e999687420e94a9afc201a14861-ddotson-tyk2-demo')>]
Now checking the AlchemicalNetwork status, we see 24 Tasks 'waiting', 1 for each Transformation:
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 0 │ │ running │ 0 │ │ waiting │ 24 │ │ error │ 0 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'waiting': 24}
But there is one more step: Tasks must also be actioned on the AlchemicalNetwork in order to be picked up by compute services:
asc.action_tasks(tasks, an_sk)
[<ScopedKey('Task-ef9ccf8712d6458582e075a5be87d06e-ddotson-tyk2-demo')>,
<ScopedKey('Task-f6a91ccc6f9e4bcbace3331a0109c9da-ddotson-tyk2-demo')>,
<ScopedKey('Task-cccc389fc316494994cb6c8f834aa006-ddotson-tyk2-demo')>,
<ScopedKey('Task-6205c49da54647a491ae11bce9105288-ddotson-tyk2-demo')>,
<ScopedKey('Task-7003ead7ec9a4ff18965067add5e80fb-ddotson-tyk2-demo')>,
<ScopedKey('Task-de763d3ed6f045f596e96373ae231cfd-ddotson-tyk2-demo')>,
<ScopedKey('Task-a2505069149d47738f2b30a47f24d860-ddotson-tyk2-demo')>,
<ScopedKey('Task-eee82adb878f4701bd61ce765abc5da6-ddotson-tyk2-demo')>,
<ScopedKey('Task-4e9d22efa6364a96854fd2e45dfe1939-ddotson-tyk2-demo')>,
<ScopedKey('Task-3fade5ba936a4d1fb967c695137aa210-ddotson-tyk2-demo')>,
<ScopedKey('Task-320fc3cb408a4e5d9cea7304cf73d7f1-ddotson-tyk2-demo')>,
<ScopedKey('Task-191cb8e4a03c4595a249cf1b7b8771db-ddotson-tyk2-demo')>,
<ScopedKey('Task-112dcba14c0342fba611c9079c0c6004-ddotson-tyk2-demo')>,
<ScopedKey('Task-1200e42c2a8a4d43b1ae39528cb48e8d-ddotson-tyk2-demo')>,
<ScopedKey('Task-35577a6ccebc44f58a129f0ad061d7b7-ddotson-tyk2-demo')>,
<ScopedKey('Task-859860913fa54db68d28ed2e2ec7e494-ddotson-tyk2-demo')>,
<ScopedKey('Task-0d67a73e932f47a094e9c1111a0833bc-ddotson-tyk2-demo')>,
<ScopedKey('Task-5e85cac6182b48258d572d037b52bc71-ddotson-tyk2-demo')>,
<ScopedKey('Task-6ec4facb5ba54ebaa9905bacb0e39a5c-ddotson-tyk2-demo')>,
<ScopedKey('Task-91e70ed59c084d3388df9009ab1553e2-ddotson-tyk2-demo')>,
<ScopedKey('Task-21b0cba98a594da4a15f02c5ab20ed62-ddotson-tyk2-demo')>,
<ScopedKey('Task-ce83f523df324508a6e6cfbce1c34393-ddotson-tyk2-demo')>,
<ScopedKey('Task-acbad8588c5b4f0c962f347785277c94-ddotson-tyk2-demo')>,
<ScopedKey('Task-e8d59e999687420e94a9afc201a14861-ddotson-tyk2-demo')>]
The reason for this extra step is that a given Transformation may be a member of multiple AlchemicalNetworks in its Scope. Computational effort is independently allocated on each AlchemicalNetwork through its actioned Tasks, and these can differ between AlchemicalNetworks in the same Scope. In the future, automated Strategys will operate on individual AlchemicalNetworks, and this ability to independently choose where to apply effort is critical for this.
After actioning our Tasks, we can see that some have been picked up and are now 'running':
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 0 │ │ running │ 0 │ │ waiting │ 24 │ │ error │ 0 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'waiting': 24}
Later on, some are now 'complete':
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 4 │ │ running │ 12 │ │ waiting │ 8 │ │ error │ 0 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'waiting': 8, 'complete': 4, 'running': 12}
…and some might have errored:
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 14 │ │ running │ 9 │ │ waiting │ 0 │ │ error │ 1 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'error': 1, 'running': 9, 'complete': 14}
…and at some point, none are waiting or running; they are either complete or errored:
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 21 │ │ running │ 0 │ │ waiting │ 0 │ │ error │ 3 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'error': 3, 'complete': 21}
Dealing with errors¶
Inevitably, some of your Tasks will encounter problems in execution, either random or systematic errors. When this happens, the Task status will be set to 'error'.
We’ll get all the Tasks in this AlchemicalNetwork with status 'error':
failed_tasks = asc.get_network_tasks(an_sk, status='error')
failed_tasks
[<ScopedKey('Task-ef9ccf8712d6458582e075a5be87d06e-ddotson-tyk2-demo')>,
<ScopedKey('Task-f6a91ccc6f9e4bcbace3331a0109c9da-ddotson-tyk2-demo')>,
<ScopedKey('Task-ce83f523df324508a6e6cfbce1c34393-ddotson-tyk2-demo')>]
We’ll choose one, and get all the failed results for it:
failures = asc.get_task_failures(failed_tasks[0])
failures
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-15fdd9ff3f68cde67e706ffa87f6a9b7-ddotson-tyk2-demo/failures/ProtocolDAGResultRef-2bc722e46be751f9ec1a770028f393f7-ddotson-tyk2-demo "HTTP/1.1 200 OK"
[<ProtocolDAGResult-4dd969e50a25508075044c4b9a61b7fe>]
Each of these is a ProtocolDAGResult, but with at least one ProtocolUnitFailure:
failures[0]
<ProtocolDAGResult-4dd969e50a25508075044c4b9a61b7fe>
failures[0].protocol_unit_failures
[ProtocolUnitFailure(lig_jmc_27 to lig_jmc_28 repeat 0 generation 0),
ProtocolUnitFailure(lig_jmc_27 to lig_jmc_28 repeat 0 generation 0),
ProtocolUnitFailure(lig_jmc_27 to lig_jmc_28 repeat 0 generation 0),
ProtocolUnitFailure(lig_jmc_27 to lig_jmc_28 repeat 0 generation 0)]
We can introspect the traceback on any of the ProtocolUnitFailures present in the ProtocolDAGResult to understand what went wrong here:
print(failures[0].protocol_unit_failures[0].traceback)
Traceback (most recent call last):
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/gufe/protocols/protocolunit.py", line 320, in execute
outputs = self._execute(context, **inputs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/equil_rfe_methods.py", line 1129, in _execute
outputs = self.run(scratch_basepath=ctx.scratch,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/equil_rfe_methods.py", line 998, in run
sampler.setup(
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/_rfe_utils/multistate.py", line 121, in setup
minimize(compound_thermostate_copy, sampler_state,
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/_rfe_utils/multistate.py", line 295, in minimize
context, integrator = dummy_cache.get_context(
^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openmmtools/cache.py", line 770, in get_context
context = thermodynamic_state.create_context(integrator, self.platform)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openmmtools/states.py", line 1179, in create_context
return openmm.Context(system, integrator, platform)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openmm/openmm.py", line 8037, in __init__
_openmm.Context_swiginit(self, _openmm.new_Context(*args))
^^^^^^^^^^^^^^^^^^^^^^^^^^
openmm.OpenMMException: Error initializing CUDA: CUDA_ERROR_MPS_CONNECTION_FAILED (805) at /home/conda/feedstock_root/build_artifacts/openmm_1721257909416/work/platforms/cuda/src/CudaContext.cpp:91
In this case, it like it may be due to a configuration issue on a compute resource our Tasks are landing on.
If we wish to re-run an 'error'ed Task, we can do so by setting its status back to 'waiting':
asc.set_tasks_status(failed_tasks[:1], 'waiting')
INFO: HTTP Request: POST https://api.alchemiscale.org/bulk/tasks/status/set "HTTP/1.1 200 OK"
[<ScopedKey('Task-ef9ccf8712d6458582e075a5be87d06e-ddotson-tyk2-demo')>]
Gathering and analyzing results¶
We can gather results for each Transformation with:
results = dict()
for tf_sk in asc.get_network_transformations(an_sk):
results[str(tf_sk)] = asc.get_transformation_results(tf_sk, visualize=False)
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-82e90d15b483a813a6501073e3964930-ddotson-tyk2-demo/results/ProtocolDAGResultRef-01a2e46b4c33a9c52aab71b8c01ba5cd-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-73164716e2ad3c2e410fce413fa41be4-ddotson-tyk2-demo/results/ProtocolDAGResultRef-caabc7fbbe8e9bae4ad0da4f5db7cb61-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-e3022bfd607e772136f5c262875a8e78-ddotson-tyk2-demo/results/ProtocolDAGResultRef-b6b00c75b4f66e546623a98d42b165c8-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-6064e55fed6e90b90b727b299462cbad-ddotson-tyk2-demo/results/ProtocolDAGResultRef-f87f2f38715e0e03c30b574fe3c0e01d-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-28453f83bb03486a08afc096a8dc5298-ddotson-tyk2-demo/results/ProtocolDAGResultRef-d206475f1ae25c4818964579eabca7b0-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-eebbe29feefc097086749081844d8adc-ddotson-tyk2-demo/results/ProtocolDAGResultRef-2a9f4953c83460388614a4cb872bc601-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-2a345646fa5b3900aba8b74935fc4ec6-ddotson-tyk2-demo/results/ProtocolDAGResultRef-821c1a7ccc7740a28f4dd3d553b2f035-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-f8fb4a09881bce886dcb55a88be49016-ddotson-tyk2-demo/results/ProtocolDAGResultRef-8295104d5c32f10bdd86b69675bb7494-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-874c39cca17041ceb7a91f142104e8eb-ddotson-tyk2-demo/results/ProtocolDAGResultRef-2d7d1eed2ee1b96e67c95385c14b220c-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-b4857c06e9564cf3c2451f436f7af410-ddotson-tyk2-demo/results/ProtocolDAGResultRef-fdea668674894a4ebc21cd51a9ca8e02-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-958a3d55d581ec276b72394435df4bd0-ddotson-tyk2-demo/results/ProtocolDAGResultRef-19825c22d09bf2a8d70e47e5cac34829-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-22a2847c2a28584956e90d549fc6bcc6-ddotson-tyk2-demo/results/ProtocolDAGResultRef-83f986401353591db65c9ed401af260b-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-d5cdd9f4f4edd380361129082a01360b-ddotson-tyk2-demo/results/ProtocolDAGResultRef-57a8c3097273d30ff5b0350943883824-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-20fd18a057e4affd5c78cad41a81d3b4-ddotson-tyk2-demo/results/ProtocolDAGResultRef-2798447525f7139b158919d08faadb95-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-44cf523e9d0e53674fadc59c9b9110b7-ddotson-tyk2-demo/results/ProtocolDAGResultRef-f56c970e363dbe075f2dda74386055a8-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-856f85d23d9b8b84916d5b83ebfa1773-ddotson-tyk2-demo/results/ProtocolDAGResultRef-f94dec1b5cc13aff45aba30eb1647c3c-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-ef57a778259ba88b9fb4c23cc2bd78f4-ddotson-tyk2-demo/results/ProtocolDAGResultRef-474df4f73b5bdb270a43361d8651dfa4-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-c33480c9f0dfd80852541205f46914a7-ddotson-tyk2-demo/results/ProtocolDAGResultRef-943e801e0c7774f3c962ca8127150ace-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-4bfeb3edd7f0d47ec52f533a5cf435a0-ddotson-tyk2-demo/results/ProtocolDAGResultRef-0d56400cc2c06fe557f170b872c9a465-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-06eac733376ab1ff566248040fe37e01-ddotson-tyk2-demo/results/ProtocolDAGResultRef-cc1df82df45efd8fba44f029a28094b6-ddotson-tyk2-demo "HTTP/1.1 200 OK"
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-76de739d87d53981abfff4ab9ce071e4-ddotson-tyk2-demo/results/ProtocolDAGResultRef-810f1f48475da9a8cbb7ed364848cdf2-ddotson-tyk2-demo "HTTP/1.1 200 OK"
results
{'Transformation-15fdd9ff3f68cde67e706ffa87f6a9b7-ddotson-tyk2-demo': None,
'Transformation-152928d4d69298dcedf2dd622073a95b-ddotson-tyk2-demo': None,
'Transformation-82e90d15b483a813a6501073e3964930-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-948bfdff0b72d5ef4f67adc85c2634c6>,
'Transformation-73164716e2ad3c2e410fce413fa41be4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-df5750fe54c2c580cfd73fa84ce840f1>,
'Transformation-e3022bfd607e772136f5c262875a8e78-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-337e4d8e80363a149a0d04f2b1ac2446>,
'Transformation-6064e55fed6e90b90b727b299462cbad-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-cf56bc1ab11ffae35e0dd045d6fe29dc>,
'Transformation-28453f83bb03486a08afc096a8dc5298-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-8226e5d9a3fa772e8923f309174d4671>,
'Transformation-eebbe29feefc097086749081844d8adc-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-f91bea18368378e55ace96616d6ec6f1>,
'Transformation-2a345646fa5b3900aba8b74935fc4ec6-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-fa766a1fc145b6700102b9fdff3d0db3>,
'Transformation-f8fb4a09881bce886dcb55a88be49016-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-a7a764961d833469c7da71645ddcc7cd>,
'Transformation-874c39cca17041ceb7a91f142104e8eb-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-c7576997d867d44c08191252879de95b>,
'Transformation-b4857c06e9564cf3c2451f436f7af410-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-4c344a156e08ebb132223c1b86dd3e45>,
'Transformation-958a3d55d581ec276b72394435df4bd0-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-0dff0aaf2d05bb1245a49961851d462c>,
'Transformation-22a2847c2a28584956e90d549fc6bcc6-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-6e452d9c821c26e75d2dc8bae917088d>,
'Transformation-d5cdd9f4f4edd380361129082a01360b-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-aad5127a208dc1f5dc23a3eca56b63da>,
'Transformation-20fd18a057e4affd5c78cad41a81d3b4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-d0d4d4749d8e0bab09409247b1f40d44>,
'Transformation-44cf523e9d0e53674fadc59c9b9110b7-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-11ede0f561d2d787f0ce7d49d825cea1>,
'Transformation-856f85d23d9b8b84916d5b83ebfa1773-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-be791dff8238a59df7d08e6973844fbb>,
'Transformation-ef57a778259ba88b9fb4c23cc2bd78f4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-626117d35fc898d368179c703b79ae44>,
'Transformation-c33480c9f0dfd80852541205f46914a7-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-6948f63f1a437cfffa21869cc8c59205>,
'Transformation-4bfeb3edd7f0d47ec52f533a5cf435a0-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-8a44fd759abe0b2ee5e5e1b3354712b1>,
'Transformation-8e74110792d4c68bb5a772733d490401-ddotson-tyk2-demo': None,
'Transformation-06eac733376ab1ff566248040fe37e01-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-b8c7ed6500e0de9a622001aa69793dbc>,
'Transformation-76de739d87d53981abfff4ab9ce071e4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-49834c8fefc70456a6bcf909bc86d2dd>}
results['Transformation-82e90d15b483a813a6501073e3964930-ddotson-tyk2-demo'].get_estimate()
20.937595533249475 kilocalorie_per_mole
results['Transformation-82e90d15b483a813a6501073e3964930-ddotson-tyk2-demo'].get_uncertainty()
0.0 kilocalorie_per_mole
In this case, we have only a single ProtocolDAGResult for each Transformation (since we created and actioned only 1 Task for each), and so the uncertainty (standard deviation between replicate simulations) given for this ProtocolResult is 0.0. The RelativeHybridTopologyProtocol combines ProtocolDAGResult values statistically, reducing the uncertainty but not increasing convergence with additional Tasks
By contrast other Protocols, such as the feflow NonEquilibriumCyclingProtocol, will improve convergence with more Tasks on a given Transformation; in the case of the NonEquilibriumCyclingProtocol, the non-equilibrium work values for each ProtocolDAGResult are combined together and fed to BAR to produce a single estimate with its own uncertainty.
Let’s create and action 2 additional Tasks for each Transformation:
tasks = []
for tf_sk in asc.get_network_transformations(an_sk):
tasks.extend(asc.create_tasks(tf_sk, count=2))
asc.action_tasks(tasks, an_sk)
[<ScopedKey('Task-a3bb16b506354cc1920c60931fc3b6e7-ddotson-tyk2-demo')>,
<ScopedKey('Task-5519d7feced84701bcd4caa4edf280de-ddotson-tyk2-demo')>,
<ScopedKey('Task-111a3cff46504f3abe702a16b5b164c1-ddotson-tyk2-demo')>,
<ScopedKey('Task-32c8103f4d624fcca7985557caeae893-ddotson-tyk2-demo')>,
<ScopedKey('Task-0cff85cbbcef4ae49336957c4bbbb285-ddotson-tyk2-demo')>,
<ScopedKey('Task-0a749b1896db48bf80210690b0d9ac6f-ddotson-tyk2-demo')>,
<ScopedKey('Task-88a7d011602c4ce39c5ff94cbb7a98ae-ddotson-tyk2-demo')>,
<ScopedKey('Task-1555231f3cd24d8db03a8e8a327adbba-ddotson-tyk2-demo')>,
<ScopedKey('Task-974c0495a7624e1aa2996a76c80c32fb-ddotson-tyk2-demo')>,
<ScopedKey('Task-3c291b7e214a4dba8ee403f939af516c-ddotson-tyk2-demo')>,
<ScopedKey('Task-7ff6c66b32644095bd881b3b9c723e3a-ddotson-tyk2-demo')>,
<ScopedKey('Task-8872ddfbd2534af390543d769b8587a9-ddotson-tyk2-demo')>,
<ScopedKey('Task-c8db674e068241cabb0f95ad03dfbb68-ddotson-tyk2-demo')>,
<ScopedKey('Task-0d72db080df544068687f3c4cfe0c4b4-ddotson-tyk2-demo')>,
<ScopedKey('Task-77457ab4507d43b996eb004856af7a1b-ddotson-tyk2-demo')>,
<ScopedKey('Task-c1f28d1a77de45bfb0acbd287f343fac-ddotson-tyk2-demo')>,
<ScopedKey('Task-61bee1dabd264ef88d558acec688f5fc-ddotson-tyk2-demo')>,
<ScopedKey('Task-c0c2070639064571b205f56eb7a695ab-ddotson-tyk2-demo')>,
<ScopedKey('Task-4aeb583ab91d40388eee8e8f13b6cc1e-ddotson-tyk2-demo')>,
<ScopedKey('Task-a0537c24e96c4cc4950f3e37b807a879-ddotson-tyk2-demo')>,
<ScopedKey('Task-1213ffa058ec4395bc96ba779f46609f-ddotson-tyk2-demo')>,
<ScopedKey('Task-7ac8ed14ac0f427da03689a87c9cf44c-ddotson-tyk2-demo')>,
<ScopedKey('Task-4b02f8758a934223b71395eaa8b26db6-ddotson-tyk2-demo')>,
<ScopedKey('Task-c5fd174873f5433cb23da737fa4642aa-ddotson-tyk2-demo')>,
<ScopedKey('Task-c3e1dc4d32634258acda412544ae778e-ddotson-tyk2-demo')>,
<ScopedKey('Task-ed05fb435acd4fca8661684a92cbe246-ddotson-tyk2-demo')>,
<ScopedKey('Task-134d60f5433e4e0babb1a81e24bdd63a-ddotson-tyk2-demo')>,
<ScopedKey('Task-19154647d7264f729cef865c47c55122-ddotson-tyk2-demo')>,
<ScopedKey('Task-a318bc8e78dd4625a6028ae0801a998c-ddotson-tyk2-demo')>,
<ScopedKey('Task-5312ba1a78784f909f5f4a767fc6fc49-ddotson-tyk2-demo')>,
<ScopedKey('Task-bb556e64e1a24bf8998b23b3f78a7ece-ddotson-tyk2-demo')>,
<ScopedKey('Task-ad7a64014b704fa685254441bf16e5ef-ddotson-tyk2-demo')>,
<ScopedKey('Task-aeee5cf01ef0494e8cb51c2e86324cf6-ddotson-tyk2-demo')>,
<ScopedKey('Task-789fbed60ab342dea8d071e2d925279f-ddotson-tyk2-demo')>,
<ScopedKey('Task-50033fbe0ffe4ee185566164cb33af65-ddotson-tyk2-demo')>,
<ScopedKey('Task-87a6955b36be404293b9491f476f7045-ddotson-tyk2-demo')>,
<ScopedKey('Task-8dc6bc5c62f044258c5b129a8c4d570f-ddotson-tyk2-demo')>,
<ScopedKey('Task-3c387ea6205447d786f9c1431591bef9-ddotson-tyk2-demo')>,
<ScopedKey('Task-8e6815ffb28b447988722bca9ddc8039-ddotson-tyk2-demo')>,
<ScopedKey('Task-5779831c1fcd4bdc819aefc608890940-ddotson-tyk2-demo')>,
<ScopedKey('Task-1b17a99a7acf4fea92879158e87b82ff-ddotson-tyk2-demo')>,
<ScopedKey('Task-37b15e41fba948c286cf73d70eff627f-ddotson-tyk2-demo')>,
<ScopedKey('Task-dea84bf82e354bb39c9b2d70fc91709d-ddotson-tyk2-demo')>,
<ScopedKey('Task-f1a16bc51f934dd19736cbf34e76804f-ddotson-tyk2-demo')>,
<ScopedKey('Task-1c8808c0d25d40e3bb02ceaa71132402-ddotson-tyk2-demo')>,
<ScopedKey('Task-fa2d45744f214c09a53246456d7dd27e-ddotson-tyk2-demo')>,
<ScopedKey('Task-52a5e061bec9413e912eaafaa1967456-ddotson-tyk2-demo')>,
<ScopedKey('Task-7d1c27206cd44bf296bf46f8f5267a2b-ddotson-tyk2-demo')>]
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 21 │ │ running │ 0 │ │ waiting │ 49 │ │ error │ 2 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'waiting': 49, 'error': 2, 'complete': 21}
As we watch our Tasks 'complete' or 'error', we can introspect and deal with errors as before:
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 46 │ │ running │ 24 │ │ waiting │ 0 │ │ error │ 2 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'complete': 46, 'running': 24, 'error': 2}
asc.get_network_tasks(an_sk, status='error')
[<ScopedKey('Task-f6a91ccc6f9e4bcbace3331a0109c9da-ddotson-tyk2-demo')>,
<ScopedKey('Task-ce83f523df324508a6e6cfbce1c34393-ddotson-tyk2-demo')>]
print(asc.get_task_failures(asc.get_network_tasks(an_sk, status='error')[0])[0].protocol_unit_failures[0].traceback)
INFO: HTTP Request: GET https://api.alchemiscale.org/transformations/Transformation-152928d4d69298dcedf2dd622073a95b-ddotson-tyk2-demo/failures/ProtocolDAGResultRef-08d10b9130b46d11e70cdf91343e686e-ddotson-tyk2-demo "HTTP/1.1 200 OK"
Traceback (most recent call last):
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/gufe/protocols/protocolunit.py", line 320, in execute
outputs = self._execute(context, **inputs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/equil_rfe_methods.py", line 1129, in _execute
outputs = self.run(scratch_basepath=ctx.scratch,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/equil_rfe_methods.py", line 998, in run
sampler.setup(
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/_rfe_utils/multistate.py", line 121, in setup
minimize(compound_thermostate_copy, sampler_state,
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openfe/protocols/openmm_rfe/_rfe_utils/multistate.py", line 295, in minimize
context, integrator = dummy_cache.get_context(
^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openmmtools/cache.py", line 770, in get_context
context = thermodynamic_state.create_context(integrator, self.platform)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openmmtools/states.py", line 1179, in create_context
return openmm.Context(system, integrator, platform)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/lila/home/dotson/mambaforge/envs/alchemiscalemiscale-compute-ddotson-v0.5.0-2024.08.15/lib/python3.12/site-packages/openmm/openmm.py", line 8037, in __init__
_openmm.Context_swiginit(self, _openmm.new_Context(*args))
^^^^^^^^^^^^^^^^^^^^^^^^^^
openmm.OpenMMException: Error initializing CUDA: CUDA_ERROR_MPS_CONNECTION_FAILED (805) at /home/conda/feedstock_root/build_artifacts/openmm_1721257909416/work/platforms/cuda/src/CudaContext.cpp:91
asc.set_tasks_status(asc.get_network_tasks(an_sk, status='error'), 'waiting')
INFO: HTTP Request: POST https://api.alchemiscale.org/bulk/tasks/status/set "HTTP/1.1 200 OK"
[<ScopedKey('Task-f6a91ccc6f9e4bcbace3331a0109c9da-ddotson-tyk2-demo')>,
<ScopedKey('Task-ce83f523df324508a6e6cfbce1c34393-ddotson-tyk2-demo')>]
asc.get_network_status(an_sk)
AlchemicalNetwork-6f22642592f789c4e7f918412ca947c5-ddotson-tyk2-demo ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ status ┃ count ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ complete │ 72 │ │ running │ 0 │ │ waiting │ 0 │ │ error │ 0 │ │ invalid │ 0 │ │ deleted │ 0 │ └──────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────┘
{'complete': 72}
Finally, with more sampling, we can pull our results again, this time with more ProtocolDAGResults included in the ProtocolResult for our Transformations:
results = dict()
for tf_sk in asc.get_network_transformations(an_sk):
results[str(tf_sk)] = asc.get_transformation_results(tf_sk, visualize=False)
results
{'Transformation-15fdd9ff3f68cde67e706ffa87f6a9b7-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-394e0443ec86ffc0c16c8e00c82e9b12>,
'Transformation-152928d4d69298dcedf2dd622073a95b-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-d116d14de4407acbeaac9f29587719e9>,
'Transformation-82e90d15b483a813a6501073e3964930-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-dcd71f35d39c66af70ed89aa7a279000>,
'Transformation-73164716e2ad3c2e410fce413fa41be4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-ad94c858081e514585b43f65c5a68f18>,
'Transformation-e3022bfd607e772136f5c262875a8e78-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-18c2551684f6b806fe8b6a9ece095009>,
'Transformation-6064e55fed6e90b90b727b299462cbad-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-d55f3325367e5a7f0a211077df9ccc7e>,
'Transformation-28453f83bb03486a08afc096a8dc5298-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-09120f442dacd96143ee12c6d166ef79>,
'Transformation-eebbe29feefc097086749081844d8adc-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-fe22eaa7d7ebd7a5d23adae41e98e1f6>,
'Transformation-2a345646fa5b3900aba8b74935fc4ec6-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-ed6c30bbaeb8124335f6080b648f78da>,
'Transformation-f8fb4a09881bce886dcb55a88be49016-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-6365aa29ece7fbc051405ba330d52328>,
'Transformation-874c39cca17041ceb7a91f142104e8eb-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-1d47e63b207e04a4960741d14d80480c>,
'Transformation-b4857c06e9564cf3c2451f436f7af410-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-febae6aa0f84830b17f89913b20b3e31>,
'Transformation-958a3d55d581ec276b72394435df4bd0-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-501a8034121f24d44579bc6edd1ba167>,
'Transformation-22a2847c2a28584956e90d549fc6bcc6-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-ef58753a1b7c856a2831287b22a11fa2>,
'Transformation-d5cdd9f4f4edd380361129082a01360b-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-0172ce2b7010333b1c341e7d14eb3a81>,
'Transformation-20fd18a057e4affd5c78cad41a81d3b4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-32536721cc9db99a4929e743344fb5ff>,
'Transformation-44cf523e9d0e53674fadc59c9b9110b7-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-761f181afdb778b924eaaf8f096aff42>,
'Transformation-856f85d23d9b8b84916d5b83ebfa1773-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-ed7a2dd2de1cb55ecd338eee10989819>,
'Transformation-ef57a778259ba88b9fb4c23cc2bd78f4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-889ff92aa5ec0aa090eca47f055ecd0e>,
'Transformation-c33480c9f0dfd80852541205f46914a7-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-bb3df5d7f37c55d3f4eddd4fa553d0db>,
'Transformation-4bfeb3edd7f0d47ec52f533a5cf435a0-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-09615a191ec58d54947dd7c974c9da3d>,
'Transformation-8e74110792d4c68bb5a772733d490401-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-2e0e7ef5064dc92335ad8baf158a2052>,
'Transformation-06eac733376ab1ff566248040fe37e01-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-836ae5c50d7bdeaa4273fa9aa0835226>,
'Transformation-76de739d87d53981abfff4ab9ce071e4-ddotson-tyk2-demo': <RelativeHybridTopologyProtocolResult-2dbee9e919321cdfac39967b797f27de>}
Advanced example: comparing results to experimental values with cinnabar¶
Source notebook: https://github.com/IAlibay/alchemiscale-utilities/blob/main/notebooks/TYK2_MST/Alchemiscale_submit_and_analyze.ipynb
We can gather up results for each Transformation as above and put into a form suitable for use with cinnabar. This will allow us to compare to experimental values for the binding free energies of these ligands. We use experimental values listed in ligands.yml, originally sourced from openforcefield/protein-ligand-benchmark:
%cat ligands.yml
lig_ejm_31:
measurement:
comment: Table 4, entry 31
doi: 10.1016/j.ejmech.2013.03.070
error: 0.029
type: ki
unit: uM
value: 0.096
name: lig_ejm_31
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C([H])([H])[H])[H])[H])Cl)[H]'
lig_ejm_42:
measurement:
comment: Table 4, entry 42
doi: 10.1016/j.ejmech.2013.03.070
error: 0.019
type: ki
unit: uM
value: 0.064
name: lig_ejm_42
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C([H])([H])C([H])([H])[H])[H])[H])Cl)[H]'
lig_ejm_43:
measurement:
comment: Table 4, entry 43
doi: 10.1016/j.ejmech.2013.03.070
error: 0.25
type: ki
unit: uM
value: 0.84
name: lig_ejm_43
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C([H])(C([H])([H])[H])C([H])([H])[H])[H])[H])Cl)[H]'
lig_ejm_45:
measurement:
comment: Table 4, entry 45
doi: 10.1016/j.ejmech.2013.03.070
error: 0.028
type: ki
unit: uM
value: 0.094
name: lig_ejm_45
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C([H])([H])C3(C(C3([H])[H])([H])[H])[H])[H])[H])Cl)[H]'
lig_ejm_46:
measurement:
comment: Table 4, entry 46
doi: 10.1016/j.ejmech.2013.03.070
error: 0.0014
type: ki
unit: uM
value: 0.0048
name: lig_ejm_46
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C3(C(C3([H])[H])([H])[H])[H])[H])[H])Cl)[H]'
lig_ejm_47:
measurement:
comment: Table 4, entry 47
doi: 10.1016/j.ejmech.2013.03.070
error: 0.022
type: ki
unit: uM
value: 0.074
name: lig_ejm_47
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C3(C(C(C3([H])[H])([H])[H])([H])[H])[H])[H])[H])Cl)[H]'
lig_ejm_48:
measurement:
comment: Table 4, entry 48
doi: 10.1016/j.ejmech.2013.03.070
error: 0.072
type: ki
unit: uM
value: 0.24
name: lig_ejm_48
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C3(C(C(C(C3([H])[H])([H])[H])([H])[H])([H])[H])[H])[H])[H])Cl)[H]'
lig_ejm_50:
measurement:
comment: Table 4, entry 50
doi: 10.1016/j.ejmech.2013.03.070
error: 0.075
type: ki
unit: uM
value: 0.25
name: lig_ejm_50
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)C([H])([H])O[H])[H])[H])Cl)[H]'
lig_ejm_54:
measurement:
comment: Table 4, entry 54
doi: 10.1016/j.ejmech.2013.03.070
error: 0.0054
type: ki
unit: uM
value: 0.018
name: lig_ejm_54
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)N([H])C([H])([H])C([H])([H])[H])[H])[H])Cl)[H]'
lig_ejm_55:
measurement:
comment: Table 4, entry 55
doi: 10.1016/j.ejmech.2013.03.070
error: 0.051
type: ki
unit: uM
value: 0.17
name: lig_ejm_55
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)OC([H])([H])[H])[H])[H])Cl)[H]'
lig_jmc_23:
measurement:
comment: Table 2, entry 23; values for four different enantionmers listed, this
one is for cis
doi: 10.1021/jm400266t
error: 0.75
type: ki
unit: nM
value: 2.5
name: lig_jmc_23
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)[C@@]3([C@](C3([H])[H])([H])F)[H])[H])[H])Cl)[H]'
lig_jmc_27:
measurement:
comment: Table 2, entry 27; cis enantiomer
doi: 10.1021/jm400266t
error: 1.5
type: ki
unit: nM
value: 5.1
name: lig_jmc_27
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)[C@@]3([C@](C3([H])[H])([H])Cl)[H])[H])[H])Cl)[H]'
lig_jmc_28:
measurement:
comment: Table 2, entry 28; cis enantiomer
doi: 10.1021/jm400266t
error: 2.6
type: ki
unit: nM
value: 8.5
name: lig_jmc_28
smiles: '[H]c1c(c(c(c(c1[H])Cl)C(=O)N([H])c2c(c(nc(c2[H])N([H])C(=O)[C@@]3([C@](C3([H])[H])([H])C([H])([H])[H])[H])[H])[H])Cl)[H]'
from pathlib import Path
# Next we define a convenience method to check what type of transformation we are looking at
def _scan_components(system):
comps = system.components.values()
if any([isinstance(comp, openfe.ProteinComponent) for comp in comps]):
return "complex"
elif any([isinstance(comp, openfe.SolventComponent) for comp in comps]):
return "solvent"
else:
return "vacuum"
# Next we create a results dictionary and scan our network edges to accumulate all
# the free energy results and their uncertainty
results_dir = Path('results')
results_dir.mkdir(parents=True, exist_ok=True)
results = dict()
for tf_sk in asc.get_network_transformations(an_sk):
transformation = asc.get_transformation(tf_sk)
result = asc.get_transformation_results(tf_sk)
if result is None:
continue
runtype = _scan_components(transformation.stateA)
mapping = transformation.mapping[0]
nameA = mapping.componentA.name
nameB = mapping.componentB.name
# store in accumulator
if f"{nameA}_{nameB}" in results.keys():
results[f"{nameA}_{nameB}"][runtype] = result
else:
results[f"{nameA}_{nameB}"] = {runtype: result, 'molA': nameA, 'molB': nameB}
# output individual results to a separate `.dat` file for future use
filename = results_dir / f"{nameA}_{nameB}.{runtype}.results.dat"
output = f"{result.get_estimate()},{result.get_uncertainty()}"
with open(filename, 'w') as f:
f.write(output)
Writing out a cinnabar input CSV file¶
Since this is a known benchmark system, we have experimental values for each of our ligands. We can combine them with our results to create a cinnabar input CSV file.
# Here we create a dictionary of experimental values from an input protein-ligand benchmark ligands.yml
# First load the yaml data
import yaml
with open('ligands.yml') as stream:
exp_data = yaml.safe_load(stream)
# Define a method for converting between Ki to estimated DG
from openff.units import unit
import math
def ki_to_dg(
ki: unit.Quantity, uncertainty: unit.Quantity,
temperature: unit.Quantity = 298.15 * unit.kelvin
) -> tuple[unit.Quantity, unit.Quantity]:
"""
Convenience method to convert a Ki w/ a given uncertainty to an
experimental estimate of the binding free energy.
Parameters
----------
ki : unit.Quantity
Experimental Ki value (e.g. 5 * unit.nanomolar)
uncertainty : unit.Quantity
Experimental error. Note: returns 0 if =< 0 * unit.nanomolar.
temperature : unit.Quantity
Experimental temperature. Default: 298.15 * unit.kelvin.
Returns
-------
DG : unit.Quantity
Gibbs binding free energy.
dDG : unit.Quantity
Error in binding free energy.
"""
if ki > 1e-15 * unit.nanomolar:
DG = (unit.molar_gas_constant * temperature.to(unit.kelvin)
* math.log(ki / unit.molar)).to(unit.kilocalorie_per_mole)
else:
raise ValueError("negative Ki values are not supported")
if uncertainty > 0 * unit.molar:
dDG = (unit.molar_gas_constant * temperature.to(unit.kelvin)
* uncertainty / ki).to(unit.kilocalorie_per_mole)
else:
dDG = 0 * unit.kilocalorie_per_mole
return DG, dDG
from openff.units import unit
exp_values = {}
for lig in exp_data:
exp_units = unit(exp_data[lig]['measurement']['unit'])
exp_values[lig] = {}
DG, dDG = ki_to_dg(exp_data[lig]['measurement']['value'] * exp_units,
exp_data[lig]['measurement']['error'] * exp_units)
exp_values[lig]['value'] = DG
exp_values[lig]['error'] = dDG
import numpy as np
# write out the cinnabar input file
with open('cinnabar_input.csv', 'w') as f:
f.write("# Experimental block\n")
f.write("# Ligand, expt_DDG, expt_dDDG\n")
for entry in exp_values:
f.write(f"{entry},{exp_values[entry]['value'].m:.2f},{exp_values[entry]['error'].m:.2f}\n")
f.write('\n')
f.write('# Calculated block\n')
f.write('# Ligand1,Ligand2,calc_DDG,calc_dDDG(MBAR),calc_dDDG(additional)\n')
for entry in results:
estimate = (results[entry]['complex'].get_estimate()
- results[entry]['solvent'].get_estimate())
err = np.sqrt(results[entry]['complex'].get_uncertainty()**2
+ results[entry]['solvent'].get_uncertainty()**2)
molA = results[entry]['molA']
molB = results[entry]['molB']
f.write(f"{molA},{molB},{estimate.m:.2f},0,{err.m:.2f}\n")
Plotting out results using cinnabar components¶
Note: the cinnabar API will change in its next release.
import cinnabar
from cinnabar import plotting as cinnabar_plotting
from cinnabar import FEMap, femap
%matplotlib inline
Generating a Cinnabar FEMap and plotting out the network¶
First let’s load the data into cinnabar and draw out the network of free energy results:
fe = FEMap.from_csv('cinnabar_input.csv')
fe.generate_absolute_values() # Get MLE generated estimates of the absolute values
fe.draw_graph()
Plotting out the relative free energy results¶
Next we can go ahead and plot out the relative free energy results:
# note you can pass the filename argument to write things out to file
cinnabar_plotting.plot_DDGs(fe.to_legacy_graph(), figsize=5, xy_lim=[5, -5])
Plotting out the absolute free energy results¶
Finally let’s go ahead and plot out the MLE derived absolute free energies:
# note you can pass the filename argument to write to file
cinnabar_plotting.plot_DGs(fe.to_legacy_graph(), figsize=5, xy_lim=[5, -5])
We can also shift our free energies by the average experimental value to have DGs on the same scale as experiment:
data = femap.read_csv('cinnabar_input.csv')
exp_DG_sum = sum([data['Experimental'][i].DG for i in data['Experimental'].keys()])
shift = exp_DG_sum / len(data['Experimental'].keys())
cinnabar_plotting.plot_DGs(fe.to_legacy_graph(), figsize=5, shift=shift.m)
/home/david/micromamba/envs/alchemiscale-client-user-guide-2024.09.26/lib/python3.12/site-packages/cinnabar/femap.py:35: UserWarning: Assuming kcal/mol units on measurements
warnings.warn("Assuming kcal/mol units on measurements")
