From fb37bb6fe0f8a05202e81495d8e3c3c1935a42e6 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 24 Oct 2025 00:20:19 +0100 Subject: [PATCH 1/7] Add net charge corrections for ABFEs --- openfe/protocols/openmm_afe/base.py | 70 +++-- .../openmm_afe/equil_afe_settings.py | 17 ++ .../openmm_afe/equil_binding_afe_method.py | 261 ++++++++++++------ .../openmm_abfe/test_abfe_protocol.py | 8 +- 4 files changed, 251 insertions(+), 105 deletions(-) diff --git a/openfe/protocols/openmm_afe/base.py b/openfe/protocols/openmm_afe/base.py index d5a7b8f90..6bd650641 100644 --- a/openfe/protocols/openmm_afe/base.py +++ b/openfe/protocols/openmm_afe/base.py @@ -598,14 +598,28 @@ def _get_lambda_schedule( return lambdas + def _get_alchemical_ion( + self, + alchem_comps: dict[str, list[Component]], + comp_resids: dict[Component, npt.NDArray], + omm_topology: app.Topology, + positions: openmm.unit.Quantity, + settings: dict[str, SettingsBaseModel], + ) -> int | None: + """ + Placeholder method to find alchemical ions if necessary. + """ + return None + def _add_restraints( self, system: openmm.System, - topology: GlobalParameterState, + topology: app.Topology, positions: openmm.unit.Quantity, alchem_comps: dict[str, list[Component]], comp_resids: dict[Component, npt.NDArray], settings: dict[str, SettingsBaseModel], + alchem_ion: int | None, ) -> tuple[ Optional[GlobalParameterState], Optional[Quantity], @@ -623,7 +637,8 @@ def _get_alchemical_system( system: openmm.System, comp_resids: dict[Component, npt.NDArray], alchem_comps: dict[str, list[Component]], - ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, list[int]]: + alchem_ion: int | None, + ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, list[list[int]]]: """ Get an alchemically modified system and its associated factory @@ -644,7 +659,7 @@ def _get_alchemical_system( Factory for creating an alchemically modified system. alchemical_system : openmm.System Alchemically modified system - alchemical_indices : list[int] + alchemical_indices : list[list[int]] A list of atom indices for the alchemically modified species in the system. @@ -652,14 +667,21 @@ def _get_alchemical_system( ---- * Add support for all alchemical factory options """ - alchemical_indices = self._get_alchemical_indices(topology, comp_resids, alchem_comps) + alchemical_indices = [self._get_alchemical_indices(topology, comp_resids, alchem_comps)] - alchemical_region = AlchemicalRegion( - alchemical_atoms=alchemical_indices, - ) + if alchem_ion is not None: + alchemical_indices.append(alchem_ion) + + alchemical_regions = [] + + for inds in alchemical_indices: + alchemical_regions.append(AlchemicalRegion(alchemical_atoms=inds)) alchemical_factory = AbsoluteAlchemicalFactory() - alchemical_system = alchemical_factory.create_alchemical_system(system, alchemical_region) + alchemical_system = alchemical_factory.create_alchemical_system( + reference_system=system, + alchemical_regions=alchemical_regions + ) return alchemical_factory, alchemical_system, alchemical_indices @@ -1127,7 +1149,16 @@ def run( # 5. Get lambdas lambdas = self._get_lambda_schedule(settings) - # 6. Add restraints + # 6. Get alchemical ions + alchem_ion = self._get_alchemical_ion( + alchem_comps=alchem_comps, + comp_resids=comp_resids, + omm_topology=omm_topology, + positions=positions, + settings=settings, + ) + + # 7. Add restraints # Note: when no restraint is applied, restrained_omm_system == omm_system ( restraint_parameter_state, @@ -1141,14 +1172,19 @@ def run( alchem_comps, comp_resids, settings, + alchem_ion, ) - # 7. Get alchemical system + # 8. Get alchemical system alchem_factory, alchem_system, alchem_indices = self._get_alchemical_system( - omm_topology, restrained_omm_system, comp_resids, alchem_comps + omm_topology, + restrained_omm_system, + comp_resids, + alchem_comps, + alchem_ion, ) - # 8. Get compound and sampler states + # 9. Get compound and sampler states sampler_states, cmp_states = self._get_states( alchem_system, positions, @@ -1159,7 +1195,7 @@ def run( restraint_parameter_state, ) - # 9. Create the multistate reporter & create PDB + # 10. Create the multistate reporter & create PDB reporter = self._get_reporter( omm_topology, positions, @@ -1169,18 +1205,18 @@ def run( # Wrap in try/finally to avoid memory leak issues try: - # 10. Get context caches + # 11. Get context caches energy_ctx_cache, sampler_ctx_cache = self._get_ctx_caches( settings["forcefield_settings"], settings["engine_settings"] ) - # 11. Get integrator + # 12. Get integrator integrator = self._get_integrator( settings["integrator_settings"], settings["simulation_settings"], ) - # 12. Get sampler + # 13. Get sampler sampler = self._get_sampler( integrator, reporter, @@ -1192,7 +1228,7 @@ def run( sampler_ctx_cache, ) - # 13. Run simulation + # 14. Run simulation unit_result_dict = self._run_simulation( sampler, reporter, settings, standard_state_corr, dry ) diff --git a/openfe/protocols/openmm_afe/equil_afe_settings.py b/openfe/protocols/openmm_afe/equil_afe_settings.py index ba20db4dc..e0c472aec 100644 --- a/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -15,12 +15,14 @@ * Add support for restraints """ +from openff.units import unit as offunit from gufe.settings import ( SettingsBaseModel, OpenMMSystemGeneratorFFSettings, ThermoSettings, ) +from gufe.settings.typing import NanometerQuantity from openfe.protocols.openmm_utils.omm_settings import ( MultiStateSimulationSettings, BaseSolvationSettings, @@ -49,6 +51,21 @@ class AlchemicalSettings(SettingsBaseModel): """ +class ABFEAlchemicalSettings(AlchemicalSettings): + """ + Absolute binding free energy alchemical settings. + """ + explicit_charge_correction: bool = True + """ + Whether or not to use explicit charge correction using + a co-alchemical ion. + """ + alchemical_ion_min_distance: NanometerQuantity = 1.2 * offunit.nanometer + """ + The minimum distance to search for a co-alchemical ion. + """ + + class LambdaSettings(SettingsBaseModel): """Lambda schedule settings. diff --git a/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 9eca0d7f2..342707007 100644 --- a/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -47,7 +47,7 @@ from openfe.due import Doi, due from openfe.protocols.openmm_afe.equil_afe_settings import ( AbsoluteBindingSettings, - AlchemicalSettings, + ABFEAlchemicalSettings, BoreschRestraintSettings, IntegratorSettings, LambdaSettings, @@ -63,6 +63,7 @@ from openfe.protocols.openmm_utils import settings_validation, system_validation from openfe.protocols.restraint_utils import geometry from openfe.protocols.restraint_utils.geometry.boresch import BoreschRestraintGeometry +from openfe.protocols.restraint_utils.geometry.utils import FindHostAtoms from openfe.protocols.restraint_utils.openmm import omm_restraints from openfe.protocols.restraint_utils.openmm.omm_restraints import BoreschRestraint from openff.units import unit as offunit @@ -551,8 +552,10 @@ def _default_settings(cls): thermo_settings=settings.ThermoSettings( temperature=298.15 * offunit.kelvin, pressure=1 * offunit.bar, + ph=None, + redox_potential=None, ), - alchemical_settings=AlchemicalSettings(), + alchemical_settings=ABFEAlchemicalSettings(), solvent_lambda_settings=LambdaSettings( lambda_elec=[ 0.0, 0.25, 0.5, 0.75, 1.0, @@ -911,7 +914,171 @@ def _gather( return repeats -class AbsoluteBindingComplexUnit(BaseAbsoluteUnit): +def _get_mda_universe( + topology: omm_topology, + positions: ommunit.Quantity, + trajectory: Optional[pathlib.Path], +) -> mda.Universe: + """ + Helper method to get a Universe from an openmm Topology, + and either an input trajectory or a set of positions. + + Parameters + ---------- + topology : openmm.app.Topology + An OpenMM Topology that defines the System. + positions: openmm.unit.Quantity + The System's current positions. + Used if a trajectory file is None or is not a file. + trajectory: pathlib.Path + A Path to a trajectory file to read positions from. + + Returns + ------- + mda.Universe + An MDAnalysis Universe of the System. + """ + from MDAnalysis.coordinates.memory import MemoryReader + + # If the trajectory file doesn't exist, then we use positions + if trajectory is not None and trajectory.is_file(): + return mda.Universe( + topology, + trajectory, + topology_format="OPENMMTOPOLOGY", + ) + else: + # Positions is an openmm Quantity in nm we need + # to convert to angstroms + return mda.Universe( + topology, + np.array(positions._value) * 10, + topology_format="OPENMMTOPOLOGY", + trajectory_format=MemoryReader, + ) + + +def _get_idxs_from_residxs( + topology: omm_topology, + residxs: list[int], +) -> list[int]: + """ + Helper method to get the a list of atom indices which belong to a list + of residues. + + Parameters + ---------- + topology : openmm.app.Topology + An OpenMM Topology that defines the System. + residxs : list[int] + A list of residue numbers who's atoms we should get atom indices. + + Returns + ------- + atom_ids : list[int] + A list of atom indices. + + TODO + ---- + * Check how this works when we deal with virtual sites. + """ + atom_ids = [] + + for r in topology.residues(): + if r.index in residxs: + atom_ids.extend([at.index for at in r.atoms()]) + + return atom_ids + + +class AbsoluteBindingUnitMixin: + """ + Mixin for common class methods between Units + """ + def _get_alchemical_ion( + self, + alchem_comps: dict[str, list[Component]], + comp_resids: dict[Component, npt.NDArray], + omm_topology: omm_topology, + positions: ommunit.Quantity, + settings: dict[str, SettingsBaseModel], + ) -> int | None: + """ + Find a suitable alchemical ion for a net charge transformation. + + Parameters + ---------- + alchem_comps: dict[str, list[Component]] + A dictionary with a list of alchemical components + in both state A and B. + comp_resids: dict[Component, npt.NDArray] + A dictionary keyed by each Component in the System + which contains arrays with the residue indices that is contained + by that Component. + omm_topology : openmm.app.Topology + The OpenMM Topology of the system. + positions : openmm.unit.Quantity + The positions of the system. + settings : dict[str, SettingsBaseModel] + A dictionary of settings that defines how to find and set + the restraint. + + Returns + ------- + int + The index of the alchemical ion. + + """ + IONS = { + -1: ['Na', 'K', 'Li', 'Cs', 'Rb'], + 1: ['Cl', 'Br', 'F', 'I'] + } + + total_charge = alchem_comps["stateA"][0].total_charge + + # Don't add an alchemical ion if we have zero net charge + # or we didn't request it. + if total_charge == 0 or not settings['alchemical_settings'].explicit_charge_correction: + return None + + # For the ABFE protocol you'll never hit this, because we deal with it + # during validation. But just adding it in case it ends up in a subclass + if abs(total_charge) > 1: + errmsg = "Cannot handle net charge correction on charges greater than one" + raise ValueError(errmsg) + + univ = _get_mda_universe( + omm_topology, + positions, + self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, + ) + + counter_ions = [ + at.element in IONS[total_charge] + for at in univ.atoms + ] + + # get the alchemical atoms + residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]]) + alchem_idxs = _get_idxs_from_residxs(topology=omm_topology, residxs=residxs) + alchem_atoms = univ.atoms[alchem_idxs] + + # Re-using a utility from the restraints utilities + # TODO: rename this class! + atom_finder = FindHostAtoms( + host_atoms=counter_ions, + guest_atoms=alchem_atoms, + min_search_distance=settings['alchemical_settings'].alchemical_ion_min_distance, + max_search_distance=999 * offunit.nm # TODO: change to None? + ) + + atom_finder.run() + + # Just use the first one that comes back ok + return atom_finder.results.host_idxs[0] + + +class AbsoluteBindingComplexUnit(BaseAbsoluteUnit, AbsoluteBindingUnitMixin): """ Protocol Unit for the complex phase of an absolute binding free energy """ @@ -986,83 +1153,6 @@ def _handle_settings(self) -> dict[str, SettingsBaseModel]: return settings - @staticmethod - def _get_mda_universe( - topology: omm_topology, - positions: ommunit.Quantity, - trajectory: Optional[pathlib.Path], - ) -> mda.Universe: - """ - Helper method to get a Universe from an openmm Topology, - and either an input trajectory or a set of positions. - - Parameters - ---------- - topology : openmm.app.Topology - An OpenMM Topology that defines the System. - positions: openmm.unit.Quantity - The System's current positions. - Used if a trajectory file is None or is not a file. - trajectory: pathlib.Path - A Path to a trajectory file to read positions from. - - Returns - ------- - mda.Universe - An MDAnalysis Universe of the System. - """ - from MDAnalysis.coordinates.memory import MemoryReader - - # If the trajectory file doesn't exist, then we use positions - if trajectory is not None and trajectory.is_file(): - return mda.Universe( - topology, - trajectory, - topology_format="OPENMMTOPOLOGY", - ) - else: - # Positions is an openmm Quantity in nm we need - # to convert to angstroms - return mda.Universe( - topology, - np.array(positions._value) * 10, - topology_format="OPENMMTOPOLOGY", - trajectory_format=MemoryReader, - ) - - @staticmethod - def _get_idxs_from_residxs( - topology: omm_topology, - residxs: list[int], - ) -> list[int]: - """ - Helper method to get the a list of atom indices which belong to a list - of residues. - - Parameters - ---------- - topology : openmm.app.Topology - An OpenMM Topology that defines the System. - residxs : list[int] - A list of residue numbers who's atoms we should get atom indices. - - Returns - ------- - atom_ids : list[int] - A list of atom indices. - - TODO - ---- - * Check how this works when we deal with virtual sites. - """ - atom_ids = [] - - for r in topology.residues(): - if r.index in residxs: - atom_ids.extend([at.index for at in r.atoms()]) - - return atom_ids - @staticmethod def _get_boresch_restraint( universe: mda.Universe, @@ -1127,6 +1217,7 @@ def _add_restraints( alchem_comps: dict[str, list[Component]], comp_resids: dict[Component, npt.NDArray], settings: dict[str, SettingsBaseModel], + alchem_ion: int | None, ) -> tuple[ GlobalParameterState, Quantity, @@ -1159,6 +1250,8 @@ def _add_restraints( settings : dict[str, SettingsBaseModel] A dictionary of settings that defines how to find and set the restraint. + alchem_ions : int | None + The index of an alchemical ion, if it exists. Returns ------- @@ -1190,7 +1283,7 @@ def _add_restraints( residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]]) # get the alchemicical atom ids - guest_atom_ids = self._get_idxs_from_residxs(topology, residxs) + guest_atom_ids = _get_idxs_from_residxs(topology, residxs) # Now get the host idxs # We assume this is everything but the alchemical component @@ -1199,13 +1292,13 @@ def _add_restraints( exclude_comps = [alchem_comps["stateA"]] + solv_comps residxs = np.concatenate([v for i, v in comp_resids.items() if i not in exclude_comps]) - host_atom_ids = self._get_idxs_from_residxs(topology, residxs) + host_atom_ids = _get_idxs_from_residxs(topology, residxs) # Finally create an MDAnalysis Universe # We try to pass the equilibration production file path through # In some cases (debugging / dry runs) this won't be available # so we'll default to using input positions. - univ = self._get_mda_universe( + univ = _get_mda_universe( topology, positions, self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, @@ -1259,7 +1352,7 @@ def _add_restraints( ) -class AbsoluteBindingSolventUnit(BaseAbsoluteUnit): +class AbsoluteBindingSolventUnit(BaseAbsoluteUnit, AbsoluteBindingUnitMixin): """ Protocol Unit for the solvent phase of an absolute binding free energy """ diff --git a/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index c86773be1..0ac644e2c 100644 --- a/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -376,7 +376,7 @@ def test_complex_dry_run(self, complex_units, settings, tmpdir): expected_indices = [ i + self.num_complex_atoms for i in range(self.num_solvent_atoms) ] - assert expected_indices == data["alchem_indices"] + assert expected_indices == data["alchem_indices"][0] # Check the non-alchemical system self._assert_expected_nonalchemical_forces(data["system"], settings) @@ -392,7 +392,7 @@ def test_complex_dry_run(self, complex_units, settings, tmpdir): assert pdb.n_atoms == self.num_all_not_water # Check energies - alchem_region = AlchemicalRegion(alchemical_atoms=data["alchem_indices"]) + alchem_region = AlchemicalRegion(alchemical_atoms=data["alchem_indices"][0]) self._test_energies( reference_system=data["system"], alchemical_system=data["alchem_system"], @@ -415,7 +415,7 @@ def test_solvent_dry_run(self, solvent_units, settings, tmpdir): # Check the alchemical indices expected_indices = [i for i in range(self.num_solvent_atoms)] - assert expected_indices == data["alchem_indices"] + assert expected_indices == data["alchem_indices"][0] # Check the non-alchemical system self._assert_expected_nonalchemical_forces(data["system"], settings) @@ -431,7 +431,7 @@ def test_solvent_dry_run(self, solvent_units, settings, tmpdir): assert pdb.n_atoms == self.num_solvent_atoms # Check energies - alchem_region = AlchemicalRegion(alchemical_atoms=data["alchem_indices"]) + alchem_region = AlchemicalRegion(alchemical_atoms=data["alchem_indices"][0]) self._test_energies( reference_system=data["system"], From e35894f15839d8b1542e3a608277a12044655b6a Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 24 Oct 2025 00:21:45 +0100 Subject: [PATCH 2/7] alter validation --- openfe/protocols/openmm_afe/equil_binding_afe_method.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 342707007..4c15dc9e5 100644 --- a/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -683,10 +683,11 @@ def _validate_endstates( raise ValueError(errmsg) # Check that the state A unique isn't charged - if diff[0][0].total_charge != 0: + if abs(diff[0][0].total_charge) > 1: errmsg = ( "Charged alchemical molecules are not currently " - "supported for solvation free energies. " + "supported for absolute binding free energies with " + "absolute net charge greater than 1. " f"Molecule total charge: {diff[0][0].total_charge}." ) raise ValueError(errmsg) From e5b884341cbcaf4b0168a70d40929dd43c3e9004 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 24 Oct 2025 00:31:30 +0100 Subject: [PATCH 3/7] Add validation test --- .../openmm_afe/equil_binding_afe_method.py | 6 +-- .../data/openmm_rfe/charged_benzenes.sdf | 39 +++++++++++++++++++ .../openmm_abfe/test_abfe_validation.py | 26 +++++++++++-- 3 files changed, 65 insertions(+), 6 deletions(-) diff --git a/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 4c15dc9e5..9de518e11 100644 --- a/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -685,9 +685,9 @@ def _validate_endstates( # Check that the state A unique isn't charged if abs(diff[0][0].total_charge) > 1: errmsg = ( - "Charged alchemical molecules are not currently " - "supported for absolute binding free energies with " - "absolute net charge greater than 1. " + "Charged alchemical molecules with absolute net " + "charge greater than 1 are not currently " + "supported with absolute binding free energies. " f"Molecule total charge: {diff[0][0].total_charge}." ) raise ValueError(errmsg) diff --git a/openfe/tests/data/openmm_rfe/charged_benzenes.sdf b/openfe/tests/data/openmm_rfe/charged_benzenes.sdf index ea2735835..6ceb3ee02 100644 --- a/openfe/tests/data/openmm_rfe/charged_benzenes.sdf +++ b/openfe/tests/data/openmm_rfe/charged_benzenes.sdf @@ -98,3 +98,42 @@ benzoic_acid 9 14 1 0 0 0 0 M END $$$$ +phthalic_acid + PyMOL3.1 3D 0 + + 16 16 0 0 0 0 0 0 0 0999 V2000 + 1.4045 -0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.7445 -0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.2607 -1.0833 0.0000 O 0 5 0 0 0 0 0 0 0 0 0 0 + 0.7022 1.2164 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7023 1.2164 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.4265 1.0159 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.4045 -0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7023 -1.2164 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.7023 -1.2164 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.2540 2.1720 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.5079 -0.0000 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.2540 -2.1719 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.2540 -2.1720 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.3723 2.3768 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.5685 2.2821 0.0000 O 0 5 0 0 0 0 0 0 0 0 0 0 + 0.8335 3.4754 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 0 0 0 + 1 4 2 0 0 0 0 + 1 9 1 0 0 0 0 + 2 3 1 0 0 0 0 + 2 6 2 0 0 0 0 + 4 5 1 0 0 0 0 + 5 7 2 0 0 0 0 + 5 10 1 0 0 0 0 + 7 8 1 0 0 0 0 + 7 11 1 0 0 0 0 + 8 9 2 0 0 0 0 + 8 12 1 0 0 0 0 + 9 13 1 0 0 0 0 + 14 15 1 0 0 0 0 + 14 16 2 0 0 0 0 + 4 14 1 0 0 0 0 +M END +$$$$ + diff --git a/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py b/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py index 47f1502b1..2e351120d 100644 --- a/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py +++ b/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py @@ -243,7 +243,7 @@ def test_validate_endstates_unique_stateB(benzene_modifications, T4_protein_comp AbsoluteBindingProtocol._validate_endstates(stateA, stateB) -def test_charged_endstate(charged_benzene_modifications, T4_protein_component): +def test_charged_endstate_no_error(charged_benzene_modifications, T4_protein_component): stateA = ChemicalSystem( { @@ -260,8 +260,28 @@ def test_charged_endstate(charged_benzene_modifications, T4_protein_component): } ) - errmsg = "Charged alchemical molecules are not currently supported" - with pytest.raises(ValueError, match=errmsg): + # No error should be raised + AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + + +def test_charged_endstate_error(charged_benzene_modifications, T4_protein_component): + + stateA = ChemicalSystem( + { + "benzene": charged_benzene_modifications["phthalic_acid"], + "protein": T4_protein_component, + "solvent": SolventComponent(), + } + ) + + stateB = ChemicalSystem( + { + "protein": T4_protein_component, + "solvent": SolventComponent(), + } + ) + + with pytest.raises(ValueError, match="net charge greater than 1"): AbsoluteBindingProtocol._validate_endstates(stateA, stateB) From cef9cb750a6eedcd79507165222449bf8c30a892 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 24 Oct 2025 00:46:17 +0100 Subject: [PATCH 4/7] Add a warning if you're trying to do a net charge transformation --- .../openmm_afe/equil_binding_afe_method.py | 21 ++++++++++++------- .../openmm_abfe/test_abfe_validation.py | 9 +++++--- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 9de518e11..2e4c5b1c4 100644 --- a/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -683,14 +683,19 @@ def _validate_endstates( raise ValueError(errmsg) # Check that the state A unique isn't charged - if abs(diff[0][0].total_charge) > 1: - errmsg = ( - "Charged alchemical molecules with absolute net " - "charge greater than 1 are not currently " - "supported with absolute binding free energies. " - f"Molecule total charge: {diff[0][0].total_charge}." - ) - raise ValueError(errmsg) + if abs(diff[0][0].total_charge) >= 1: + + if abs(diff[0][0].total_charge) > 1: + errmsg = ( + "Charged alchemical molecules with absolute net " + "charge greater than 1 are not currently " + "supported with absolute binding free energies. " + f"Molecule total charge: {diff[0][0].total_charge}." + ) + raise ValueError(errmsg) + + wmsg = "Net charge transformations are experimental and error prone." + warnings.warn(wmsg) # If there are any alchemical Components in state B if len(diff[1]) > 0: diff --git a/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py b/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py index 2e351120d..f3f452374 100644 --- a/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py +++ b/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py @@ -243,7 +243,7 @@ def test_validate_endstates_unique_stateB(benzene_modifications, T4_protein_comp AbsoluteBindingProtocol._validate_endstates(stateA, stateB) -def test_charged_endstate_no_error(charged_benzene_modifications, T4_protein_component): +def test_charged_endstate_warn(charged_benzene_modifications, T4_protein_component): stateA = ChemicalSystem( { @@ -260,8 +260,11 @@ def test_charged_endstate_no_error(charged_benzene_modifications, T4_protein_com } ) - # No error should be raised - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + # No error should be raised but an warning should + wmsg = "Net charge transformations are experimental and error prone." + + with pytest.warns(UserWarning, match=wmsg): + AbsoluteBindingProtocol._validate_endstates(stateA, stateB) def test_charged_endstate_error(charged_benzene_modifications, T4_protein_component): From 345da59ebaa3d92e8930dbb0bafedf1654f74587 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 24 Oct 2025 01:54:25 +0100 Subject: [PATCH 5/7] Add solvent restraint --- openfe/protocols/openmm_afe/base.py | 2 +- .../openmm_afe/equil_afe_settings.py | 10 +- .../openmm_afe/equil_binding_afe_method.py | 123 +++++++++++++++++- 3 files changed, 128 insertions(+), 7 deletions(-) diff --git a/openfe/protocols/openmm_afe/base.py b/openfe/protocols/openmm_afe/base.py index 6bd650641..712214b31 100644 --- a/openfe/protocols/openmm_afe/base.py +++ b/openfe/protocols/openmm_afe/base.py @@ -623,7 +623,7 @@ def _add_restraints( ) -> tuple[ Optional[GlobalParameterState], Optional[Quantity], - Optional[openmm.System], + openmm.System, Optional[geometry.BaseRestraintGeometry], ]: """ diff --git a/openfe/protocols/openmm_afe/equil_afe_settings.py b/openfe/protocols/openmm_afe/equil_afe_settings.py index e0c472aec..a8e664af2 100644 --- a/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -37,6 +37,7 @@ from openfe.protocols.restraint_utils.settings import ( BaseRestraintSettings, BoreschRestraintSettings, + SpringConstantLinearQuantity, ) import numpy as np @@ -60,10 +61,15 @@ class ABFEAlchemicalSettings(AlchemicalSettings): Whether or not to use explicit charge correction using a co-alchemical ion. """ - alchemical_ion_min_distance: NanometerQuantity = 1.2 * offunit.nanometer + alchemical_ion_min_distance: NanometerQuantity = 1.0 * offunit.nanometer """ The minimum distance to search for a co-alchemical ion. """ + alchemical_ion_solvent_spring_constant: SpringConstantLinearQuantity = 1000.0 * offunit.kilojoule_per_mole / offunit.nm**2 + """ + The spring constant holding the ion away from the alchemical solute + in the solvent leg. + """ class LambdaSettings(SettingsBaseModel): @@ -339,7 +345,7 @@ def must_be_positive(cls, v): """Settings for solvating the system in the complex.""" # Alchemical settings - alchemical_settings: AlchemicalSettings + alchemical_settings: ABFEAlchemicalSettings """ Alchemical protocol settings. """ diff --git a/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 2e4c5b1c4..bcaec26cf 100644 --- a/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -34,6 +34,7 @@ import gufe import MDAnalysis as mda +from MDAnalysis.lib.distances import calc_bonds import numpy as np import numpy.typing as npt from gufe import ( @@ -63,13 +64,14 @@ from openfe.protocols.openmm_utils import settings_validation, system_validation from openfe.protocols.restraint_utils import geometry from openfe.protocols.restraint_utils.geometry.boresch import BoreschRestraintGeometry -from openfe.protocols.restraint_utils.geometry.utils import FindHostAtoms +from openfe.protocols.restraint_utils.geometry.utils import FindHostAtoms, get_central_atom_idx from openfe.protocols.restraint_utils.openmm import omm_restraints from openfe.protocols.restraint_utils.openmm.omm_restraints import BoreschRestraint +from openfe.protocols.restraint_utils.openmm.omm_forces import add_force_in_separate_group from openff.units import unit as offunit from openff.units import Quantity from openff.units.openmm import to_openmm -from openmm import System +from openmm import System, HarmonicBondForce from openmm import unit as ommunit from openmm.app import Topology as omm_topology from openmmtools import multistate @@ -1075,10 +1077,17 @@ def _get_alchemical_ion( host_atoms=counter_ions, guest_atoms=alchem_atoms, min_search_distance=settings['alchemical_settings'].alchemical_ion_min_distance, - max_search_distance=999 * offunit.nm # TODO: change to None? + # set max distance to just above solvent padding to avoid picking + # an ion more than half a box distance away + max_search_distance=settings['solvation_settings'].solvent_padding + 0.1 * offunit.nanometer, ) - atom_finder.run() + # only run on the final frame + atom_finder.run(frames=[-1]) + + if len(atom_finder.results.host_idxs) == 0: + errmsg = "No suitable alchemical ion was found" + raise ValueError(errmsg) # Just use the first one that comes back ok return atom_finder.results.host_idxs[0] @@ -1430,3 +1439,109 @@ def _handle_settings(self) -> dict[str, SettingsBaseModel]: settings["output_settings"] = prot_settings.solvent_output_settings return settings + + def _add_restraints( + self, + system: System, + topology: omm_topology, + positions: ommunit.Quantity, + alchem_comps: dict[str, list[Component]], + comp_resids: dict[Component, npt.NDArray], + settings: dict[str, SettingsBaseModel], + alchem_ion: int | None, + ) -> tuple[ + GlobalParameterState | None, + Quantity | None, + System, + geometry.HostGuestRestraintGeometry | None, + ]: + """ + Find and add restraints to the OpenMM System. + + Notes + ----- + Currently, only Boresch-like restraints are supported. + + Parameters + ---------- + system : openmm.System + The System to add the restraint to. + topology : openmm.app.Topology + An OpenMM Topology that defines the System. + positions: openmm.unit.Quantity + The System's current positions. + Used if a trajectory file isn't found. + alchem_comps: dict[str, list[Component]] + A dictionary with a list of alchemical components + in both state A and B. + comp_resids: dict[Component, npt.NDArray] + A dictionary keyed by each Component in the System + which contains arrays with the residue indices that is contained + by that Component. + settings : dict[str, SettingsBaseModel] + A dictionary of settings that defines how to find and set + the restraint. + alchem_ions : int | None + The index of an alchemical ion, if it exists. + + Returns + ------- + restraint_parameter_state : RestraintParameterState + A RestraintParameterState object that defines the control + parameter for the restraint. + correction : openff.units.Quantity + The standard state correction for the restraint. + system : openmm.System + A copy of the System with the restraint added. + rest_geom : geometry.HostGuestRestraintGeometry + The restraint Geometry object. + """ + if alchem_ion is None: + return None, None, system, None + + if self.verbose: + self.logger.info("Generating restraints") + + univ = _get_mda_universe( + topology, + positions, + self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, + ) + + # alchemical ion atom + alchem_ion = univ.atoms[alchem_ion] + + # get the alchemical ligand atoms + ligand_rdmol = alchem_comps["stateA"][0].to_rdkit() + residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]]) + ligand_alchem_idxs = _get_idxs_from_residxs(topology=omm_topology, residxs=residxs) + ligand_central_atom_idx = ligand_alchem_idxs[get_central_atom_idx(ligand_rdmol)] + ligand_central_atom = univ.atoms[ligand_central_atom_idx] + + # go to the final frame + univ.trajectory[-1] + + distance = float( + calc_bonds( + alchem_ion.position, + ligand_central_atom.position, + box=univ.dimensions + ) + ) + + spring_constant = to_openmm( + settings['alchemical_settings'].alchemical_ion_solvent_spring_constant + ) + + force = HarmonicBondForce() + force.addBond( + ligand_central_atom_idx, + alchem_ion, + distance * ommunit.angstrom, + spring_constant, + ) + + force.setName("ion_restraint") + add_force_in_separate_group(system, force) + + return None, None, system, None From 2631f721d0bad1d09f442cc9f269835fb8976369 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Sat, 25 Oct 2025 23:46:00 +0100 Subject: [PATCH 6/7] Add custom restraint defn --- .../openmm_afe/equil_afe_settings.py | 5 +++ .../openmm_afe/equil_binding_afe_method.py | 36 +++++++++++++------ .../openmm_abfe/test_abfe_protocol.py | 5 +++ 3 files changed, 36 insertions(+), 10 deletions(-) diff --git a/openfe/protocols/openmm_afe/equil_afe_settings.py b/openfe/protocols/openmm_afe/equil_afe_settings.py index a8e664af2..57cafd20e 100644 --- a/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -72,6 +72,11 @@ class ABFEAlchemicalSettings(AlchemicalSettings): """ +class ABFERestraintSettings(BoreschRestraintSettings): + host_restraint_ids: list[int] | None = None + guest_restraint_ids: list[int] | None = None + + class LambdaSettings(SettingsBaseModel): """Lambda schedule settings. diff --git a/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/openfe/protocols/openmm_afe/equil_binding_afe_method.py index bcaec26cf..3be0559c3 100644 --- a/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -50,6 +50,7 @@ AbsoluteBindingSettings, ABFEAlchemicalSettings, BoreschRestraintSettings, + ABFERestraintSettings, IntegratorSettings, LambdaSettings, ABFEPreEquilOutputSettings, @@ -596,7 +597,7 @@ def _default_settings(cls): solvent_solvation_settings=OpenMMSolvationSettings(), engine_settings=OpenMMEngineSettings(), integrator_settings=IntegratorSettings(), - restraint_settings=BoreschRestraintSettings(), + restraint_settings=ABFERestraintSettings(), solvent_equil_simulation_settings=MDSimulationSettings( equilibration_length_nvt=0.1 * offunit.nanosecond, equilibration_length=0.2 * offunit.nanosecond, @@ -660,8 +661,9 @@ def _validate_endstates( If the alchemical species is charged. """ if not (stateA.contains(ProteinComponent) and stateB.contains(ProteinComponent)): - errmsg = "No ProteinComponent found" - raise ValueError(errmsg) + if not (stateA.contains(SmallMoleculeComponent) and stateB.contains(SmallMoleculeComponent)): + errmsg = "No suitable host component found" + raise ValueError(errmsg) if not (stateA.contains(SolventComponent) and stateB.contains(SolventComponent)): errmsg = "No SolventComponent found" @@ -1176,6 +1178,8 @@ def _get_boresch_restraint( host_atom_ids: list[int], temperature: Quantity, settings: BoreschRestraintSettings, + guest_restraint_idxs: list[int] | None, + host_restraint_idxs: list[int] | None, ) -> tuple[BoreschRestraintGeometry, BoreschRestraint]: """ Get a Boresch-like restraint Geometry and OpenMM restraint force @@ -1195,6 +1199,14 @@ def _get_boresch_restraint( The temperature of the simulation where the restraint will be added. settings : BoreschRestraintSettings Settings on how the Boresch-like restraint should be defined. + guest_restraint_idxs : list[int] | None + A list indices that define the four guest restraint atoms. This is + a user override, if not None, these will be used instead of searching + for them. + host_restraint_idxs : list[int] | None + A list of indices that define the four host restraint atoms. This is + a user override, if not None, these will be used instead of searching + for them. Returns ------- @@ -1211,6 +1223,8 @@ def _get_boresch_restraint( guest_rdmol=guest_rdmol, guest_idxs=guest_atom_ids, host_idxs=host_atom_ids, + guest_restraint_atoms_idxs=guest_restraint_idxs, + host_restraint_atoms_idxs=host_restraint_idxs, host_selection=settings.host_selection, anchor_finding_strategy=settings.anchor_finding_strategy, dssp_filter=settings.dssp_filter, @@ -1319,14 +1333,16 @@ def _add_restraints( self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, ) - if isinstance(settings["restraint_settings"], BoreschRestraintSettings): + if isinstance(settings["restraint_settings"], ABFERestraintSettings): rest_geom, restraint = self._get_boresch_restraint( - univ, - guest_rdmol, - guest_atom_ids, - host_atom_ids, - settings["thermo_settings"].temperature, - settings["restraint_settings"], + universe=univ, + guest_rdmol=guest_rdmol, + guest_atom_ids=guest_atom_ids, + host_atom_ids=host_atom_ids, + temperature=settings["thermo_settings"].temperature, + settings=settings["restraint_settings"], + guest_restraint_idxs=settings["restraint_settings"].guest_restraint_ids, + host_restraint_idxs=settings["restraint_settings"].host_restraint_ids, ) else: # TODO turn this into a direction for different restraint types supported? diff --git a/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index 0ac644e2c..1cdf316b8 100644 --- a/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -540,5 +540,10 @@ def assign_fictitious_charges(offmol): ommunit.elementary_charge ) + # Alchemical sytstem should be 0 charge in the standard parameters + # and all charge in the offset + c, s, e = alchem_system_nbf.getParticleParameters(index) + assert pytest.approx(0.0) == c.value_in_unit(ommunit.elementary_charge) + offsets = alchem_system_nbf.getParticleParameterOffset(i) assert pytest.approx(prop_chgs[i]) == offsets[2] From d55d90cd8b9d19eff320c6796c23fb7060f67876 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Sun, 26 Oct 2025 20:23:15 +0000 Subject: [PATCH 7/7] Try to make this work properly --- openfe/protocols/openmm_afe/base.py | 92 +++- .../openmm_afe/equil_binding_afe_method.py | 31 +- .../tests/data/host_guest/OA/OA_charged.sdf | 422 ++++++++++++++++++ openfe/tests/data/host_guest/OA/__init__.py | 0 .../data/host_guest/OA/guests_charged.sdf | 414 +++++++++++++++++ openfe/tests/data/host_guest/__init__.py | 0 .../tests/protocols/openmm_abfe/conftest.py | 21 + .../openmm_abfe/test_abfe_protocol.py | 228 +++++++++- 8 files changed, 1172 insertions(+), 36 deletions(-) create mode 100644 openfe/tests/data/host_guest/OA/OA_charged.sdf create mode 100644 openfe/tests/data/host_guest/OA/__init__.py create mode 100644 openfe/tests/data/host_guest/OA/guests_charged.sdf create mode 100644 openfe/tests/data/host_guest/__init__.py diff --git a/openfe/protocols/openmm_afe/base.py b/openfe/protocols/openmm_afe/base.py index 712214b31..f0508336d 100644 --- a/openfe/protocols/openmm_afe/base.py +++ b/openfe/protocols/openmm_afe/base.py @@ -40,7 +40,6 @@ from openmmtools.alchemy import ( AlchemicalRegion, AbsoluteAlchemicalFactory, - AlchemicalState, ) from typing import Optional from openmm import app @@ -83,6 +82,50 @@ logger = logging.getLogger(__name__) +class SingleRegionAlchemicalState(GlobalParameterState): + class _LambdaParameter(GlobalParameterState.GlobalParameter): + """A global parameter in the interval [0, 1] with standard value 1.""" + + def __init__(self, parameter_name): + super().__init__(parameter_name, standard_value=1.0, validator=self.lambda_validator) + + @staticmethod + def lambda_validator(self, instance, parameter_value): + if parameter_value is None: + return parameter_value + if not (0.0 <= parameter_value <= 1.0): + raise ValueError('{} must be between 0 and 1.'.format(self.parameter_name)) + return float(parameter_value) + + lambda_sterics_A = _LambdaParameter("lambda_sterics_A") + lambda_electrostatics_A = _LambdaParameter("lambda_electrostatics_A") + lambda_bonds_A = _LambdaParameter("lambda_bonds_A") + lambda_angles_A = _LambdaParameter("lambda_angles_A") + lambda_torsions_A = _LambdaParameter("lambda_torsions_A") + + +class TwoRegionAlchemicalState(SingleRegionAlchemicalState): + class _LambdaParameter(GlobalParameterState.GlobalParameter): + """A global parameter in the interval [0, 1] with standard value 1.""" + + def __init__(self, parameter_name): + super().__init__(parameter_name, standard_value=1.0, validator=self.lambda_validator) + + @staticmethod + def lambda_validator(self, instance, parameter_value): + if parameter_value is None: + return parameter_value + if not (0.0 <= parameter_value <= 1.0): + raise ValueError('{} must be between 0 and 1.'.format(self.parameter_name)) + return float(parameter_value) + + lambda_sterics_B = _LambdaParameter("lambda_sterics_B") + lambda_electrostatics_B = _LambdaParameter("lambda_electrostatics_B") + lambda_bonds_B = _LambdaParameter("lambda_bonds_B") + lambda_angles_B = _LambdaParameter("lambda_angles_B") + lambda_torsions_B = _LambdaParameter("lambda_torsions_B") + + class BaseAbsoluteUnit(gufe.ProtocolUnit): """ Base class for ligand absolute free energy transformations. @@ -638,7 +681,7 @@ def _get_alchemical_system( comp_resids: dict[Component, npt.NDArray], alchem_comps: dict[str, list[Component]], alchem_ion: int | None, - ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, list[list[int]]]: + ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, dict[str, list[int]]]: """ Get an alchemically modified system and its associated factory @@ -659,7 +702,7 @@ def _get_alchemical_system( Factory for creating an alchemically modified system. alchemical_system : openmm.System Alchemically modified system - alchemical_indices : list[list[int]] + alchemical_indices : dict[str, list[int]] A list of atom indices for the alchemically modified species in the system. @@ -667,15 +710,22 @@ def _get_alchemical_system( ---- * Add support for all alchemical factory options """ - alchemical_indices = [self._get_alchemical_indices(topology, comp_resids, alchem_comps)] + alchemical_indices = { + 'A': self._get_alchemical_indices(topology, comp_resids, alchem_comps) + } if alchem_ion is not None: - alchemical_indices.append(alchem_ion) + alchemical_indices['B'] = [alchem_ion] alchemical_regions = [] - for inds in alchemical_indices: - alchemical_regions.append(AlchemicalRegion(alchemical_atoms=inds)) + for region in alchemical_indices: + alchemical_regions.append( + AlchemicalRegion( + alchemical_atoms=alchemical_indices[region], + name=region, + ) + ) alchemical_factory = AbsoluteAlchemicalFactory() alchemical_system = alchemical_factory.create_alchemical_system( @@ -694,6 +744,7 @@ def _get_states( lambdas: dict[str, list[float]], solvent_comp: Optional[SolventComponent], restraint_state: Optional[GlobalParameterState], + alchemical_indices: dict[str, list[int]], ) -> tuple[list[SamplerState], list[ThermodynamicState]]: """ Get a list of sampler and thermodynmic states from an @@ -715,6 +766,9 @@ def _get_states( The solvent component of the system, if there is one. restraint_state : Optional[GlobalParameterState] The restraint parameter control state, if there is one. + alchemical_indices : dict[str, list[int]] + Dictionary of the alchemical indices for each alchemical + region in the system. Returns ------- @@ -724,7 +778,17 @@ def _get_states( A list of ThermodynamicState for each replica in the system. """ # Fetch an alchemical state - alchemical_state = AlchemicalState.from_system(alchemical_system) + if len(alchemical_indices.keys()) == 1: + alchemical_state = SingleRegionAlchemicalState.from_system( + alchemical_system + ) + elif len(alchemical_indices.keys()) == 2: + alchemical_state = TwoRegionAlchemicalState.from_system( + alchemical_system + ) + else: + errmsg = "more than two regions are not supported" + raise ValueError(errmsg) # Set up the system constants temperature = settings["thermo_settings"].temperature @@ -736,7 +800,16 @@ def _get_states( constants["pressure"] = ensure_quantity(pressure, "openmm") # Get the thermodynamic parameter protocol - param_protocol = copy.deepcopy(lambdas) + param_protocol = {} + + def _add_lambdas_to_protocol(protocol, lambdas, region_name): + protocol[f"lambda_electrostatics_{region_name}"] = lambdas["lambda_electrostatics"] + protocol[f"lambda_sterics_{region_name}"] = lambdas["lambda_sterics"] + + param_protocol["lambda_restraints"] = lambdas["lambda_restraints"] + _add_lambdas_to_protocol(param_protocol, lambdas, "A") + if len(alchemical_indices.keys()) == 2: + _add_lambdas_to_protocol(param_protocol, lambdas, "B") # Get the composable states if restraint_state is not None: @@ -1193,6 +1266,7 @@ def run( lambdas, solv_comp, restraint_parameter_state, + alchem_indices, ) # 10. Create the multistate reporter & create PDB diff --git a/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 3be0559c3..5c3824a3b 100644 --- a/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -30,6 +30,7 @@ import uuid import warnings from collections import defaultdict +from copy import deepcopy from typing import Any, Iterable, Optional, Union import gufe @@ -1063,10 +1064,12 @@ def _get_alchemical_ion( self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, ) - counter_ions = [ - at.element in IONS[total_charge] + counter_ions_idxs = [ + at.ix for at in univ.atoms + if at.element in IONS[total_charge] ] + counter_ions = univ.atoms[counter_ions_idxs] # get the alchemical atoms residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]]) @@ -1095,7 +1098,7 @@ def _get_alchemical_ion( return atom_finder.results.host_idxs[0] -class AbsoluteBindingComplexUnit(BaseAbsoluteUnit, AbsoluteBindingUnitMixin): +class AbsoluteBindingComplexUnit(AbsoluteBindingUnitMixin, BaseAbsoluteUnit): """ Protocol Unit for the complex phase of an absolute binding free energy """ @@ -1383,7 +1386,7 @@ def _add_restraints( ) -class AbsoluteBindingSolventUnit(BaseAbsoluteUnit, AbsoluteBindingUnitMixin): +class AbsoluteBindingSolventUnit(AbsoluteBindingUnitMixin, BaseAbsoluteUnit): """ Protocol Unit for the solvent phase of an absolute binding free energy """ @@ -1515,6 +1518,8 @@ def _add_restraints( if alchem_ion is None: return None, None, system, None + restrained_system = deepcopy(system) + if self.verbose: self.logger.info("Generating restraints") @@ -1525,22 +1530,22 @@ def _add_restraints( ) # alchemical ion atom - alchem_ion = univ.atoms[alchem_ion] + alchem_ion_ag = univ.atoms[alchem_ion] # get the alchemical ligand atoms ligand_rdmol = alchem_comps["stateA"][0].to_rdkit() residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]]) - ligand_alchem_idxs = _get_idxs_from_residxs(topology=omm_topology, residxs=residxs) - ligand_central_atom_idx = ligand_alchem_idxs[get_central_atom_idx(ligand_rdmol)] - ligand_central_atom = univ.atoms[ligand_central_atom_idx] + ligand_alchem_idxs = _get_idxs_from_residxs(topology=topology, residxs=residxs) + ligand_central_atom = ligand_alchem_idxs[get_central_atom_idx(ligand_rdmol)] + ligand_central_atom_ag = univ.atoms[ligand_central_atom] # go to the final frame univ.trajectory[-1] distance = float( calc_bonds( - alchem_ion.position, - ligand_central_atom.position, + alchem_ion_ag.position, + ligand_central_atom_ag.position, box=univ.dimensions ) ) @@ -1551,13 +1556,13 @@ def _add_restraints( force = HarmonicBondForce() force.addBond( - ligand_central_atom_idx, + ligand_central_atom, alchem_ion, distance * ommunit.angstrom, spring_constant, ) force.setName("ion_restraint") - add_force_in_separate_group(system, force) + add_force_in_separate_group(restrained_system, force) - return None, None, system, None + return None, None, restrained_system, None diff --git a/openfe/tests/data/host_guest/OA/OA_charged.sdf b/openfe/tests/data/host_guest/OA/OA_charged.sdf new file mode 100644 index 000000000..c49fe096f --- /dev/null +++ b/openfe/tests/data/host_guest/OA/OA_charged.sdf @@ -0,0 +1,422 @@ + + RDKit 3D + +184204 0 0 0 0 0 0 0 0999 V2000 + 1.9780 -2.9130 2.5560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.7890 -0.8910 1.7980 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.0950 -1.7370 3.3190 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.8090 -3.1130 1.4400 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.7500 -2.1260 1.0940 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.8730 -0.6510 2.8620 C 0 0 0 0 0 0 0 0 0 0 0 0 + 4.5480 0.2040 1.3130 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.9840 0.5840 3.7170 C 0 0 1 0 0 0 0 0 0 0 0 0 + 3.9320 1.3630 0.6550 C 0 0 1 0 0 0 0 0 0 0 0 0 + 3.6580 2.3950 1.6720 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.3640 2.5270 2.1880 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.0150 1.6640 3.2190 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.7110 1.8220 3.7590 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.5080 3.4420 1.6300 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.1860 2.7330 3.1850 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2210 3.5300 2.1100 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6470 4.3030 1.4010 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.6530 2.7550 3.6810 C 0 0 1 0 0 0 0 0 0 0 0 0 + -2.5740 1.8700 2.8770 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.1320 2.4120 1.6520 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5380 3.6600 0.4250 C 0 0 1 0 0 0 0 0 0 0 0 0 + -2.7830 3.6460 1.1110 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.5440 -4.1280 0.4970 O 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8190 -3.8400 2.8300 C 0 0 2 0 0 0 0 0 0 0 0 0 + -0.4320 -3.4440 2.0490 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2820 -4.7530 -0.0200 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6540 -3.9910 0.7700 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.9750 1.6730 0.8700 C 0 0 0 0 0 0 0 0 0 0 0 0 + -4.2640 0.3670 1.1370 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8020 -3.7270 0.0330 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.7750 -2.9180 0.5540 C 0 0 0 0 0 0 0 0 0 0 0 0 + -4.8390 -0.4100 0.1310 O 0 0 0 0 0 0 0 0 0 0 0 0 + -3.8440 -2.5740 -0.2240 O 0 0 0 0 0 0 0 0 0 0 0 0 + -3.9300 -1.1640 -0.7330 C 0 0 1 0 0 0 0 0 0 0 0 0 + -1.4680 -2.7260 2.6070 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.6720 -2.4640 1.9030 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.9270 0.5100 3.1980 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.9200 -1.7400 2.6170 C 0 0 2 0 0 0 0 0 0 0 0 0 + -3.7670 -0.2670 2.3280 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.4480 -3.9480 -0.4130 C 0 0 2 0 0 0 0 0 0 0 0 0 + -1.7230 2.5680 5.1830 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.1160 2.7320 5.7950 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.0940 2.4680 7.3220 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.1470 1.9470 7.8190 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.9710 0.3330 5.2440 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.4740 1.5040 6.0680 C 0 0 0 0 0 0 0 0 0 0 0 0 + 4.9650 1.5110 5.9710 C 0 0 0 0 0 0 0 0 0 0 0 0 + 5.6550 0.8410 6.7490 O 0 0 0 0 0 0 0 0 0 0 0 0 + 0.6190 -4.2620 4.3570 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4420 -5.3440 4.4880 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2540 -6.6130 4.0330 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.0480 -7.1760 4.7900 O 0 0 0 0 0 0 0 0 0 0 0 0 + -3.9730 -2.2770 4.1800 C 0 0 0 0 0 0 0 0 0 0 0 0 + -5.1780 -1.5800 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + -5.1500 -1.9220 6.3420 C 0 0 0 0 0 0 0 0 0 0 0 0 + -4.2720 -1.6220 7.0980 O 0 0 0 0 0 0 0 0 0 0 0 0 + 4.8530 1.8910 -0.4090 C 0 0 0 0 0 0 0 0 0 0 0 0 + 6.4190 2.7220 -2.5260 C 0 0 0 0 0 0 0 0 0 0 0 0 + 4.8060 3.1880 -0.8380 C 0 0 0 0 0 0 0 0 0 0 0 0 + 5.7010 1.0210 -1.0690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 6.5160 1.4180 -2.1280 C 0 0 0 0 0 0 0 0 0 0 0 0 + 5.5390 3.5760 -1.8830 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.9240 -4.1690 -1.8430 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.5740 -4.1720 -4.5420 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.0450 -4.6830 -2.7060 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.1640 -3.6710 -2.3470 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.5120 -3.7000 -3.6730 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.3590 -4.6800 -4.0870 C 0 0 0 0 0 0 0 0 0 0 0 0 + -4.3640 -1.2630 -2.1710 C 0 0 0 0 0 0 0 0 0 0 0 0 + -5.1030 -1.2900 -4.8080 C 0 0 0 0 0 0 0 0 0 0 0 0 + -5.1920 -0.2470 -2.6580 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.9340 -2.3460 -3.0260 C 0 0 0 0 0 0 0 0 0 0 0 0 + -4.2790 -2.3050 -4.3980 C 0 0 0 0 0 0 0 0 0 0 0 0 + -5.5090 -0.2170 -4.0010 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5850 4.4050 -0.9270 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5140 5.2850 -3.5240 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.7380 4.3570 -1.7190 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.3570 4.8830 -1.4230 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.3650 5.2850 -2.7420 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.6800 4.7550 -2.9980 C 0 0 0 0 0 0 0 0 0 0 0 0 + 4.6560 -3.1320 -4.2310 O 0 0 0 0 0 0 0 0 0 0 0 0 + -3.9210 -3.3160 -5.2690 O 0 0 0 0 0 0 0 0 0 0 0 0 + -6.1530 0.7740 -4.5970 O 0 0 0 0 0 0 0 0 0 0 0 0 + 7.3710 0.5580 -2.7490 O 0 0 0 0 0 0 0 0 0 0 0 0 + 5.5640 4.8870 -2.3690 O 0 0 0 0 0 0 0 0 0 0 0 0 + -3.8000 4.8130 -3.8840 O 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8470 5.5530 -3.3740 O 0 0 0 0 0 0 0 0 0 0 0 0 + -3.1250 -4.3820 -4.8430 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.3880 -6.3310 -3.9960 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.7660 -4.2420 -5.0990 C 0 0 0 0 0 0 0 0 0 0 0 0 + -3.6200 -5.4880 -4.0900 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.7470 -6.4710 -3.6320 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8840 -5.2010 -4.6980 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4360 -5.0770 -4.9660 O 0 0 0 0 0 0 0 0 0 0 0 0 + 7.0140 -0.7940 -2.8120 C 0 0 0 0 0 0 0 0 0 0 0 0 + 6.6030 -3.4800 -2.7760 C 0 0 0 0 0 0 0 0 0 0 0 0 + 7.8680 -1.6190 -2.0560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 5.8570 -1.2990 -3.4370 C 0 0 0 0 0 0 0 0 0 0 0 0 + 5.6910 -2.6750 -3.4480 C 0 0 0 0 0 0 0 0 0 0 0 0 + 7.6440 -2.9510 -2.0130 C 0 0 0 0 0 0 0 0 0 0 0 0 + -6.1080 2.0390 -4.1700 C 0 0 0 0 0 0 0 0 0 0 0 0 + -6.0420 4.5660 -3.0460 C 0 0 0 0 0 0 0 0 0 0 0 0 + -7.2560 2.5070 -3.5050 C 0 0 0 0 0 0 0 0 0 0 0 0 + -4.9680 2.7640 -4.2470 C 0 0 0 0 0 0 0 0 0 0 0 0 + -4.9020 4.0350 -3.6650 C 0 0 0 0 0 0 0 0 0 0 0 0 + -7.2160 3.7920 -2.9240 C 0 0 0 0 0 0 0 0 0 0 0 0 + -8.3420 4.2640 -2.0140 C 0 0 0 0 0 0 0 0 0 0 0 0 + -9.2730 3.5770 -1.6770 O 0 0 0 0 0 0 0 0 0 0 0 0 + 3.5430 9.3010 -1.4370 C 0 0 0 0 0 0 0 0 0 0 0 0 + 4.6040 9.6750 -1.0440 O 0 0 0 0 0 0 0 0 0 0 0 0 + 8.5220 -3.8700 -1.0650 C 0 0 0 0 0 0 0 0 0 0 0 0 + 8.5050 -5.0240 -1.2770 O 0 0 0 0 0 0 0 0 0 0 0 0 + -4.0100 2.7470 7.9930 O 0 0 0 0 0 0 0 0 0 0 0 0 + -6.1860 -2.4440 6.7280 O 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0950 -6.9830 2.9000 O 0 0 0 0 0 0 0 0 0 0 0 0 + -3.1550 -7.5970 -2.7690 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.3410 -8.3000 -2.1910 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.5090 9.9980 -1.3550 O 0 0 0 0 0 0 0 0 0 0 0 0 + -8.3240 5.4290 -1.6400 O 0 0 0 0 0 0 0 0 0 0 0 0 + -4.3750 -7.7630 -2.5820 O 0 0 0 0 0 0 0 0 0 0 0 0 + 9.0540 -3.4610 -0.0610 O 0 0 0 0 0 0 0 0 0 0 0 0 + 5.5220 2.2790 5.1950 O 0 0 0 0 0 0 0 0 0 0 0 0 + 4.4110 5.6700 -2.4190 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.1750 7.2800 -2.2470 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.1540 5.1300 -2.8380 C 0 0 0 0 0 0 0 0 0 0 0 0 + 4.5390 7.0130 -1.9880 C 0 0 0 0 0 0 0 0 0 0 0 0 + 3.3890 7.8860 -1.8580 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.0370 5.9900 -2.7670 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.3140 -1.5380 4.0470 H 0 0 0 0 0 0 0 0 0 0 0 0 + 4.3950 -2.2610 0.2290 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.9410 1.0910 3.5770 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.9520 1.0690 0.2500 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4670 1.0980 4.5330 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.8430 4.0220 0.7740 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.9970 3.7840 3.5520 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.2860 2.6080 0.2230 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.1220 -4.7900 2.3850 H 0 0 0 0 0 0 0 0 0 0 0 0 + -4.3700 2.1250 -0.0360 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8910 -4.1970 -0.9440 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.8730 -0.8670 -0.6830 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5000 -2.4000 3.6440 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.5320 0.0400 4.0950 H 0 0 0 0 0 0 0 0 0 0 0 0 + -4.8830 -1.9290 2.1380 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.0410 -2.9290 -0.3430 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.2640 1.6140 5.4480 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.0850 3.3470 5.6060 H 0 0 0 0 0 0 0 0 0 0 0 0 + -3.4140 3.7700 5.6300 H 0 0 0 0 0 0 0 0 0 0 0 0 + -3.7740 1.9840 5.3470 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.4060 -0.6560 5.4020 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.9960 0.0970 5.6760 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.0110 2.4500 5.7790 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.2700 1.1450 7.0800 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.5540 -4.5780 4.8230 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2700 -3.3680 4.8770 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6790 -5.3400 5.5540 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.3680 -5.2640 3.9150 H 0 0 0 0 0 0 0 0 0 0 0 0 + -4.1990 -3.3400 4.2770 H 0 0 0 0 0 0 0 0 0 0 0 0 + -3.1350 -2.0140 4.8300 H 0 0 0 0 0 0 0 0 0 0 0 0 + -5.0920 -0.4960 4.7620 H 0 0 0 0 0 0 0 0 0 0 0 0 + -6.1050 -1.8690 4.3570 H 0 0 0 0 0 0 0 0 0 0 0 0 + 6.9260 3.1600 -3.3810 H 0 0 0 0 0 0 0 0 0 0 0 0 + 4.0220 3.8410 -0.4650 H 0 0 0 0 0 0 0 0 0 0 0 0 + 5.8180 0.0340 -0.6280 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.7470 -4.1360 -5.6140 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0270 -4.9510 -2.4310 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.8240 -3.2270 -1.6060 H 0 0 0 0 0 0 0 0 0 0 0 0 + -5.3290 -1.2620 -5.8720 H 0 0 0 0 0 0 0 0 0 0 0 0 + -5.4360 0.6030 -2.0260 H 0 0 0 0 0 0 0 0 0 0 0 0 + -3.2850 -3.1120 -2.6070 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5800 5.8250 -4.4650 H 0 0 0 0 0 0 0 0 0 0 0 0 + -3.6250 3.8900 -1.2980 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.5490 4.8370 -0.8240 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6830 -7.0780 -3.6390 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.4480 -3.3280 -5.5950 H 0 0 0 0 0 0 0 0 0 0 0 0 + -4.6980 -5.5820 -3.9820 H 0 0 0 0 0 0 0 0 0 0 0 0 + 6.2860 -4.5190 -2.7310 H 0 0 0 0 0 0 0 0 0 0 0 0 + 8.6710 -1.1640 -1.4820 H 0 0 0 0 0 0 0 0 0 0 0 0 + 5.0750 -0.6430 -3.8120 H 0 0 0 0 0 0 0 0 0 0 0 0 + -6.1300 5.5770 -2.6560 H 0 0 0 0 0 0 0 0 0 0 0 0 + -8.1780 1.9330 -3.5290 H 0 0 0 0 0 0 0 0 0 0 0 0 + -4.1080 2.2910 -4.7150 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.2730 7.8750 -2.3700 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.0630 4.0900 -3.1410 H 0 0 0 0 0 0 0 0 0 0 0 0 + 5.4990 7.4650 -1.7500 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 3 2 0 + 1 4 1 0 + 4 5 2 0 + 2 5 1 0 + 3 6 1 0 + 2 6 2 0 + 2 7 1 0 + 6 8 1 0 + 7 9 1 0 + 9 10 1 0 + 10 11 1 0 + 11 12 2 0 + 8 12 1 0 + 12 13 1 0 + 11 14 1 0 + 13 15 2 0 + 15 16 1 0 + 14 16 2 0 + 16 17 1 0 + 15 18 1 0 + 18 19 1 0 + 19 20 2 0 + 17 21 1 0 + 21 22 1 0 + 20 22 1 0 + 4 23 1 0 + 1 24 1 0 + 24 25 1 0 + 26 27 1 0 + 25 27 2 0 + 20 28 1 0 + 28 29 2 0 + 27 30 1 0 + 30 31 2 0 + 29 32 1 0 + 31 33 1 0 + 32 34 1 0 + 33 34 1 0 + 25 35 1 0 + 35 36 2 0 + 31 36 1 0 + 19 37 1 0 + 36 38 1 0 + 38 39 1 0 + 37 39 2 0 + 29 39 1 0 + 26 40 1 0 + 23 40 1 0 + 18 41 1 0 + 41 42 1 0 + 42 43 1 0 + 43 44 1 0 + 8 45 1 0 + 45 46 1 0 + 46 47 1 0 + 47 48 1 0 + 24 49 1 0 + 49 50 1 0 + 50 51 1 0 + 51 52 1 0 + 38 53 1 0 + 53 54 1 0 + 54 55 1 0 + 55 56 2 0 + 9 57 1 0 + 57 59 2 0 + 57 60 1 0 + 60 61 2 0 + 58 61 1 0 + 59 62 1 0 + 58 62 2 0 + 40 63 1 0 + 63 65 2 0 + 63 66 1 0 + 64 67 1 0 + 66 67 2 0 + 65 68 1 0 + 64 68 2 0 + 34 69 1 0 + 69 71 2 0 + 69 72 1 0 + 70 73 1 0 + 72 73 2 0 + 70 74 2 0 + 71 74 1 0 + 21 75 1 0 + 75 77 2 0 + 75 78 1 0 + 78 79 2 0 + 76 79 1 0 + 77 80 1 0 + 76 80 2 0 + 67 81 1 0 + 73 82 1 0 + 74 83 1 0 + 61 84 1 0 + 62 85 1 0 + 80 86 1 0 + 79 87 1 0 + 82 88 1 0 + 88 90 2 0 + 88 91 1 0 + 91 92 2 0 + 89 92 1 0 + 89 93 2 0 + 90 93 1 0 + 93 94 1 0 + 68 94 1 0 + 84 95 1 0 + 95 97 2 0 + 95 98 1 0 + 98 99 2 0 + 96 99 1 0 + 81 99 1 0 + 97100 1 0 + 96100 2 0 + 83101 1 0 +101103 2 0 +101104 1 0 +104105 2 0 +102105 1 0 + 86105 1 0 +103106 1 0 +102106 2 0 +106107 1 0 +107108 2 0 +109110 2 0 +100111 1 0 +111112 2 0 + 43113 2 0 + 55114 1 0 + 51115 2 0 + 92116 1 0 +116117 2 0 +109118 1 0 +107119 1 0 +116120 1 0 +111121 1 0 + 47122 2 0 + 85123 1 0 +123125 2 0 +123126 1 0 +124127 1 0 +109127 1 0 +126127 2 0 +125128 1 0 +124128 2 0 + 87128 1 0 + 3129 1 0 + 5130 1 0 + 8131 1 1 + 9132 1 1 + 13133 1 0 + 14134 1 0 + 18135 1 6 + 21136 1 1 + 24137 1 6 + 28138 1 0 + 30139 1 0 + 34140 1 1 + 35141 1 0 + 37142 1 0 + 38143 1 6 + 40144 1 1 + 41145 1 0 + 41146 1 0 + 42147 1 0 + 42148 1 0 + 45149 1 0 + 45150 1 0 + 46151 1 0 + 46152 1 0 + 49153 1 0 + 49154 1 0 + 50155 1 0 + 50156 1 0 + 53157 1 0 + 53158 1 0 + 54159 1 0 + 54160 1 0 + 58161 1 0 + 59162 1 0 + 60163 1 0 + 64164 1 0 + 65165 1 0 + 66166 1 0 + 70167 1 0 + 71168 1 0 + 72169 1 0 + 76170 1 0 + 77171 1 0 + 78172 1 0 + 89173 1 0 + 90174 1 0 + 91175 1 0 + 96176 1 0 + 97177 1 0 + 98178 1 0 +102179 1 0 +103180 1 0 +104181 1 0 +124182 1 0 +125183 1 0 +126184 1 0 +M CHG 8 44 -1 48 -1 52 -1 114 -1 118 -1 119 -1 120 -1 121 -1 +M END + +> + + +> +-0.091408472103269203 0.10657664056381454 -0.067896451396138771 0.10657664056381454 -0.16556027118602526 -0.091408472103269203 -0.3436575490657402 -0.02315850887933503 0.31904197152218094 +-0.3436575490657402 0.10657664056381454 -0.091408531707913979 -0.067896451396138771 -0.16556027118602526 -0.091408472103269203 0.10657664056381454 -0.3436575490657402 -0.02315850887933503 +-0.091408531707913979 0.10657664056381454 0.31904197152218094 -0.3436575490657402 -0.3436575490657402 -0.02315850887933503 -0.091408531707913979 -0.3436575490657402 0.10657664056381454 +-0.16556027118602526 0.10657664056381454 -0.16556027118602526 0.10657664056381454 -0.3436575490657402 -0.3436575490657402 0.31904197152218094 -0.067896451396138771 -0.091408472103269203 +-0.067896451396138771 -0.02315850887933503 -0.091408472103269203 0.31904197152218094 -0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777 +-0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777 -0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777 +-0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777 -0.091567691010625465 -0.16220025305190813 -0.14929413382449877 -0.14929413382449877 0.1104655456121849 +0.1104655456121849 -0.091567691010625465 -0.16220025305190813 -0.14929410402217638 -0.14929413382449877 0.1104655456121849 0.1104655456121849 -0.091567691010625465 -0.16220025305190813 +-0.14929413382449877 -0.14929413382449877 0.1104655456121849 0.1104655456121849 -0.091567691010625465 -0.16220025305190813 -0.14929413382449877 -0.14929413382449877 0.1104655456121849 +0.1104655456121849 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 0.085777629571764366 +-0.11220825795570145 -0.18432724062839281 -0.11220825795570145 -0.082448165458829506 0.085777629571764366 -0.26024511043468246 0.085777629571764366 -0.11220825795570145 +-0.11220825795570145 -0.18432724062839281 0.085777629571764366 -0.082448165458829506 0.085777629571764366 -0.11220825795570145 -0.11220825795570145 -0.18432724062839281 +0.085777629571764366 -0.082448165458829506 0.89783174212536088 -0.81537824456134567 0.89783174212536088 -0.81537824456134567 0.89783174212536088 -0.81537824456134567 -0.84737205092349777 +-0.84737205092349777 -0.84737205092349777 0.89783174212536088 -0.81537824456134567 -0.81537824456134567 -0.81537824456134567 -0.81537824456134567 -0.81537824456134567 -0.84737205092349777 +0.085777629571764366 -0.11220825795570145 -0.1843272704307152 -0.11220825795570145 -0.082448165458829506 0.085777629571764366 0.12293871160110702 0.1574754160221504 0.077845622258989708 +0.048985436963646309 0.12293871160110702 0.1574754160221504 0.077845622258989708 0.048985436963646309 0.077845622258989708 0.1574754160221504 0.1574754160221504 0.048985436963646309 +0.12293871160110702 0.12293871160110702 0.077845622258989708 0.048985436963646309 0.066841882127134697 0.066841882127134697 0.022772260415165321 0.022772260415165321 0.066841882127134697 +0.066841882127134697 0.022772260415165321 0.022772260415165321 0.066841882127134697 0.066841882127134697 0.022772260415165321 0.022772260415165321 0.066841882127134697 0.066841882127134697 +0.022772260415165321 0.022772260415165321 0.15592389221748579 0.14880136068424452 0.14880136068424452 0.15592389221748579 0.14880136068424452 0.14880136068424452 0.15592389221748579 +0.14880136068424452 0.14880136068424452 0.15592389221748579 0.14880136068424452 0.14880136068424452 0.1685743671234535 0.13864818628391493 0.16857433732113111 0.16857433732113111 +0.16857433732113111 0.13864818628391493 0.16857433732113111 0.16857433732113111 0.13864818628391493 0.16857433732113111 0.13864812667927015 0.16857439692577589 + +$$$$ diff --git a/openfe/tests/data/host_guest/OA/__init__.py b/openfe/tests/data/host_guest/OA/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/openfe/tests/data/host_guest/OA/guests_charged.sdf b/openfe/tests/data/host_guest/OA/guests_charged.sdf new file mode 100644 index 000000000..b1f356172 --- /dev/null +++ b/openfe/tests/data/host_guest/OA/guests_charged.sdf @@ -0,0 +1,414 @@ +OA-G0 + RDKit 3D + + 20 20 0 0 0 0 0 0 0 0999 V2000 + -0.7955 1.2297 -3.7854 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.3091 -0.3847 -0.7486 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.5096 0.2717 0.3768 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2838 -0.8860 -1.7613 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.9597 0.1420 -0.0138 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.9546 -0.0222 -1.5321 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8697 1.3180 -2.2595 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2049 0.6341 -4.2709 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.7693 1.7767 -4.3800 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.0136 0.3154 -1.2099 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.8934 -1.2223 -0.3497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.7004 -0.2524 1.3210 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.7925 1.3196 0.5218 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.6733 -0.7466 -2.7812 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0732 -1.9566 -1.6356 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.4456 -0.7024 0.4948 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5023 1.0547 0.2785 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8449 -0.5532 -1.8532 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0170 1.8572 -1.9078 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.7461 1.9209 -1.9970 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 7 1 0 + 1 8 1 0 + 1 9 2 0 + 2 3 1 0 + 2 4 1 0 + 3 5 1 0 + 4 6 1 0 + 5 6 1 0 + 6 7 1 0 + 2 10 1 0 + 2 11 1 0 + 3 12 1 0 + 3 13 1 0 + 4 14 1 0 + 4 15 1 0 + 5 16 1 0 + 5 17 1 0 + 6 18 1 0 + 7 19 1 0 + 7 20 1 0 +M CHG 1 8 -1 +M END + +> +-7.9565659999999996 + +> +OA-G0 + +> +0.90913933776319023 -0.080904421582818034 -0.080904421582818034 -0.08442279435694218 -0.08442279435694218 -0.062293932959437373 -0.19467459358274936 -0.846026298776269 -0.846026298776269 +0.025389452278614045 0.025389452278614045 0.025389452278614045 0.025389452278614045 0.043140637502074239 0.043140637502074239 0.043140637502074239 0.043140637502074239 0.053682352229952809 +0.021366753429174424 0.021366753429174424 + +$$$$ +OA-G1 + RDKit 3D + + 17 16 0 0 0 0 0 0 0 0999 V2000 + -1.0006 -0.7469 -3.8833 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.0084 -0.5611 -3.0020 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8433 -0.3922 -5.3688 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.0381 0.0211 0.8063 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.1259 -0.8917 -1.5409 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.1069 0.3450 -0.6725 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2681 0.0857 -5.7303 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8736 -0.6292 -6.0659 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.9447 -1.1634 -3.5343 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.9430 -0.1434 -3.3206 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.0312 -0.3847 1.0235 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0983 0.9255 1.4073 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.7104 -0.7126 1.1224 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.1300 -1.2920 -1.3561 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.5854 -1.6783 -1.2604 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.1106 0.7483 -0.8545 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6092 1.1326 -0.9379 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 2 0 + 1 3 1 0 + 2 5 1 0 + 3 7 1 0 + 3 8 2 0 + 4 6 1 0 + 5 6 1 0 + 1 9 1 0 + 2 10 1 0 + 4 11 1 0 + 4 12 1 0 + 4 13 1 0 + 5 14 1 0 + 5 15 1 0 + 6 16 1 0 + 6 17 1 0 +M CHG 1 7 -1 +M END + +> +-6.8312020000000002 + +> +OA-G1 + +> +-0.17128288493875196 -0.23447102056268385 0.90111100925680465 -0.10154167310718228 -0.040073477937018168 -0.082682838006054651 -0.84050458417657548 -0.84050458417657548 0.07679569973226856 +0.11721654488321613 0.027137529762352213 0.027137529762352213 0.027137529762352213 0.031102622585261568 0.031102622585261568 0.036159987287486303 0.036159987287486303 + +$$$$ +OA-G2 + RDKit 3D + + 25 25 0 0 0 0 0 0 0 0999 V2000 + 0.2241 -1.1357 -4.2804 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4652 0.1472 -4.5856 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.7894 -0.5542 -1.5717 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.1825 0.4275 -5.8319 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.5566 -0.2851 -1.1171 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4939 -1.5520 -3.0259 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0275 1.3042 -3.7403 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.0051 0.9021 -2.6793 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6591 -0.4247 -1.9881 C 0 0 2 0 0 0 0 0 0 0 0 0 + 0.3129 0.1732 0.2902 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.6132 0.0157 -6.8813 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.2758 1.0443 -5.6932 O 0 0 0 0 0 0 0 0 0 0 0 0 + 0.5476 -1.9343 -4.9434 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.6613 -0.4502 -0.9346 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.9554 -0.8848 -2.5917 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.4888 -1.9348 -3.3003 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0542 -2.4021 -2.5903 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8919 1.7520 -3.2369 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4008 2.0735 -4.3938 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.0802 1.7062 -1.9245 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.0002 0.8451 -3.1539 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5025 -0.6880 -1.3241 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.3537 1.2659 0.3463 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6741 -0.1560 0.6311 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.0569 -0.2247 0.9886 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 2 0 + 1 6 1 0 + 2 4 1 0 + 2 7 1 0 + 3 5 2 0 + 4 11 1 0 + 4 12 2 0 + 5 9 1 0 + 5 10 1 0 + 6 9 1 0 + 7 8 1 0 + 8 9 1 0 + 1 13 1 0 + 3 14 1 0 + 3 15 1 0 + 6 16 1 0 + 6 17 1 0 + 7 18 1 0 + 7 19 1 0 + 8 20 1 0 + 8 21 1 0 + 9 22 1 1 + 10 23 1 0 + 10 24 1 0 + 10 25 1 0 +M CHG 1 11 -1 +M END + +> +-9.8733339999999998 + +> +OA-G2 + +> +-0.24234159156680107 -0.13431934878230095 -0.236714898198843 0.91454404726624494 -0.1058359132707119 -0.0076737576723098751 -0.029078564196825026 -0.066287456601858141 +-0.033445252627134325 -0.06594990059733391 -0.83899009093642229 -0.83899009093642229 0.11854273959994316 0.10526807740330696 0.10526807740330696 0.021749406531453134 0.021749406531453134 +0.047736430019140241 0.047736430019140241 0.03281161695718765 0.03281161695718765 0.046814215779304502 0.034864944815635679 0.034864944815635679 0.034864911288022993 + +$$$$ +OA-G3 + RDKit 3D + + 17 16 0 0 0 0 0 0 0 0999 V2000 + -0.2421 -0.0718 0.3818 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8393 -0.2077 -0.3937 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.3051 0.2551 -3.3541 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.1850 0.7201 -1.5211 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.1018 -0.6816 -3.2022 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.2330 -0.0262 -2.8578 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5879 0.9108 -2.3102 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8592 0.2770 -4.4873 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4542 -0.7664 1.1869 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.9353 0.7462 0.2188 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.5207 -1.0316 -0.1977 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.1638 1.1681 -1.3121 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4650 1.5450 -1.5737 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.3602 -1.3853 -2.4008 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0074 -1.2732 -4.1197 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.4879 0.6810 -3.6576 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.0316 -0.7782 -2.8399 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 2 0 + 2 4 1 0 + 3 5 1 0 + 3 7 1 0 + 3 8 2 0 + 4 6 1 0 + 5 6 1 0 + 1 9 1 0 + 1 10 1 0 + 2 11 1 0 + 4 12 1 0 + 4 13 1 0 + 5 14 1 0 + 5 15 1 0 + 6 16 1 0 + 6 17 1 0 +M CHG 1 7 -1 +M END + +> +-7.7345790000000001 + +> +OA-G3 + +> +-0.25688242101494002 -0.14571022176567247 0.90786136007484264 -0.056502479402457964 -0.21427863025489977 -0.055603332160150301 -0.84741937303367787 -0.84741937303367787 0.10201528734144043 +0.10201528734144043 0.10772500997957062 0.038299766095245588 0.038299766095245588 0.016157646389568552 0.016157646389568552 0.047642030479276884 0.047642030479276884 + +$$$$ +OA-G4 + RDKit 3D + + 29 28 0 0 0 0 0 0 0 0999 V2000 + -0.3023 -0.8943 -1.9056 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4373 -0.5293 -0.8413 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.0556 1.4485 -5.1381 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.1104 -0.2112 0.5212 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.9350 -0.4263 -0.9512 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.9787 1.0573 -5.6003 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.1827 -1.2245 -3.2898 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.0475 1.5230 -4.0753 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8663 -0.8574 -4.3441 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.2784 0.6305 -4.3059 C 0 0 2 0 0 0 0 0 0 0 0 0 + 2.0098 2.2581 -4.9567 O 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8803 0.6125 -6.0669 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.3802 -0.9841 -1.7762 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.1488 -1.1197 1.1320 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.1238 0.1959 0.4987 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.5344 0.5154 1.0259 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.2639 0.5732 -1.2544 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.3111 -1.1387 -1.6928 H 0 0 0 0 0 0 0 0 0 0 0 0 + 2.4012 -0.6500 0.0139 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.3241 0.9347 -6.4681 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.2613 2.1149 -5.5589 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.8879 0.4712 -5.7716 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.3887 -2.3007 -3.3170 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.1330 -0.7307 -3.5143 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.3900 1.2984 -3.0936 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.3600 2.5753 -4.0472 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4864 -1.1087 -5.3434 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.7539 -1.4826 -4.1850 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.9857 0.7778 -3.4799 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 2 0 + 1 7 1 0 + 2 4 1 0 + 2 5 1 0 + 3 8 1 0 + 3 11 1 0 + 3 12 2 0 + 6 10 1 0 + 7 9 1 0 + 8 10 1 0 + 9 10 1 0 + 1 13 1 0 + 4 14 1 0 + 4 15 1 0 + 4 16 1 0 + 5 17 1 0 + 5 18 1 0 + 5 19 1 0 + 6 20 1 0 + 6 21 1 0 + 6 22 1 0 + 7 23 1 0 + 7 24 1 0 + 8 25 1 0 + 8 26 1 0 + 9 27 1 0 + 9 28 1 0 + 10 29 1 1 +M CHG 1 11 -1 +M END + +> +-9.1737099999999998 + +> +OA-G4 + +> +-0.16766302782142983 -0.12061256664837229 0.90602863830482139 -0.070390922306426643 -0.070390922306426643 -0.078514670310863136 -0.042349432884105323 -0.19354121881569253 +-0.064572726665385841 -0.05895461264098513 -0.84640865522468911 -0.84640865522468911 0.11097088706647527 0.034675578297726037 0.034675578297726037 0.034675578297726037 0.034675578297726037 +0.034675578297726037 0.034675578297726037 0.024976786868325596 0.024976786868325596 0.024976786868325596 0.042367330463281991 0.042367330463281991 0.018510187687031155 0.018510187687031155 +0.037587809087387444 0.037587809087387444 0.062893400611034753 + +$$$$ +OA-G5 + RDKit 3D + + 17 16 0 0 0 0 0 0 0 0999 V2000 + 0.4721 -0.3596 0.1049 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.6140 0.0086 -0.5838 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.3138 -0.7087 -3.7837 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.6475 1.3143 -4.0310 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.5825 0.9091 -1.7860 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.2344 0.2559 -3.0098 C 0 0 2 0 0 0 0 0 0 0 0 0 + 0.5382 -1.3387 -3.0968 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.5250 -0.7473 -5.0305 O 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4022 -1.0124 0.9678 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.4547 -0.0073 -0.1892 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.5854 -0.3514 -0.2552 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7758 1.8832 -4.3717 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.3607 2.0256 -3.5996 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.1200 0.8627 -4.9107 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.4463 1.2159 -2.0129 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.1278 1.8263 -1.5296 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.1320 -0.2920 -2.6927 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 2 0 + 2 5 1 0 + 3 6 1 0 + 3 7 1 0 + 3 8 2 0 + 4 6 1 0 + 5 6 1 0 + 1 9 1 0 + 1 10 1 0 + 2 11 1 0 + 4 12 1 0 + 4 13 1 0 + 4 14 1 0 + 5 15 1 0 + 5 16 1 0 + 6 17 1 1 +M CHG 1 7 -1 +M END + +> +-7.6388220000000002 + +> +OA-G5 + +> +-0.2979469792369534 -0.10777685506378903 0.91644908027613869 -0.067875248325221682 -0.053967733812682772 -0.20121777517830625 -0.84354089660679588 -0.84354089660679588 0.094638693201191282 +0.094638693201191282 0.11095144884551272 0.018116577583200792 0.018116577583200792 0.018116577583200792 0.05721881898010478 0.05721881898010478 0.030401098596699098 + +$$$$ +OA-G6 + RDKit 3D + + 19 18 0 0 0 0 0 0 0 0999 V2000 + -1.3428 0.3977 -4.1555 C 0 0 0 0 0 0 0 0 0 0 0 0 + 1.3494 0.8562 -1.6206 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.1465 0.0031 0.1950 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4127 -0.7040 -3.6398 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4163 -0.9172 -2.1271 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.0498 0.3297 -1.2978 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.4251 0.5122 -5.4098 O 0 0 0 0 0 0 0 0 0 0 0 0 + -1.9062 1.0746 -3.2485 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.1072 0.0897 -1.4250 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.6012 1.7234 -1.0027 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.4304 1.1499 -2.6726 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.5555 -0.7928 0.4655 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.1538 -0.3402 0.4542 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0812 0.8805 0.8099 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.6080 -0.5200 -3.9934 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7392 -1.6346 -4.1217 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.2862 -1.7273 -1.8917 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.4072 -1.2739 -1.8180 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7712 1.1284 -1.5083 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 4 1 0 + 1 7 1 0 + 1 8 2 0 + 2 6 1 0 + 3 6 1 0 + 4 5 1 0 + 5 6 1 0 + 2 9 1 0 + 2 10 1 0 + 2 11 1 0 + 3 12 1 0 + 3 13 1 0 + 3 14 1 0 + 4 15 1 0 + 4 16 1 0 + 5 17 1 0 + 5 18 1 0 + 6 19 1 0 +M CHG 1 7 -1 +M END + +> +-7.5684560000000003 + +> +OA-G6 + +> +0.90856356252180903 -0.083272322228080342 -0.083272322228080342 -0.21121091317189367 -0.054996575375920849 -0.083306490590697835 -0.84895735155594976 -0.84895735155594976 +0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022712499687546177 0.022712499687546177 0.034720391819351597 +0.034720391819351597 0.0582043871675667 + +$$$$ diff --git a/openfe/tests/data/host_guest/__init__.py b/openfe/tests/data/host_guest/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/openfe/tests/protocols/openmm_abfe/conftest.py b/openfe/tests/protocols/openmm_abfe/conftest.py index 2bca3b1f9..7a356e0cb 100644 --- a/openfe/tests/protocols/openmm_abfe/conftest.py +++ b/openfe/tests/protocols/openmm_abfe/conftest.py @@ -1,5 +1,7 @@ # This code is part of OpenFE and is licensed under the MIT license. # For details, see https://github.com/OpenFreeEnergy/openfe +from importlib import resources +from rdkit import Chem import gufe import pytest from openfe.protocols import openmm_afe @@ -31,3 +33,22 @@ def benzene_complex_dag(benzene_modifications, T4_protein_component): ) return protocol.create(stateA=stateA, stateB=stateB, mapping=None) + + +@pytest.fixture(scope="session") +def guests_OA(): + files = {} + with resources.as_file(resources.files("openfe.tests.data.host_guest.OA")) as d: + fn = str(d / "guests_charged.sdf") + supp = Chem.SDMolSupplier(str(fn), removeHs=False) + for rdmol in supp: + files[rdmol.GetProp("_Name")] = gufe.SmallMoleculeComponent(rdmol) + return files + + +@pytest.fixture(scope="session") +def host_OA(): + with resources.as_file(resources.files("openfe.tests.data.host_guest.OA")) as d: + fn = str(d / "OA_charged.sdf") + rdmol = [m for m in Chem.SDMolSupplier(str(fn), removeHs=False)][0] + return gufe.SmallMoleculeComponent(rdmol) diff --git a/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index 1d9d328df..2a9b53c41 100644 --- a/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -1,6 +1,5 @@ # This code is part of OpenFE and is licensed under the MIT license. # For details, see https://github.com/OpenFreeEnergy/openfe -from importlib import resources from math import sqrt from unittest import mock @@ -10,18 +9,15 @@ import openfe import pytest from numpy.testing import assert_allclose -from openmmtools.alchemy import ( - AlchemicalRegion, - AlchemicalState, - AbsoluteAlchemicalFactory, -) +from openmmtools.alchemy import AlchemicalRegion from openmmtools.tests.test_alchemy import ( compare_system_energies, check_noninteracting_energy_components, check_interacting_energy_components, + check_multi_noninteracting_energy_components, + check_multi_interacting_energy_components, ) from openfe import ChemicalSystem, SolventComponent, SmallMoleculeComponent -from openfe.protocols.openmm_utils.omm_settings import OpenMMSolvationSettings from openfe.protocols import openmm_afe from openfe.protocols.openmm_afe import ( AbsoluteBindingComplexUnit, @@ -330,14 +326,14 @@ def _test_energies(reference_system, alchemical_system, alchemical_regions, posi positions=positions, ) - check_noninteracting_energy_components( + check_multi_noninteracting_energy_components( reference_system=reference_system, alchemical_system=alchemical_system, alchemical_regions=alchemical_regions, positions=positions, ) - check_interacting_energy_components( + check_multi_interacting_energy_components( reference_system=reference_system, alchemical_system=alchemical_system, alchemical_regions=alchemical_regions, @@ -359,7 +355,7 @@ def test_complex_dry_run(self, complex_units, settings, tmpdir): # Check the alchemical indices expected_indices = [i + self.num_complex_atoms for i in range(self.num_solvent_atoms)] - assert expected_indices == data["alchem_indices"][0] + assert expected_indices == data["alchem_indices"]["A"] # Check the non-alchemical system self._assert_expected_nonalchemical_forces(data["system"], settings) @@ -375,11 +371,11 @@ def test_complex_dry_run(self, complex_units, settings, tmpdir): assert pdb.n_atoms == self.num_all_not_water # Check energies - alchem_region = AlchemicalRegion(alchemical_atoms=data["alchem_indices"][0]) + alchem_regions = [AlchemicalRegion(alchemical_atoms=data["alchem_indices"]["A"], name="A")] self._test_energies( reference_system=data["system"], alchemical_system=data["alchem_system"], - alchemical_regions=alchem_region, + alchemical_regions=alchem_regions, positions=data["positions"], ) @@ -398,7 +394,7 @@ def test_solvent_dry_run(self, solvent_units, settings, tmpdir): # Check the alchemical indices expected_indices = [i for i in range(self.num_solvent_atoms)] - assert expected_indices == data["alchem_indices"][0] + assert expected_indices == data["alchem_indices"]["A"] # Check the non-alchemical system self._assert_expected_nonalchemical_forces(data["system"], settings) @@ -414,7 +410,7 @@ def test_solvent_dry_run(self, solvent_units, settings, tmpdir): assert pdb.n_atoms == self.num_solvent_atoms # Check energies - alchem_region = AlchemicalRegion(alchemical_atoms=data["alchem_indices"][0]) + alchem_region = [AlchemicalRegion(alchemical_atoms=data["alchem_indices"]["A"], name="A")] self._test_energies( reference_system=data["system"], @@ -453,6 +449,210 @@ def settings(self): return s +class TestOAHostGuestNetCharge(TestT4LysozymeDryRun): + """ + OA-G0 host guest complex, with a net charged ligand and a + user defined restraint. + """ + solvent = SolventComponent(ion_concentration=0.15 * offunit.molar) + num_all_not_water = 219 + num_complex_atoms = 192 + num_ligand_atoms = 20 + num_solvent_atoms = 25 + + @pytest.fixture(scope="class") + def settings(self): + s = openmm_afe.AbsoluteBindingProtocol.default_settings() + s.protocol_repeats = 1 + s.engine_settings.compute_platform = "cpu" + s.forcefield_settings.nonbonded_cutoff = 1.2 * offunit.nanometer + s.complex_output_settings.output_indices = "not water" + s.restraint_settings.guest_restraint_ids = [1, 2, 4] + s.restraint_settings.host_restraint_ids = [1, 4, 3] + s.complex_solvation_settings.solvent_padding = 1.5 * offunit.nanometer + s.alchemical_settings.alchemical_ion_min_distance = 1.0 * offunit.nanometer + return s + + @pytest.fixture(scope="class") + def dag(self, protocol, guests_OA, host_OA): + stateA = ChemicalSystem( + { + "ligand": guests_OA["OA-G0"], + "host": host_OA, + "solvent": self.solvent, + } + ) + + stateB = ChemicalSystem( + { + "host": host_OA, + "solvent": self.solvent, + } + ) + + return protocol.create( + stateA=stateA, + stateB=stateB, + mapping=None, + ) + + def _assert_expected_alchemical_forces(self, system, complexed: bool, settings): + """ + Assert the forces expected in the alchemical system. + """ + self._assert_force_num(system, NonbondedForce, 1) + self._assert_force_num(system, CustomNonbondedForce, 4) + self._assert_force_num(system, CustomBondForce, 4) + if complexed: + self._assert_force_num(system, HarmonicBondForce, 1) + else: + self._assert_force_num(system, HarmonicBondForce, 2) + self._assert_force_num(system, HarmonicAngleForce, 1) + self._assert_force_num(system, PeriodicTorsionForce, 1) + self._assert_force_num(system, MonteCarloBarostat, 1) + + if complexed: + self._assert_force_num(system, CustomCompoundBondForce, 1) + + assert len(system.getForces()) == 14 + + # Check the nonbonded force has the right contents + nonbond = [f for f in system.getForces() if isinstance(f, NonbondedForce)] + assert len(nonbond) == 1 + assert nonbond[0].getNonbondedMethod() == NonbondedForce.PME + assert ( + from_openmm(nonbond[0].getCutoffDistance()) + == settings.forcefield_settings.nonbonded_cutoff + ) + + # Check the barostat made it all the way through + barostat = [f for f in system.getForces() if isinstance(f, MonteCarloBarostat)] + assert len(barostat) == 1 + assert barostat[0].getFrequency() == int(settings.integrator_settings.barostat_frequency.m) + assert barostat[0].getDefaultPressure() == to_openmm(settings.thermo_settings.pressure) + assert barostat[0].getDefaultTemperature() == to_openmm( + settings.thermo_settings.temperature + ) + + @staticmethod + def _test_energies(reference_system, alchemical_system, alchemical_regions, positions): + #compare_system_energies( + # reference_system=reference_system, + # alchemical_system=alchemical_system, + # alchemical_regions=alchemical_regions, + # positions=positions, + #) + + check_multi_noninteracting_energy_components( + reference_system=reference_system, + alchemical_system=alchemical_system, + alchemical_regions=alchemical_regions, + positions=positions, + ) + + check_multi_interacting_energy_components( + reference_system=reference_system, + alchemical_system=alchemical_system, + alchemical_regions=alchemical_regions, + positions=positions, + ) + + def test_complex_dry_run(self, complex_units, settings, tmpdir): + with tmpdir.as_cwd(): + data = complex_units[0].run(dry=True, verbose=True)["debug"] + + # Check the sampler + self._verify_sampler(data["sampler"], complexed=True, settings=settings) + + # Check the alchemical system + self._assert_expected_alchemical_forces( + data["alchem_system"], complexed=True, settings=settings + ) + self._test_dodecahedron_vectors(data["alchem_system"]) + + # Check the alchemical indices + # since the host is an SMC, it's guest then host + expected_indices = [i for i in range(self.num_ligand_atoms)] + assert expected_indices == data["alchem_indices"]["A"] + + # check there's only one alchemical ion + assert len(data["alchem_indices"]["B"]) == 1 + + # Check the non-alchemical system + self._assert_expected_nonalchemical_forces(data["system"], settings) + self._test_dodecahedron_vectors(data["system"]) + # Check the box vectors haven't changed (they shouldn't have because we didn't do MD) + assert_allclose( + from_openmm(data["alchem_system"].getDefaultPeriodicBoxVectors()), + from_openmm(data["system"].getDefaultPeriodicBoxVectors()), + ) + + # Check the PDB + pdb = mdt.load_pdb("alchemical_system.pdb") + assert pdb.n_atoms == self.num_all_not_water + + # Check energies + alchem_regions = [ + AlchemicalRegion(alchemical_atoms=entry) + for entry in data["alchem_indices"] + ] + + assert len(alchem_regions) == 2 # should be 2 with alchemical ion + + self._test_energies( + reference_system=data["system"], + alchemical_system=data["alchem_system"], + alchemical_regions=alchem_regions, + positions=data["positions"], + ) + + def test_solvent_dry_run(self, solvent_units, settings, tmpdir): + with tmpdir.as_cwd(): + data = solvent_units[0].run(dry=True, verbose=True)["debug"] + + # Check the sampler + self._verify_sampler(data["sampler"], complexed=False, settings=settings) + + # Check the alchemical system + self._assert_expected_alchemical_forces( + data["alchem_system"], complexed=False, settings=settings + ) + self._test_dodecahedron_vectors(data["alchem_system"]) + + # Check the alchemical indices + expected_indices = [i for i in range(self.num_ligand_atoms)] + assert expected_indices == data["alchem_indices"]["A"] + + # Check the non-alchemical system + self._assert_expected_nonalchemical_forces(data["system"], settings) + self._test_dodecahedron_vectors(data["system"]) + # Check the box vectors haven't changed (they shouldn't have because we didn't do MD) + assert_allclose( + from_openmm(data["alchem_system"].getDefaultPeriodicBoxVectors()), + from_openmm(data["system"].getDefaultPeriodicBoxVectors()), + ) + + # Check the PDB + pdb = mdt.load_pdb("alchemical_system.pdb") + assert pdb.n_atoms == self.num_solvent_atoms + + # Check energies + assert len(data["alchem_indices"]) == 2 + + alchem_regions = [] + for k, v in data["alchem_indices"].items(): + alchem_regions.append( + AlchemicalRegion(alchemical_atoms=v, name=k) + ) + + self._test_energies( + reference_system=data["restrained_system"], + alchemical_system=data["alchem_system"], + alchemical_regions=alchem_regions, + positions=data["positions"], + ) + + def test_user_charges(benzene_modifications, T4_protein_component, tmpdir): s = openmm_afe.AbsoluteBindingProtocol.default_settings() s.protocol_repeats = 1