From e47b90a37228b985d0ec39e8ea67dc2cec989cba Mon Sep 17 00:00:00 2001 From: IAlibay Date: Tue, 16 Jun 2026 11:32:33 +0100 Subject: [PATCH 01/31] enable user defined restraints in abfe calculations --- src/openfe/protocols/openmm_afe/abfe_units.py | 2 + .../protocols/restraint_utils/settings.py | 40 ++++++++++--------- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 343cbd2e8..f8be7618e 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -254,6 +254,8 @@ def _get_boresch_restraint( guest_rdmol=guest_rdmol, guest_idxs=guest_atom_ids, host_idxs=host_atom_ids, + guest_restraint_atoms_idxs=settings.guest_restraint_ids, + host_restraint_atoms_idxs=settings.host_restraint_ids, host_selection=settings.host_selection, anchor_finding_strategy=settings.anchor_finding_strategy, dssp_filter=settings.dssp_filter, diff --git a/src/openfe/protocols/restraint_utils/settings.py b/src/openfe/protocols/restraint_utils/settings.py index 231b6c141..c7e6aa7cf 100644 --- a/src/openfe/protocols/restraint_utils/settings.py +++ b/src/openfe/protocols/restraint_utils/settings.py @@ -174,17 +174,18 @@ class BoreschRestraintSettings(BaseRestraintSettings): Boresch-like restraint search parameter. The maximum distance between any host atom and the guest G0 atom. Must be in units compatible with nanometer. """ - # TODO: re-enable this (Issue #1555) - # host_atoms: Optional[list[int]] = None - # """ - # The indices of the host component atoms to restrain. - # If defined, these will override any automatic selection. - # """ - # guest_atoms: Optional[list[int]] = None - # """ - # The indices of the guest component atoms to restraint. - # If defined, these will override any automatic selection. - # """ + host_atoms: Optional[tuple[int, int, int]] = None + """ + The indices of the host component atoms to restrain. + The entries define the H0, H1, and H2 atoms in order. + If defined, these will override any automatic selection. + """ + guest_atoms: Optional[tuple[int, int, int]] = None + """ + The indices of the guest component atoms to restraint. + The entries define the G0, G1, and G2 atoms in order. + If defined, these will override any automatic selection. + """ anchor_finding_strategy: Literal["multi-residue", "bonded"] = "bonded" """ The Boresch atom picking strategy to use. @@ -194,10 +195,13 @@ class BoreschRestraintSettings(BaseRestraintSettings): * `multi-residue`: pick host atoms which can span multiple residues. """ - -# @field_validator("guest_atoms", "host_atoms") -# def positive_idxs_list(cls, v): -# if v is not None and any([i < 0 for i in v]): -# errmsg = "negative indices passed" -# raise ValueError(errmsg) -# return v + @field_validator("guest_atoms", "host_atoms") + def positive_idxs_three_tuple(cls, v): + if v is not None: + if len(v) != 3: + errmsg = "``guest_atoms`` and ``host_atoms`` must contain three elements." + raise ValueError(errmsg) + if any([i < 0 for i in v]): + errmsg = "``guest_atoms`` and ``host_atoms`` cannot have negative indices." + raise ValueError(errmsg) + return v From a453b53a4e713d3d06644a050a73f02df5630fe2 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Tue, 16 Jun 2026 11:48:07 +0100 Subject: [PATCH 02/31] Add something specific to abfe restraints --- .../openmm_afe/equil_afe_settings.py | 25 +++++++++++++++++++ .../protocols/restraint_utils/settings.py | 23 ----------------- 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index 0941b6784..d17329617 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -213,6 +213,31 @@ def must_be_all(cls, v): return v +class ABFEBoreschRestraintSettings(BoreschRestraintSettings): + host_atoms: Optional[tuple[int, int, int]] = None + """ + The indices of the host component atoms to restrain. + The entries define the H0, H1, and H2 atoms in order. + If defined, these will override any automatic selection. + """ + guest_atoms: Optional[tuple[int, int, int]] = None + """ + The indices of the guest component atoms to restraint. + The entries define the G0, G1, and G2 atoms in order. + If defined, these will override any automatic selection. + """ + @field_validator("guest_atoms", "host_atoms") + def positive_idxs_three_tuple(cls, v): + if v is not None: + if len(v) != 3: + errmsg = "``guest_atoms`` and ``host_atoms`` must contain three elements." + raise ValueError(errmsg) + if any([i < 0 for i in v]): + errmsg = "``guest_atoms`` and ``host_atoms`` cannot have negative indices." + raise ValueError(errmsg) + return v + + # This subclasses from SettingsBaseModel as it has vacuum_forcefield and # solvent_forcefield fields, not just a single forcefield_settings field class AbsoluteSolvationSettings(SettingsBaseModel): diff --git a/src/openfe/protocols/restraint_utils/settings.py b/src/openfe/protocols/restraint_utils/settings.py index c7e6aa7cf..8ab8dcd86 100644 --- a/src/openfe/protocols/restraint_utils/settings.py +++ b/src/openfe/protocols/restraint_utils/settings.py @@ -174,18 +174,6 @@ class BoreschRestraintSettings(BaseRestraintSettings): Boresch-like restraint search parameter. The maximum distance between any host atom and the guest G0 atom. Must be in units compatible with nanometer. """ - host_atoms: Optional[tuple[int, int, int]] = None - """ - The indices of the host component atoms to restrain. - The entries define the H0, H1, and H2 atoms in order. - If defined, these will override any automatic selection. - """ - guest_atoms: Optional[tuple[int, int, int]] = None - """ - The indices of the guest component atoms to restraint. - The entries define the G0, G1, and G2 atoms in order. - If defined, these will override any automatic selection. - """ anchor_finding_strategy: Literal["multi-residue", "bonded"] = "bonded" """ The Boresch atom picking strategy to use. @@ -194,14 +182,3 @@ class BoreschRestraintSettings(BaseRestraintSettings): * `bonded`: pick host atoms that are bonded to each other. * `multi-residue`: pick host atoms which can span multiple residues. """ - - @field_validator("guest_atoms", "host_atoms") - def positive_idxs_three_tuple(cls, v): - if v is not None: - if len(v) != 3: - errmsg = "``guest_atoms`` and ``host_atoms`` must contain three elements." - raise ValueError(errmsg) - if any([i < 0 for i in v]): - errmsg = "``guest_atoms`` and ``host_atoms`` cannot have negative indices." - raise ValueError(errmsg) - return v From 12c509f42f020bebed4f67f9ab01f3d9a543e746 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Tue, 16 Jun 2026 11:53:45 +0100 Subject: [PATCH 03/31] fix some typing --- src/openfe/protocols/openmm_afe/abfe_units.py | 8 ++++---- src/openfe/protocols/openmm_afe/equil_afe_settings.py | 4 ++-- .../protocols/openmm_afe/equil_binding_afe_method.py | 4 ++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index f8be7618e..cdd5ee947 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -26,7 +26,7 @@ from rdkit import Chem from openfe.protocols.openmm_afe.equil_afe_settings import ( - BoreschRestraintSettings, + ABFEBoreschRestraintSettings, SettingsBaseModel, ) from openfe.protocols.openmm_utils import system_validation @@ -218,7 +218,7 @@ def _get_boresch_restraint( guest_atom_ids: list[int], host_atom_ids: list[int], temperature: Quantity, - settings: BoreschRestraintSettings, + settings: ABFEBoreschRestraintSettings, ) -> tuple[BoreschRestraintGeometry, BoreschRestraint]: """ Get a Boresch-like restraint Geometry and OpenMM restraint force @@ -236,7 +236,7 @@ def _get_boresch_restraint( A list of atom indices defining the host molecules in the universe. temperature : openff.units.Quantity The temperature of the simulation where the restraint will be added. - settings : BoreschRestraintSettings + settings : ABFEBoreschRestraintSettings Settings on how the Boresch-like restraint should be defined. Returns @@ -357,7 +357,7 @@ def _add_restraints( self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, ) - if isinstance(settings["restraint_settings"], BoreschRestraintSettings): + if isinstance(settings["restraint_settings"], ABFEBoreschRestraintSettings): rest_geom, restraint = self._get_boresch_restraint( univ, guest_rdmol, diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index d17329617..fb6a3faee 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -214,13 +214,13 @@ def must_be_all(cls, v): class ABFEBoreschRestraintSettings(BoreschRestraintSettings): - host_atoms: Optional[tuple[int, int, int]] = None + host_atoms: tuple[int, int, int] | None = None """ The indices of the host component atoms to restrain. The entries define the H0, H1, and H2 atoms in order. If defined, these will override any automatic selection. """ - guest_atoms: Optional[tuple[int, int, int]] = None + guest_atoms: tuple[int, int, int] | None = None """ The indices of the guest component atoms to restraint. The entries define the G0, G1, and G2 atoms in order. diff --git a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 06dcbaed0..9048ca1ab 100644 --- a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -46,7 +46,7 @@ ABFEPreEquilOutputSettings, AbsoluteBindingSettings, AlchemicalSettings, - BoreschRestraintSettings, + ABFEBoreschRestraintSettings, IntegratorSettings, LambdaSettings, MDSimulationSettings, @@ -174,7 +174,7 @@ def _default_settings(cls): engine_settings=OpenMMEngineSettings(), solvent_integrator_settings=IntegratorSettings(), complex_integrator_settings=IntegratorSettings(), - restraint_settings=BoreschRestraintSettings(), + restraint_settings=ABFEBoreschRestraintSettings(), solvent_equil_simulation_settings=MDSimulationSettings( equilibration_length_nvt=0.1 * offunit.nanosecond, equilibration_length=0.2 * offunit.nanosecond, From f5933279baedc720dab351bf90e40ebdf923a5c1 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Tue, 16 Jun 2026 12:14:18 +0100 Subject: [PATCH 04/31] isolate restraints only for ABFEs for now --- src/openfe/protocols/openmm_afe/equil_afe_settings.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index fb6a3faee..602f9f6cb 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -214,19 +214,19 @@ def must_be_all(cls, v): class ABFEBoreschRestraintSettings(BoreschRestraintSettings): - host_atoms: tuple[int, int, int] | None = None + host_restraint_ids: tuple[int, int, int] | None = None """ The indices of the host component atoms to restrain. The entries define the H0, H1, and H2 atoms in order. If defined, these will override any automatic selection. """ - guest_atoms: tuple[int, int, int] | None = None + guest_restraint_ids: tuple[int, int, int] | None = None """ The indices of the guest component atoms to restraint. The entries define the G0, G1, and G2 atoms in order. If defined, these will override any automatic selection. """ - @field_validator("guest_atoms", "host_atoms") + @field_validator("guest_restraint_ids", "host_restraint_ids") def positive_idxs_three_tuple(cls, v): if v is not None: if len(v) != 3: From 505f4b701296f40b579d7c1807a9028ba356b975 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Tue, 16 Jun 2026 13:34:51 +0100 Subject: [PATCH 05/31] Add settings test --- src/openfe/protocols/openmm_afe/equil_afe_settings.py | 3 --- .../tests/protocols/openmm_abfe/test_abfe_settings.py | 10 ++++++++++ 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index 602f9f6cb..ee8a154d6 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -229,9 +229,6 @@ class ABFEBoreschRestraintSettings(BoreschRestraintSettings): @field_validator("guest_restraint_ids", "host_restraint_ids") def positive_idxs_three_tuple(cls, v): if v is not None: - if len(v) != 3: - errmsg = "``guest_atoms`` and ``host_atoms`` must contain three elements." - raise ValueError(errmsg) if any([i < 0 for i in v]): errmsg = "``guest_atoms`` and ``host_atoms`` cannot have negative indices." raise ValueError(errmsg) diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py index fad97e6ad..013c775ff 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py @@ -7,6 +7,7 @@ AbsoluteBindingProtocol, AbsoluteBindingSettings, ) +from openfe.protocols.openmm_afe.equil_afe_settings import ABFEBoreschRestraintSettings @pytest.fixture() @@ -70,6 +71,15 @@ def test_monotonic_lambda_windows(val, default_settings): lambda_settings.lambda_restraints = val["restraints"] +@pytest.mark.parametrize("parameter", ["host_restraint_ids", "guest_restraint_ids"]) +def test_boresch_restraints_restraint_negative_ids(parameter): + setting = ABFEBoreschRestraintSettings() + + errmsg = "``guest_atoms`` and ``host_atoms`` cannot have negative indices." + with pytest.raises(ValueError, match=errmsg): + setattr(setting, parameter, [1, 2, -3]) + + def test_equil_not_all_complex(default_settings): with pytest.raises(ValueError, match="output_indices must be all"): default_settings.complex_equil_output_settings.output_indices = "not water" From 6ddcba14a353ca48ff2a84736ba86c73633d53c5 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 13:16:42 +0100 Subject: [PATCH 06/31] Add test for user defined restraints --- src/openfe/protocols/openmm_afe/abfe_units.py | 4 +- src/openfe/tests/conftest.py | 11 + .../data/benzene_modifications_am1bcc.sdf | 295 ++++++++++++++++++ .../openmm_abfe/test_abfe_protocol.py | 44 +++ 4 files changed, 352 insertions(+), 2 deletions(-) create mode 100644 src/openfe/tests/data/benzene_modifications_am1bcc.sdf diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index cdd5ee947..72f8352ad 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -254,8 +254,8 @@ def _get_boresch_restraint( guest_rdmol=guest_rdmol, guest_idxs=guest_atom_ids, host_idxs=host_atom_ids, - guest_restraint_atoms_idxs=settings.guest_restraint_ids, - host_restraint_atoms_idxs=settings.host_restraint_ids, + guest_restraint_atoms_idxs=list(settings.guest_restraint_ids), + host_restraint_atoms_idxs=list(settings.host_restraint_ids), host_selection=settings.host_selection, anchor_finding_strategy=settings.anchor_finding_strategy, dssp_filter=settings.dssp_filter, diff --git a/src/openfe/tests/conftest.py b/src/openfe/tests/conftest.py index fcd41412f..3e9046b12 100644 --- a/src/openfe/tests/conftest.py +++ b/src/openfe/tests/conftest.py @@ -233,6 +233,17 @@ def benzene_modifications(): return files +@pytest.fixture(scope="session") +def benzene_modifications_am1bcc(): + files = {} + with resources.as_file(resources.files("openfe.tests.data")) as d: + fn = str(d / "benzene_modifications_am1bcc.sdf") + supp = Chem.SDMolSupplier(str(fn), removeHs=False) + for rdmol in supp: + files[rdmol.GetProp("_Name")] = SmallMoleculeComponent(rdmol) + return files + + @pytest.fixture(scope="session") def charged_benzene_modifications(): files = {} diff --git a/src/openfe/tests/data/benzene_modifications_am1bcc.sdf b/src/openfe/tests/data/benzene_modifications_am1bcc.sdf new file mode 100644 index 000000000..206c82d61 --- /dev/null +++ b/src/openfe/tests/data/benzene_modifications_am1bcc.sdf @@ -0,0 +1,295 @@ +benzene + RDKit 3D + + 12 12 0 0 0 0 0 0 0 0999 V2000 + 25.9780 5.3270 4.7790 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.3950 5.0740 3.4990 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.3400 5.8600 2.9020 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8370 6.9210 3.5690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.4200 7.1960 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.4980 6.3790 5.4690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.2298 4.6859 5.2451 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9676 4.2351 2.9497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.6890 5.6311 1.8951 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.5730 7.5660 3.0889 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8209 8.0598 5.3863 H 0 0 0 0 0 0 0 0 0 0 0 0 + 26.1874 6.5720 6.4958 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 2 0 + 1 6 1 0 + 1 7 1 0 + 2 3 1 0 + 2 8 1 0 + 3 4 2 0 + 3 9 1 0 + 4 5 1 0 + 4 10 1 0 + 5 6 2 0 + 5 11 1 0 + 6 12 1 0 +M END + +> +benzene + +> +-0.13 -0.13 -0.13 -0.13 -0.13 -0.13 0.13 0.13 0.13 0.13 0.13 0.13 + +$$$$ +toluene + RDKit 3D + + 15 15 0 0 0 0 0 0 0 0999 V2000 + 28.9072 8.7434 5.1220 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.1966 8.1433 6.6393 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.9864 8.4164 5.6052 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.2579 9.2269 5.5838 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9780 5.3270 4.7790 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.3950 5.0740 3.4990 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.3400 5.8600 2.9020 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8370 6.9210 3.5690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.4200 7.1960 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.4980 6.3790 5.4690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.2298 4.6859 5.2451 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9676 4.2351 2.9497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.6890 5.6311 1.8951 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.5730 7.5660 3.0889 H 0 0 0 0 0 0 0 0 0 0 0 0 + 26.1874 6.5720 6.4958 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 3 1 0 + 2 3 1 0 + 3 4 1 0 + 5 6 2 0 + 5 10 1 0 + 5 11 1 0 + 6 7 1 0 + 6 12 1 0 + 7 8 2 0 + 7 13 1 0 + 8 9 1 0 + 8 14 1 0 + 3 9 1 0 + 9 10 2 0 + 10 15 1 0 +M END + +> +toluene + +> +0.044033066666666669 0.044033066666666669 -0.053799933333333334 0.044033066666666669 -0.12699993333333334 -0.13499993333333335 -0.12699993333333334 -0.13099993333333335 +-0.077299933333333321 -0.13099993333333335 0.13000006666666666 0.13000006666666666 0.13000006666666666 0.13000006666666666 0.13000006666666666 + +$$$$ +phenol + RDKit 3D + + 13 13 0 0 0 0 0 0 0 0999 V2000 + 25.9780 5.3270 4.7790 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.3950 5.0740 3.4990 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.3400 5.8600 2.9020 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8370 6.9210 3.5690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.4200 7.1960 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.4980 6.3790 5.4690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.2298 4.6859 5.2451 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9676 4.2351 2.9497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.6890 5.6311 1.8951 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.5730 7.5660 3.0889 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.1311 8.0887 6.4624 H 0 0 0 0 0 0 0 0 0 0 0 0 + 26.1874 6.5720 6.4958 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.9460 8.3293 5.5517 O 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 2 0 + 1 6 1 0 + 1 7 1 0 + 2 3 1 0 + 2 8 1 0 + 3 4 2 0 + 3 9 1 0 + 4 5 1 0 + 4 10 1 0 + 5 6 2 0 + 5 13 1 0 + 6 12 1 0 + 11 13 1 0 +M END + +> +phenol + +> +-0.094423076923076915 -0.16592307692307692 -0.094423076923076915 -0.18492307692307691 0.12317692307692309 -0.18492307692307691 0.13307692307692309 0.13307692307692309 0.13307692307692309 +0.14157692307692307 0.41807692307692307 0.14157692307692307 -0.4990230769230769 + +$$$$ +benzonitrile + RDKit 3D + + 13 13 0 0 0 0 0 0 0 0999 V2000 + 28.5559 9.5700 6.2831 N 0 0 0 0 0 0 0 0 0 0 0 0 + 27.9981 8.4043 5.5824 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9780 5.3270 4.7790 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.3950 5.0740 3.4990 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.3400 5.8600 2.9020 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8370 6.9210 3.5690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.4200 7.1960 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.4980 6.3790 5.4690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.2298 4.6859 5.2451 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9676 4.2351 2.9497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.6890 5.6311 1.8951 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.5730 7.5660 3.0889 H 0 0 0 0 0 0 0 0 0 0 0 0 + 26.1874 6.5720 6.4958 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 3 0 + 3 4 2 0 + 3 8 1 0 + 3 9 1 0 + 4 5 1 0 + 4 10 1 0 + 5 6 2 0 + 5 11 1 0 + 6 7 1 0 + 6 12 1 0 + 2 7 1 0 + 7 8 2 0 + 8 13 1 0 +M END + +> +benzonitrile + +> +-0.36380000000000001 0.23380000000000001 -0.13500000000000001 -0.107 -0.13500000000000001 -0.090999999999999998 -0.019000000000000003 -0.090999999999999998 0.14000000000000001 +0.13800000000000001 0.14000000000000001 0.14499999999999999 0.14499999999999999 + +$$$$ +anisole + RDKit 3D + + 16 16 0 0 0 0 0 0 0 0999 V2000 + 29.2873 8.8784 4.9226 C 0 0 0 0 0 0 0 0 0 0 0 0 + 29.5502 9.7990 5.4437 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.0548 8.3459 5.4720 O 0 0 0 0 0 0 0 0 0 0 0 0 + 30.0866 8.1484 5.0502 H 0 0 0 0 0 0 0 0 0 0 0 0 + 29.1525 9.0868 3.8612 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9780 5.3270 4.7790 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.3950 5.0740 3.4990 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.3400 5.8600 2.9020 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8370 6.9210 3.5690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.4200 7.1960 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.4980 6.3790 5.4690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.2298 4.6859 5.2451 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9676 4.2351 2.9497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.6890 5.6311 1.8951 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.5730 7.5660 3.0889 H 0 0 0 0 0 0 0 0 0 0 0 0 + 26.1874 6.5720 6.4958 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 + 1 4 1 0 + 1 5 1 0 + 1 3 1 0 + 6 7 2 0 + 6 11 1 0 + 6 12 1 0 + 7 8 1 0 + 7 13 1 0 + 8 9 2 0 + 8 14 1 0 + 9 10 1 0 + 9 15 1 0 + 3 10 1 0 + 10 11 2 0 + 11 16 1 0 +M END + +> +anisole + +> +0.113825 0.043825000000000003 -0.32877500000000004 0.043825000000000003 0.043825000000000003 -0.097875000000000004 -0.16487500000000002 -0.097875000000000004 -0.17987500000000001 0.123225 +-0.17987500000000001 0.13212499999999999 0.13212499999999999 0.13212499999999999 0.14212499999999997 0.14212499999999997 + +$$$$ +benzaldehyde + RDKit 3D + + 14 14 0 0 0 0 0 0 0 0999 V2000 + 29.2079 8.8492 4.9632 O 0 0 0 0 0 0 0 0 0 0 0 0 + 27.5482 8.8691 6.4597 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.9981 8.4043 5.5824 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9780 5.3270 4.7790 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.3950 5.0740 3.4990 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.3400 5.8600 2.9020 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8370 6.9210 3.5690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.4200 7.1960 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.4980 6.3790 5.4690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.2298 4.6859 5.2451 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9676 4.2351 2.9497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.6890 5.6311 1.8951 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.5730 7.5660 3.0889 H 0 0 0 0 0 0 0 0 0 0 0 0 + 26.1874 6.5720 6.4958 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 3 2 0 + 2 3 1 0 + 4 5 2 0 + 4 9 1 0 + 4 10 1 0 + 5 6 1 0 + 5 11 1 0 + 6 7 2 0 + 6 12 1 0 + 7 8 1 0 + 7 13 1 0 + 3 8 1 0 + 8 9 2 0 + 9 14 1 0 +M END + +> +benzaldehyde + +> +-0.52817142857142862 -0.0028714285714285795 0.5754285714285714 -0.14507142857142857 -0.098071428571428587 -0.14507142857142857 -0.078071428571428583 -0.19767142857142858 +-0.078071428571428583 0.13742857142857143 0.13492857142857143 0.13742857142857143 0.14392857142857141 0.14392857142857141 + +$$$$ +styrene + RDKit 3D + + 16 16 0 0 0 0 0 0 0 0999 V2000 + 29.2873 8.8784 4.9226 C 0 0 0 0 0 0 0 0 0 0 0 0 + 29.6609 8.3486 4.0463 H 0 0 0 0 0 0 0 0 0 0 0 0 + 29.8344 9.7353 5.3157 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.5365 8.8812 6.4825 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.9864 8.4164 5.6052 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9780 5.3270 4.7790 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.3950 5.0740 3.4990 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.3400 5.8600 2.9020 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.8370 6.9210 3.5690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 27.4200 7.1960 4.8560 C 0 0 0 0 0 0 0 0 0 0 0 0 + 26.4980 6.3790 5.4690 C 0 0 0 0 0 0 0 0 0 0 0 0 + 25.2298 4.6859 5.2451 H 0 0 0 0 0 0 0 0 0 0 0 0 + 25.9676 4.2351 2.9497 H 0 0 0 0 0 0 0 0 0 0 0 0 + 27.6890 5.6311 1.8951 H 0 0 0 0 0 0 0 0 0 0 0 0 + 28.5730 7.5660 3.0889 H 0 0 0 0 0 0 0 0 0 0 0 0 + 26.1874 6.5720 6.4958 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 + 1 3 1 0 + 1 5 2 0 + 4 5 1 0 + 6 7 2 0 + 6 11 1 0 + 6 12 1 0 + 7 8 1 0 + 7 13 1 0 + 8 9 2 0 + 8 14 1 0 + 9 10 1 0 + 9 15 1 0 + 5 10 1 0 + 10 11 2 0 + 11 16 1 0 +M END + +> +styrene + +> +-0.20899999999999999 0.11349999999999999 0.11349999999999999 0.123 -0.1132 -0.13100000000000001 -0.127 -0.13100000000000001 -0.11849999999999999 -0.057800000000000004 -0.11849999999999999 +0.13100000000000001 0.13100000000000001 0.13100000000000001 0.13150000000000001 0.13150000000000001 + +$$$$ diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index 6ba27ca2b..1186f313f 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -41,6 +41,7 @@ from openfe.protocols.openmm_afe import ( AbsoluteBindingProtocol, ) +from openfe.protocols.restraint_utils.geometry import BoreschRestraintGeometry from .utils import UNIT_TYPES, _get_units @@ -614,6 +615,49 @@ def settings(self): return s +def test_user_restraint(benzene_modifications_am1bcc, T4_protein_component, tmp_path): + s = openmm_afe.AbsoluteBindingProtocol.default_settings() + s.protocol_repeats = 1 + s.engine_settings.compute_platform = "cpu" + s.restraint_settings.guest_restraint_ids = [0, 1, 2] + # Ca and C from VAL 87, and N from TYR 88 + s.restraint_settings.host_restraint_ids = [1383, 1384, 1398] + + protocol = openmm_afe.AbsoluteBindingProtocol(settings=s) + + stateA = gufe.ChemicalSystem( + { + "protein": T4_protein_component, + "benzene": benzene_modifications_am1bcc["benzene"], + "solvent": gufe.SolventComponent(), + } + ) + + stateB = gufe.ChemicalSystem( + { + "protein": T4_protein_component, + "solvent": gufe.SolventComponent(), + } + ) + + dag = protocol.create(stateA=stateA, stateB=stateB, mapping=None) + + complex_setup_units = _get_units(dag.protocol_units, UNIT_TYPES["complex"]["setup"]) + + results = complex_setup_units[0].run(dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path) + + geom = BoreschRestraintGeometry.model_validate(results['restraint_geometry']) + # This should be C1, C2, and C3 on the benzene + assert geom.guest_atoms == [2613, 2614, 2615] + assert geom.host_atoms == [1383, 1384, 1398] + assert pytest.approx(geom.r_aA0.to("nanometer").m, rel=1e-4) == 0.510798 + assert pytest.approx(geom.theta_A0.to("radians").m, rel=1e-4) == 1.20278 + assert pytest.approx(geom.theta_B0.to("radians").m, rel=1e-4) == 1.25705 + assert pytest.approx(geom.phi_A0.to("radians").m, rel=1e-4) == 0.86035 + assert pytest.approx(geom.phi_B0.to("radians").m, rel=1e-4) == 1.59444 + assert pytest.approx(geom.phi_C0.to("radians").m, rel=1e-4) == 2.92365 + + def test_user_charges(benzene_modifications, T4_protein_component, tmp_path): s = openmm_afe.AbsoluteBindingProtocol.default_settings() s.protocol_repeats = 1 From 801220a8047f1661e7029e70cbe23d18c5a3d0b6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 17 Jun 2026 12:17:24 +0000 Subject: [PATCH 07/31] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/openfe/protocols/openmm_afe/equil_afe_settings.py | 1 + src/openfe/protocols/openmm_afe/equil_binding_afe_method.py | 2 +- .../tests/protocols/openmm_abfe/test_abfe_protocol.py | 6 ++++-- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index ee8a154d6..5853b4f76 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -226,6 +226,7 @@ class ABFEBoreschRestraintSettings(BoreschRestraintSettings): The entries define the G0, G1, and G2 atoms in order. If defined, these will override any automatic selection. """ + @field_validator("guest_restraint_ids", "host_restraint_ids") def positive_idxs_three_tuple(cls, v): if v is not None: diff --git a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 9048ca1ab..d62a21ac2 100644 --- a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -43,10 +43,10 @@ from openfe.due import Doi, due from openfe.protocols.openmm_afe.equil_afe_settings import ( + ABFEBoreschRestraintSettings, ABFEPreEquilOutputSettings, AbsoluteBindingSettings, AlchemicalSettings, - ABFEBoreschRestraintSettings, IntegratorSettings, LambdaSettings, MDSimulationSettings, diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index 1186f313f..29139c285 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -644,9 +644,11 @@ def test_user_restraint(benzene_modifications_am1bcc, T4_protein_component, tmp_ complex_setup_units = _get_units(dag.protocol_units, UNIT_TYPES["complex"]["setup"]) - results = complex_setup_units[0].run(dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path) + results = complex_setup_units[0].run( + dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path + ) - geom = BoreschRestraintGeometry.model_validate(results['restraint_geometry']) + geom = BoreschRestraintGeometry.model_validate(results["restraint_geometry"]) # This should be C1, C2, and C3 on the benzene assert geom.guest_atoms == [2613, 2614, 2615] assert geom.host_atoms == [1383, 1384, 1398] From b3de6fd0fe0fdd7478d6244814b1c6962d9b6ea9 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 13:18:47 +0100 Subject: [PATCH 08/31] Add a news item --- news/abfe_user_dfn_restraints.rst | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 news/abfe_user_dfn_restraints.rst diff --git a/news/abfe_user_dfn_restraints.rst b/news/abfe_user_dfn_restraints.rst new file mode 100644 index 000000000..65c1d3b75 --- /dev/null +++ b/news/abfe_user_dfn_restraints.rst @@ -0,0 +1,24 @@ +**Added:** + +* Added support for user-defined Boresch restraints in the + ABFE Protocol (`PR #2019 `_). + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* From 1fb2991b85108379a6510bd64ac303a97228e76d Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 13:43:00 +0100 Subject: [PATCH 09/31] Add support for smcs as host molecules --- .../openmm_afe/equil_binding_afe_method.py | 8 ++++-- .../openmm_abfe/test_abfe_validation.py | 26 +++++++++++++++++-- 2 files changed, 30 insertions(+), 4 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py index d62a21ac2..d7b14cea2 100644 --- a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -278,8 +278,12 @@ 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) + # Check if there is a suitable SmallMoleculeComponent that could + # be the host molecule instead. + smcs = stateA.get_components_of_type(SmallMoleculeComponent) + if all(smc in stateA.component_diff(stateB)[0] for smc in smcs): + errmsg = "No suitable host molecule found" + raise ValueError(errmsg) if not (stateA.contains(SolventComponent) and stateB.contains(SolventComponent)): errmsg = "No SolventComponent found" diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py index 378702434..5d0a7844a 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py @@ -121,16 +121,38 @@ def test_validate_no_protcomp( stateB = ChemicalSystem( { - "benzene": benzene_modifications["benzene"], "solvent": SolventComponent(), } ) - errmsg = "No ProteinComponent found" + errmsg = "No suitable host molecule found" with pytest.raises(ValueError, match=errmsg): AbsoluteBindingProtocol._validate_endstates(stateA, stateB) +def test_validate_smc_host( + benzene_modifications, +): + """ + Should pass if there's another smc present that could be a host. + """ + stateA = ChemicalSystem( + { + "benzene": benzene_modifications["benzene"], + "toluene": benzene_modifications["toluene"], + "solvent": SolventComponent()} + ) + + stateB = ChemicalSystem( + { + "toluene": benzene_modifications["toluene"], + "solvent": SolventComponent(), + } + ) + + AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + + def test_validate_endstates_nosolvcomp_stateA(benzene_modifications, T4_protein_component): stateA = ChemicalSystem( { From f56d0ae09460389ea2803c7379d1e16753d35a3d Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 13:46:09 +0100 Subject: [PATCH 10/31] Add news item --- news/abfe_host_smc.rst | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 news/abfe_host_smc.rst diff --git a/news/abfe_host_smc.rst b/news/abfe_host_smc.rst new file mode 100644 index 000000000..b7c31a871 --- /dev/null +++ b/news/abfe_host_smc.rst @@ -0,0 +1,25 @@ +**Added:** + +* Add support for SmallMoleculeComponents to act + as host in the ``AbsoluteBindingProtocol`` + (`PR 2020 `_). + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* From 3f08aaf5012e0a4dc66101dd97e9d13f275f5ee6 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 13:50:06 +0100 Subject: [PATCH 11/31] assign None if needed --- src/openfe/protocols/openmm_afe/abfe_units.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 72f8352ad..7714137b5 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -254,8 +254,8 @@ def _get_boresch_restraint( guest_rdmol=guest_rdmol, guest_idxs=guest_atom_ids, host_idxs=host_atom_ids, - guest_restraint_atoms_idxs=list(settings.guest_restraint_ids), - host_restraint_atoms_idxs=list(settings.host_restraint_ids), + guest_restraint_atoms_idxs=list(settings.guest_restraint_ids) if settings.guest_restraint_ids is not None else None, + host_restraint_atoms_idxs=list(settings.host_restraint_ids) if settings.host_restraint_ids is not None else None, host_selection=settings.host_selection, anchor_finding_strategy=settings.anchor_finding_strategy, dssp_filter=settings.dssp_filter, From 44931cb9616415aadf0364a03d07e6ad0371ee87 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 17 Jun 2026 13:11:48 +0000 Subject: [PATCH 12/31] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/openfe/protocols/openmm_afe/abfe_units.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 7714137b5..3969db0cd 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -254,8 +254,12 @@ def _get_boresch_restraint( guest_rdmol=guest_rdmol, guest_idxs=guest_atom_ids, host_idxs=host_atom_ids, - guest_restraint_atoms_idxs=list(settings.guest_restraint_ids) if settings.guest_restraint_ids is not None else None, - host_restraint_atoms_idxs=list(settings.host_restraint_ids) if settings.host_restraint_ids is not None else None, + guest_restraint_atoms_idxs=list(settings.guest_restraint_ids) + if settings.guest_restraint_ids is not None + else None, + host_restraint_atoms_idxs=list(settings.host_restraint_ids) + if settings.host_restraint_ids is not None + else None, host_selection=settings.host_selection, anchor_finding_strategy=settings.anchor_finding_strategy, dssp_filter=settings.dssp_filter, From 21f9e9bb4c715b224e1ab9e15a3e2e2570541ded Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 17 Jun 2026 13:13:05 +0000 Subject: [PATCH 13/31] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py index 5d0a7844a..e8f7108c3 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py @@ -140,7 +140,8 @@ def test_validate_smc_host( { "benzene": benzene_modifications["benzene"], "toluene": benzene_modifications["toluene"], - "solvent": SolventComponent()} + "solvent": SolventComponent(), + } ) stateB = ChemicalSystem( From 0e6c1e05f958e1ad06594b66a669c750c03a31bf Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 14:22:26 +0100 Subject: [PATCH 14/31] Add settings to control alchemical ion --- .../openmm_afe/equil_afe_settings.py | 26 ++++++++++++++++++- .../openmm_afe/equil_binding_afe_method.py | 4 +-- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index 5853b4f76..17b8c0fe7 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -17,11 +17,15 @@ """ import numpy as np +from openff.units import unit as offunit from gufe.settings import ( OpenMMSystemGeneratorFFSettings, SettingsBaseModel, ThermoSettings, ) +from gufe.settings.typing import ( + NanometerQuantity, +) from pydantic import field_validator from openfe.protocols.openmm_utils.omm_settings import ( @@ -38,6 +42,7 @@ from openfe.protocols.restraint_utils.settings import ( BaseRestraintSettings, BoreschRestraintSettings, + SpringConstantLinearQuantity, ) @@ -87,6 +92,25 @@ class AlchemicalSettings(SettingsBaseModel): """ +class ABFEAlchemicalSettings(AlchemicalSettings): + # Dev note: we make a separate class for ABFEs so that SepTop and AHFE can + # use the parent class without having ABFE-specific net charge 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.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): """Lambda schedule settings. @@ -383,7 +407,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/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py index d7b14cea2..6f2d6db72 100644 --- a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -46,7 +46,7 @@ ABFEBoreschRestraintSettings, ABFEPreEquilOutputSettings, AbsoluteBindingSettings, - AlchemicalSettings, + ABFEAlchemicalSettings, IntegratorSettings, LambdaSettings, MDSimulationSettings, @@ -134,7 +134,7 @@ def _default_settings(cls): temperature=298.15 * offunit.kelvin, pressure=1 * offunit.bar, ), - alchemical_settings=AlchemicalSettings(), + alchemical_settings=ABFEAlchemicalSettings(), solvent_lambda_settings=LambdaSettings( lambda_elec=[ 0.0, 0.25, 0.5, 0.75, 1.0, From da65918181073591eb5978e04807a876081efca2 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 21:04:50 +0100 Subject: [PATCH 15/31] First pass at adding net charge support for ABFEs --- news/net_charge_abfe.rst | 25 + src/openfe/protocols/openmm_afe/abfe_units.py | 503 +++++++++++++++--- .../protocols/openmm_afe/base_afe_units.py | 149 ++++-- .../protocols/openmm_septop/base_units.py | 4 +- src/openfe/protocols/openmm_septop/utils.py | 85 --- src/openfe/protocols/openmm_utils/states.py | 125 +++++ 6 files changed, 681 insertions(+), 210 deletions(-) create mode 100644 news/net_charge_abfe.rst delete mode 100644 src/openfe/protocols/openmm_septop/utils.py create mode 100644 src/openfe/protocols/openmm_utils/states.py diff --git a/news/net_charge_abfe.rst b/news/net_charge_abfe.rst new file mode 100644 index 000000000..7adbffaaa --- /dev/null +++ b/news/net_charge_abfe.rst @@ -0,0 +1,25 @@ +**Added:** + +* Add support for transformations involving net + charges in the ABFE Protocol + (`PR #2020 `_). + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 3969db0cd..3190c0a8e 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -31,6 +31,7 @@ ) from openfe.protocols.openmm_utils import system_validation from openfe.protocols.restraint_utils import geometry +from openfe.protocols.restraint_utils.geometry.utils import FindHostAtoms, get_central_atom_idx from openfe.protocols.restraint_utils.geometry.boresch import BoreschRestraintGeometry from openfe.protocols.restraint_utils.openmm import omm_restraints from openfe.protocols.restraint_utils.openmm.omm_restraints import BoreschRestraint @@ -44,6 +45,296 @@ logger = logging.getLogger(__name__) +def _get_mda_universe( + topology: omm_topology, + positions: ommunit.Quantity | None, + trajectory: pathlib.Path | None, +) -> 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 | None + The System's current positions. + Used if a trajectory file is None or is not a file. + trajectory: pathlib.Path | None + 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: + if positions is None: + raise ValueError("No positions to create the Universe with") + + # 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: Iterable[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 : Iterable[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 + + +def _get_minimum_image_distance(box_dimensions: npt.NDArray) -> Quantity: + """ + Get the minimum image distance using the minimum perpendicular width + of the triclinic vectors. + + Parameters + ---------- + box_dimensions : npt.NDArray + The box dimensions as obtained from an MDAnalysis Universe. + + Returns + ------- + openff.units.Quantity + The minimum perpendicular width in units of Angstrom. + + Acknowledgements + ---------------- + Originally contributed by Bendict Tan (@aqemia-benedict-tan). + """ + from MDAnalysis.lib import mdamath + + box_vectors = mdamath.triclinic_vectors(box_dimensions) + + # Calculate the volume based on the scalar triple product + volume = mdamath.stp(box_vectors[0], box_vectors[1], box_vectors[2]) + + # Now calculate the perpendicular widths using perp_width_i = Volume / Area_of_face_i + # Where Area_of_face_i is |b_{i+1} × b_{i+2}|. + areas = np.cross(bvecs[[1, 2, 0]], bvecs[[2, 0, 1]]) + perp_widths = vol / np.linalg.norm(face_normals, axis=1) + + return perp_widths.min() * offunit.angstrom + + +def _find_most_common_ions( + openmm_topology: app.Topology, + openmm_system: System, + target_charge: int, +) -> list[int] | None: + """ + Get the most common ions of a given net charge in a system. + + Parameters + ---------- + openmm_topology : openmm.app.Topology + The Topology of the OpenMM System. + openmm_system : openmm.System + The OpenMM System. + target_charge : int + The charge the ion should have. + + Returns + ------- + list[int] | None + If present, the list of indices matching the most common ion. + + Notes + ----- + This is similar to what is done in ``_get_ion_parameters`` in + :mod:`openfe.protocols.openmm_rfe._rfe_utils.topologyhelpers`. + """ + from collections import Counter, defaultdict + + nbf = [i for i in openmm_system.getForces() if isinstance(i, NonbondedForce)][0] + + ion_counts: Counter = Counter() + ion_atom_indices: dict[str, int] = defaultdict(list) + + for residue in topology.residues(): + atoms = list(residue.atoms()) + + # We are only interested in single atom counterions + if len(atoms) != 1: + continue + + charge, _, _ = nbf.getParticleParameters(atoms[0].index) + charge_val = charge.value_in_unit(omm_unit.elementary_charge) + + if np.isclose(charge_val, target_charge, atol=0.01): + ion_counts[residue.name] += 1 + ion_atom_indices[residue.name].append(atoms[0].index) + + if ion_counts: + best_resname = ion_counts.most_common(1)[0][0] + return ion_atom_indices[best_resname] + else: + return None + + +class ABFESetupUnitMixin: + """ + Mixin for common class methods between Units + """ + def _get_alchemical_ions( + self, + alchem_comps: dict[str, list[Component]], + comp_resids: dict[Component, npt.NDArray], + openmm_topology: omm_topology, + openmm_system: System, + positions: ommunit.Quantity, + settings: dict[str, SettingsBaseModel], + dry: bool, + ) -> list[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. + openmm_topology : openmm.app.Topology + The OpenMM Topology of the system. + openmm_system : openmm.System + The OpenMM System to work on. + 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. + dry: bool + ``True`` if we are dry-running. + + Returns + ------- + list[int] + The indices of the alchemical ions. + + """ + 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 + + # TODO: For now, let's stick with a single -1/+1 case, but we should expand to more + if abs(total_charge) > 1: + errmsg = "Cannot handle net charge correction on charges greater than one" + raise ValueError(errmsg) + + # Get the indices of the most common ion type that can act as a counterion + ion_indices = _find_most_common_ions(openmm_topology, openmm_system, -total_charge) + + if ion_indices is None: + errmsg = "No suitable ions could be found to act as counterion in the system" + raise ValueError(errmsg) + + univ = _get_mda_universe( + omm_topology, + positions, + self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, + ) + + # get an atomgroup of the possible alchemical ions + ion_atomgroup = univ.atoms[ion_indices] + + # 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_atomgroup = univ.atoms[alchem_idxs] + + # Get the maximum distance we can use to find ions + univ.trajectory[-1] # use the box dimensions from the last frame + box = univ.dimensions + + if box is None or np.all(np.isinfinite(box)) or np.any(box[:3] <= 0.0): + # If it's not a dry simulation then error out + if not dry: + errmsg = ( + f"Invalid box for co-alchemical ion search: {box}" + ) + raise ValueError(errmsg) + + # For a dry execution, just assign a super high value + max_search_distance= 999 * offunit.nanometer + + # Set the max search distance to half the smallest perpendicular width + # with a 1 Angstrom padding + max_search_distance = (_get_minimum_image_distance(box) * 0.5) - 1 * offunit.angstrom + + # Re-using a utility from the restraints utilities + # TODO: rename this class! + atom_finder = FindHostAtoms( + host_atoms=ions_atomgroup, + guest_atoms=alchem_atomgroup, + min_search_distance=settings['alchemical_settings'].alchemical_ion_min_distance, + max_search_distance=max_search_distance, + ) + + # 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] + + class ComplexComponentsMixin: def _get_components(self): """ @@ -123,7 +414,7 @@ def _get_settings(self) -> dict[str, SettingsBaseModel]: return settings -class ABFEComplexSetupUnit(ComplexComponentsMixin, ComplexSettingsMixin, BaseAbsoluteSetupUnit): +class ABFEComplexSetupUnit(ABFESetupUnitMixin, ComplexComponentsMixin, ComplexSettingsMixin, BaseAbsoluteSetupUnit): """ Setup unit for the complex phase of absolute binding free energy transformations. @@ -131,86 +422,6 @@ class ABFEComplexSetupUnit(ComplexComponentsMixin, ComplexSettingsMixin, BaseAbs simtype = "complex" - @staticmethod - def _get_mda_universe( - topology: omm_topology, - positions: ommunit.Quantity | None, - trajectory: pathlib.Path | None, - ) -> 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 | None - The System's current positions. - Used if a trajectory file is None or is not a file. - trajectory: pathlib.Path | None - 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: - if positions is None: - raise ValueError("No positions to create the Universe with") - - # 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: Iterable[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 : Iterable[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, @@ -281,6 +492,7 @@ def _add_restraints( alchem_comps: dict[str, list[Component]], comp_resids: dict[Component, npt.NDArray], settings: dict[str, SettingsBaseModel], + alchemical_ions: list[int] | None, ) -> tuple[ Quantity, System, @@ -312,6 +524,8 @@ def _add_restraints( settings : dict[str, SettingsBaseModel] A dictionary of settings that defines how to find and set the restraint. + alchemical_ions : list[int] | None + The alchemical ion indices, if they exist. Returns ------- @@ -321,6 +535,10 @@ def _add_restraints( A copy of the System with the restraint added. rest_geom : geometry.HostGuestRestraintGeometry The restraint Geometry object. + + TODO + ---- + Add a restraint between the alchemical ion and the guest molecule? """ if self.verbose: self.logger.info("Generating restraints") @@ -340,7 +558,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 @@ -349,13 +567,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, @@ -389,7 +607,7 @@ def _add_restraints( restraint.add_force( thermodynamic_state, rest_geom, - controlling_parameter_name="lambda_restraints", + controlling_parameter_name="lambda_restraints_A", ) # Get the standard state correction as a unit.Quantity correction = restraint.get_standard_state_correction( @@ -497,7 +715,7 @@ def _get_settings(self) -> dict[str, SettingsBaseModel]: return settings -class ABFESolventSetupUnit(SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteSetupUnit): +class ABFESolventSetupUnit(ABFESetupUnitMixin, SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteSetupUnit): """ Setup unit for the solvent phase of absolute binding free energy transformations. @@ -505,6 +723,119 @@ class ABFESolventSetupUnit(SolventComponentsMixin, SolventSettingsMixin, BaseAbs simtype = "solvent" + 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], + alchemical_ions: list[int] | None, + ) -> tuple[ + 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. + alchemical_ions : list[int] | None + The alchemical ion indices, if they exist. + + Returns + ------- + correction : openff.units.Quantity | None + The standard state correction for the restraint. + system : openmm.System + A copy of the System with the restraint added. + rest_geom : geometry.HostGuestRestraintGeometry | None + The restraint Geometry object. + + TODO + ---- + Expand to support restraining multiple ions. + """ + if alchemical_ions is None: + return None, None, system, None + + if len(alchemical_ions) > 1: + errmsg = "Currentlyl cannot handle more than one alchemical ion" + raise ValueError(errmsg) + + restrained_system = deepcopy(system) + + if self.verbose: + self.logger.info("Generating restraints for alchemical ions") + + universe = _get_mda_universe( + topology, + positions, + self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, + ) + + # alchemical ion atom atomgroup + alchem_ion_ag = universe.atoms[alchemical_ions] + + # 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=topology, residxs=residxs) + ligand_central_atom = ligand_alchem_idxs[get_central_atom_idx(ligand_rdmol)] + ligand_central_atom_ag = universe.atoms[ligand_central_atom] + + # Get the ligand-ion distance based on the final frame + universe.trajectory[-1] + + distance = float( + calc_bonds( + alchem_ion_ag.position, + ligand_central_atom_ag.position, + box=universe.dimensions, + ) + ) + + spring_constant = to_openmm( + settings["alchemical_settings"].alchemical_ion_solvent_spring_constant + ) + + force = HarmonicBondForce() + + force.addBond( + ligand_central_atom, + alchemical_ions[0], + distance * ommunit.angstrom, + spring_constant, + ) + + force.setName("ion_restraint") + add_force_in_separate_group(restrained_system, force) + + return None, None, restrained_system, None + class ABFESolventSimUnit( SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteMultiStateSimulationUnit diff --git a/src/openfe/protocols/openmm_afe/base_afe_units.py b/src/openfe/protocols/openmm_afe/base_afe_units.py index 676ba4c51..7dfe2c375 100644 --- a/src/openfe/protocols/openmm_afe/base_afe_units.py +++ b/src/openfe/protocols/openmm_afe/base_afe_units.py @@ -46,7 +46,6 @@ from openmmtools.alchemy import ( AbsoluteAlchemicalFactory, AlchemicalRegion, - AlchemicalState, ) from openmmtools.states import ( GlobalParameterState, @@ -85,6 +84,10 @@ make_vec3_box, serialize, ) +from openfe.protocols.openmm_utils.states import ( + SingleRegionAlchemicalState, + DualRegionAlchemicalState, +) from openfe.protocols.restraint_utils import geometry from openfe.protocols.restraint_utils.openmm import omm_restraints from openfe.utils import log_system_probe, without_oechem_backend @@ -534,6 +537,25 @@ def _get_omm_objects( return topology, system, positions, comp_resids + def _get_alchemical_ions( + self, + alchemical_components: dict[str, list[Component]], + comp_resids: dict[Component, npt.NDArray], + openmm_topology: app.Topology, + openmm_system: openmm.System, + positions: openmm.unit.Quantity, + settings: dict[str, SettingsBaseModel], + dry: bool, + ) -> list[int] | None: + """ + Placeholder method for finding alchemical ions. + + Returns + ------- + None + """ + None + def _add_restraints( self, system: openmm.System, @@ -542,9 +564,10 @@ def _add_restraints( alchem_comps: dict[str, list[Component]], comp_resids: dict[Component, npt.NDArray], settings: dict[str, SettingsBaseModel], + alchemical_ions : list[int] | None, ) -> tuple[ Quantity | None, - openmm.System | None, + openmm.System, geometry.BaseRestraintGeometry | None, ]: """ @@ -557,7 +580,8 @@ def _get_alchemical_system( topology: app.Topology, system: openmm.System, comp_resids: dict[Component, npt.NDArray], - alchem_comps: dict[str, list[Component]], + alchemical_components: dict[str, list[Component]], + alchemical_ions: list[int] | None, alchemical_settings: AlchemicalSettings, ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, list[int]]: """ @@ -571,8 +595,10 @@ def _get_alchemical_system( System to alchemically modify. comp_resids : dict[str, npt.NDArray] A dictionary of residues for each component in the System. - alchem_comps : dict[str, list[Component]] + alchemical_components : dict[str, list[Component]] A dictionary of alchemical components for each end state. + alchemical_ions : list[int] | None + List of indices for any alchemical ions, there are any. alchemical_settings : AlchemicalSettings Settings controlling how the alchemical system is built. @@ -582,29 +608,41 @@ def _get_alchemical_system( Factory for creating an alchemically modified system. alchemical_system : openmm.System Alchemically modified system - alchemical_indices : list[int] - A list of atom indices for the alchemically modified - species in the system. + alchemical_indices : dict[str, list[int]] + A dictionary containing a list of atom indices + for each independent alchemically modified species in the system. TODO ---- * Add support for all alchemical factory options """ - alchemical_indices = self._get_alchemical_indices(topology, comp_resids, alchem_comps) - - alchemical_region = AlchemicalRegion( - alchemical_atoms=alchemical_indices, - softcore_alpha=alchemical_settings.softcore_alpha, - annihilate_electrostatics=True, - annihilate_sterics=alchemical_settings.annihilate_sterics, - softcore_a=alchemical_settings.softcore_a, - softcore_b=alchemical_settings.softcore_b, - softcore_c=alchemical_settings.softcore_c, - softcore_beta=0.0, - softcore_d=1.0, - softcore_e=1.0, - softcore_f=2.0, - ) + # first region is the alchemically changing species + alchemical_indices = { + 'A': self._get_alchemical_indices(topology, comp_resids, alchemical_components) + } + # second region is any alchemical ions + if alchemical_ions is not None: + alchemical_indices['B'] = alchemical_ions + + alchemical_regions = [] + + for region in alchemical_indices: + alchemical_regions.append( + AlchemicalRegion( + alchemical_atoms=alchemical_indices[region], + name=region, + softcore_alpha=alchemical_settings.softcore_alpha, + annihilate_electrostatics=True, + annihilate_sterics=alchemical_settings.annihilate_sterics, + softcore_a=alchemical_settings.softcore_a, + softcore_b=alchemical_settings.softcore_b, + softcore_c=alchemical_settings.softcore_c, + softcore_beta=0.0, + softcore_d=1.0, + softcore_e=1.0, + softcore_f=2.0, + ) + ) alchemical_factory = AbsoluteAlchemicalFactory( consistent_exceptions=False, @@ -614,7 +652,10 @@ def _get_alchemical_system( disable_alchemical_dispersion_correction=alchemical_settings.disable_alchemical_dispersion_correction, split_alchemical_forces=True, ) - 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 @@ -715,6 +756,17 @@ def run( omm_system, omm_topology, positions, settings, dry ) + # Get alchemical ions, if needed / allowed + alchem_ion = self._get_alchemical_ions( + alchemical_components=alchem_comps, + comp_resids=comp_resids, + openmm_topology=omm_topology, + openmm_system=omm_system, + positions=positions, + settings=settings, + dry=dry, + ) + # Add restraints # Note: when no restraint is applied, restrained_omm_system == omm_system ( @@ -728,6 +780,7 @@ def run( alchem_comps, comp_resids, settings, + alchem_ion, ) # Get alchemical system @@ -735,7 +788,8 @@ def run( topology=omm_topology, system=restrained_omm_system, comp_resids=comp_resids, - alchem_comps=alchem_comps, + alchemical_components=alchem_comps, + alchemical_ions=alchem_ions, alchemical_settings=settings["alchemical_settings"], ) @@ -767,6 +821,7 @@ def run( "pdb_structure": pdb_structure, "selection_indices": selection_indices, "box_vectors": from_openmm(box_vectors), + "alchemical_indices": alchem_indices, } if standard_state_corr is not None: @@ -786,7 +841,6 @@ def run( "standard_system": omm_system, "restrained_system": restrained_omm_system, "alchem_system": alchem_system, - "alchem_indices": alchem_indices, "alchem_factory": alchem_factory, "debug_positions": positions, } @@ -919,6 +973,7 @@ def _get_states( lambdas: dict[str, list[float]], solvent_component: BaseSolventComponent | None, alchemically_restrained: bool, + alchemical_indices: dict[str, list[int]], ) -> tuple[list[SamplerState], list[ThermodynamicState]]: """ Get a list of sampler and thermodynmic states from an @@ -941,6 +996,9 @@ def _get_states( alchemically_restrained : bool Whether or not the system requires a control parameter for any alchemical restraints. + alchemical_indices : dict[str, list[int]] + Dictionary of alchemical indices for each alchemical + region in the system. Returns ------- @@ -950,7 +1008,13 @@ 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) == 1: + alchemical_state = SingleAlchemicalState.from_system(alchemical_system) + elif len(alchemical_indices) == 2: + alchemical_state = DoubleAclhemicalState.from_system(alchemical_system) + else: + errmsg = "More than two alchemical regions are not supported" + raise ValueError(errmsg) # Set up the system constants temperature = thermodynamic_settings.temperature @@ -962,24 +1026,29 @@ def _get_states( constants["pressure"] = ensure_quantity(pressure, "openmm") # Get the thermodynamic parameter protocol - param_protocol = copy.deepcopy(lambdas) + # We populate the protocol from lambdas based on the number of regions + # we have. + def _add_lambdas_to_protocol(protocol, lambdas, region_name): + for key in ["lambda_electrostatics", "lambda_sterics"]: + protocol[f"{key}_{region_name}"] = copy.deepcopy(lambdas[key]) - # Get the composable states - if alchemically_restrained: - restraint_state = omm_restraints.RestraintParameterState(lambda_restraints=1.0) - composable_states = [alchemical_state, restraint_state] - else: - composable_states = [alchemical_state] + param_protocol = {} + # Always add for the first region + _add_lambdas_to_protocol(param_protocol, lambdas, "A") + + # If we have two regions (i.e. an alchemical ion) add + if len(alchemical_indices) == 2: + _add_lambdas_to_protocol(param_protocol, lambdas, "A") - # In this case we also don't have a restraint being controlled - # so we drop it from the protocol - param_protocol.pop("lambda_restraints", None) + # Only the first region (ligand) can be restrained + if alchemically_restrained: + param_protocol["lambda_restraints_A"] = copy.deepcopy(lambdas["lambda_restraints"]) cmp_states = create_thermodynamic_state_protocol( alchemical_system, protocol=param_protocol, constants=constants, - composable_states=composable_states, + composable_states=[alchemical_state], ) sampler_state = SamplerState(positions=positions) @@ -1355,6 +1424,7 @@ def run( positions: openmm.unit.Quantity, box_vectors: Quantity, selection_indices: npt.NDArray, + alchemical_indices: dict[str, list[int]], alchemical_restraints: bool, dry: bool = False, verbose: bool = True, @@ -1374,6 +1444,8 @@ def run( The box vectors of the System. selection_indices : npt.NDArray Indices of the System particles to write to file. + alchemical_indices : dict[str, list[int]] + Dictionary of alchemical indices for each alchemical region. alchemical_restraints: bool, Whether or not the system has alchemical restraints. dry: bool @@ -1434,6 +1506,7 @@ def run( lambdas=lambdas, solvent_component=solv_comp, alchemically_restrained=alchemical_restraints, + alchemical_indices=alchemical_indices, ) # Get the integrator @@ -1528,6 +1601,7 @@ def _execute( system = deserialize(setup_results.outputs["system"]) positions = to_openmm(np.load(setup_results.outputs["positions"]) * offunit.nanometer) selection_indices = setup_results.outputs["selection_indices"] + alchemical_indices = setup_rersults.outputs["alchemical_indices"] box_vectors = setup_results.outputs["box_vectors"] if setup_results.outputs["restraint_geometry"] is not None: @@ -1540,6 +1614,7 @@ def _execute( positions=positions, box_vectors=box_vectors, selection_indices=selection_indices, + alchemical_indices=alchemical_indices, alchemical_restraints=alchemical_restraints, scratch_basepath=ctx.scratch, shared_basepath=ctx.shared, diff --git a/src/openfe/protocols/openmm_septop/base_units.py b/src/openfe/protocols/openmm_septop/base_units.py index bf19f8dd5..f64b84ca4 100644 --- a/src/openfe/protocols/openmm_septop/base_units.py +++ b/src/openfe/protocols/openmm_septop/base_units.py @@ -86,7 +86,7 @@ system_validation, ) from ..openmm_utils.mdtraj_utils import mdtraj_from_openmm -from .utils import SepTopParameterState +from openfe.protocols.openmm_utils.states import DualRegionAlchemicalState logger = logging.getLogger(__name__) @@ -879,7 +879,7 @@ def _get_states( cmp_states : list[ThermodynamicState] A list of ThermodynamicState for each replica in the system. """ - alchemical_state = SepTopParameterState.from_system(alchemical_system) + alchemical_state = DualRegionAlchemicalState.from_system(alchemical_system) # Set up the system constants temperature = settings["thermo_settings"].temperature diff --git a/src/openfe/protocols/openmm_septop/utils.py b/src/openfe/protocols/openmm_septop/utils.py deleted file mode 100644 index 38d7f1d34..000000000 --- a/src/openfe/protocols/openmm_septop/utils.py +++ /dev/null @@ -1,85 +0,0 @@ -from openmmtools import states -from openmmtools.states import GlobalParameterState - - -class SepTopParameterState(GlobalParameterState): - """ - Composable state to control lambda parameters for two ligands. - See :class:`openmmtools.states.GlobalParameterState` for more details. - Parameters - ---------- - parameters_name_suffix : Optional[str] - If specified, the state will control a modified version of the parameter - ``lambda_restraints_{parameters_name_suffix}` instead of just - ``lambda_restraints``. - lambda_sterics_A : Optional[float] - The value for the vdW interactions for ligand A. - If defined, must be between 0 and 1. - lambda_electrosterics_A : Optional[float] - The value for the electrostatics interactions for ligand A. - If defined, must be between 0 and 1. - lambda_restraints_A : Optional[float] - The strength of the restraint for ligand A. - If defined, must be between 0 and 1. - lambda_bonds_A : Optional[float] - The value for modifying bonds for ligand A. - If defined, must be between 0 and 1. - lambda_angles_A : Optional[float] - The value for modifying angles for ligand A. - If defined, must be between 0 and 1. - lambda_dihedrals_A : Optional[float] - The value for modifying dihedrals for ligand A. - If defined, must be between 0 and 1. - lambda_sterics_B : Optional[float] - The value for the vdW interactions for ligand B. - If defined, must be between 0 and 1. - lambda_electrosterics_B : Optional[float] - The value for the electrostatics interactions for ligand B. - If defined, must be between 0 and 1. - lambda_restraints_B : Optional[float] - The strength of the restraint for ligand B. - If defined, must be between 0 and 1. - lambda_bonds_B : Optional[float] - The value for modifying bonds for ligand B. - If defined, must be between 0 and 1. - lambda_angles_B : Optional[float] - The value for modifying angles for ligand B. - If defined, must be between 0 and 1. - lambda_dihedrals_B : Optional[float] - The value for modifying dihedrals for ligand B. - If defined, must be between 0 and 1. - """ - - class _LambdaParameter(states.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 parameters for ligand A - lambda_sterics_A = _LambdaParameter("lambda_sterics_A") - lambda_electrostatics_A = _LambdaParameter("lambda_electrostatics_A") - lambda_restraints_A = _LambdaParameter("lambda_restraints_A") - lambda_bonds_A = _LambdaParameter("lambda_bonds_A") - lambda_angles_A = _LambdaParameter("lambda_angles_A") - lambda_torsions_A = _LambdaParameter("lambda_torsions_A") - - # Lambda parameters for ligand B - lambda_sterics_B = _LambdaParameter("lambda_sterics_B") - lambda_electrostatics_B = _LambdaParameter("lambda_electrostatics_B") - lambda_restraints_B = _LambdaParameter("lambda_restraints_B") - lambda_bonds_B = _LambdaParameter("lambda_bonds_B") - lambda_angles_B = _LambdaParameter("lambda_angles_B") - lambda_torsions_B = _LambdaParameter("lambda_torsions_B") - - # # Restraints solvent - # lambda_restraints = _LambdaParameter('lambda_restraints') diff --git a/src/openfe/protocols/openmm_utils/states.py b/src/openfe/protocols/openmm_utils/states.py new file mode 100644 index 000000000..1d07fd8da --- /dev/null +++ b/src/openfe/protocols/openmm_utils/states.py @@ -0,0 +1,125 @@ +# This code is part of OpenFE and is licensed under the MIT license. +# For details, see https://github.com/OpenFreeEnergy/openfe +# +# Acknowledgements: +# This module derives heavily from the AlchemicalState class +# in openmmtools.alchemy (https://github.com/choderalab/openmmtools). + +from openmmtools.states import 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) + + +class SingleRegionAlchemicalState(GlobalParameterState): + """ + Composable state to control lambda parameters for a single + alchemical molecule / region (``ligand_A``). + + Parameters + ---------- + parameters_name_suffix : str | None + If specified, the state will control a modified version of the parameter + ``lambda_restraints_{parameters_name_suffix}` instead of just + ``lambda_restraints``. + lambda_sterics_A : float | None + Control parameter for the vdW interactions for ligand A. + If defined, must be between 0 and 1. + lambda_electrostatics_A : float | None + Control parameter for the electrostatics interactions for ligand A. + If defined, must be between 0 and 1. + lambda_bonds_A : float | None + Control parameter for alchemically modified bonds for ligand A. + If defined, must be between 0 and 1. + lambda_angles_A : float | None + Control parameter for alchemically modified angles for ligand A. + If defined, must be between 0 and 1. + lambda_torsions_A : float | None + Control parameter for alchemically modified dihedrals for ligand A. + If defined, must be between 0 and 1. + lambda_restraints_A : float | None + Control parameter for alchmemically modified restraints for ligand A. + + See Also + -------- + :class:`openmmtools.states.GlobalParameterState` + :class:`openfe.protocols.restraint_utils.geometry.DualRegionAlchemicalState` + """ + 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") + lambda_restraints_A = _LambdaParameter("lambda_restraints_A") + + +class DualRegionAlchemicalState(SingleRegionAlchemicalState): + """ + Composable state to control lambda parameters for a system + with two alchemical molecules / regions (``ligand_A`` and ``ligand_B``). + + Parameters + ---------- + parameters_name_suffix : str | None + If specified, the state will control a modified version of the parameter + ``lambda_restraints_A_{parameters_name_suffix}` instead of just + ``lambda_restraints_A``. + lambda_sterics_A : float | None + Control parameter for the vdW interactions for ligand A. + If defined, must be between 0 and 1. + lambda_electrostatics_A : float | None + Control parameter for the electrostatics interactions for ligand A. + If defined, must be between 0 and 1. + lambda_bonds_A : float | None + Control parameter for alchemically modified bonds for ligand A. + If defined, must be between 0 and 1. + lambda_angles_A : float | None + Control parameter for alchemically modified angles for ligand A. + If defined, must be between 0 and 1. + lambda_torsions_A : float | None + Control parameter for alchemically modified dihedrals for ligand A. + If defined, must be between 0 and 1. + lambda_restraints_A : float | None + Control parameter for alchmemically modified restraints for ligand A. + lambda_sterics_B : float | None + Control parameter for the vdW interactions for ligand B. + If defined, must be between 0 and 1. + lambda_electrostatics_B : float | None + Control parameter for the electrostatics interactions for ligand B. + If defined, must be between 0 and 1. + lambda_bonds_B : float | None + Control parameter for alchemically modified bonds for ligand B. + If defined, must be between 0 and 1. + lambda_angles_B : float | None + Control parameter for alchemically modified angles for ligand B. + If defined, must be between 0 and 1. + lambda_torsions_B : float | None + Control parameter for alchemically modified dihedrals for ligand B. + If defined, must be between 0 and 1. + lambda_restraints_B : float | None + Control parameter for alchmemically modified restraints for ligand B. + + See Also + -------- + :class:`openmmtools.states.GlobalParameterState` + :class:`openfe.protocols.restraint_utils.geometry.SingleRegionAlchemicalState` + """ + 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") + lambda_restraints_B = _LambdaParameter("lambda_restraints_B") From db7567dd703a857d935aae39d0042d970bf3de23 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 21:19:19 +0100 Subject: [PATCH 16/31] various bug fixes --- src/openfe/protocols/openmm_afe/abfe_units.py | 20 +++++++++---------- .../protocols/openmm_afe/base_afe_units.py | 10 +++++----- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 3190c0a8e..8dbbd2a2a 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -152,8 +152,8 @@ def _get_minimum_image_distance(box_dimensions: npt.NDArray) -> Quantity: volume = mdamath.stp(box_vectors[0], box_vectors[1], box_vectors[2]) # Now calculate the perpendicular widths using perp_width_i = Volume / Area_of_face_i - # Where Area_of_face_i is |b_{i+1} × b_{i+2}|. - areas = np.cross(bvecs[[1, 2, 0]], bvecs[[2, 0, 1]]) + # Where Area_of_face_i is |box_vectors_{i+1} × box_vectors_{i+2}|. + areas = np.cross(box_vectors[[1, 2, 0]], box_vectors[[2, 0, 1]]) perp_widths = vol / np.linalg.norm(face_normals, axis=1) return perp_widths.min() * offunit.angstrom @@ -193,7 +193,7 @@ def _find_most_common_ions( ion_counts: Counter = Counter() ion_atom_indices: dict[str, int] = defaultdict(list) - for residue in topology.residues(): + for residue in openmm_topology.residues(): atoms = list(residue.atoms()) # We are only interested in single atom counterions @@ -289,7 +289,7 @@ def _get_alchemical_ions( ) # get an atomgroup of the possible alchemical ions - ion_atomgroup = univ.atoms[ion_indices] + ions_atomgroup = univ.atoms[ion_indices] # get the alchemical atoms residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]]) @@ -310,10 +310,10 @@ def _get_alchemical_ions( # For a dry execution, just assign a super high value max_search_distance= 999 * offunit.nanometer - - # Set the max search distance to half the smallest perpendicular width - # with a 1 Angstrom padding - max_search_distance = (_get_minimum_image_distance(box) * 0.5) - 1 * offunit.angstrom + else: + # Set the max search distance to half the smallest perpendicular width + # with a 1 Angstrom padding + max_search_distance = (_get_minimum_image_distance(box) * 0.5) - 1 * offunit.angstrom # Re-using a utility from the restraints utilities # TODO: rename this class! @@ -780,7 +780,7 @@ def _add_restraints( Expand to support restraining multiple ions. """ if alchemical_ions is None: - return None, None, system, None + return None, system, None if len(alchemical_ions) > 1: errmsg = "Currentlyl cannot handle more than one alchemical ion" @@ -834,7 +834,7 @@ def _add_restraints( force.setName("ion_restraint") add_force_in_separate_group(restrained_system, force) - return None, None, restrained_system, None + return None, restrained_system, None class ABFESolventSimUnit( diff --git a/src/openfe/protocols/openmm_afe/base_afe_units.py b/src/openfe/protocols/openmm_afe/base_afe_units.py index 7dfe2c375..9056dc746 100644 --- a/src/openfe/protocols/openmm_afe/base_afe_units.py +++ b/src/openfe/protocols/openmm_afe/base_afe_units.py @@ -757,7 +757,7 @@ def run( ) # Get alchemical ions, if needed / allowed - alchem_ion = self._get_alchemical_ions( + alchem_ions = self._get_alchemical_ions( alchemical_components=alchem_comps, comp_resids=comp_resids, openmm_topology=omm_topology, @@ -1009,9 +1009,9 @@ def _get_states( """ # Fetch an alchemical state if len(alchemical_indices) == 1: - alchemical_state = SingleAlchemicalState.from_system(alchemical_system) + alchemical_state = SingleRegionAlchemicalState.from_system(alchemical_system) elif len(alchemical_indices) == 2: - alchemical_state = DoubleAclhemicalState.from_system(alchemical_system) + alchemical_state = DoubleRegionAclhemicalState.from_system(alchemical_system) else: errmsg = "More than two alchemical regions are not supported" raise ValueError(errmsg) @@ -1038,7 +1038,7 @@ def _add_lambdas_to_protocol(protocol, lambdas, region_name): # If we have two regions (i.e. an alchemical ion) add if len(alchemical_indices) == 2: - _add_lambdas_to_protocol(param_protocol, lambdas, "A") + _add_lambdas_to_protocol(param_protocol, lambdas, "B") # Only the first region (ligand) can be restrained if alchemically_restrained: @@ -1601,7 +1601,7 @@ def _execute( system = deserialize(setup_results.outputs["system"]) positions = to_openmm(np.load(setup_results.outputs["positions"]) * offunit.nanometer) selection_indices = setup_results.outputs["selection_indices"] - alchemical_indices = setup_rersults.outputs["alchemical_indices"] + alchemical_indices = setup_results.outputs["alchemical_indices"] box_vectors = setup_results.outputs["box_vectors"] if setup_results.outputs["restraint_geometry"] is not None: From 543dd7f2914c50701974f721b4c204a4813ac13b Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 21:23:26 +0100 Subject: [PATCH 17/31] some bugfixes --- src/openfe/protocols/openmm_afe/abfe_units.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 8dbbd2a2a..a703cf77d 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -154,7 +154,7 @@ def _get_minimum_image_distance(box_dimensions: npt.NDArray) -> Quantity: # Now calculate the perpendicular widths using perp_width_i = Volume / Area_of_face_i # Where Area_of_face_i is |box_vectors_{i+1} × box_vectors_{i+2}|. areas = np.cross(box_vectors[[1, 2, 0]], box_vectors[[2, 0, 1]]) - perp_widths = vol / np.linalg.norm(face_normals, axis=1) + perp_widths = volume / np.linalg.norm(areas, axis=1) return perp_widths.min() * offunit.angstrom @@ -258,11 +258,6 @@ def _get_alchemical_ions( The indices of the alchemical ions. """ - 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 From 6d2ced5ebaef0d08c9c5b349c8a9d2c32914bc98 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 21:28:47 +0100 Subject: [PATCH 18/31] more bug fixes --- src/openfe/protocols/openmm_afe/abfe_units.py | 8 ++++---- src/openfe/protocols/openmm_afe/base_afe_units.py | 6 +++--- src/openfe/protocols/openmm_utils/states.py | 12 ++++++------ 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index a703cf77d..dd2d65c3a 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -278,7 +278,7 @@ def _get_alchemical_ions( raise ValueError(errmsg) univ = _get_mda_universe( - omm_topology, + openmm_topology, positions, self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, ) @@ -288,7 +288,7 @@ def _get_alchemical_ions( # 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_idxs = _get_idxs_from_residxs(topology=openmm_topology, residxs=residxs) alchem_atomgroup = univ.atoms[alchem_idxs] # Get the maximum distance we can use to find ions @@ -327,7 +327,7 @@ def _get_alchemical_ions( raise ValueError(errmsg) # Just use the first one that comes back ok - return atom_finder.results.host_idxs[0] + return [atom_finder.results.host_idxs[0]] class ComplexComponentsMixin: @@ -778,7 +778,7 @@ def _add_restraints( return None, system, None if len(alchemical_ions) > 1: - errmsg = "Currentlyl cannot handle more than one alchemical ion" + errmsg = "Currently cannot handle more than one alchemical ion" raise ValueError(errmsg) restrained_system = deepcopy(system) diff --git a/src/openfe/protocols/openmm_afe/base_afe_units.py b/src/openfe/protocols/openmm_afe/base_afe_units.py index 9056dc746..bb3a3513e 100644 --- a/src/openfe/protocols/openmm_afe/base_afe_units.py +++ b/src/openfe/protocols/openmm_afe/base_afe_units.py @@ -554,7 +554,7 @@ def _get_alchemical_ions( ------- None """ - None + return None def _add_restraints( self, @@ -780,7 +780,7 @@ def run( alchem_comps, comp_resids, settings, - alchem_ion, + alchem_ions, ) # Get alchemical system @@ -1011,7 +1011,7 @@ def _get_states( if len(alchemical_indices) == 1: alchemical_state = SingleRegionAlchemicalState.from_system(alchemical_system) elif len(alchemical_indices) == 2: - alchemical_state = DoubleRegionAclhemicalState.from_system(alchemical_system) + alchemical_state = DualRegionAlchemicalState.from_system(alchemical_system) else: errmsg = "More than two alchemical regions are not supported" raise ValueError(errmsg) diff --git a/src/openfe/protocols/openmm_utils/states.py b/src/openfe/protocols/openmm_utils/states.py index 1d07fd8da..2edf6edf4 100644 --- a/src/openfe/protocols/openmm_utils/states.py +++ b/src/openfe/protocols/openmm_utils/states.py @@ -10,7 +10,7 @@ class _LambdaParameter(GlobalParameterState.GlobalParameter): """ - A global parameter in the interval [0,, 1] with standard value 1. + 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) @@ -51,12 +51,12 @@ class SingleRegionAlchemicalState(GlobalParameterState): Control parameter for alchemically modified dihedrals for ligand A. If defined, must be between 0 and 1. lambda_restraints_A : float | None - Control parameter for alchmemically modified restraints for ligand A. + Control parameter for alchemically modified restraints for ligand A. See Also -------- :class:`openmmtools.states.GlobalParameterState` - :class:`openfe.protocols.restraint_utils.geometry.DualRegionAlchemicalState` + :class:`openfe.protocols.openmm_utils.states.DualRegionAlchemicalState` """ lambda_sterics_A = _LambdaParameter("lambda_sterics_A") lambda_electrostatics_A = _LambdaParameter("lambda_electrostatics_A") @@ -93,7 +93,7 @@ class DualRegionAlchemicalState(SingleRegionAlchemicalState): Control parameter for alchemically modified dihedrals for ligand A. If defined, must be between 0 and 1. lambda_restraints_A : float | None - Control parameter for alchmemically modified restraints for ligand A. + Control parameter for alchemically modified restraints for ligand A. lambda_sterics_B : float | None Control parameter for the vdW interactions for ligand B. If defined, must be between 0 and 1. @@ -110,12 +110,12 @@ class DualRegionAlchemicalState(SingleRegionAlchemicalState): Control parameter for alchemically modified dihedrals for ligand B. If defined, must be between 0 and 1. lambda_restraints_B : float | None - Control parameter for alchmemically modified restraints for ligand B. + Control parameter for alchemically modified restraints for ligand B. See Also -------- :class:`openmmtools.states.GlobalParameterState` - :class:`openfe.protocols.restraint_utils.geometry.SingleRegionAlchemicalState` + :class:`openfe.protocols.openmm_utils.states.SingleRegionAlchemicalState` """ lambda_sterics_B = _LambdaParameter("lambda_sterics_B") lambda_electrostatics_B = _LambdaParameter("lambda_electrostatics_B") From fcfe3b10ee07c89f32e23f5f65111f12a94bc020 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 23:16:51 +0100 Subject: [PATCH 19/31] Add fixes for most tests --- src/openfe/protocols/openmm_afe/abfe_units.py | 10 ++--- .../openmm_afe/equil_afe_settings.py | 11 +++++ .../openmm_abfe/test_abfe_energies.py | 39 ++++++------------ .../openmm_abfe/test_abfe_protocol.py | 40 ++++++++++++++----- .../openmm_abfe/test_abfe_protocol_results.py | 2 + .../openmm_ahfe/test_ahfe_protocol.py | 9 ++++- .../openmm_ahfe/test_ahfe_protocol_results.py | 2 + .../protocols/openmm_ahfe/test_ahfe_resume.py | 9 ++++- 8 files changed, 77 insertions(+), 45 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index dd2d65c3a..048ef4156 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -160,7 +160,7 @@ def _get_minimum_image_distance(box_dimensions: npt.NDArray) -> Quantity: def _find_most_common_ions( - openmm_topology: app.Topology, + openmm_topology: omm_topology, openmm_system: System, target_charge: int, ) -> list[int] | None: @@ -220,7 +220,7 @@ class ABFESetupUnitMixin: """ def _get_alchemical_ions( self, - alchem_comps: dict[str, list[Component]], + alchemical_components: dict[str, list[Component]], comp_resids: dict[Component, npt.NDArray], openmm_topology: omm_topology, openmm_system: System, @@ -233,7 +233,7 @@ def _get_alchemical_ions( Parameters ---------- - alchem_comps: dict[str, list[Component]] + alchemical_components: dict[str, list[Component]] A dictionary with a list of alchemical components in both state A and B. comp_resids: dict[Component, npt.NDArray] @@ -258,7 +258,7 @@ def _get_alchemical_ions( The indices of the alchemical ions. """ - total_charge = alchem_comps["stateA"][0].total_charge + total_charge = alchemical_components["stateA"][0].total_charge # Don't add an alchemical ion if we have zero net charge # or we didn't request it. @@ -287,7 +287,7 @@ def _get_alchemical_ions( ions_atomgroup = univ.atoms[ion_indices] # get the alchemical atoms - residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]]) + residxs = np.concatenate([comp_resids[key] for key in alchemical_components["stateA"]]) alchem_idxs = _get_idxs_from_residxs(topology=openmm_topology, residxs=residxs) alchem_atomgroup = univ.atoms[alchem_idxs] diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index 17b8c0fe7..687cbdc6b 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -411,6 +411,17 @@ def must_be_positive(cls, v): """ Alchemical protocol settings. """ + + @field_validator("alchemical_settings", mode="before") + @classmethod + def coerce_alchemical_settings(cls, v): + """ + Migration from previous default to the new one + """ + if isinstance(v, AlchemicalSettings) and not isinstance(v, ABFEAlchemicalSettings): + return ABFEAlchemicalSettings(**v.model_dump()) + return v + complex_lambda_settings: LambdaSettings """ Settings for controlling the complex transformation leg diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py index 14539d968..a6ac366f8 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py @@ -13,28 +13,17 @@ from openmmtools.alchemy import ( AbsoluteAlchemicalFactory, AlchemicalRegion, - AlchemicalState, ) from openfe.protocols import openmm_afe from openfe.protocols.openmm_afe.abfe_units import ( ABFEComplexSetupUnit, ) +from openfe.protocols.openmm_utils.states import SingleRegionAlchemicalState from openfe.protocols.openmm_utils.omm_settings import OpenMMSolvationSettings from openfe.protocols.openmm_utils.serialization import deserialize -class AlchemStateRest(AlchemicalState): - """ - A modified AlchemicalState for testing. - - Note: we don't need this in the main protocol since we use composable - thermodynamic states. - """ - - lambda_restraints = AlchemicalState._LambdaParameter("lambda_restraints") - - def get_alchemical_energy_components(alchemical_system, alchemical_state, positions, platform): """Compute potential energy of the alchemical system by Force. @@ -166,13 +155,11 @@ def get_energy_components( platform = Platform.getPlatformByName("Reference") alchemical_region = AlchemicalRegion(alchemical_atoms=indices) - alchemical_state = AlchemStateRest.from_system( - system, parameters_name_suffix=alchemical_region.name - ) + alchemical_state = SingleRegionAlchemicalState.from_system(system) - alchemical_state.lambda_sterics = lambda_sterics - alchemical_state.lambda_electrostatics = lambda_electrostatics - alchemical_state.lambda_restraints = lambda_restraints + alchemical_state.lambda_sterics_A = lambda_sterics + alchemical_state.lambda_electrostatics_A = lambda_electrostatics + alchemical_state.lambda_restraints_A = lambda_restraints return get_alchemical_energy_components( system, alchemical_state, positions, platform=platform @@ -182,7 +169,7 @@ def get_energy_components( def test_energies_regression(self, lambda_val, t4_reference_system, t4_validation_data): energies_ref = self.get_energy_components( t4_reference_system, - t4_validation_data["alchem_indices"], + t4_validation_data["alchemical_indices"], t4_validation_data["debug_positions"], lambda_val, lambda_val, @@ -191,7 +178,7 @@ def test_energies_regression(self, lambda_val, t4_reference_system, t4_validatio energies_val = self.get_energy_components( t4_validation_data["alchem_system"], - t4_validation_data["alchem_indices"], + t4_validation_data["alchemical_indices"], t4_validation_data["debug_positions"], lambda_val, lambda_val, @@ -223,7 +210,7 @@ def assert_energies(actual, expected, nonbonded_lower: bool): # lambda 1 on all energies = self.get_energy_components( t4_validation_data["alchem_system"], - t4_validation_data["alchem_indices"], + t4_validation_data["alchemical_indices"], t4_validation_data["debug_positions"], lambda_sterics=1.0, lambda_electrostatics=1.0, @@ -236,7 +223,7 @@ def assert_energies(actual, expected, nonbonded_lower: bool): energies = self.get_energy_components( t4_validation_data["alchem_system"], - t4_validation_data["alchem_indices"], + t4_validation_data["alchemical_indices"], t4_validation_data["debug_positions"], lambda_sterics=1.0, lambda_electrostatics=1.0, @@ -249,7 +236,7 @@ def assert_energies(actual, expected, nonbonded_lower: bool): energies = self.get_energy_components( t4_validation_data["alchem_system"], - t4_validation_data["alchem_indices"], + t4_validation_data["alchemical_indices"], t4_validation_data["debug_positions"], lambda_sterics=1.0, lambda_electrostatics=0.0, @@ -261,16 +248,16 @@ def assert_energies(actual, expected, nonbonded_lower: bool): # turn off sterics expected = copy.deepcopy(energies) - expected["alchemically modified NonbondedForce for non-alchemical/alchemical sterics"] = ( + expected["alchemically modified NonbondedForce for non-alchemical/alchemical sterics for region A"] = ( 0 * ommunit.kilojoule_per_mole ) expected[ - "alchemically modified BondForce for non-alchemical/alchemical sterics exceptions" + "alchemically modified BondForce for non-alchemical/alchemical sterics exceptions for region A" ] = 0 * ommunit.kilojoule_per_mole energies = self.get_energy_components( t4_validation_data["alchem_system"], - t4_validation_data["alchem_indices"], + t4_validation_data["alchemical_indices"], t4_validation_data["debug_positions"], lambda_sterics=0.0, lambda_electrostatics=0.0, diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index 29139c285..37db3f2a4 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -30,8 +30,8 @@ from openmmtools.multistate.multistatesampler import MultiStateSampler from openmmtools.tests import test_alchemy from openmmtools.tests.test_alchemy import ( - check_interacting_energy_components, - check_noninteracting_energy_components, + check_multi_interacting_energy_components, + check_multi_noninteracting_energy_components, compare_system_energies, ) @@ -41,6 +41,7 @@ from openfe.protocols.openmm_afe import ( AbsoluteBindingProtocol, ) +from openfe.protocols.openmm_afe.abfe_units import _get_mda_universe from openfe.protocols.restraint_utils.geometry import BoreschRestraintGeometry from .utils import UNIT_TYPES, _get_units @@ -169,7 +170,7 @@ def test_mda_universe_error(): when calling the mda Universe getter. """ with pytest.raises(ValueError, match="No positions to create"): - _ = openmm_afe.ABFEComplexSetupUnit._get_mda_universe( + _ = _get_mda_universe( topology="foo", positions=None, trajectory=None ) @@ -437,14 +438,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, @@ -462,6 +463,7 @@ def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings, system=setup_results["alchem_system"], positions=setup_results["debug_positions"], selection_indices=setup_results["selection_indices"], + alchemical_indices=setup_results["alchemical_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=True, dry=True, @@ -483,7 +485,7 @@ def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings, expected_indices = [ i + self.num_protein_component_atoms - 1 for i in range(self.num_ligand_atoms) ] - assert expected_indices == setup_results["alchem_indices"] + assert expected_indices == setup_results["alchemical_indices"]["A"] # Check the non-alchemical system assert setup_results["standard_system"].getNumParticles() == self.expected_complex_particles @@ -504,11 +506,19 @@ def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings, assert pdb.n_atoms == self.num_all_not_water # Check energies - alchem_region = AlchemicalRegion(alchemical_atoms=setup_results["alchem_indices"]) + alchemical_regions = [] + for key in setup_results["alchemical_indices"]: + alchemical_regions.append( + AlchemicalRegion( + alchemical_atoms=setup_results["alchemical_indices"][key], + name=key, + ) + ) + self._test_energies( reference_system=setup_results["standard_system"], alchemical_system=setup_results["alchem_system"], - alchemical_regions=alchem_region, + alchemical_regions=alchemical_regions, positions=setup_results["debug_positions"], ) @@ -523,6 +533,7 @@ def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings, system=setup_results["alchem_system"], positions=setup_results["debug_positions"], selection_indices=setup_results["selection_indices"], + alchemical_indices=setup_results["alchemical_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, dry=True, @@ -545,7 +556,7 @@ def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings, # Check the alchemical indices expected_indices = [i for i in range(self.num_ligand_atoms)] - assert expected_indices == setup_results["alchem_indices"] + assert expected_indices == setup_results["alchemical_indices"]["A"] # Check the non-alchemical system assert ( @@ -568,12 +579,19 @@ def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings, assert pdb.n_atoms == self.num_ligand_atoms # Check energies - alchem_region = AlchemicalRegion(alchemical_atoms=setup_results["alchem_indices"]) + alchemical_regions = [] + for key in setup_results["alchemical_indices"]: + alchemical_regions.append( + AlchemicalRegion( + alchemical_atoms=setup_results["alchemical_indices"][key], + name=key, + ) + ) self._test_energies( reference_system=setup_results["standard_system"], alchemical_system=setup_results["alchem_system"], - alchemical_regions=alchem_region, + alchemical_regions=alchemical_regions, positions=setup_results["debug_positions"], ) diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py index 2c63ea819..9586cc8c7 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py @@ -29,6 +29,7 @@ def patcher(): "positions": Path("positions.npy"), "pdb_structure": Path("hybrid_system.pdb"), "selection_indices": np.zeros(100), + "alchemical_indices": {"A": list(range(10))}, "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm, "standard_state_correction": 0 * offunit.kilocalorie_per_mole, "restraint_geometry": None, @@ -44,6 +45,7 @@ def patcher(): "positions": Path("positions.npy"), "pdb_structure": Path("hybrid_system.pdb"), "selection_indices": np.zeros(100), + "alchemical_indices": {"A": list(range(10))}, "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm, "standard_state_correction": 0 * offunit.kilocalorie_per_mole, "restraint_geometry": True, diff --git a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py index a5193a530..44cec8f4e 100644 --- a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py @@ -230,6 +230,7 @@ def test_setup_dry_sim_vac_benzene(benzene_system, method, protocol_dry_settings selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], dry=True, ) @@ -398,6 +399,7 @@ def test_setup_solv_benzene(benzene_system, protocol_dry_settings, tmp_path): selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path, @@ -451,6 +453,7 @@ def test_dry_run_vsite_fail(benzene_system, tmp_path, protocol_dry_settings): selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path, @@ -498,6 +501,7 @@ def test_setup_dry_sim_solv_benzene_tip4p(benzene_system, protocol_dry_settings, selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path, @@ -608,7 +612,7 @@ def assign_fictitious_charges(offmol): nonbond = [f for f in system.getForces() if isinstance(f, CustomNonbondedForce)] assert len(nonbond) == 4 - custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics"][0] + custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][0] # loop through the 12 benzene atoms for i in range(12): @@ -675,7 +679,7 @@ def test_dry_run_charge_backends( nonbond = [f for f in system.getForces() if isinstance(f, CustomNonbondedForce)] assert len(nonbond) == 4 - custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics"][0] + custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][0] charges = [] for i in range(system.getNumParticles()): @@ -756,6 +760,7 @@ def test_dry_run_vacuum_write_frequency( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path, diff --git a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py index c8ca84020..87cd90aa1 100644 --- a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py +++ b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py @@ -52,6 +52,7 @@ def patcher(): "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm, "standard_state_correction": 0 * offunit.kilocalorie_per_mole, "restraint_geometry": None, + "alchemical_indices": [1, 2, 3, 4, 5, 6], "gufe_version": gufe.__version__, "openfe_version": openfe.__version__, "openmm_version": openmm.__version__, @@ -64,6 +65,7 @@ def patcher(): "positions": Path("positions.npy"), "pdb_structure": Path("hybrid_system.pdb"), "selection_indices": np.zeros(100), + "alchemical_indices": [1, 2, 3, 4, 5, 6], "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm, "standard_state_correction": 0 * offunit.kilocalorie_per_mole, "restraint_geometry": None, diff --git a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py index e62767f3d..2ebd8cbbe 100644 --- a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py +++ b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py @@ -228,6 +228,7 @@ def test_resume( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], scratch_basepath=tmp_path, shared_basepath=tmp_path, ) @@ -303,6 +304,7 @@ def test_resume_fail_particles( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], scratch_basepath=tmp_path, shared_basepath=tmp_path, ) @@ -345,6 +347,7 @@ def test_resume_fail_constraints( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], scratch_basepath=tmp_path, shared_basepath=tmp_path, ) @@ -385,6 +388,7 @@ def test_resume_fail_forces( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], scratch_basepath=tmp_path, shared_basepath=tmp_path, ) @@ -441,6 +445,7 @@ def test_resume_differ_barostat( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], scratch_basepath=tmp_path, shared_basepath=tmp_path, ) @@ -486,7 +491,7 @@ def test_resume_differ_forces( # Now add a fake force new_force = openmm.NonbondedForce() new_force.setNonbondedMethod(openmm.NonbondedForce.PME) - new_force.addGlobalParameter("lambda_electrostatics", 1.0) + new_force.addGlobalParameter("lambda_electrostatics_A", 1.0) fake_system.addForce(new_force) @@ -500,6 +505,7 @@ def test_resume_differ_forces( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], scratch_basepath=tmp_path, shared_basepath=tmp_path, dry=True, @@ -546,6 +552,7 @@ def test_resume_bad_files( selection_indices=setup_results["selection_indices"], box_vectors=setup_results["box_vectors"], alchemical_restraints=False, + alchemical_indices=setup_results["alchemical_indices"], scratch_basepath=tmp_path, shared_basepath=tmp_path, ) From be5ed1e06b823f10e14cb8aab032ebbd6bb2a8b3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 17 Jun 2026 22:17:42 +0000 Subject: [PATCH 20/31] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/openfe/protocols/openmm_afe/abfe_units.py | 25 +++++++++++-------- .../protocols/openmm_afe/base_afe_units.py | 8 +++--- .../openmm_afe/equil_afe_settings.py | 6 +++-- .../openmm_afe/equil_binding_afe_method.py | 2 +- .../protocols/openmm_septop/base_units.py | 2 +- src/openfe/protocols/openmm_utils/states.py | 7 ++++-- .../openmm_abfe/test_abfe_energies.py | 8 +++--- .../openmm_abfe/test_abfe_protocol.py | 4 +-- .../openmm_ahfe/test_ahfe_protocol.py | 8 ++++-- 9 files changed, 40 insertions(+), 30 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 048ef4156..66021d821 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -31,8 +31,8 @@ ) from openfe.protocols.openmm_utils import system_validation from openfe.protocols.restraint_utils import geometry -from openfe.protocols.restraint_utils.geometry.utils import FindHostAtoms, get_central_atom_idx from openfe.protocols.restraint_utils.geometry.boresch import BoreschRestraintGeometry +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 @@ -147,10 +147,10 @@ def _get_minimum_image_distance(box_dimensions: npt.NDArray) -> Quantity: from MDAnalysis.lib import mdamath box_vectors = mdamath.triclinic_vectors(box_dimensions) - + # Calculate the volume based on the scalar triple product volume = mdamath.stp(box_vectors[0], box_vectors[1], box_vectors[2]) - + # Now calculate the perpendicular widths using perp_width_i = Volume / Area_of_face_i # Where Area_of_face_i is |box_vectors_{i+1} × box_vectors_{i+2}|. areas = np.cross(box_vectors[[1, 2, 0]], box_vectors[[2, 0, 1]]) @@ -218,6 +218,7 @@ class ABFESetupUnitMixin: """ Mixin for common class methods between Units """ + def _get_alchemical_ions( self, alchemical_components: dict[str, list[Component]], @@ -262,7 +263,7 @@ def _get_alchemical_ions( # 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: + if total_charge == 0 or not settings["alchemical_settings"].explicit_charge_correction: return None # TODO: For now, let's stick with a single -1/+1 case, but we should expand to more @@ -298,13 +299,11 @@ def _get_alchemical_ions( if box is None or np.all(np.isinfinite(box)) or np.any(box[:3] <= 0.0): # If it's not a dry simulation then error out if not dry: - errmsg = ( - f"Invalid box for co-alchemical ion search: {box}" - ) + errmsg = f"Invalid box for co-alchemical ion search: {box}" raise ValueError(errmsg) # For a dry execution, just assign a super high value - max_search_distance= 999 * offunit.nanometer + max_search_distance = 999 * offunit.nanometer else: # Set the max search distance to half the smallest perpendicular width # with a 1 Angstrom padding @@ -315,7 +314,7 @@ def _get_alchemical_ions( atom_finder = FindHostAtoms( host_atoms=ions_atomgroup, guest_atoms=alchem_atomgroup, - min_search_distance=settings['alchemical_settings'].alchemical_ion_min_distance, + min_search_distance=settings["alchemical_settings"].alchemical_ion_min_distance, max_search_distance=max_search_distance, ) @@ -409,7 +408,9 @@ def _get_settings(self) -> dict[str, SettingsBaseModel]: return settings -class ABFEComplexSetupUnit(ABFESetupUnitMixin, ComplexComponentsMixin, ComplexSettingsMixin, BaseAbsoluteSetupUnit): +class ABFEComplexSetupUnit( + ABFESetupUnitMixin, ComplexComponentsMixin, ComplexSettingsMixin, BaseAbsoluteSetupUnit +): """ Setup unit for the complex phase of absolute binding free energy transformations. @@ -710,7 +711,9 @@ def _get_settings(self) -> dict[str, SettingsBaseModel]: return settings -class ABFESolventSetupUnit(ABFESetupUnitMixin, SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteSetupUnit): +class ABFESolventSetupUnit( + ABFESetupUnitMixin, SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteSetupUnit +): """ Setup unit for the solvent phase of absolute binding free energy transformations. diff --git a/src/openfe/protocols/openmm_afe/base_afe_units.py b/src/openfe/protocols/openmm_afe/base_afe_units.py index bb3a3513e..a221afef4 100644 --- a/src/openfe/protocols/openmm_afe/base_afe_units.py +++ b/src/openfe/protocols/openmm_afe/base_afe_units.py @@ -85,8 +85,8 @@ serialize, ) from openfe.protocols.openmm_utils.states import ( - SingleRegionAlchemicalState, DualRegionAlchemicalState, + SingleRegionAlchemicalState, ) from openfe.protocols.restraint_utils import geometry from openfe.protocols.restraint_utils.openmm import omm_restraints @@ -564,7 +564,7 @@ def _add_restraints( alchem_comps: dict[str, list[Component]], comp_resids: dict[Component, npt.NDArray], settings: dict[str, SettingsBaseModel], - alchemical_ions : list[int] | None, + alchemical_ions: list[int] | None, ) -> tuple[ Quantity | None, openmm.System, @@ -618,11 +618,11 @@ def _get_alchemical_system( """ # first region is the alchemically changing species alchemical_indices = { - 'A': self._get_alchemical_indices(topology, comp_resids, alchemical_components) + "A": self._get_alchemical_indices(topology, comp_resids, alchemical_components) } # second region is any alchemical ions if alchemical_ions is not None: - alchemical_indices['B'] = alchemical_ions + alchemical_indices["B"] = alchemical_ions alchemical_regions = [] diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py index 687cbdc6b..dcdbdc002 100644 --- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py +++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py @@ -17,7 +17,6 @@ """ import numpy as np -from openff.units import unit as offunit from gufe.settings import ( OpenMMSystemGeneratorFFSettings, SettingsBaseModel, @@ -26,6 +25,7 @@ from gufe.settings.typing import ( NanometerQuantity, ) +from openff.units import unit as offunit from pydantic import field_validator from openfe.protocols.openmm_utils.omm_settings import ( @@ -104,7 +104,9 @@ class ABFEAlchemicalSettings(AlchemicalSettings): """ 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 + 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. diff --git a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py index 6f2d6db72..da0d5ccb1 100644 --- a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -43,10 +43,10 @@ from openfe.due import Doi, due from openfe.protocols.openmm_afe.equil_afe_settings import ( + ABFEAlchemicalSettings, ABFEBoreschRestraintSettings, ABFEPreEquilOutputSettings, AbsoluteBindingSettings, - ABFEAlchemicalSettings, IntegratorSettings, LambdaSettings, MDSimulationSettings, diff --git a/src/openfe/protocols/openmm_septop/base_units.py b/src/openfe/protocols/openmm_septop/base_units.py index f64b84ca4..173499ba0 100644 --- a/src/openfe/protocols/openmm_septop/base_units.py +++ b/src/openfe/protocols/openmm_septop/base_units.py @@ -76,6 +76,7 @@ from openfe.protocols.openmm_utils import omm_compute from openfe.protocols.openmm_utils.omm_settings import MultiStateAnalysisSettings, SettingsBaseModel from openfe.protocols.openmm_utils.serialization import deserialize +from openfe.protocols.openmm_utils.states import DualRegionAlchemicalState from openfe.utils import log_system_probe, without_oechem_backend from ..openmm_utils import ( @@ -86,7 +87,6 @@ system_validation, ) from ..openmm_utils.mdtraj_utils import mdtraj_from_openmm -from openfe.protocols.openmm_utils.states import DualRegionAlchemicalState logger = logging.getLogger(__name__) diff --git a/src/openfe/protocols/openmm_utils/states.py b/src/openfe/protocols/openmm_utils/states.py index 2edf6edf4..c79bfacd9 100644 --- a/src/openfe/protocols/openmm_utils/states.py +++ b/src/openfe/protocols/openmm_utils/states.py @@ -1,6 +1,6 @@ # This code is part of OpenFE and is licensed under the MIT license. # For details, see https://github.com/OpenFreeEnergy/openfe -# +# # Acknowledgements: # This module derives heavily from the AlchemicalState class # in openmmtools.alchemy (https://github.com/choderalab/openmmtools). @@ -12,6 +12,7 @@ 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) @@ -20,7 +21,7 @@ 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)) + raise ValueError("{} must be between 0 and 1.".format(self.parameter_name)) return float(parameter_value) @@ -58,6 +59,7 @@ class SingleRegionAlchemicalState(GlobalParameterState): :class:`openmmtools.states.GlobalParameterState` :class:`openfe.protocols.openmm_utils.states.DualRegionAlchemicalState` """ + lambda_sterics_A = _LambdaParameter("lambda_sterics_A") lambda_electrostatics_A = _LambdaParameter("lambda_electrostatics_A") lambda_bonds_A = _LambdaParameter("lambda_bonds_A") @@ -117,6 +119,7 @@ class DualRegionAlchemicalState(SingleRegionAlchemicalState): :class:`openmmtools.states.GlobalParameterState` :class:`openfe.protocols.openmm_utils.states.SingleRegionAlchemicalState` """ + lambda_sterics_B = _LambdaParameter("lambda_sterics_B") lambda_electrostatics_B = _LambdaParameter("lambda_electrostatics_B") lambda_bonds_B = _LambdaParameter("lambda_bonds_B") diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py index a6ac366f8..22e487d09 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py @@ -19,9 +19,9 @@ from openfe.protocols.openmm_afe.abfe_units import ( ABFEComplexSetupUnit, ) -from openfe.protocols.openmm_utils.states import SingleRegionAlchemicalState from openfe.protocols.openmm_utils.omm_settings import OpenMMSolvationSettings from openfe.protocols.openmm_utils.serialization import deserialize +from openfe.protocols.openmm_utils.states import SingleRegionAlchemicalState def get_alchemical_energy_components(alchemical_system, alchemical_state, positions, platform): @@ -248,9 +248,9 @@ def assert_energies(actual, expected, nonbonded_lower: bool): # turn off sterics expected = copy.deepcopy(energies) - expected["alchemically modified NonbondedForce for non-alchemical/alchemical sterics for region A"] = ( - 0 * ommunit.kilojoule_per_mole - ) + expected[ + "alchemically modified NonbondedForce for non-alchemical/alchemical sterics for region A" + ] = 0 * ommunit.kilojoule_per_mole expected[ "alchemically modified BondForce for non-alchemical/alchemical sterics exceptions for region A" ] = 0 * ommunit.kilojoule_per_mole diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index 37db3f2a4..c91abf856 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -170,9 +170,7 @@ def test_mda_universe_error(): when calling the mda Universe getter. """ with pytest.raises(ValueError, match="No positions to create"): - _ = _get_mda_universe( - topology="foo", positions=None, trajectory=None - ) + _ = _get_mda_universe(topology="foo", positions=None, trajectory=None) class TestT4LysozymeDryRun: diff --git a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py index 44cec8f4e..2a047ead0 100644 --- a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py @@ -612,7 +612,9 @@ def assign_fictitious_charges(offmol): nonbond = [f for f in system.getForces() if isinstance(f, CustomNonbondedForce)] assert len(nonbond) == 4 - custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][0] + custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][ + 0 + ] # loop through the 12 benzene atoms for i in range(12): @@ -679,7 +681,9 @@ def test_dry_run_charge_backends( nonbond = [f for f in system.getForces() if isinstance(f, CustomNonbondedForce)] assert len(nonbond) == 4 - custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][0] + custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][ + 0 + ] charges = [] for i in range(system.getNumParticles()): From 6930b72382e678fc326050debeaf77638daf8233 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 23:21:13 +0100 Subject: [PATCH 21/31] add missing imports --- src/openfe/protocols/openmm_afe/abfe_units.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 048ef4156..95929ae24 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -9,6 +9,7 @@ import logging import pathlib from collections.abc import Iterable +from copy import deepcopy import MDAnalysis as mda import numpy as np @@ -17,9 +18,10 @@ SolventComponent, ) from gufe.components import Component, SolvatedPDBComponent +from openff.units import unit and 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.states import ThermodynamicState @@ -35,6 +37,7 @@ from openfe.protocols.restraint_utils.geometry.boresch import BoreschRestraintGeometry 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 .base_afe_units import ( BaseAbsoluteMultiStateAnalysisUnit, @@ -201,7 +204,7 @@ def _find_most_common_ions( continue charge, _, _ = nbf.getParticleParameters(atoms[0].index) - charge_val = charge.value_in_unit(omm_unit.elementary_charge) + charge_val = charge.value_in_unit(ommunit.elementary_charge) if np.isclose(charge_val, target_charge, atol=0.01): ion_counts[residue.name] += 1 From 98b8101b91bc3f405169ac39ef649bdae8a8556c Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 23:22:27 +0100 Subject: [PATCH 22/31] fix typo --- src/openfe/protocols/openmm_afe/abfe_units.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 6ff0185a8..e4a9aaf1a 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -18,7 +18,7 @@ SolventComponent, ) from gufe.components import Component, SolvatedPDBComponent -from openff.units import unit and offunit +from openff.units import unit as offunit from openff.units import Quantity from openff.units.openmm import to_openmm from openmm import System, HarmonicBondForce From c059b6cee1ba6d03a55780a7443c112abe3e2fbd Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 23:23:45 +0100 Subject: [PATCH 23/31] more missing imports --- src/openfe/protocols/openmm_afe/abfe_units.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index e4a9aaf1a..a54c77dfa 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -12,6 +12,7 @@ from copy import deepcopy import MDAnalysis as mda +from MDAnalysis.lib.distances import calc_bonds import numpy as np import numpy.typing as npt from gufe import ( @@ -21,7 +22,11 @@ from openff.units import unit as offunit from openff.units import Quantity from openff.units.openmm import to_openmm -from openmm import System, HarmonicBondForce +from openmm import ( + System, + HarmonicBondForce, + NonbondedForce, +) from openmm import unit as ommunit from openmm.app import Topology as omm_topology from openmmtools.states import ThermodynamicState From fdb5679be374c7a2cf6ee91c1165c250897bbfe4 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 17 Jun 2026 22:24:33 +0000 Subject: [PATCH 24/31] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/openfe/protocols/openmm_afe/abfe_units.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index a54c77dfa..cc5143bf0 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -12,20 +12,20 @@ from copy import deepcopy import MDAnalysis as mda -from MDAnalysis.lib.distances import calc_bonds import numpy as np import numpy.typing as npt from gufe import ( SolventComponent, ) from gufe.components import Component, SolvatedPDBComponent -from openff.units import unit as offunit +from MDAnalysis.lib.distances import calc_bonds from openff.units import Quantity +from openff.units import unit as offunit from openff.units.openmm import to_openmm from openmm import ( - System, HarmonicBondForce, NonbondedForce, + System, ) from openmm import unit as ommunit from openmm.app import Topology as omm_topology @@ -41,8 +41,8 @@ from openfe.protocols.restraint_utils.geometry.boresch import BoreschRestraintGeometry 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 openfe.protocols.restraint_utils.openmm.omm_restraints import BoreschRestraint from .base_afe_units import ( BaseAbsoluteMultiStateAnalysisUnit, From 0e7bc94c8b71cf0dc4ce5a9d925508a043dc45aa Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 23:33:56 +0100 Subject: [PATCH 25/31] fix various mypy issues --- src/openfe/protocols/openmm_afe/abfe_units.py | 8 ++++---- src/openfe/protocols/openmm_afe/base_afe_units.py | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index a54c77dfa..8473d969b 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -199,7 +199,7 @@ def _find_most_common_ions( nbf = [i for i in openmm_system.getForces() if isinstance(i, NonbondedForce)][0] ion_counts: Counter = Counter() - ion_atom_indices: dict[str, int] = defaultdict(list) + ion_atom_indices: dict[str, list[int]] = defaultdict(list) # type: ignore[arg-type] for residue in openmm_topology.residues(): atoms = list(residue.atoms()) @@ -289,7 +289,7 @@ def _get_alchemical_ions( univ = _get_mda_universe( openmm_topology, positions, - self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, + self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename, # type: ignore[attr-defined] ) # get an atomgroup of the possible alchemical ions @@ -304,7 +304,7 @@ def _get_alchemical_ions( univ.trajectory[-1] # use the box dimensions from the last frame box = univ.dimensions - if box is None or np.all(np.isinfinite(box)) or np.any(box[:3] <= 0.0): + if box is None or np.all(np.isinf(box)) or np.any(box[:3] <= 0.0): # If it's not a dry simulation then error out if not dry: errmsg = f"Invalid box for co-alchemical ion search: {box}" @@ -315,7 +315,7 @@ def _get_alchemical_ions( else: # Set the max search distance to half the smallest perpendicular width # with a 1 Angstrom padding - max_search_distance = (_get_minimum_image_distance(box) * 0.5) - 1 * offunit.angstrom + max_search_distance = (_get_minimum_image_distance(box) * 0.5) - 1 * offunit.angstrom # type: ignore[operator] # Re-using a utility from the restraints utilities # TODO: rename this class! diff --git a/src/openfe/protocols/openmm_afe/base_afe_units.py b/src/openfe/protocols/openmm_afe/base_afe_units.py index a221afef4..5a5379d47 100644 --- a/src/openfe/protocols/openmm_afe/base_afe_units.py +++ b/src/openfe/protocols/openmm_afe/base_afe_units.py @@ -583,7 +583,7 @@ def _get_alchemical_system( alchemical_components: dict[str, list[Component]], alchemical_ions: list[int] | None, alchemical_settings: AlchemicalSettings, - ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, list[int]]: + ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, dict[str, list[int]]]: """ Get an alchemically modified system and its associated factory @@ -1032,7 +1032,7 @@ def _add_lambdas_to_protocol(protocol, lambdas, region_name): for key in ["lambda_electrostatics", "lambda_sterics"]: protocol[f"{key}_{region_name}"] = copy.deepcopy(lambdas[key]) - param_protocol = {} + param_protocol: dict[str, list[float]] = {} # Always add for the first region _add_lambdas_to_protocol(param_protocol, lambdas, "A") From b2e3ef06f96ca1ef63dce68749a63c9f0ce61940 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 23:36:59 +0100 Subject: [PATCH 26/31] remove unused ignore --- src/openfe/protocols/openmm_afe/abfe_units.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 869438c54..0d067c97d 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -199,7 +199,7 @@ def _find_most_common_ions( nbf = [i for i in openmm_system.getForces() if isinstance(i, NonbondedForce)][0] ion_counts: Counter = Counter() - ion_atom_indices: dict[str, list[int]] = defaultdict(list) # type: ignore[arg-type] + ion_atom_indices: dict[str, list[int]] = defaultdict(list) for residue in openmm_topology.residues(): atoms = list(residue.atoms()) From 977b3859b206fe1a5890b27f86bc38525a3d49ac Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 17 Jun 2026 23:56:28 +0100 Subject: [PATCH 27/31] Slightly improve user charge tests --- .../tests/protocols/openmm_abfe/test_abfe_protocol.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index c91abf856..d7ed153f7 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -744,6 +744,12 @@ def assign_fictitious_charges(offmol): c, s, e = system_nbf.getParticleParameters(index) assert pytest.approx(prop_chgs[i]) == c.value_in_unit(ommunit.elementary_charge) + + # Alchemical system should be 0 charge in 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 88a014cc73af786089445d31f5e82d299942dd4e Mon Sep 17 00:00:00 2001 From: IAlibay Date: Thu, 18 Jun 2026 00:27:42 +0100 Subject: [PATCH 28/31] Modify validation checks to allow for net charge ligands --- .../openmm_afe/equil_binding_afe_method.py | 39 +++++++++++++---- .../data/openmm_rfe/charged_benzenes.sdf | 38 ++++++++++++++++ .../openmm_abfe/test_abfe_validation.py | 43 ++++++++++++++----- 3 files changed, 100 insertions(+), 20 deletions(-) diff --git a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py index da0d5ccb1..e68dd2699 100644 --- a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py +++ b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py @@ -254,6 +254,7 @@ def _adaptive_settings( def _validate_endstates( stateA: ChemicalSystem, stateB: ChemicalSystem, + charge_correction: bool, ) -> None: """ A binding transformation is defined (in terms of gufe components) @@ -266,6 +267,8 @@ def _validate_endstates( The chemical system of end state A stateB : ChemicalSystem The chemical system of end state B + charge_correction : bool + If we are using a charge correction scheme. Raises ------ @@ -275,7 +278,10 @@ def _validate_endstates( If stateA has more than one unique Component. If the stateA unique Component is not a SmallMoleculeComponent. If stateB contains any unique Components. - If the alchemical species is charged. + If the alchemical species is has greater than one unit charge. + UserWarning + If the alchemical species has a net charge but ``charge_correction`` + is ``False``. """ if not (stateA.contains(ProteinComponent) and stateB.contains(ProteinComponent)): # Check if there is a suitable SmallMoleculeComponent that could @@ -306,14 +312,25 @@ def _validate_endstates( ) raise ValueError(errmsg) - # Check that the state A unique isn't charged + # Check that the state A unique total charge if diff[0][0].total_charge != 0: - errmsg = ( - "Charged alchemical molecules are not currently " - "supported for solvation free energies. " - f"Molecule total charge: {diff[0][0].total_charge}." - ) - raise ValueError(errmsg) + # Error if the total charge is greater than 1 + if abs(diff[0][0].total_charge) > 1: + errmsg = ( + "Alchemical molecules with a formal charge of greater than 1 " + "are not currently supported in binding free energy calculations. " + f"Molecule total charge: {diff[0][0].total_charge}." + ) + raise ValueError(errmsg) + + # Warn if we aren't using a charge correction + if not charge_correction: + wmsg = ( + "Ligand has a net charge but no charge correction scheme has been requested. " + "Please note that you will need to apply your own correction scheme " + "to account for finite-size effects." + ) + warnings.warn(wmsg) # If there are any alchemical Components in state B if len(diff[1]) > 0: @@ -405,7 +422,11 @@ def _validate( # Validate the end states & alchemical components system_validation.validate_chemical_system(stateA) system_validation.validate_chemical_system(stateB) - self._validate_endstates(stateA, stateB) + self._validate_endstates( + stateA, + stateB, + self.settings.alchemical_settings.explicit_charge_correction, + ) # Validate the complex lambda schedule self._validate_lambda_schedule( diff --git a/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf b/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf index ea2735835..18a65a49c 100644 --- a/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf +++ b/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf @@ -98,3 +98,41 @@ 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/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py index e8f7108c3..eb2a4d79d 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py @@ -127,7 +127,7 @@ def test_validate_no_protcomp( errmsg = "No suitable host molecule found" with pytest.raises(ValueError, match=errmsg): - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) def test_validate_smc_host( @@ -151,7 +151,7 @@ def test_validate_smc_host( } ) - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) def test_validate_endstates_nosolvcomp_stateA(benzene_modifications, T4_protein_component): @@ -171,7 +171,7 @@ def test_validate_endstates_nosolvcomp_stateA(benzene_modifications, T4_protein_ ) with pytest.raises(ValueError, match="No SolventComponent found"): - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) def test_validate_endstates_nosolvcomp_stateB(benzene_modifications, T4_protein_component): @@ -191,7 +191,7 @@ def test_validate_endstates_nosolvcomp_stateB(benzene_modifications, T4_protein_ ) with pytest.raises(ValueError, match="No SolventComponent found"): - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) def test_validate_endstates_multiple_uniqueA(benzene_modifications, T4_protein_component): @@ -212,7 +212,7 @@ def test_validate_endstates_multiple_uniqueA(benzene_modifications, T4_protein_c ) with pytest.raises(ValueError, match="Only one alchemical species"): - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) def test_validate_solvent_endstates_solvent_dissapearing( @@ -237,7 +237,7 @@ def test_validate_solvent_endstates_solvent_dissapearing( errmsg = "Only disappearing small molecule components" with pytest.raises(ValueError, match=errmsg): - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) def test_validate_endstates_unique_stateB(benzene_modifications, T4_protein_component): @@ -258,13 +258,13 @@ def test_validate_endstates_unique_stateB(benzene_modifications, T4_protein_comp ) with pytest.raises(ValueError, match="appearing in state B"): - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) -def test_charged_endstate(charged_benzene_modifications, T4_protein_component): +def test_highly_charged_endstate_error(charged_benzene_modifications, T4_protein_component): stateA = ChemicalSystem( { - "benzene": charged_benzene_modifications["benzoic_acid"], + "benzene": charged_benzene_modifications["phthalic_acid"], "protein": T4_protein_component, "solvent": SolventComponent(), } @@ -277,9 +277,30 @@ def test_charged_endstate(charged_benzene_modifications, T4_protein_component): } ) - errmsg = "Charged alchemical molecules are not currently supported" + errmsg = "Alchemical molecules with a formal charge of greater than 1" with pytest.raises(ValueError, match=errmsg): - AbsoluteBindingProtocol._validate_endstates(stateA, stateB) + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True) + + +def test_charge_endstate_warning(charged_benzene_modifications, T4_protein_component): + stateA = ChemicalSystem( + { + "benzene": charged_benzene_modifications["benzoic_acid"], + "protein": T4_protein_component, + "solvent": SolventComponent(), + } + ) + + stateB = ChemicalSystem( + { + "protein": T4_protein_component, + "solvent": SolventComponent(), + } + ) + + msg = "Ligand has a net charge but no charge correction scheme has been requested." + with pytest.warns(UserWarning, match=msg): + AbsoluteBindingProtocol._validate_endstates(stateA, stateB, False) def test_validate_fail_extends(benzene_modifications, T4_protein_component, default_settings): From b1a788d34e0aad4f6fd2d5fab983497243c0c8c8 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Thu, 18 Jun 2026 02:33:57 +0100 Subject: [PATCH 29/31] small fix --- src/openfe/protocols/openmm_afe/abfe_units.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 0d067c97d..2ae6c36b8 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -818,7 +818,7 @@ def _add_restraints( distance = float( calc_bonds( - alchem_ion_ag.position, + alchem_ion_ag.atoms[0].position, ligand_central_atom_ag.position, box=universe.dimensions, ) From 8d28d963c6cff8783c123f17a1b2884bda6030c4 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Thu, 18 Jun 2026 14:40:10 +0100 Subject: [PATCH 30/31] Some fixes for SepTop tests --- src/openfe/tests/protocols/openmm_septop/test_septop_slow.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py b/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py index c003432a0..022b57b6d 100644 --- a/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py +++ b/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py @@ -18,7 +18,7 @@ SepTopProtocol, SepTopSolventSetupUnit, ) -from openfe.protocols.openmm_septop.utils import SepTopParameterState +from openfe.protocols.openmm_utils.states import DualRegionAlchemicalState from openfe.protocols.openmm_utils.serialization import deserialize @@ -30,7 +30,7 @@ def default_settings(): def compare_energies(alchemical_system, positions): - alchemical_state = SepTopParameterState.from_system(alchemical_system) + alchemical_state = DualRegionAlchemicalState.from_system(alchemical_system) from openmmtools.alchemy import AbsoluteAlchemicalFactory From 1de4be685c9de89c846117c503d4568ff1aba216 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Mon, 22 Jun 2026 08:55:07 +0100 Subject: [PATCH 31/31] Temporarily disable harmonic bond addition in solvent --- src/openfe/protocols/openmm_afe/abfe_units.py | 7 +- src/openfe/tests/conftest.py | 19 + .../host_guest/OA_guests_nagl_charged.sdf | 414 +++++++++++++++++ .../data/host_guest/OA_host_nagl_charges.sdf | 422 ++++++++++++++++++ src/openfe/tests/data/host_guest/__init__.py | 0 .../openmm_abfe/test_abfe_protocol.py | 254 ++++++++++- 6 files changed, 1110 insertions(+), 6 deletions(-) create mode 100644 src/openfe/tests/data/host_guest/OA_guests_nagl_charged.sdf create mode 100644 src/openfe/tests/data/host_guest/OA_host_nagl_charges.sdf create mode 100644 src/openfe/tests/data/host_guest/__init__.py diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py index 2ae6c36b8..d3a475e2e 100644 --- a/src/openfe/protocols/openmm_afe/abfe_units.py +++ b/src/openfe/protocols/openmm_afe/abfe_units.py @@ -310,8 +310,8 @@ def _get_alchemical_ions( errmsg = f"Invalid box for co-alchemical ion search: {box}" raise ValueError(errmsg) - # For a dry execution, just assign a super high value - max_search_distance = 999 * offunit.nanometer + # For a dry execution, just assign a reasonably high value + max_search_distance = 1.5 * offunit.nanometer else: # Set the max search distance to half the smallest perpendicular width # with a 1 Angstrom padding @@ -838,7 +838,8 @@ def _add_restraints( ) force.setName("ion_restraint") - add_force_in_separate_group(restrained_system, force) + # TODO: Temporarily disabled this for testing + # add_force_in_separate_group(restrained_system, force) return None, restrained_system, None diff --git a/src/openfe/tests/conftest.py b/src/openfe/tests/conftest.py index 3e9046b12..10e7047b0 100644 --- a/src/openfe/tests/conftest.py +++ b/src/openfe/tests/conftest.py @@ -255,6 +255,25 @@ def charged_benzene_modifications(): return files +@pytest.fixture(scope="session") +def OA_guests_charged() -> dict[str, SmallMoleculeComponent]: + files = {} + with resources.as_file(resources.files("openfe.tests.data.host_guest")) as d: + fn = str(d / "OA_guests_nagl_charged.sdf") + supp = Chem.SDMolSupplier(str(fn), removeHs=False) + for rdmol in supp: + files[rdmol.GetProp("_Name")] = SmallMoleculeComponent(rdmol) + return files + + +@pytest.fixture(scope="session") +def OA_host_charged() -> SmallMoleculeComponent: + with resources.as_file(resources.files("openfe.tests.data.host_guest")) as d: + fn = str(d / "OA_host_nagl_charges.sdf") + mol = [m for m in Chem.SDMolSupplier(str(fn), removeHs=False)][0] + yield SmallMoleculeComponent(mol) + + @pytest.fixture(scope="session") def T4L_septop_reference_xml(): with resources.as_file(resources.files("openfe.tests.data.openmm_septop")) as d: diff --git a/src/openfe/tests/data/host_guest/OA_guests_nagl_charged.sdf b/src/openfe/tests/data/host_guest/OA_guests_nagl_charged.sdf new file mode 100644 index 000000000..b1f356172 --- /dev/null +++ b/src/openfe/tests/data/host_guest/OA_guests_nagl_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/src/openfe/tests/data/host_guest/OA_host_nagl_charges.sdf b/src/openfe/tests/data/host_guest/OA_host_nagl_charges.sdf new file mode 100644 index 000000000..c49fe096f --- /dev/null +++ b/src/openfe/tests/data/host_guest/OA_host_nagl_charges.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/src/openfe/tests/data/host_guest/__init__.py b/src/openfe/tests/data/host_guest/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py index d7ed153f7..ceb421e1c 100644 --- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py +++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py @@ -176,7 +176,7 @@ def test_mda_universe_error(): class TestT4LysozymeDryRun: solvent = SolventComponent(ion_concentration=0 * offunit.molar) num_all_not_water = 2634 # 9 counterions to neutralize - num_protein_component_atoms = 2614 + num_host_component_atoms = 2614 # No ions num_ligand_atoms = 12 expected_complex_particles = 32607 @@ -481,7 +481,7 @@ def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings, # Check the alchemical indices expected_indices = [ - i + self.num_protein_component_atoms - 1 for i in range(self.num_ligand_atoms) + i + self.num_host_component_atoms - 1 for i in range(self.num_ligand_atoms) ] assert expected_indices == setup_results["alchemical_indices"]["A"] @@ -631,6 +631,254 @@ def settings(self): return s +class TestHostGuestCharged(TestT4LysozymeDryRun): + """ + A host guest system with user defined Boresch restraints + and a net charge. + """ + solvent = SolventComponent(ion_concentration=0.15 * offunit.molar) + num_all_not_water = 223 + num_host_component_atoms = 192 + # No ions + num_ligand_atoms = 75 + expected_complex_particles = 6166 + expected_ligand_solvent_particles = 29910 + + barostat_by_phase = { + "complex": MonteCarloBarostat, + "solvent": MonteCarloBarostat, + } + + @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.complex_solvation_settings.solvent_padding = None + s.complex_solvation_settings.number_of_solvent_molecules = 2000 + s.solvent_solvation_settings.box_shape = "cube" + s.solvent_solvation_settings.solvent_padding = None + s.solvent_solvation_settings.number_of_solvent_molecules = 10000 + s.restraint_settings.guest_restraint_ids = [1, 2, 4] + s.restraint_settings.host_restraint_ids = [1, 4, 3] + s.alchemical_settings.alchemical_ion_min_distance = 1.0 * offunit.nanometer + return s + + + @pytest.fixture(scope="class") + def dag(self, protocol, OA_guests_charged, OA_host_charged): + stateA = ChemicalSystem( + { + "ligand": OA_guests_charged["OA-G0"], + "host": OA_host_charged, + "solvent": self.solvent, + } + ) + + stateB = ChemicalSystem( + { + "host": OA_host_charged, + "solvent": self.solvent, + } + ) + + return protocol.create( + stateA=stateA, + stateB=stateB, + mapping=None, + ) + + def _assert_expected_alchemical_forces(self, system, phase: str, settings): + """ + Assert the forces expected in the alchemical system. + """ + barostat_type = self.barostat_by_phase[phase] + self._assert_force_num(system, NonbondedForce, 1) + self._assert_force_num(system, CustomNonbondedForce, 4) + self._assert_force_num(system, CustomBondForce, 4) + self._assert_force_num(system, HarmonicAngleForce, 1) + self._assert_force_num(system, PeriodicTorsionForce, 1) + self._assert_force_num(system, barostat_type, 1) + + if phase == "complex": + self._assert_force_num(system, HarmonicBondForce, 1) + self._assert_force_num(system, CustomCompoundBondForce, 1) + else: + self._assert_force_num(system, HarmonicBondForce, 2) + + 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, barostat_type)] + assert len(barostat) == 1 + expected_frequency = int( + ( + settings.complex_integrator_settings + if phase == "complex" + else settings.solvent_integrator_settings + ).barostat_frequency.m + ) + assert barostat[0].getFrequency() == expected_frequency + assert barostat[0].getDefaultPressure() == to_openmm(settings.thermo_settings.pressure) + assert barostat[0].getDefaultTemperature() == to_openmm( + settings.thermo_settings.temperature + ) + + def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings, tmp_path): + setup_results = complex_setup_units[0].run( + dry=True, + verbose=True, + scratch_basepath=tmp_path, + shared_basepath=tmp_path, + ) + sim_results = complex_sim_units[0].run( + system=setup_results["alchem_system"], + positions=setup_results["debug_positions"], + selection_indices=setup_results["selection_indices"], + alchemical_indices=setup_results["alchemical_indices"], + box_vectors=setup_results["box_vectors"], + alchemical_restraints=True, + dry=True, + scratch_basepath=tmp_path, + shared_basepath=tmp_path, + ) + + # Check the sampler + self._verify_sampler(sim_results["sampler"], phase="complex", settings=settings) + + # Check the alchemical system + assert setup_results["alchem_system"].getNumParticles() == self.expected_complex_particles + self._assert_expected_alchemical_forces( + setup_results["alchem_system"], phase="complex", settings=settings + ) + self._check_box_vectors(setup_results["alchem_system"]) + + # # Check the alchemical indices + # expected_indices = [ + # i + self.num_host_component_atoms - 1 for i in range(self.num_ligand_atoms) + # ] + # assert expected_indices == setup_results["alchemical_indices"]["A"] + + # Check the non-alchemical system + assert setup_results["standard_system"].getNumParticles() == self.expected_complex_particles + self._assert_expected_nonalchemical_forces( + setup_results["standard_system"], + "complex", + settings=settings, + ) + self._check_box_vectors(setup_results["standard_system"]) + # Check the box vectors haven't changed (they shouldn't have because we didn't do MD) + assert_allclose( + from_openmm(setup_results["alchem_system"].getDefaultPeriodicBoxVectors()), + from_openmm(setup_results["standard_system"].getDefaultPeriodicBoxVectors()), + ) + + # Check the PDB + pdb = mdt.load_pdb(setup_results["pdb_structure"]) + assert pdb.n_atoms == self.num_all_not_water + + # Check energies + alchemical_regions = [] + for key in setup_results["alchemical_indices"]: + alchemical_regions.append( + AlchemicalRegion( + alchemical_atoms=setup_results["alchemical_indices"][key], + name=key, + ) + ) + + self._test_energies( + reference_system=setup_results["standard_system"], + alchemical_system=setup_results["alchem_system"], + alchemical_regions=alchemical_regions, + positions=setup_results["debug_positions"], + ) + + def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings, tmp_path): + setup_results = solvent_setup_units[0].run( + dry=True, + verbose=True, + scratch_basepath=tmp_path, + shared_basepath=tmp_path, + ) + sim_results = solvent_sim_units[0].run( + system=setup_results["alchem_system"], + positions=setup_results["debug_positions"], + selection_indices=setup_results["selection_indices"], + alchemical_indices=setup_results["alchemical_indices"], + box_vectors=setup_results["box_vectors"], + alchemical_restraints=False, + dry=True, + scratch_basepath=tmp_path, + shared_basepath=tmp_path, + ) + + # Check the sampler + self._verify_sampler(sim_results["sampler"], phase="solvent", settings=settings) + + # Check the alchemical system + assert ( + setup_results["alchem_system"].getNumParticles() + == self.expected_ligand_solvent_particles + ) + self._assert_expected_alchemical_forces( + setup_results["alchem_system"], phase="solvent", settings=settings + ) + self._test_cubic_vectors(setup_results["alchem_system"]) + + # # Check the alchemical indices + # expected_indices = [i for i in range(self.num_ligand_atoms)] + # assert expected_indices == setup_results["alchemical_indices"]["A"] + + # Check the non-alchemical system + assert ( + setup_results["standard_system"].getNumParticles() + == self.expected_ligand_solvent_particles + ) + self._assert_expected_nonalchemical_forces( + setup_results["standard_system"], phase="solvent", settings=settings + ) + self._test_cubic_vectors(setup_results["standard_system"]) + + # Check the box vectors haven't changed (they shouldn't have because we didn't do MD) + assert_allclose( + from_openmm(setup_results["alchem_system"].getDefaultPeriodicBoxVectors()), + from_openmm(setup_results["standard_system"].getDefaultPeriodicBoxVectors()), + ) + + # Check the PDB + pdb = mdt.load_pdb(setup_results["pdb_structure"]) + assert pdb.n_atoms == self.num_ligand_atoms + + # Check energies + alchemical_regions = [] + for key in setup_results["alchemical_indices"]: + alchemical_regions.append( + AlchemicalRegion( + alchemical_atoms=setup_results["alchemical_indices"][key], + name=key, + ) + ) + + self._test_energies( + reference_system=setup_results["standard_system"], + alchemical_system=setup_results["alchem_system"], + alchemical_regions=alchemical_regions, + positions=setup_results["debug_positions"], + ) + + def test_user_restraint(benzene_modifications_am1bcc, T4_protein_component, tmp_path): s = openmm_afe.AbsoluteBindingProtocol.default_settings() s.protocol_repeats = 1 @@ -767,7 +1015,7 @@ class TestA2AMembraneDryRun(TestT4LysozymeDryRun): solvent = SolventComponent(ion_concentration=0 * offunit.molar) num_all_not_water = 16080 - num_protein_component_atoms = 39391 + num_host_component_atoms = 39391 expected_complex_particles = 39426 expected_ligand_solvent_particles = 3036