From f7a4be9cb3af8114f5a3d5c39c658c2d2750807c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Mon, 2 Oct 2023 11:54:10 -0400 Subject: [PATCH 01/12] WIP -- neq cycling openfe rebase --- feflow/protocols/nonequilibrium_cycling.py | 36 ++++++++++++++++------ feflow/settings/neq_cycling.py | 7 ++++- 2 files changed, 33 insertions(+), 10 deletions(-) diff --git a/feflow/protocols/nonequilibrium_cycling.py b/feflow/protocols/nonequilibrium_cycling.py index 4f15882..717699b 100644 --- a/feflow/protocols/nonequilibrium_cycling.py +++ b/feflow/protocols/nonequilibrium_cycling.py @@ -16,6 +16,7 @@ ProtocolDAGResult, ) +from openfe.protocols.openmm_utils import system_creation from perses.app.relative_setup import RelativeFEPSetup from perses.app.setup_relative_calculation import get_openmm_platform from perses.annihilation.relative import HybridTopologyFactory @@ -180,7 +181,6 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): import openmm import openmm.unit as openmm_unit from openmmtools.integrators import PeriodicNonequilibriumIntegrator - from perses.utils.openeye import generate_unique_atom_names # Setting up logging to file in shared filesystem file_logger = logging.getLogger("neq-cycling") @@ -210,13 +210,6 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): else: receptor_top, receptor_pos = None, None - # Get ligands cheminformatics molecules - ligand_a = ligand_a.to_openff().to_openeye() - ligand_b = ligand_b.to_openff().to_openeye() - # Generating unique atom names for ligands -- openmmforcefields needs them - ligand_a = generate_unique_atom_names(ligand_a) - ligand_b = generate_unique_atom_names(ligand_b) - # Get solvent parameters from component if solvent_a: ion_concentration = solvent_a.ion_concentration.to_openmm() @@ -225,8 +218,33 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): else: ion_concentration, positive_ion, negative_ion = None, None, None - # Get settings + ## Up to this point we have protein top/pos, ligand components and solvent attributes in OpenMM units + + + # Get settings for system generator + forcefield_settings = settings.forcefield_settings thermodynamic_settings = settings.thermo_settings + system_settings = settings.system_settings + + # handle cache for system generator + if settings.forcefield_cache is not None: + ffcache = ctx.shared / settings.forcefield_cache + else: + ffcache = None + + system_generator = system_creation.get_system_generator( + forcefield_settings=forcefield_settings, + thermo_settings=thermodynamic_settings, + system_settings=system_settings, + cache=ffcache, + has_solvent=solvent_a is not None, + ) + + + + + # Get settings + phase = self._detect_phase(state_a, state_b) traj_save_frequency = settings.traj_save_frequency work_save_frequency = settings.work_save_frequency # Note: this is divisor of traj save freq. diff --git a/feflow/settings/neq_cycling.py b/feflow/settings/neq_cycling.py index 0b547b6..1ae1281 100644 --- a/feflow/settings/neq_cycling.py +++ b/feflow/settings/neq_cycling.py @@ -4,10 +4,11 @@ This module implements the objects that will be needed to run relative binding free energy calculations using perses. """ - +from typing import Optional from gufe.settings import Settings from openff.units import unit from pydantic import root_validator +from openfe.protocols.openmm_utils.omm_settings import SystemSettings # Default settings for the lambda functions x = 'lambda' @@ -44,6 +45,10 @@ class NonEquilibriumCyclingSettings(Settings): class Config: arbitrary_types_allowed = True + # System Settings (from openfe) + system_settings: SystemSettings + forcefield_cache: Optional[str] = "db.json" # TODO: Remove once it has been integrated with openfe settings + # Lambda settings lambda_functions = DEFAULT_ALCHEMICAL_FUNCTIONS From 39b41120a1795e25e72ad74eb82e31c8258f20d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 3 Oct 2023 23:54:45 -0400 Subject: [PATCH 02/12] WIP -- setup from openfe --- feflow/protocols/nonequilibrium_cycling.py | 86 +++++++++++++++++++++- 1 file changed, 85 insertions(+), 1 deletion(-) diff --git a/feflow/protocols/nonequilibrium_cycling.py b/feflow/protocols/nonequilibrium_cycling.py index 717699b..3acb29f 100644 --- a/feflow/protocols/nonequilibrium_cycling.py +++ b/feflow/protocols/nonequilibrium_cycling.py @@ -22,7 +22,7 @@ from perses.annihilation.relative import HybridTopologyFactory from openff.units import unit -from openff.units.openmm import to_openmm +from openff.units.openmm import to_openmm, from_openmm # Specific instance of logger for this module # logger = logging.getLogger(__name__) @@ -180,7 +180,10 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): import numpy as np import openmm import openmm.unit as openmm_unit + from openff.units.openmm import ensure_quantity from openmmtools.integrators import PeriodicNonequilibriumIntegrator + from gufe.components import SmallMoleculeComponent + from openfe.protocols.openmm_rfe import _rfe_utils # Setting up logging to file in shared filesystem file_logger = logging.getLogger("neq-cycling") @@ -221,6 +224,8 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): ## Up to this point we have protein top/pos, ligand components and solvent attributes in OpenMM units + #### START OF COPY/PASTE OPENFE #### + # Get settings for system generator forcefield_settings = settings.forcefield_settings thermodynamic_settings = settings.thermo_settings @@ -240,6 +245,85 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): has_solvent=solvent_a is not None, ) + # b. force the creation of parameters + # This is necessary because we need to have the FF generated ahead of + # solvating the system. + # Note: by default this is cached to ctx.shared/db.json so shouldn't + # incur too large a cost + self.logger.info("Parameterizing molecules") + small_mols_a = [] + for comp in state_a.components.values(): + if isinstance(comp, SmallMoleculeComponent): + small_mols_a.append(comp) + + for comp in small_mols_a: + offmol = comp.to_openff() + system_generator.create_system(offmol.to_topology().to_openmm(), + molecules=[offmol]) + if comp == mapping.componentA: + molB = mapping.componentB.to_openff() + system_generator.create_system(molB.to_topology().to_openmm(), + molecules=[molB]) + + # c. get OpenMM Modeller + a dictionary of resids for each component + solvation_settings = settings.solvation_settings + stateA_modeller, comp_resids = system_creation.get_omm_modeller( + protein_comp=receptor_a, + solvent_comp=solvent_a, + small_mols=ligand_a, + omm_forcefield=system_generator.forcefield, + solvent_settings=solvation_settings, + ) + + # d. get topology & positions + # Note: roundtrip positions to remove vec3 issues + stateA_topology = stateA_modeller.getTopology() + stateA_positions = to_openmm( + from_openmm(stateA_modeller.getPositions()) + ) + + # e. create the stateA System + stateA_system = system_generator.create_system( + stateA_modeller.topology, + molecules=[s.to_openff() for s in small_mols_a], + ) + + # 2. Get stateB system + # a. get the topology + stateB_topology, stateB_alchem_resids = _rfe_utils.topologyhelpers.combined_topology( + stateA_topology, + mapping.componentB.to_openff().to_topology().to_openmm(), + exclude_resids=comp_resids[mapping.componentA], + ) + + # b. get a list of small molecules for stateB + off_mols_stateB = [mapping.componentB.to_openff(), ] + for comp in small_mols_a: + if comp != mapping.componentA: + off_mols_stateB.append(comp.to_openff()) + + stateB_system = system_generator.create_system( + stateB_topology, + molecules=off_mols_stateB, + ) + + # c. Define correspondence mappings between the two systems + ligand_mappings = _rfe_utils.topologyhelpers.get_system_mappings( + mapping.componentA_to_componentB, + stateA_system, stateA_topology, comp_resids[mapping.componentA], + stateB_system, stateB_topology, stateB_alchem_resids, + # These are non-optional settings for this method + fix_constraints=True, + ) + + # d. Finally get the positions + stateB_positions = _rfe_utils.topologyhelpers.set_and_check_new_positions( + ligand_mappings, stateA_topology, stateB_topology, + old_positions=ensure_quantity(stateA_positions, 'openmm'), + insert_positions=ensure_quantity(mapping.componentB.to_openff().conformers[0], 'openmm'), + ) + + ####### TO THIS PART IS COPY/PASTE FROM OPENFE ######### From c49c82f78d409a23c3d289655e1a71b7b9d368aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 3 Oct 2023 23:54:52 -0400 Subject: [PATCH 03/12] Solvation settings --- feflow/settings/neq_cycling.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/feflow/settings/neq_cycling.py b/feflow/settings/neq_cycling.py index 1ae1281..1860d5a 100644 --- a/feflow/settings/neq_cycling.py +++ b/feflow/settings/neq_cycling.py @@ -8,7 +8,7 @@ from gufe.settings import Settings from openff.units import unit from pydantic import root_validator -from openfe.protocols.openmm_utils.omm_settings import SystemSettings +from openfe.protocols.openmm_utils.omm_settings import SystemSettings, SolvationSettings # Default settings for the lambda functions x = 'lambda' @@ -49,6 +49,9 @@ class Config: system_settings: SystemSettings forcefield_cache: Optional[str] = "db.json" # TODO: Remove once it has been integrated with openfe settings + # Solvation settings + solvation_settings: SolvationSettings + # Lambda settings lambda_functions = DEFAULT_ALCHEMICAL_FUNCTIONS @@ -68,7 +71,7 @@ class Config: work_save_frequency: int = 500 atom_selection_expression: str = "not water" - # Number of cycles to run + # Number of replicates to run (1 cycle/replicate) num_replicates: int = 1 @root_validator From 856f5c0b80ce2ff1c557ea05b2af7101771cb1a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Wed, 4 Oct 2023 19:17:02 -0400 Subject: [PATCH 04/12] Using openfe -- settings and HTF --- feflow/protocols/nonequilibrium_cycling.py | 153 ++++++++++----------- feflow/settings/neq_cycling.py | 4 +- feflow/tests/conftest.py | 2 +- 3 files changed, 73 insertions(+), 86 deletions(-) diff --git a/feflow/protocols/nonequilibrium_cycling.py b/feflow/protocols/nonequilibrium_cycling.py index 3acb29f..9e84157 100644 --- a/feflow/protocols/nonequilibrium_cycling.py +++ b/feflow/protocols/nonequilibrium_cycling.py @@ -6,7 +6,7 @@ import logging import time -from gufe.settings import Settings, OpenMMSystemGeneratorFFSettings, ThermoSettings +from gufe.settings import Settings from gufe.chemicalsystem import ChemicalSystem from gufe.mapping import ComponentMapping from gufe.protocols import ( @@ -17,9 +17,7 @@ ) from openfe.protocols.openmm_utils import system_creation -from perses.app.relative_setup import RelativeFEPSetup from perses.app.setup_relative_calculation import get_openmm_platform -from perses.annihilation.relative import HybridTopologyFactory from openff.units import unit from openff.units.openmm import to_openmm, from_openmm @@ -260,14 +258,14 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): offmol = comp.to_openff() system_generator.create_system(offmol.to_topology().to_openmm(), molecules=[offmol]) - if comp == mapping.componentA: - molB = mapping.componentB.to_openff() - system_generator.create_system(molB.to_topology().to_openmm(), - molecules=[molB]) + if comp == ligand_a: + mol_b = ligand_b.to_openff() + system_generator.create_system(mol_b.to_topology().to_openmm(), + molecules=[mol_b]) # c. get OpenMM Modeller + a dictionary of resids for each component solvation_settings = settings.solvation_settings - stateA_modeller, comp_resids = system_creation.get_omm_modeller( + state_a_modeller, comp_resids = system_creation.get_omm_modeller( protein_comp=receptor_a, solvent_comp=solvent_a, small_mols=ligand_a, @@ -277,98 +275,76 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): # d. get topology & positions # Note: roundtrip positions to remove vec3 issues - stateA_topology = stateA_modeller.getTopology() - stateA_positions = to_openmm( - from_openmm(stateA_modeller.getPositions()) + state_a_topology = state_a_modeller.getTopology() + state_a_positions = to_openmm( + from_openmm(state_a_modeller.getPositions()) ) # e. create the stateA System - stateA_system = system_generator.create_system( - stateA_modeller.topology, + state_a_system = system_generator.create_system( + state_a_modeller.topology, molecules=[s.to_openff() for s in small_mols_a], ) # 2. Get stateB system # a. get the topology - stateB_topology, stateB_alchem_resids = _rfe_utils.topologyhelpers.combined_topology( - stateA_topology, - mapping.componentB.to_openff().to_topology().to_openmm(), - exclude_resids=comp_resids[mapping.componentA], + state_b_topology, state_b_alchem_resids = _rfe_utils.topologyhelpers.combined_topology( + state_a_topology, + ligand_b.to_openff().to_topology().to_openmm(), + exclude_resids=comp_resids[ligand_a], ) # b. get a list of small molecules for stateB - off_mols_stateB = [mapping.componentB.to_openff(), ] + off_mols_state_b = [ligand_b.to_openff(), ] for comp in small_mols_a: - if comp != mapping.componentA: - off_mols_stateB.append(comp.to_openff()) + if comp != ligand_a: + off_mols_state_b.append(comp.to_openff()) - stateB_system = system_generator.create_system( - stateB_topology, - molecules=off_mols_stateB, + state_b_system = system_generator.create_system( + state_b_topology, + molecules=off_mols_state_b, ) # c. Define correspondence mappings between the two systems ligand_mappings = _rfe_utils.topologyhelpers.get_system_mappings( - mapping.componentA_to_componentB, - stateA_system, stateA_topology, comp_resids[mapping.componentA], - stateB_system, stateB_topology, stateB_alchem_resids, + mapping.get("ligandmapping").componentA_to_componentB, + state_a_system, state_a_topology, comp_resids[ligand_a], + state_b_system, state_b_topology, state_b_alchem_resids, # These are non-optional settings for this method fix_constraints=True, ) # d. Finally get the positions - stateB_positions = _rfe_utils.topologyhelpers.set_and_check_new_positions( - ligand_mappings, stateA_topology, stateB_topology, - old_positions=ensure_quantity(stateA_positions, 'openmm'), - insert_positions=ensure_quantity(mapping.componentB.to_openff().conformers[0], 'openmm'), + state_b_positions = _rfe_utils.topologyhelpers.set_and_check_new_positions( + ligand_mappings, state_a_topology, state_b_topology, + old_positions=ensure_quantity(state_a_positions, 'openmm'), + insert_positions=ensure_quantity(ligand_b.to_openff().conformers[0], 'openmm'), ) - ####### TO THIS PART IS COPY/PASTE FROM OPENFE ######### - - - - # Get settings + # Get alchemical settings + alchemical_settings = settings.alchemical_settings + + # Now we can create the HTF from the previous objects + hybrid_factory = _rfe_utils.relative.HybridTopologyFactory( + state_a_system, state_a_positions, state_a_topology, + state_b_system, state_b_positions, state_b_topology, + old_to_new_atom_map=ligand_mappings['old_to_new_atom_map'], + old_to_new_core_atom_map=ligand_mappings['old_to_new_core_atom_map'], + use_dispersion_correction=alchemical_settings.use_dispersion_correction, + softcore_alpha=alchemical_settings.softcore_alpha, + softcore_LJ_v2=alchemical_settings.softcore_LJ_v2, + softcore_LJ_v2_alpha=alchemical_settings.softcore_alpha, + interpolate_old_and_new_14s=alchemical_settings.interpolate_old_and_new_14s, + flatten_torsions=alchemical_settings.flatten_torsions, + ) - phase = self._detect_phase(state_a, state_b) + ####### UP TO THIS PART IS COPY/PASTE FROM OPENFE ######### traj_save_frequency = settings.traj_save_frequency work_save_frequency = settings.work_save_frequency # Note: this is divisor of traj save freq. selection_expression = settings.atom_selection_expression - # Get the ligand mapping from ComponentMapping object - # NOTE: perses to date has a different "directionality" sense in terms of the mapping, - # see perses.rjmc.topology_proposal.propose docstring for detailed information. - ligand_mapping = mapping['ligand'].componentB_to_componentA - - # Setup relative FE calculation - fe_setup = RelativeFEPSetup( - old_ligand=ligand_a, - new_ligand=ligand_b, - receptor=receptor_top, - receptor_positions=receptor_pos, - forcefield_files=settings.forcefield_settings.forcefields, - small_molecule_forcefield=settings.forcefield_settings.small_molecule_forcefield, - phases=[phase], - transformation_atom_map=ligand_mapping, # Handle atom mapping between systems - ionic_strength=ion_concentration, - positive_ion=positive_ion, - negative_ion=negative_ion, - ) - - topology_proposals = fe_setup.topology_proposals - old_positions = fe_setup.old_positions - new_positions = fe_setup.new_positions - - # Generate Hybrid Topology Factory - Generic HTF - htf = HybridTopologyFactory( - topology_proposal=topology_proposals[phase], - current_positions=old_positions[phase], - new_positions=new_positions[phase], - softcore_LJ_v2=settings.softcore_LJ_v2, - interpolate_old_and_new_14s=settings.interpolate_old_and_new_14s, - ) - - system = htf.hybrid_system - positions = htf.hybrid_positions + system = hybrid_factory.hybrid_system + positions = hybrid_factory.hybrid_positions # Set up integrator temperature = to_openmm(thermodynamic_settings.temperature) @@ -415,13 +391,14 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): # Save positions if step % traj_save_frequency == 0: file_logger.debug(f"coarse step: {step}: saving trajectory (freq {traj_save_frequency})") - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_eq_old.append(initial_positions) forward_eq_new.append(final_positions) # Make sure trajectories are stored at the end of the eq loop file_logger.debug(f"coarse step: {step}: saving trajectory (freq {traj_save_frequency})") - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_eq_old.append(initial_positions) forward_eq_new.append(final_positions) @@ -438,12 +415,13 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): integrator.step(work_save_frequency) forward_works.append(integrator.get_protocol_work(dimensionless=True)) if fwd_step % traj_save_frequency == 0: - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_neq_old.append(initial_positions) forward_neq_new.append(final_positions) # Make sure trajectories are stored at the end of the neq loop - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_neq_old.append(initial_positions) forward_neq_new.append(final_positions) @@ -456,12 +434,13 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): for step in range(coarse_eq_steps): integrator.step(work_save_frequency) if step % traj_save_frequency == 0: - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) reverse_eq_new.append(initial_positions) # TODO: Maybe better naming not old/new but initial/final reverse_eq_old.append(final_positions) # Make sure trajectories are stored at the end of the eq loop - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) reverse_eq_old.append(initial_positions) reverse_eq_new.append(final_positions) @@ -476,12 +455,13 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): integrator.step(work_save_frequency) reverse_works.append(integrator.get_protocol_work(dimensionless=True)) if rev_step % traj_save_frequency == 0: - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) reverse_neq_old.append(initial_positions) reverse_neq_new.append(final_positions) # Make sure trajectories are stored at the end of the neq loop - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_eq_old.append(initial_positions) forward_eq_new.append(final_positions) @@ -501,6 +481,7 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): file_logger.info(f"replicate_{self.name} Estimated performance: {estimated_performance} ns/day") # Serialize works + phase = self._detect_phase(state_a, state_b) # infer phase from systems and components forward_work_path = ctx.shared / f"forward_{phase}_{self.name}.npy" reverse_work_path = ctx.shared / f"reverse_{phase}_{self.name}.npy" with open(forward_work_path, 'wb') as out_file: @@ -717,11 +698,17 @@ def __init__(self, settings: Settings): @classmethod def _default_settings(cls): - from perses.protocols.settings import NonEquilibriumCyclingSettings + from feflow.settings.neq_cycling import NonEquilibriumCyclingSettings + from gufe.settings import OpenMMSystemGeneratorFFSettings, ThermoSettings + from openfe.protocols.openmm_utils.omm_settings import SystemSettings, SolvationSettings + from openfe.protocols.openmm_rfe.equil_rfe_settings import AlchemicalSettings return NonEquilibriumCyclingSettings( - forcefield_settings=OpenMMSystemGeneratorFFSettings(), - thermo_settings=ThermoSettings(temperature=300 * unit.kelvin), - ) + forcefield_settings=OpenMMSystemGeneratorFFSettings(), + thermo_settings=ThermoSettings(temperature=300 * unit.kelvin), + system_settings=SystemSettings(), + solvation_settings=SolvationSettings(), + alchemical_settings=AlchemicalSettings(), + ) # NOTE: create method should be really fast, since it would be running in the work units not the clients!! def _create( diff --git a/feflow/settings/neq_cycling.py b/feflow/settings/neq_cycling.py index 1860d5a..e2c97d3 100644 --- a/feflow/settings/neq_cycling.py +++ b/feflow/settings/neq_cycling.py @@ -9,6 +9,7 @@ from openff.units import unit from pydantic import root_validator from openfe.protocols.openmm_utils.omm_settings import SystemSettings, SolvationSettings +from openfe.protocols.openmm_rfe.equil_rfe_settings import AlchemicalSettings # Default settings for the lambda functions x = 'lambda' @@ -56,8 +57,7 @@ class Config: lambda_functions = DEFAULT_ALCHEMICAL_FUNCTIONS # alchemical settings - softcore_LJ_v2 = True - interpolate_old_and_new_14s = False + alchemical_settings: AlchemicalSettings # NEQ integration settings timestep = 4.0 * unit.femtoseconds diff --git a/feflow/tests/conftest.py b/feflow/tests/conftest.py index 0f10b4a..32aee40 100644 --- a/feflow/tests/conftest.py +++ b/feflow/tests/conftest.py @@ -69,7 +69,7 @@ def toluene_solvent_system(toluene, solvent_comp): @pytest.fixture def short_settings(): from openff.units import unit - from perses.protocols import NonEquilibriumCyclingProtocol + from feflow.protocols import NonEquilibriumCyclingProtocol settings = NonEquilibriumCyclingProtocol.default_settings() From f0fcd6233f970354ff506aa3a34542f35a7c935d Mon Sep 17 00:00:00 2001 From: pulidoi <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 10 Oct 2023 13:51:11 -0400 Subject: [PATCH 05/12] Setting up CI --- .github/workflows/ci.yaml | 86 +++++++++++++++++++++++++++++++ devtools/conda-envs/test_env.yaml | 10 ---- pyproject.toml | 6 +-- 3 files changed, 89 insertions(+), 13 deletions(-) create mode 100644 .github/workflows/ci.yaml diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml new file mode 100644 index 0000000..b855a12 --- /dev/null +++ b/.github/workflows/ci.yaml @@ -0,0 +1,86 @@ +name: "CI" +on: + pull_request: + branches: + - main + push: + branches: + - main + schedule: + # At 07:00 UTC on Monday and Thursday. + - cron: "0 7 * * *" + release: + types: + - published + +concurrency: + group: "${{ github.workflow }}-${{ github.ref }}" + cancel-in-progress: true + +defaults: + run: + shell: bash -l {0} + +jobs: + tests: + runs-on: ${{ matrix.os }}-latest + name: "๐Ÿ’ป-${{matrix.os }} ๐Ÿ-${{ matrix.python-version }} ๐Ÿ—ƒ๏ธ${{ matrix.pydantic-version }}" + strategy: + fail-fast: false + matrix: + os: ["ubuntu"] + pydantic-version: [">1"] + python-version: + - "3.9" + - "3.10" + - "3.11" + include: + - os: "macos" + python-version: "3.11" + pydantic-version: ">1" + + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + + - name: "Setup Micromamba" + uses: mamba-org/setup-micromamba@v1 + with: + environment-file: devtools/conda-envs/test_env.yaml + environment-name: feflow-env + cache-environment: true + cache-downloads: true + create-args: >- + python=${{ matrix.python-version }} + init-shell: bash + + - name: "Install" + run: python -m pip install --no-deps -e . + + - name: "Test imports" + run: | + # if we add more to this, consider changing to for + env vars + python -Ic "import feflow; print(feflow.__version__)" + + - name: "Environment Information" + run: | + micromamba info + micromamba list + + - name: "Run tests" + env: + # Set the OFE_SLOW_TESTS to True if running a Cron job + OFE_SLOW_TESTS: ${{ fromJSON('{"false":"false","true":"true"}')[github.event_name != 'pull_request'] }} + run: | + pytest -n auto -v --cov=feflow --cov-report=xml --durations=10 + + - name: codecov + if: ${{ github.repository == 'choderalab/feflow' + && github.event_name == 'pull_request' }} + uses: codecov/codecov-action@v3 + with: + token: ${{ secrets.CODECOV_TOKEN }} + file: coverage.xml + fail_ci_if_error: False + verbose: True diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 0849b26..41cdbb6 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -2,26 +2,16 @@ name: feflow-test channels: - conda-forge - defaults - - openeye # TODO: Remove once we don't depend on openeye dependencies: # Base depends - openfe # TODO: Remove once we don't depend on openfe - python - pip - # Perses depends (TODO: Remove once we don't depend on perses) - - openmoltools - - cloudpathlib - - dask - - openeye-toolkits - # Testing - pytest - pytest-cov - codecov - # Pip-only installs - - pip: - - git+https://github.com/choderalab/perses.git@protocol-neqcyc # TODO: Remove once it only depends on openmmtools/openfe # - codecov diff --git a/pyproject.toml b/pyproject.toml index 683fa9c..c81476a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,7 +10,7 @@ description = "Recipes and protocols for molecular free energy calculations usin dynamic = ["version"] readme = "README.md" authors = [ - { name = "Ivรกn Pulido", email = "ivan.pulido@choderalab.org" } + { name = "ChoderaLab", email = "ivan.pulido@choderalab.org" } ] license = { text = "MIT" } # See https://pypi.org/classifiers/ @@ -67,7 +67,7 @@ feflow = [ ] [tool.versioningit] -default-version = "1+unknown" +default-version = "0.1.0+unknown" [tool.versioningit.format] distance = "{base_version}+{distance}.{vcs}{rev}" @@ -79,7 +79,7 @@ distance-dirty = "{base_version}+{distance}.{vcs}{rev}.dirty" method = "git" # <- The method name # Parameters to pass to the method: match = ["*"] -default-tag = "1.0.0" +default-tag = "0.1.0" [tool.versioningit.write] file = "feflow/_version.py" From e5c2481484603538348e8f604f4dd2e234a48af3 Mon Sep 17 00:00:00 2001 From: pulidoi <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 10 Oct 2023 13:51:22 -0400 Subject: [PATCH 06/12] Use feflow in tests --- feflow/tests/conftest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/feflow/tests/conftest.py b/feflow/tests/conftest.py index 32aee40..d9473c2 100644 --- a/feflow/tests/conftest.py +++ b/feflow/tests/conftest.py @@ -94,7 +94,7 @@ def short_settings_gpu(short_settings): @pytest.fixture def short_settings_multiple_cycles(): from openff.units import unit - from perses.protocols import NonEquilibriumCyclingProtocol + from feflow.protocols import NonEquilibriumCyclingProtocol settings = NonEquilibriumCyclingProtocol.default_settings() From b9e367e10158f8793884d1342e2abdc6669d0838 Mon Sep 17 00:00:00 2001 From: pulidoi <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 10 Oct 2023 13:54:20 -0400 Subject: [PATCH 07/12] Finishing rebase --- feflow/protocols/nonequilibrium_cycling.py | 14 +-- feflow/utils/hybrid_topology.py | 104 +++++++++++++++++++++ 2 files changed, 112 insertions(+), 6 deletions(-) create mode 100644 feflow/utils/hybrid_topology.py diff --git a/feflow/protocols/nonequilibrium_cycling.py b/feflow/protocols/nonequilibrium_cycling.py index 9e84157..f3cbe23 100644 --- a/feflow/protocols/nonequilibrium_cycling.py +++ b/feflow/protocols/nonequilibrium_cycling.py @@ -16,8 +16,9 @@ ProtocolDAGResult, ) +# TODO: Remove/change when things get migrated to openmmtools or feflow from openfe.protocols.openmm_utils import system_creation -from perses.app.setup_relative_calculation import get_openmm_platform +from openfe.protocols.openmm_rfe._rfe_utils.compute import get_openmm_platform from openff.units import unit from openff.units.openmm import to_openmm, from_openmm @@ -131,8 +132,8 @@ def extract_positions(context, hybrid_topology_factory, atom_selection_exp="not positions = context.getState(getPositions=True).getPositions(asNumpy=True) # Get topology from HTF - indices for initial and final topologies in hybrid topology - initial_indices = np.asarray(hybrid_topology_factory.old_atom_indices) - final_indices = np.asarray(hybrid_topology_factory.new_atom_indices) + initial_indices = np.asarray(hybrid_topology_factory.initial_atom_indices) + final_indices = np.asarray(hybrid_topology_factory.final_atom_indices) hybrid_topology = hybrid_topology_factory.hybrid_topology selection = atom_selection_exp md_trajectory = md.Trajectory(xyz=positions, topology=hybrid_topology) @@ -182,6 +183,7 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): from openmmtools.integrators import PeriodicNonequilibriumIntegrator from gufe.components import SmallMoleculeComponent from openfe.protocols.openmm_rfe import _rfe_utils + from feflow.utils.hybrid_topology import HybridTopologyFactoryModded as HybridTopologyFactory # Setting up logging to file in shared filesystem file_logger = logging.getLogger("neq-cycling") @@ -268,7 +270,7 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): state_a_modeller, comp_resids = system_creation.get_omm_modeller( protein_comp=receptor_a, solvent_comp=solvent_a, - small_mols=ligand_a, + small_mols=small_mols_a, omm_forcefield=system_generator.forcefield, solvent_settings=solvation_settings, ) @@ -307,7 +309,7 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): # c. Define correspondence mappings between the two systems ligand_mappings = _rfe_utils.topologyhelpers.get_system_mappings( - mapping.get("ligandmapping").componentA_to_componentB, + mapping["ligand"].componentA_to_componentB, state_a_system, state_a_topology, comp_resids[ligand_a], state_b_system, state_b_topology, state_b_alchem_resids, # These are non-optional settings for this method @@ -325,7 +327,7 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): alchemical_settings = settings.alchemical_settings # Now we can create the HTF from the previous objects - hybrid_factory = _rfe_utils.relative.HybridTopologyFactory( + hybrid_factory = HybridTopologyFactory( state_a_system, state_a_positions, state_a_topology, state_b_system, state_b_positions, state_b_topology, old_to_new_atom_map=ligand_mappings['old_to_new_atom_map'], diff --git a/feflow/utils/hybrid_topology.py b/feflow/utils/hybrid_topology.py new file mode 100644 index 0000000..cce0d3d --- /dev/null +++ b/feflow/utils/hybrid_topology.py @@ -0,0 +1,104 @@ +from openfe.protocols.openmm_rfe._rfe_utils.relative import HybridTopologyFactory + + +# TODO: This is an utility class. To be deleted when we migrate/extend the base HybridTopologyFactory +class HybridTopologyFactoryModded(HybridTopologyFactory): + """ + Utility class that extends the base HybridTopologyFactory class with properties for + getting the indices from initial and final states. + """ + + def __init__(self, + old_system, old_positions, old_topology, + new_system, new_positions, new_topology, + old_to_new_atom_map, old_to_new_core_atom_map, + use_dispersion_correction=False, + softcore_alpha=0.5, + softcore_LJ_v2=True, + softcore_LJ_v2_alpha=0.85, + interpolate_old_and_new_14s=False, + flatten_torsions=False, + **kwargs): + """ + Initialize the Hybrid topology factory. + + Parameters + ---------- + old_system : openmm.System + OpenMM system defining the "old" (i.e. starting) state. + old_positions : [n,3] np.ndarray of float + The positions of the "old system". + old_topology : openmm.Topology + OpenMM topology defining the "old" state. + new_system: opemm.System + OpenMM system defining the "new" (i.e. end) state. + new_positions : [m,3] np.ndarray of float + The positions of the "new system" + new_topology : openmm.Topology + OpenMM topology defining the "new" state. + old_to_new_atom_map : dict of int : int + Dictionary of corresponding atoms between the old and new systems. + Unique atoms are not included in this atom map. + old_to_new_core_atom_map : dict of int : int + Dictionary of corresponding atoms between the alchemical "core + atoms" (i.e. residues which are changing) between the old and + new systems. + use_dispersion_correction : bool, default False + Whether to use the long range correction in the custom sterics + force. This can be very expensive for NCMC. + softcore_alpha: float, default None + "alpha" parameter of softcore sterics, default 0.5. + softcore_LJ_v2 : bool, default True + Implement the softcore LJ as defined by Gapsys et al. JCTC 2012. + softcore_LJ_v2_alpha : float, default 0.85 + Softcore alpha parameter for LJ v2 + interpolate_old_and_new_14s : bool, default False + Whether to turn off interactions for new exceptions (not just + 1,4s) at lambda = 0 and old exceptions at lambda = 1; if False, + they are present in the nonbonded force. + flatten_torsions : bool, default False + If True, torsion terms involving `unique_new_atoms` will be + scaled such that at lambda=0,1, the torsion term is turned off/on + respectively. The opposite is true for `unique_old_atoms`. + """ + super().__init__(old_system, old_positions, old_topology, + new_system, new_positions, new_topology, + old_to_new_atom_map, old_to_new_core_atom_map, + use_dispersion_correction=False, + softcore_alpha=0.5, + softcore_LJ_v2=True, + softcore_LJ_v2_alpha=0.85, + interpolate_old_and_new_14s=False, + flatten_torsions=False, + **kwargs) + + # TODO: We need to refactor for the init to use these properties and have an attribute with the indices + @property + def initial_atom_indices(self): + """ + Indices of atoms in hybrid topology for atoms in the old topology. + + Returns + ------- + hybrid_indices: list[int] + Indices of old atoms in hybrid topology + + """ + n_atoms_old = self._old_system.getNumParticles() + hybrid_indices = [self._old_to_hybrid_map[idx] for idx in range(n_atoms_old)] + return hybrid_indices + + @property + def final_atom_indices(self): + """ + Indices of atoms in hybrid topology for atoms in the new topology. + + Returns + ------- + hybrid_indices: list[int] + Indices of new atoms in hybrid topology + + """ + n_atoms_new = self._new_system.getNumParticles() + hybrid_indices = [self._new_to_hybrid_map[idx] for idx in range(n_atoms_new)] + return hybrid_indices From b7426febc811219c2439771e5c4af32cd771080a Mon Sep 17 00:00:00 2001 From: pulidoi <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 10 Oct 2023 14:03:13 -0400 Subject: [PATCH 08/12] version variable --- feflow/__init__.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/feflow/__init__.py b/feflow/__init__.py index 9aef5c1..3add673 100644 --- a/feflow/__init__.py +++ b/feflow/__init__.py @@ -1,6 +1,5 @@ -"""Recipes and protocols for molecular free energy calculations using the openmmtools/perses and Open Free Energy toolkits""" - -# Add imports here -from .feflow import * +"""Recipes and protocols for molecular free energy calculations using the +openmmtools/perses and Open Free Energy toolkits""" from importlib.metadata import version +__version__ = version("openfe") From 0243d921980001cb3fd03c65c33346e51c5fc2b1 Mon Sep 17 00:00:00 2001 From: pulidoi <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 10 Oct 2023 14:10:02 -0400 Subject: [PATCH 09/12] Fix CI --- devtools/conda-envs/test_env.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 41cdbb6..2902edb 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -11,6 +11,7 @@ dependencies: # Testing - pytest - pytest-cov + - pytest-xdist - codecov # - codecov From b9afe790345c0b6bd0a6a2d3227d7f0f37b9426f Mon Sep 17 00:00:00 2001 From: pulidoi <2949729+ijpulidos@users.noreply.github.com> Date: Tue, 10 Oct 2023 16:45:22 -0400 Subject: [PATCH 10/12] cleanup setup htf --- feflow/protocols/nonequilibrium_cycling.py | 23 +--------------------- 1 file changed, 1 insertion(+), 22 deletions(-) diff --git a/feflow/protocols/nonequilibrium_cycling.py b/feflow/protocols/nonequilibrium_cycling.py index f3cbe23..d625c49 100644 --- a/feflow/protocols/nonequilibrium_cycling.py +++ b/feflow/protocols/nonequilibrium_cycling.py @@ -205,27 +205,6 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): solvent_a = state_a.components.get("solvent") # solvent_b = state_b.components.get("solvent") # Should not be needed - - # Check first state for receptor if not get receptor from second one - if receptor_a: - receptor_top = receptor_a.to_openmm_topology() - receptor_pos = receptor_a.to_openmm_positions() - else: - receptor_top, receptor_pos = None, None - - # Get solvent parameters from component - if solvent_a: - ion_concentration = solvent_a.ion_concentration.to_openmm() - positive_ion = solvent_a.positive_ion - negative_ion = solvent_a.negative_ion - else: - ion_concentration, positive_ion, negative_ion = None, None, None - - ## Up to this point we have protein top/pos, ligand components and solvent attributes in OpenMM units - - - #### START OF COPY/PASTE OPENFE #### - # Get settings for system generator forcefield_settings = settings.forcefield_settings thermodynamic_settings = settings.thermo_settings @@ -339,8 +318,8 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): interpolate_old_and_new_14s=alchemical_settings.interpolate_old_and_new_14s, flatten_torsions=alchemical_settings.flatten_torsions, ) + ####### END OF SETUP ######### - ####### UP TO THIS PART IS COPY/PASTE FROM OPENFE ######### traj_save_frequency = settings.traj_save_frequency work_save_frequency = settings.work_save_frequency # Note: this is divisor of traj save freq. selection_expression = settings.atom_selection_expression From 1fd47dd654061611b8b66764c23030c409a4a3fa Mon Sep 17 00:00:00 2001 From: David Dotson Date: Thu, 12 Oct 2023 16:39:24 -0700 Subject: [PATCH 11/12] Renamed neq_cycling settings module for consistency elsewhere --- feflow/protocols/nonequilibrium_cycling.py | 2 +- feflow/settings/{neq_cycling.py => nonequilibrium_cycling.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename feflow/settings/{neq_cycling.py => nonequilibrium_cycling.py} (100%) diff --git a/feflow/protocols/nonequilibrium_cycling.py b/feflow/protocols/nonequilibrium_cycling.py index d625c49..6a9f8ef 100644 --- a/feflow/protocols/nonequilibrium_cycling.py +++ b/feflow/protocols/nonequilibrium_cycling.py @@ -679,7 +679,7 @@ def __init__(self, settings: Settings): @classmethod def _default_settings(cls): - from feflow.settings.neq_cycling import NonEquilibriumCyclingSettings + from feflow.settings.nonequilibrium_cycling import NonEquilibriumCyclingSettings from gufe.settings import OpenMMSystemGeneratorFFSettings, ThermoSettings from openfe.protocols.openmm_utils.omm_settings import SystemSettings, SolvationSettings from openfe.protocols.openmm_rfe.equil_rfe_settings import AlchemicalSettings diff --git a/feflow/settings/neq_cycling.py b/feflow/settings/nonequilibrium_cycling.py similarity index 100% rename from feflow/settings/neq_cycling.py rename to feflow/settings/nonequilibrium_cycling.py From 81e1b80ea762e0ea78bbde6fc5e68d774ceadb2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Thu, 12 Oct 2023 23:49:36 -0400 Subject: [PATCH 12/12] Applying review suggestions --- .github/workflows/ci.yaml | 4 +--- devtools/conda-envs/test_env.yaml | 1 + feflow/utils/hybrid_topology.py | 12 ++++++------ 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index b855a12..e41b20e 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -24,12 +24,11 @@ defaults: jobs: tests: runs-on: ${{ matrix.os }}-latest - name: "๐Ÿ’ป-${{matrix.os }} ๐Ÿ-${{ matrix.python-version }} ๐Ÿ—ƒ๏ธ${{ matrix.pydantic-version }}" + name: "๐Ÿ’ป-${{matrix.os }} ๐Ÿ-${{ matrix.python-version }}" strategy: fail-fast: false matrix: os: ["ubuntu"] - pydantic-version: [">1"] python-version: - "3.9" - "3.10" @@ -37,7 +36,6 @@ jobs: include: - os: "macos" python-version: "3.11" - pydantic-version: ">1" steps: - uses: actions/checkout@v4 diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 2902edb..5456c86 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -4,6 +4,7 @@ channels: - defaults dependencies: # Base depends + - gufe - openfe # TODO: Remove once we don't depend on openfe - python - pip diff --git a/feflow/utils/hybrid_topology.py b/feflow/utils/hybrid_topology.py index cce0d3d..ee73a22 100644 --- a/feflow/utils/hybrid_topology.py +++ b/feflow/utils/hybrid_topology.py @@ -64,12 +64,12 @@ def __init__(self, super().__init__(old_system, old_positions, old_topology, new_system, new_positions, new_topology, old_to_new_atom_map, old_to_new_core_atom_map, - use_dispersion_correction=False, - softcore_alpha=0.5, - softcore_LJ_v2=True, - softcore_LJ_v2_alpha=0.85, - interpolate_old_and_new_14s=False, - flatten_torsions=False, + use_dispersion_correction=use_dispersion_correction, + softcore_alpha=softcore_alpha, + softcore_LJ_v2=softcore_LJ_v2, + softcore_LJ_v2_alpha=softcore_LJ_v2_alpha, + interpolate_old_and_new_14s=interpolate_old_and_new_14s, + flatten_torsions=flatten_torsions, **kwargs) # TODO: We need to refactor for the init to use these properties and have an attribute with the indices