From 43d60073540a84f2338709932e0de9d4450c5fe0 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 5 Dec 2025 13:18:52 +0000 Subject: [PATCH 01/20] Migrate validation to Protocol._validate --- .../protocols/openmm_rfe/equil_rfe_methods.py | 434 +++++++++++------- 1 file changed, 260 insertions(+), 174 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 514237634..f272631b8 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -117,161 +117,6 @@ def _get_resname(off_mol) -> str: return names[0] -def _get_alchemical_charge_difference( - mapping: LigandAtomMapping, - nonbonded_method: str, - explicit_charge_correction: bool, - solvent_component: SolventComponent, -) -> int: - """ - Checks and returns the difference in formal charge between state A and B. - - Raises - ------ - ValueError - * If an explicit charge correction is attempted and the - nonbonded method is not PME. - * If the absolute charge difference is greater than one - and an explicit charge correction is attempted. - UserWarning - If there is any charge difference. - - Parameters - ---------- - mapping : dict[str, ComponentMapping] - Dictionary of mappings between transforming components. - nonbonded_method : str - The OpenMM nonbonded method used for the simulation. - explicit_charge_correction : bool - Whether or not to use an explicit charge correction. - solvent_component : openfe.SolventComponent - The SolventComponent of the simulation. - - Returns - ------- - int - The formal charge difference between states A and B. - This is defined as sum(charge state A) - sum(charge state B) - """ - - difference = mapping.get_alchemical_charge_difference() - - if abs(difference) > 0: - if explicit_charge_correction: - if nonbonded_method.lower() != "pme": - errmsg = "Explicit charge correction when not using PME is not currently supported." - raise ValueError(errmsg) - if abs(difference) > 1: - errmsg = ( - f"A charge difference of {difference} is observed " - "between the end states and an explicit charge " - "correction has been requested. Unfortunately " - "only absolute differences of 1 are supported." - ) - raise ValueError(errmsg) - - ion = {-1: solvent_component.positive_ion, 1: solvent_component.negative_ion}[ - difference - ] - wmsg = ( - f"A charge difference of {difference} is observed " - "between the end states. This will be addressed by " - f"transforming a water into a {ion} ion" - ) - logger.warning(wmsg) - warnings.warn(wmsg) - else: - wmsg = ( - f"A charge difference of {difference} is observed " - "between the end states. No charge correction has " - "been requested, please account for this in your " - "final results." - ) - logger.warning(wmsg) - warnings.warn(wmsg) - - return difference - - -def _validate_alchemical_components( - alchemical_components: dict[str, list[Component]], - mapping: Optional[Union[ComponentMapping, list[ComponentMapping]]], -): - """ - Checks that the alchemical components are suitable for the RFE protocol. - - Specifically we check: - 1. That all alchemical components are mapped. - 2. That all alchemical components are SmallMoleculeComponents. - 3. If the mappings involves element changes in core atoms - - Parameters - ---------- - alchemical_components : dict[str, list[Component]] - Dictionary contatining the alchemical components for - states A and B. - mapping : Optional[Union[ComponentMapping, list[ComponentMapping]]] - all mappings between transforming components. - - Raises - ------ - ValueError - * If there are more than one mapping or mapping is None - * If there are any unmapped alchemical components. - * If there are any alchemical components that are not - SmallMoleculeComponents. - UserWarning - * Mappings which involve element changes in core atoms - """ - if isinstance(mapping, ComponentMapping): - mapping = [mapping] - # Check mapping - # For now we only allow for a single mapping, this will likely change - if mapping is None or len(mapping) != 1: - errmsg = "A single LigandAtomMapping is expected for this Protocol" - raise ValueError(errmsg) - - # Check that all alchemical components are mapped & small molecules - mapped = { - "stateA": [m.componentA for m in mapping], - "stateB": [m.componentB for m in mapping], - } - - for idx in ["stateA", "stateB"]: - if len(alchemical_components[idx]) != len(mapped[idx]): - errmsg = f"missing alchemical components in {idx}" - raise ValueError(errmsg) - for comp in alchemical_components[idx]: - if comp not in mapped[idx]: - raise ValueError(f"Unmapped alchemical component {comp}") - if not isinstance(comp, SmallMoleculeComponent): # pragma: no-cover - errmsg = ( - "Transformations involving non " - "SmallMoleculeComponent species {comp} " - "are not currently supported" - ) - raise ValueError(errmsg) - - # Validate element changes in mappings - for m in mapping: - molA = m.componentA.to_rdkit() - molB = m.componentB.to_rdkit() - for i, j in m.componentA_to_componentB.items(): - atomA = molA.GetAtomWithIdx(i) - atomB = molB.GetAtomWithIdx(j) - if atomA.GetAtomicNum() != atomB.GetAtomicNum(): - wmsg = ( - f"Element change in mapping between atoms " - f"Ligand A: {i} (element {atomA.GetAtomicNum()}) and " - f"Ligand B: {j} (element {atomB.GetAtomicNum()})\n" - "No mass scaling is attempted in the hybrid topology, " - "the average mass of the two atoms will be used in the " - "simulation" - ) - logger.warning(wmsg) - warnings.warn(wmsg) # TODO: remove this once logging is fixed - - class RelativeHybridTopologyProtocolResult(gufe.ProtocolResult): """Dict-like container for the output of a RelativeHybridTopologyProtocol""" @@ -612,21 +457,204 @@ def _adaptive_settings( return protocol_settings - def _create( + @staticmethod + def _validate_endstates( + stateA: ChemicalSystem, + stateB: ChemicalSystem, + ) -> None: + """ + Validates the end states for the RFE protocol. + + Parameters + ---------- + stateA : ChemicalSystem + The chemical system of end state A. + stateB : ChemicalSystem + The chemical system of end state B. + + Raises + ------ + ValueError + * If either state contains more than one unique Component. + * If unique components are not SmallMoleculeComponents. + """ + # Get the difference in Components between each state + diff = stateA.component_diff(stateB) + + for i, entry in enumerate(diff): + state_label = "A" if i == 0 else "B" + + # Check that there is only one unique Component in each state + if len(entry) != 0: + errmsg = ( + "Only one alchemical component is allowed per end state. " + f"Found {len(entry)} in state {state_label}." + ) + raise ValueError(errmsg) + + # Check that the unique Component is a SmallMoleculeComponent + if not isinstance(entry[0], SmallMoleculeComponent): + errmsg = ( + f"Alchemical component in state {state_label} is of type " + f"{type(entry[0])}, but only SmallMoleculeComponents " + "transformations are currently supported." + ) + raise ValueError(errmsg) + + @staticmethod + def _validate_mapping( + mapping: Optional[Union[ComponentMapping, list[ComponentMapping]]], + alchemical_components: dict[str, list[Component]], + ) -> None: + """ + Validates that the provided mapping(s) are suitable for the RFE protocol. + + Parameters + ---------- + mapping : Optional[Union[ComponentMapping, list[ComponentMapping]]] + all mappings between transforming components. + alchemical_components : dict[str, list[Component]] + Dictionary contatining the alchemical components for + states A and B. + + Raises + ------ + ValueError + * If there are more than one mapping or mapping is None + * If the mapping components are not in the alchemical components. + UserWarning + * Mappings which involve element changes in core atoms + """ + # if a single mapping is provided, convert to list + if isinstance(mapping, ComponentMapping): + mapping = [mapping] + + # For now we only support a single mapping + if mapping is None or len(mapping) > 1: + errmsg = "A single LigandAtomMapping is expected for this Protocol" + raise ValueError(errmsg) + + # check that the mapping components are in the alchemical components + for m in mapping: + if m.componentA not in alchemical_components["stateA"]: + raise ValueError(f"Mapping componentA {m.componentA} not in alchemical components of stateA") + if m.componentB not in alchemical_components["stateB"]: + raise ValueError(f"Mapping componentB {m.componentB} not in alchemical components of stateB") + + # TODO: remove - this is now the default behaviour? + # Check for element changes in mappings + for m in mapping: + molA = m.componentA.to_rdkit() + molB = m.componentB.to_rdkit() + for i, j in m.componentA_to_componentB.items(): + atomA = molA.GetAtomWithIdx(i) + atomB = molB.GetAtomWithIdx(j) + if atomA.GetAtomicNum() != atomB.GetAtomicNum(): + wmsg = ( + f"Element change in mapping between atoms " + f"Ligand A: {i} (element {atomA.GetAtomicNum()}) and " + f"Ligand B: {j} (element {atomB.GetAtomicNum()})\n" + "No mass scaling is attempted in the hybrid topology, " + "the average mass of the two atoms will be used in the " + "simulation" + ) + logger.warning(wmsg) + warnings.warn(wmsg) + + @staticmethod + def _validate_charge_difference( + mapping: LigandAtomMapping, + nonbonded_method: str, + explicit_charge_correction: bool, + solvent_component: SolventComponent | None, + ): + """ + Validates the net charge difference between the two states. + + Parameters + ---------- + mapping : dict[str, ComponentMapping] + Dictionary of mappings between transforming components. + nonbonded_method : str + The OpenMM nonbonded method used for the simulation. + explicit_charge_correction : bool + Whether or not to use an explicit charge correction. + solvent_component : openfe.SolventComponent | None + The SolventComponent of the simulation. + + Raises + ------ + ValueError + * If an explicit charge correction is attempted and the + nonbonded method is not PME. + * If the absolute charge difference is greater than one + and an explicit charge correction is attempted. + UserWarning + * If there is any charge difference. + """ + difference = mapping.get_alchemical_charge_difference() + + if abs(difference) == 0: + return + + if not explicit_charge_correction: + wmsg = ( + f"A charge difference of {difference} is observed " + "between the end states. No charge correction has " + "been requested, please account for this in your " + "final results." + ) + logger.warning(wmsg) + warnings.warn(wmsg) + return + + # We implicitly check earlier that we have to have pme for a solvated + # system, so we only need to check the nonbonded method here + if nonbonded_method.lower() != "pme": + errmsg = "Explicit charge correction when not using PME is not currently supported." + raise ValueError(errmsg) + + if abs(difference) > 1: + errmsg = ( + f"A charge difference of {difference} is observed " + "between the end states and an explicit charge " + "correction has been requested. Unfortunately " + "only absolute differences of 1 are supported." + ) + raise ValueError(errmsg) + + ion = { + -1: solvent_component.positive_ion, + 1: solvent_component.negative_ion + }[difference] + + wmsg = ( + f"A charge difference of {difference} is observed " + "between the end states. This will be addressed by " + f"transforming a water into a {ion} ion" + ) + logger.info(wmsg) + + def _validate( self, stateA: ChemicalSystem, stateB: ChemicalSystem, - mapping: Optional[Union[gufe.ComponentMapping, list[gufe.ComponentMapping]]], - extends: Optional[gufe.ProtocolDAGResult] = None, - ) -> list[gufe.ProtocolUnit]: - # TODO: Extensions? + mapping: gufe.ComponentMapping | list[gufe.ComponentMapping] | None, + extends: gufe.ProtocolDAGResult | None = None, + ) -> None: + # Check we're not trying to extend if extends: - raise NotImplementedError("Can't extend simulations yet") + # This technically should be NotImplementedError + # but gufe.Protocol.validate calls `_validate` wrapped around an + # except for NotImplementedError, so we can't raise it here + raise ValueError("Can't extend simulations yet") - # Get alchemical components & validate them + mapping + # Validate the end states + self._validate_endstates(stateA, stateB) + + # Valildate the mapping alchem_comps = system_validation.get_alchemical_components(stateA, stateB) - _validate_alchemical_components(alchem_comps, mapping) - ligandmapping = mapping[0] if isinstance(mapping, list) else mapping + self._validate_mapping(mapping, alchem_comps) # Validate solvent component nonbond = self.settings.forcefield_settings.nonbonded_method @@ -638,11 +666,78 @@ def _create( # Validate protein component system_validation.validate_protein(stateA) + # Validate charge difference + # Note: validation depends on the mapping & solvent component checks + if stateA.contains(SolventComponent): + solv_comp = stateA.get_components_of_type(SolventComponent)[0] + else: + solv_comp = None + + self._validate_charge_difference( + mapping=mapping[0] if isinstance(mapping, list) else mapping, + nonbonded_method=self.settings.forcefield_settings.nonbonded_method, + explicit_charge_correction=self.settings.alchemical_settings.explicit_charge_correction, + solvent_component=solv_comp, + ) + + # Validate integrator things + settings_validation.validate_timestep( + self.settings.forcefield_settings.hydrogen_mass, + self.settings.integrator_settings.timestep, + ) + + _ = settings_validation.convert_steps_per_iteration( + simulation_settings=self.settings.simulation_settings, + integrator_settings=self.settings.integrator_settings, + ) + + _ = settings_validation.get_simsteps( + sim_length=self.settings.simulation_settings.equilibration_length, + timestep=self.settings.integrator_settings.timestep, + mc_steps=steps_per_iteration, + ) + + _ = settings_validation.get_simsteps( + sim_length=self.settings.simulation_settings.production_length, + timestep=self.settings.integrator_settings.timestep, + mc_steps=steps_per_iteration, + ) + + _ = settings_validation.convert_checkpoint_interval_to_iterations( + checkpoint_interval=self.settings.output_settings.checkpoint_interval, + time_per_iteration=self.settings.simulation_settings.time_per_iteration, + ) + + # Validate alchemical settings + # PR #125 temporarily pin lambda schedule spacing to n_replicas + if self.settings.simulation_settings.n_replicas != self.settings.lambda_settings.n_windows: + errmsg = ( + "Number of replicas in simulation_settings must equal " + "number of lambda windows in lambda_settings." + ) + raise ValueError(errmsg) + + def _create( + self, + stateA: ChemicalSystem, + stateB: ChemicalSystem, + mapping: Optional[Union[gufe.ComponentMapping, list[gufe.ComponentMapping]]], + extends: Optional[gufe.ProtocolDAGResult] = None, + ) -> list[gufe.ProtocolUnit]: + # validate inputs + self.validate(stateA=stateA, stateB=stateB, mapping=mapping, extends=extends) + + # get alchemical components and mapping + alchem_comps = system_validation.get_alchemical_components(stateA, stateB) + ligandmapping = mapping[0] if isinstance(mapping, list) else mapping + # actually create and return Units Anames = ",".join(c.name for c in alchem_comps["stateA"]) Bnames = ",".join(c.name for c in alchem_comps["stateB"]) + # our DAG has no dependencies, so just list units n_repeats = self.settings.protocol_repeats + units = [ RelativeHybridTopologyProtocolUnit( protocol=self, @@ -816,10 +911,6 @@ def run( output_settings: MultiStateOutputSettings = protocol_settings.output_settings integrator_settings: IntegratorSettings = protocol_settings.integrator_settings - # is the timestep good for the mass? - settings_validation.validate_timestep( - forcefield_settings.hydrogen_mass, integrator_settings.timestep - ) # TODO: Also validate various conversions? # Convert various time based inputs to steps/iterations steps_per_iteration = settings_validation.convert_steps_per_iteration( @@ -842,12 +933,7 @@ def run( # Get the change difference between the end states # and check if the charge correction used is appropriate - charge_difference = _get_alchemical_charge_difference( - mapping, - forcefield_settings.nonbonded_method, - alchem_settings.explicit_charge_correction, - solvent_comp, - ) + charge_difference = mapping.get_alchemical_charge_difference() # 1. Create stateA system self.logger.info("Parameterizing molecules") From 2cd56ba75a7f5b602521a94def4a72a11750dee7 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 5 Dec 2025 13:25:43 +0000 Subject: [PATCH 02/20] some fixes --- openfe/protocols/openmm_rfe/equil_rfe_methods.py | 6 +++--- .../tests/protocols/openmm_rfe/test_hybrid_top_protocol.py | 4 ---- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index f272631b8..be6c1d39a 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -485,7 +485,7 @@ def _validate_endstates( state_label = "A" if i == 0 else "B" # Check that there is only one unique Component in each state - if len(entry) != 0: + if len(entry) != 1: errmsg = ( "Only one alchemical component is allowed per end state. " f"Found {len(entry)} in state {state_label}." @@ -686,7 +686,7 @@ def _validate( self.settings.integrator_settings.timestep, ) - _ = settings_validation.convert_steps_per_iteration( + steps_per_iteration = settings_validation.convert_steps_per_iteration( simulation_settings=self.settings.simulation_settings, integrator_settings=self.settings.integrator_settings, ) @@ -710,7 +710,7 @@ def _validate( # Validate alchemical settings # PR #125 temporarily pin lambda schedule spacing to n_replicas - if self.settings.simulation_settings.n_replicas != self.settings.lambda_settings.n_windows: + if self.settings.simulation_settings.n_replicas != self.settings.lambda_settings.lambda_windows: errmsg = ( "Number of replicas in simulation_settings must equal " "number of lambda windows in lambda_settings." diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py index f5ea92cff..1f178991a 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py @@ -30,10 +30,6 @@ from openfe import setup from openfe.protocols import openmm_rfe from openfe.protocols.openmm_rfe._rfe_utils import topologyhelpers -from openfe.protocols.openmm_rfe.equil_rfe_methods import ( - _get_alchemical_charge_difference, - _validate_alchemical_components, -) from openfe.protocols.openmm_utils import omm_compute, system_creation from openfe.protocols.openmm_utils.charge_generation import ( HAS_ESPALOMA_CHARGE, From f3305622f3d147bd22464bb8726cb12551c41c4a Mon Sep 17 00:00:00 2001 From: IAlibay Date: Mon, 8 Dec 2025 13:34:41 +0000 Subject: [PATCH 03/20] move some things around --- .../protocols/openmm_rfe/equil_rfe_methods.py | 74 +++- .../openmm_rfe/test_hybrid_top_validation.py | 391 ++++++++++++++++++ 2 files changed, 445 insertions(+), 20 deletions(-) create mode 100644 openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index be6c1d39a..59e33ef52 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -635,6 +635,55 @@ def _validate_charge_difference( ) logger.info(wmsg) + @staticmethod + def _validate_simulation_settings( + simulation_settings, + integrator_settings, + output_settings, + ): + + steps_per_iteration = settings_validation.convert_steps_per_iteration( + simulation_settings=simulation_settings, + integrator_settings=integrator_settings, + ) + + _ = settings_validation.get_simsteps( + sim_length=simulation_settings.equilibration_length, + timestep=integrator_settings.timestep, + mc_steps=steps_per_iteration, + ) + + _ = settings_validation.get_simsteps( + sim_length=simulation_settings.production_length, + timestep=integrator_settings.timestep, + mc_steps=steps_per_iteration, + ) + + _ = settings_validation.convert_checkpoint_interval_to_iterations( + checkpoint_interval=output_settings.checkpoint_interval, + time_per_iteration=simulation_settings.time_per_iteration, + ) + + if output_settings.positions_write_frequency is not None: + _ = settings_validation.divmod_time_and_check( + numerator=output_settings.positions_write_frequency, + denominator=sampler_settings.time_per_iteration, + numerator_name="output settings' position_write_frequency", + denominator_name="sampler settings' time_per_iteration", + ) + + if output_settings.velocities_write_frequency is not None: + _ = settings_validation.divmod_time_and_check( + numerator=output_settings.velocities_write_frequency, + denominator=sampler_settings.time_per_iteration, + numerator_name="output settings' velocity_write_frequency", + denominator_name="sampler settings' time_per_iteration", + ) + + _, _ = settings_validation.convert_real_time_analysis_iterations( + simulation_settings=sampler_settings, + ) + def _validate( self, stateA: ChemicalSystem, @@ -686,26 +735,11 @@ def _validate( self.settings.integrator_settings.timestep, ) - steps_per_iteration = settings_validation.convert_steps_per_iteration( - simulation_settings=self.settings.simulation_settings, - integrator_settings=self.settings.integrator_settings, - ) - - _ = settings_validation.get_simsteps( - sim_length=self.settings.simulation_settings.equilibration_length, - timestep=self.settings.integrator_settings.timestep, - mc_steps=steps_per_iteration, - ) - - _ = settings_validation.get_simsteps( - sim_length=self.settings.simulation_settings.production_length, - timestep=self.settings.integrator_settings.timestep, - mc_steps=steps_per_iteration, - ) - - _ = settings_validation.convert_checkpoint_interval_to_iterations( - checkpoint_interval=self.settings.output_settings.checkpoint_interval, - time_per_iteration=self.settings.simulation_settings.time_per_iteration, + # Validate simulation & output settings + self._validate_simulation_settings( + self.settings.simulation_settings, + self.settings.integrator_settings, + self.settings.output_settings, ) # Validate alchemical settings diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py new file mode 100644 index 000000000..62c8e5d23 --- /dev/null +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -0,0 +1,391 @@ +# This code is part of OpenFE and is licensed under the MIT license. +# For details, see https://github.com/OpenFreeEnergy/openfe +import copy +import json +import sys +import xml.etree.ElementTree as ET +from importlib import resources +from math import sqrt +from pathlib import Path +from unittest import mock + +import gufe +import mdtraj as mdt +import numpy as np +import pytest +from kartograf import KartografAtomMapper +from kartograf.atom_aligner import align_mol_shape +from numpy.testing import assert_allclose +from openff.toolkit import Molecule +from openff.units import unit +from openff.units.openmm import ensure_quantity, from_openmm, to_openmm +from openmm import CustomNonbondedForce, MonteCarloBarostat, NonbondedForce, XmlSerializer, app +from openmm import unit as omm_unit +from openmmforcefields.generators import SMIRNOFFTemplateGenerator +from openmmtools.multistate.multistatesampler import MultiStateSampler +from rdkit import Chem +from rdkit.Geometry import Point3D + +import openfe +from openfe import setup +from openfe.protocols import openmm_rfe +from openfe.protocols.openmm_rfe._rfe_utils import topologyhelpers +from openfe.protocols.openmm_utils import omm_compute, system_creation +from openfe.protocols.openmm_utils.charge_generation import ( + HAS_ESPALOMA_CHARGE, + HAS_NAGL, + HAS_OPENEYE, +) + + +@pytest.fixture() +def vac_settings(): + settings = openmm_rfe.RelativeHybridTopologyProtocol.default_settings() + settings.forcefield_settings.nonbonded_method = "nocutoff" + settings.engine_settings.compute_platform = None + settings.protocol_repeats = 1 + return settings + + +@pytest.fixture() +def solv_settings(): + settings = openmm_rfe.RelativeHybridTopologyProtocol.default_settings() + settings.engine_settings.compute_platform = None + settings.protocol_repeats = 1 + return settings + + +def test_invalid_protocol_repeats(): + settings = openmm_rfe.RelativeHybridTopologyProtocol.default_settings() + with pytest.raises(ValueError, match="must be a positive value"): + settings.protocol_repeats = -1 + + +@pytest.mark.parametrize( + "mapping", + [None, [], ["A", "B"]], +) +def test_validate_alchemical_components_wrong_mappings(mapping): + with pytest.raises(ValueError, match="A single LigandAtomMapping"): + _validate_alchemical_components({"stateA": [], "stateB": []}, mapping) + + +def test_validate_alchemical_components_missing_alchem_comp(benzene_to_toluene_mapping): + alchem_comps = {"stateA": [openfe.SolventComponent()], "stateB": []} + with pytest.raises(ValueError, match="Unmapped alchemical component"): + _validate_alchemical_components(alchem_comps, benzene_to_toluene_mapping) + + +def test_hightimestep( + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + vac_settings, + tmpdir, +): + vac_settings.forcefield_settings.hydrogen_mass = 1.0 + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=vac_settings, + ) + + dag = p.create( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + ) + dag_unit = list(dag.protocol_units)[0] + + errmsg = "too large for hydrogen mass" + with tmpdir.as_cwd(): + with pytest.raises(ValueError, match=errmsg): + dag_unit.run(dry=True) + + +def test_n_replicas_not_n_windows( + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + vac_settings, + tmpdir, +): + # For PR #125 we pin such that the number of lambda windows + # equals the numbers of replicas used - TODO: remove limitation + # default lambda windows is 11 + vac_settings.simulation_settings.n_replicas = 13 + + errmsg = "Number of replicas 13 does not equal the number of lambda windows 11" + + with tmpdir.as_cwd(): + with pytest.raises(ValueError, match=errmsg): + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=vac_settings, + ) + dag = p.create( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + ) + dag_unit = list(dag.protocol_units)[0] + dag_unit.run(dry=True) + + +def test_missing_ligand(benzene_system, benzene_to_toluene_mapping): + # state B doesn't have a ligand component + stateB = openfe.ChemicalSystem({"solvent": openfe.SolventComponent()}) + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + + match_str = "missing alchemical components in stateB" + with pytest.raises(ValueError, match=match_str): + _ = p.create( + stateA=benzene_system, + stateB=stateB, + mapping=benzene_to_toluene_mapping, + ) + + +def test_vaccuum_PME_error( + benzene_vacuum_system, benzene_modifications, benzene_to_toluene_mapping +): + # state B doesn't have a solvent component (i.e. its vacuum) + stateB = openfe.ChemicalSystem({"ligand": benzene_modifications["toluene"]}) + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + errmsg = "PME cannot be used for vacuum transform" + with pytest.raises(ValueError, match=errmsg): + _ = p.create( + stateA=benzene_vacuum_system, + stateB=stateB, + mapping=benzene_to_toluene_mapping, + ) + + +def test_incompatible_solvent(benzene_system, benzene_modifications, benzene_to_toluene_mapping): + # the solvents are different + stateB = openfe.ChemicalSystem( + { + "ligand": benzene_modifications["toluene"], + "solvent": openfe.SolventComponent(positive_ion="K", negative_ion="Cl"), + } + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + # We don't have a way to map non-ligand components so for now it + # just triggers that it's not a mapped component + errmsg = "missing alchemical components in stateA" + with pytest.raises(ValueError, match=errmsg): + _ = p.create( + stateA=benzene_system, + stateB=stateB, + mapping=benzene_to_toluene_mapping, + ) + + +def test_mapping_mismatch_A(benzene_system, toluene_system, benzene_modifications): + # the atom mapping doesn't refer to the ligands in the systems + mapping = setup.LigandAtomMapping( + componentA=benzene_system.components["ligand"], + componentB=benzene_modifications["phenol"], + componentA_to_componentB=dict(), + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + errmsg = ( + r"Unmapped alchemical component " + r"SmallMoleculeComponent\(name=toluene\)" + ) + with pytest.raises(ValueError, match=errmsg): + _ = p.create( + stateA=benzene_system, + stateB=toluene_system, + mapping=mapping, + ) + + +def test_mapping_mismatch_B(benzene_system, toluene_system, benzene_modifications): + mapping = setup.LigandAtomMapping( + componentA=benzene_modifications["phenol"], + componentB=toluene_system.components["ligand"], + componentA_to_componentB=dict(), + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + errmsg = ( + r"Unmapped alchemical component " + r"SmallMoleculeComponent\(name=benzene\)" + ) + with pytest.raises(ValueError, match=errmsg): + _ = p.create( + stateA=benzene_system, + stateB=toluene_system, + mapping=mapping, + ) + + +def test_complex_mismatch(benzene_system, toluene_complex_system, benzene_to_toluene_mapping): + # only one complex + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + with pytest.raises(ValueError): + _ = p.create( + stateA=benzene_system, + stateB=toluene_complex_system, + mapping=benzene_to_toluene_mapping, + ) + + +def test_too_many_specified_mappings(benzene_system, toluene_system, benzene_to_toluene_mapping): + # mapping dict requires 'ligand' key + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + errmsg = "A single LigandAtomMapping is expected for this Protocol" + with pytest.raises(ValueError, match=errmsg): + _ = p.create( + stateA=benzene_system, + stateB=toluene_system, + mapping=[benzene_to_toluene_mapping, benzene_to_toluene_mapping], + ) + + +def test_protein_mismatch( + benzene_complex_system, toluene_complex_system, benzene_to_toluene_mapping +): + # hack one protein to be labelled differently + prot = toluene_complex_system["protein"] + alt_prot = openfe.ProteinComponent(prot.to_rdkit(), name="Mickey Mouse") + alt_toluene_complex_system = openfe.ChemicalSystem( + { + "ligand": toluene_complex_system["ligand"], + "solvent": toluene_complex_system["solvent"], + "protein": alt_prot, + } + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + with pytest.raises(ValueError): + _ = p.create( + stateA=benzene_complex_system, + stateB=alt_toluene_complex_system, + mapping=benzene_to_toluene_mapping, + ) + + +def test_element_change_warning(atom_mapping_basic_test_files): + # check a mapping with element change gets rejected early + l1 = atom_mapping_basic_test_files["2-methylnaphthalene"] + l2 = atom_mapping_basic_test_files["2-naftanol"] + + # We use the 'old' lomap defaults because the + # basic test files inputs we use aren't fully aligned + mapper = setup.LomapAtomMapper( + time=20, threed=True, max3d=1000.0, element_change=True, seed="", shift=True + ) + + mapping = next(mapper.suggest_mappings(l1, l2)) + + sys1 = openfe.ChemicalSystem( + {"ligand": l1, "solvent": openfe.SolventComponent()}, + ) + sys2 = openfe.ChemicalSystem( + {"ligand": l2, "solvent": openfe.SolventComponent()}, + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol( + settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + ) + with pytest.warns(UserWarning, match="Element change"): + _ = p.create( + stateA=sys1, + stateB=sys2, + mapping=mapping, + ) + + +@pytest.mark.parametrize( + "mapping_name,result", + [ + ["benzene_to_toluene_mapping", 0], + ["benzene_to_benzoic_mapping", 1], + ["benzene_to_aniline_mapping", -1], + ["aniline_to_benzene_mapping", 1], + ], +) +def test_get_charge_difference(mapping_name, result, request): + mapping = request.getfixturevalue(mapping_name) + if result != 0: + ion = r"Na\+" if result == -1 else r"Cl\-" + wmsg = ( + f"A charge difference of {result} is observed " + "between the end states. This will be addressed by " + f"transforming a water into a {ion} ion" + ) + with pytest.warns(UserWarning, match=wmsg): + val = _get_alchemical_charge_difference(mapping, "pme", True, openfe.SolventComponent()) + assert result == pytest.approx(val) + else: + val = _get_alchemical_charge_difference(mapping, "pme", True, openfe.SolventComponent()) + assert result == pytest.approx(val) + + +def test_get_charge_difference_no_pme(benzene_to_benzoic_mapping): + errmsg = "Explicit charge correction when not using PME" + with pytest.raises(ValueError, match=errmsg): + _get_alchemical_charge_difference( + benzene_to_benzoic_mapping, + "nocutoff", + True, + openfe.SolventComponent(), + ) + + +def test_get_charge_difference_no_corr(benzene_to_benzoic_mapping): + wmsg = ( + "A charge difference of 1 is observed between the end states. " + "No charge correction has been requested" + ) + with pytest.warns(UserWarning, match=wmsg): + _get_alchemical_charge_difference( + benzene_to_benzoic_mapping, + "pme", + False, + openfe.SolventComponent(), + ) + + +def test_greater_than_one_charge_difference_error(aniline_to_benzoic_mapping): + errmsg = "A charge difference of 2" + with pytest.raises(ValueError, match=errmsg): + _get_alchemical_charge_difference( + aniline_to_benzoic_mapping, + "pme", + True, + openfe.SolventComponent(), + ) + + +def test_get_alchemical_waters_no_waters( + benzene_solvent_openmm_system, +): + system, topology, positions = benzene_solvent_openmm_system + + errmsg = "There are no waters" + + with pytest.raises(ValueError, match=errmsg): + topologyhelpers.get_alchemical_waters( + topology, positions, charge_difference=1, distance_cutoff=3.0 * unit.nanometer + ) From 1e0153e949fbeb4a947ffb7beafd614eefaebc80 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Mon, 15 Dec 2025 10:26:58 -1000 Subject: [PATCH 04/20] add validate endstate tests --- .../openmm_rfe/test_hybrid_top_validation.py | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 62c8e5d23..50a518606 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -61,6 +61,51 @@ def test_invalid_protocol_repeats(): settings.protocol_repeats = -1 +@pytest.mark.parametrize('state', ['A', 'B']) +def test_endstate_two_alchemcomp_stateA(state, benzene_modifications): + first_state = openfe.ChemicalSystem({ + 'ligandA': benzene_modifications['benzene'], + 'ligandB': benzene_modifications['toluene'], + 'solvent': openfe.SolventComponent(), + }) + other_state = openfe.ChemicalSystem({ + 'ligandC': benzene_modifications['phenol'], + 'solvent': openfe.SolventComponent(), + }) + + if state == 'A': + args = (first_state, other_state) + else: + args = (other_state, first_state) + + with pytest.raises(ValueError, match="Only one alchemical component"): + openmm_rfe.RelativeHybridTopologyProtocol._validate_endstates( + *args + ) + +@pytest.mark.parametrize('state', ['A', 'B']) +def test_endstates_not_smc(state, benzene_modifications): + first_state = openfe.ChemicalSystem({ + 'ligand': benzene_modifications['benzene'], + 'foo': openfe.SolventComponent(), + }) + other_state = openfe.ChemicalSystem({ + 'ligand': benzene_modifications['benzene'], + 'foo': benzene_modifications['toluene'], + }) + + if state == 'A': + args = (first_state, other_state) + else: + args = (other_state, first_state) + + errmsg = "only SmallMoleculeComponents transformations" + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_endstates( + *args + ) + + @pytest.mark.parametrize( "mapping", [None, [], ["A", "B"]], From fbc455416c5b3eb9683c31515deec80b2d664da4 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Mon, 15 Dec 2025 10:52:30 -1000 Subject: [PATCH 05/20] validate mapping tests --- .../openmm_rfe/test_hybrid_top_validation.py | 84 ++++++++++++------- 1 file changed, 53 insertions(+), 31 deletions(-) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 50a518606..241be3ebd 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -106,6 +106,59 @@ def test_endstates_not_smc(state, benzene_modifications): ) +def test_validate_mapping_none_mapping(): + errmsg = "A single LigandAtomMapping is expected" + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_mapping(None, None) + + +def test_validate_mapping_multi_mapping(benzene_to_toluene_mapping): + errmsg = "A single LigandAtomMapping is expected" + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_mapping( + [benzene_to_toluene_mapping] * 2, + None + ) + + +@pytest.mark.parametrize('state', ['A', 'B']) +def test_validate_mapping_alchem_not_in(state, benzene_to_toluene_mapping): + errmsg = f"not in alchemical components of state{state}" + + if state == "A": + alchem_comps = {"stateA": [], "stateB": [benzene_to_toluene_mapping.componentB]} + else: + alchem_comps = {"stateA": [benzene_to_toluene_mapping.componentA], "stateB": []} + + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_mapping( + [benzene_to_toluene_mapping], + alchem_comps, + ) + + +def test_element_change_warning(atom_mapping_basic_test_files): + # check a mapping with element change gets rejected early + l1 = atom_mapping_basic_test_files["2-methylnaphthalene"] + l2 = atom_mapping_basic_test_files["2-naftanol"] + + # We use the 'old' lomap defaults because the + # basic test files inputs we use aren't fully aligned + mapper = setup.LomapAtomMapper( + time=20, threed=True, max3d=1000.0, element_change=True, seed="", shift=True + ) + + mapping = next(mapper.suggest_mappings(l1, l2)) + + alchem_comps = {"stateA": [l1], "stateB": [l2]} + + with pytest.warns(UserWarning, match="Element change"): + openmm_rfe.RelativeHybridTopologyProtocol._validate_mapping( + [mapping], + alchem_comps, + ) + + @pytest.mark.parametrize( "mapping", [None, [], ["A", "B"]], @@ -330,37 +383,6 @@ def test_protein_mismatch( ) -def test_element_change_warning(atom_mapping_basic_test_files): - # check a mapping with element change gets rejected early - l1 = atom_mapping_basic_test_files["2-methylnaphthalene"] - l2 = atom_mapping_basic_test_files["2-naftanol"] - - # We use the 'old' lomap defaults because the - # basic test files inputs we use aren't fully aligned - mapper = setup.LomapAtomMapper( - time=20, threed=True, max3d=1000.0, element_change=True, seed="", shift=True - ) - - mapping = next(mapper.suggest_mappings(l1, l2)) - - sys1 = openfe.ChemicalSystem( - {"ligand": l1, "solvent": openfe.SolventComponent()}, - ) - sys2 = openfe.ChemicalSystem( - {"ligand": l2, "solvent": openfe.SolventComponent()}, - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - with pytest.warns(UserWarning, match="Element change"): - _ = p.create( - stateA=sys1, - stateB=sys2, - mapping=mapping, - ) - - @pytest.mark.parametrize( "mapping_name,result", [ From c2f49d25af7ead72a8f31ad22224bf5510e6b500 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Mon, 15 Dec 2025 11:45:04 -1000 Subject: [PATCH 06/20] net charge validation tests --- .../protocols/openmm_rfe/equil_rfe_methods.py | 4 + .../openmm_rfe/test_hybrid_top_validation.py | 186 ++++++++---------- 2 files changed, 83 insertions(+), 107 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 59e33ef52..e20c7360b 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -608,6 +608,10 @@ def _validate_charge_difference( warnings.warn(wmsg) return + if solvent_component is None: + errmsg = "Cannot use eplicit charge correction without solvent" + raise ValueError(errmsg) + # We implicitly check earlier that we have to have pme for a solvated # system, so we only need to check the nonbonded method here if nonbonded_method.lower() != "pme": diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 241be3ebd..f2b64d292 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -1,5 +1,6 @@ # This code is part of OpenFE and is licensed under the MIT license. # For details, see https://github.com/OpenFreeEnergy/openfe +import logging import copy import json import sys @@ -159,19 +160,87 @@ def test_element_change_warning(atom_mapping_basic_test_files): ) +def test_charge_difference_no_corr(benzene_to_benzoic_mapping): + wmsg = ( + "A charge difference of 1 is observed between the end states. " + "No charge correction has been requested" + ) + + with pytest.warns(UserWarning, match=wmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference( + benzene_to_benzoic_mapping, + "pme", + False, + openfe.SolventComponent(), + ) + + +def test_charge_difference_no_solvent(benzene_to_benzoic_mapping): + errmsg = "Cannot use eplicit charge correction without solvent" + + with pytest.raises(ValueError, errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference( + benzene_to_benzoic_mapping, + "pme", + True, + None, + ) + + +def test_charge_difference_no_pme(benzene_to_benzoic_mapping): + errmsg = "Explicit charge correction when not using PME" + + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference( + benzene_to_benzoic_mapping, + "nocutoff", + True, + openfe.SolventComponent(), + ) + + +def test_greater_than_one_charge_difference_error(aniline_to_benzoic_mapping): + errmsg = "A charge difference of 2" + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference( + aniline_to_benzoic_mapping, + "pme", + True, + openfe.SolventComponent(), + ) + + @pytest.mark.parametrize( - "mapping", - [None, [], ["A", "B"]], + "mapping_name,result", + [ + ["benzene_to_toluene_mapping", 0], + ["benzene_to_benzoic_mapping", 1], + ["benzene_to_aniline_mapping", -1], + ["aniline_to_benzene_mapping", 1], + ], ) -def test_validate_alchemical_components_wrong_mappings(mapping): - with pytest.raises(ValueError, match="A single LigandAtomMapping"): - _validate_alchemical_components({"stateA": [], "stateB": []}, mapping) - +def test_get_charge_difference(mapping_name, result, request, caplog): + mapping = request.getfixturevalue(mapping_name) + caplog.set_level(logging.INFO) + + ion = r"Na\+" if result == -1 else r"Cl\-" + msg = ( + f"A charge difference of {result} is observed " + "between the end states. This will be addressed by " + f"transforming a water into a {ion} ion" + ) + + openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference( + mapping, + "pme", + True, + openfe.SolventComponent() + ) -def test_validate_alchemical_components_missing_alchem_comp(benzene_to_toluene_mapping): - alchem_comps = {"stateA": [openfe.SolventComponent()], "stateB": []} - with pytest.raises(ValueError, match="Unmapped alchemical component"): - _validate_alchemical_components(alchem_comps, benzene_to_toluene_mapping) + if result != 0: + assert msg in caplog.text + else: + assert msg not in caplog.text def test_hightimestep( @@ -286,78 +355,6 @@ def test_incompatible_solvent(benzene_system, benzene_modifications, benzene_to_ ) -def test_mapping_mismatch_A(benzene_system, toluene_system, benzene_modifications): - # the atom mapping doesn't refer to the ligands in the systems - mapping = setup.LigandAtomMapping( - componentA=benzene_system.components["ligand"], - componentB=benzene_modifications["phenol"], - componentA_to_componentB=dict(), - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - errmsg = ( - r"Unmapped alchemical component " - r"SmallMoleculeComponent\(name=toluene\)" - ) - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=toluene_system, - mapping=mapping, - ) - - -def test_mapping_mismatch_B(benzene_system, toluene_system, benzene_modifications): - mapping = setup.LigandAtomMapping( - componentA=benzene_modifications["phenol"], - componentB=toluene_system.components["ligand"], - componentA_to_componentB=dict(), - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - errmsg = ( - r"Unmapped alchemical component " - r"SmallMoleculeComponent\(name=benzene\)" - ) - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=toluene_system, - mapping=mapping, - ) - - -def test_complex_mismatch(benzene_system, toluene_complex_system, benzene_to_toluene_mapping): - # only one complex - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - with pytest.raises(ValueError): - _ = p.create( - stateA=benzene_system, - stateB=toluene_complex_system, - mapping=benzene_to_toluene_mapping, - ) - - -def test_too_many_specified_mappings(benzene_system, toluene_system, benzene_to_toluene_mapping): - # mapping dict requires 'ligand' key - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - errmsg = "A single LigandAtomMapping is expected for this Protocol" - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=toluene_system, - mapping=[benzene_to_toluene_mapping, benzene_to_toluene_mapping], - ) - - def test_protein_mismatch( benzene_complex_system, toluene_complex_system, benzene_to_toluene_mapping ): @@ -409,31 +406,6 @@ def test_get_charge_difference(mapping_name, result, request): assert result == pytest.approx(val) -def test_get_charge_difference_no_pme(benzene_to_benzoic_mapping): - errmsg = "Explicit charge correction when not using PME" - with pytest.raises(ValueError, match=errmsg): - _get_alchemical_charge_difference( - benzene_to_benzoic_mapping, - "nocutoff", - True, - openfe.SolventComponent(), - ) - - -def test_get_charge_difference_no_corr(benzene_to_benzoic_mapping): - wmsg = ( - "A charge difference of 1 is observed between the end states. " - "No charge correction has been requested" - ) - with pytest.warns(UserWarning, match=wmsg): - _get_alchemical_charge_difference( - benzene_to_benzoic_mapping, - "pme", - False, - openfe.SolventComponent(), - ) - - def test_greater_than_one_charge_difference_error(aniline_to_benzoic_mapping): errmsg = "A charge difference of 2" with pytest.raises(ValueError, match=errmsg): From c50f99cad7daf09272bd81033d080d10cc5193af Mon Sep 17 00:00:00 2001 From: IAlibay Date: Mon, 22 Dec 2025 10:23:47 -1000 Subject: [PATCH 07/20] more stuff --- .../openmm_rfe/test_hybrid_top_validation.py | 102 ------------------ 1 file changed, 102 deletions(-) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index f2b64d292..d0452eb9a 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -297,23 +297,6 @@ def test_n_replicas_not_n_windows( dag_unit.run(dry=True) -def test_missing_ligand(benzene_system, benzene_to_toluene_mapping): - # state B doesn't have a ligand component - stateB = openfe.ChemicalSystem({"solvent": openfe.SolventComponent()}) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - - match_str = "missing alchemical components in stateB" - with pytest.raises(ValueError, match=match_str): - _ = p.create( - stateA=benzene_system, - stateB=stateB, - mapping=benzene_to_toluene_mapping, - ) - - def test_vaccuum_PME_error( benzene_vacuum_system, benzene_modifications, benzene_to_toluene_mapping ): @@ -332,91 +315,6 @@ def test_vaccuum_PME_error( ) -def test_incompatible_solvent(benzene_system, benzene_modifications, benzene_to_toluene_mapping): - # the solvents are different - stateB = openfe.ChemicalSystem( - { - "ligand": benzene_modifications["toluene"], - "solvent": openfe.SolventComponent(positive_ion="K", negative_ion="Cl"), - } - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - # We don't have a way to map non-ligand components so for now it - # just triggers that it's not a mapped component - errmsg = "missing alchemical components in stateA" - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=stateB, - mapping=benzene_to_toluene_mapping, - ) - - -def test_protein_mismatch( - benzene_complex_system, toluene_complex_system, benzene_to_toluene_mapping -): - # hack one protein to be labelled differently - prot = toluene_complex_system["protein"] - alt_prot = openfe.ProteinComponent(prot.to_rdkit(), name="Mickey Mouse") - alt_toluene_complex_system = openfe.ChemicalSystem( - { - "ligand": toluene_complex_system["ligand"], - "solvent": toluene_complex_system["solvent"], - "protein": alt_prot, - } - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - with pytest.raises(ValueError): - _ = p.create( - stateA=benzene_complex_system, - stateB=alt_toluene_complex_system, - mapping=benzene_to_toluene_mapping, - ) - - -@pytest.mark.parametrize( - "mapping_name,result", - [ - ["benzene_to_toluene_mapping", 0], - ["benzene_to_benzoic_mapping", 1], - ["benzene_to_aniline_mapping", -1], - ["aniline_to_benzene_mapping", 1], - ], -) -def test_get_charge_difference(mapping_name, result, request): - mapping = request.getfixturevalue(mapping_name) - if result != 0: - ion = r"Na\+" if result == -1 else r"Cl\-" - wmsg = ( - f"A charge difference of {result} is observed " - "between the end states. This will be addressed by " - f"transforming a water into a {ion} ion" - ) - with pytest.warns(UserWarning, match=wmsg): - val = _get_alchemical_charge_difference(mapping, "pme", True, openfe.SolventComponent()) - assert result == pytest.approx(val) - else: - val = _get_alchemical_charge_difference(mapping, "pme", True, openfe.SolventComponent()) - assert result == pytest.approx(val) - - -def test_greater_than_one_charge_difference_error(aniline_to_benzoic_mapping): - errmsg = "A charge difference of 2" - with pytest.raises(ValueError, match=errmsg): - _get_alchemical_charge_difference( - aniline_to_benzoic_mapping, - "pme", - True, - openfe.SolventComponent(), - ) - - def test_get_alchemical_waters_no_waters( benzene_solvent_openmm_system, ): From 9e0d29be83df2a7f747d6bf385d0853c1d92589e Mon Sep 17 00:00:00 2001 From: IAlibay Date: Wed, 24 Dec 2025 00:43:36 -0500 Subject: [PATCH 08/20] remove old tests --- .../openmm_rfe/test_hybrid_top_protocol.py | 260 ------------------ 1 file changed, 260 deletions(-) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py index 1f178991a..148f32fe1 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py @@ -192,21 +192,6 @@ def test_create_independent_repeat_ids(benzene_system, toluene_system, benzene_t assert len(repeat_ids) == 6 -@pytest.mark.parametrize( - "mapping", - [None, [], ["A", "B"]], -) -def test_validate_alchemical_components_wrong_mappings(mapping): - with pytest.raises(ValueError, match="A single LigandAtomMapping"): - _validate_alchemical_components({"stateA": [], "stateB": []}, mapping) - - -def test_validate_alchemical_components_missing_alchem_comp(benzene_to_toluene_mapping): - alchem_comps = {"stateA": [openfe.SolventComponent()], "stateB": []} - with pytest.raises(ValueError, match="Unmapped alchemical component"): - _validate_alchemical_components(alchem_comps, benzene_to_toluene_mapping) - - @pytest.mark.parametrize("method", ["repex", "sams", "independent", "InDePeNdENT"]) def test_dry_run_default_vacuum( benzene_vacuum_system, @@ -989,189 +974,6 @@ def test_hightimestep( dag_unit.run(dry=True) -def test_n_replicas_not_n_windows( - benzene_vacuum_system, - toluene_vacuum_system, - benzene_to_toluene_mapping, - vac_settings, - tmpdir, -): - # For PR #125 we pin such that the number of lambda windows - # equals the numbers of replicas used - TODO: remove limitation - # default lambda windows is 11 - vac_settings.simulation_settings.n_replicas = 13 - - errmsg = "Number of replicas 13 does not equal the number of lambda windows 11" - - with tmpdir.as_cwd(): - with pytest.raises(ValueError, match=errmsg): - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=vac_settings, - ) - dag = p.create( - stateA=benzene_vacuum_system, - stateB=toluene_vacuum_system, - mapping=benzene_to_toluene_mapping, - ) - dag_unit = list(dag.protocol_units)[0] - dag_unit.run(dry=True) - - -def test_missing_ligand(benzene_system, benzene_to_toluene_mapping): - # state B doesn't have a ligand component - stateB = openfe.ChemicalSystem({"solvent": openfe.SolventComponent()}) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - - match_str = "missing alchemical components in stateB" - with pytest.raises(ValueError, match=match_str): - _ = p.create( - stateA=benzene_system, - stateB=stateB, - mapping=benzene_to_toluene_mapping, - ) - - -def test_vaccuum_PME_error( - benzene_vacuum_system, benzene_modifications, benzene_to_toluene_mapping -): - # state B doesn't have a solvent component (i.e. its vacuum) - stateB = openfe.ChemicalSystem({"ligand": benzene_modifications["toluene"]}) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - errmsg = "PME cannot be used for vacuum transform" - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_vacuum_system, - stateB=stateB, - mapping=benzene_to_toluene_mapping, - ) - - -def test_incompatible_solvent(benzene_system, benzene_modifications, benzene_to_toluene_mapping): - # the solvents are different - stateB = openfe.ChemicalSystem( - { - "ligand": benzene_modifications["toluene"], - "solvent": openfe.SolventComponent(positive_ion="K", negative_ion="Cl"), - } - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - # We don't have a way to map non-ligand components so for now it - # just triggers that it's not a mapped component - errmsg = "missing alchemical components in stateA" - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=stateB, - mapping=benzene_to_toluene_mapping, - ) - - -def test_mapping_mismatch_A(benzene_system, toluene_system, benzene_modifications): - # the atom mapping doesn't refer to the ligands in the systems - mapping = setup.LigandAtomMapping( - componentA=benzene_system.components["ligand"], - componentB=benzene_modifications["phenol"], - componentA_to_componentB=dict(), - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - errmsg = ( - r"Unmapped alchemical component " - r"SmallMoleculeComponent\(name=toluene\)" - ) - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=toluene_system, - mapping=mapping, - ) - - -def test_mapping_mismatch_B(benzene_system, toluene_system, benzene_modifications): - mapping = setup.LigandAtomMapping( - componentA=benzene_modifications["phenol"], - componentB=toluene_system.components["ligand"], - componentA_to_componentB=dict(), - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - errmsg = ( - r"Unmapped alchemical component " - r"SmallMoleculeComponent\(name=benzene\)" - ) - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=toluene_system, - mapping=mapping, - ) - - -def test_complex_mismatch(benzene_system, toluene_complex_system, benzene_to_toluene_mapping): - # only one complex - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - with pytest.raises(ValueError): - _ = p.create( - stateA=benzene_system, - stateB=toluene_complex_system, - mapping=benzene_to_toluene_mapping, - ) - - -def test_too_many_specified_mappings(benzene_system, toluene_system, benzene_to_toluene_mapping): - # mapping dict requires 'ligand' key - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - errmsg = "A single LigandAtomMapping is expected for this Protocol" - with pytest.raises(ValueError, match=errmsg): - _ = p.create( - stateA=benzene_system, - stateB=toluene_system, - mapping=[benzene_to_toluene_mapping, benzene_to_toluene_mapping], - ) - - -def test_protein_mismatch( - benzene_complex_system, toluene_complex_system, benzene_to_toluene_mapping -): - # hack one protein to be labelled differently - prot = toluene_complex_system["protein"] - alt_prot = openfe.ProteinComponent(prot.to_rdkit(), name="Mickey Mouse") - alt_toluene_complex_system = openfe.ChemicalSystem( - { - "ligand": toluene_complex_system["ligand"], - "solvent": toluene_complex_system["solvent"], - "protein": alt_prot, - } - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - with pytest.raises(ValueError): - _ = p.create( - stateA=benzene_complex_system, - stateB=alt_toluene_complex_system, - mapping=benzene_to_toluene_mapping, - ) - - def test_element_change_warning(atom_mapping_basic_test_files): # check a mapping with element change gets rejected early l1 = atom_mapping_basic_test_files["2-methylnaphthalene"] @@ -1745,68 +1547,6 @@ def test_filenotfound_replica_states(self, protocolresult): protocolresult.get_replica_states() -@pytest.mark.parametrize( - "mapping_name,result", - [ - ["benzene_to_toluene_mapping", 0], - ["benzene_to_benzoic_mapping", 1], - ["benzene_to_aniline_mapping", -1], - ["aniline_to_benzene_mapping", 1], - ], -) -def test_get_charge_difference(mapping_name, result, request): - mapping = request.getfixturevalue(mapping_name) - if result != 0: - ion = r"Na\+" if result == -1 else r"Cl\-" - wmsg = ( - f"A charge difference of {result} is observed " - "between the end states. This will be addressed by " - f"transforming a water into a {ion} ion" - ) - with pytest.warns(UserWarning, match=wmsg): - val = _get_alchemical_charge_difference(mapping, "pme", True, openfe.SolventComponent()) - assert result == pytest.approx(val) - else: - val = _get_alchemical_charge_difference(mapping, "pme", True, openfe.SolventComponent()) - assert result == pytest.approx(val) - - -def test_get_charge_difference_no_pme(benzene_to_benzoic_mapping): - errmsg = "Explicit charge correction when not using PME" - with pytest.raises(ValueError, match=errmsg): - _get_alchemical_charge_difference( - benzene_to_benzoic_mapping, - "nocutoff", - True, - openfe.SolventComponent(), - ) - - -def test_get_charge_difference_no_corr(benzene_to_benzoic_mapping): - wmsg = ( - "A charge difference of 1 is observed between the end states. " - "No charge correction has been requested" - ) - with pytest.warns(UserWarning, match=wmsg): - _get_alchemical_charge_difference( - benzene_to_benzoic_mapping, - "pme", - False, - openfe.SolventComponent(), - ) - - -def test_greater_than_one_charge_difference_error(aniline_to_benzoic_mapping): - errmsg = "A charge difference of 2" - with pytest.raises(ValueError, match=errmsg): - _get_alchemical_charge_difference( - aniline_to_benzoic_mapping, - "pme", - True, - openfe.SolventComponent(), - ) - - @pytest.fixture(scope="session") def benzene_solvent_openmm_system(benzene_modifications): smc = benzene_modifications["benzene"] From b6d5ecd315b6ca186c62494e8aed37e4e0f69cde Mon Sep 17 00:00:00 2001 From: IAlibay Date: Thu, 25 Dec 2025 19:27:36 -0500 Subject: [PATCH 09/20] Fix up the one test --- .../protocols/openmm_rfe/equil_rfe_methods.py | 30 +++++++++++++++---- 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index e20c7360b..356423922 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -641,10 +641,28 @@ def _validate_charge_difference( @staticmethod def _validate_simulation_settings( - simulation_settings, - integrator_settings, - output_settings, + simulation_settings: MultiStateSimulationSettings, + integrator_settings: IntegratorSettings, + output_settings: MultiStateOutputSettings, ): + """ + Validate various simulation settings, including but not limited to + timestep conversions, and output file write frequencies. + + Parameters + ---------- + simulation_settings : MultiStateSimulationSettings + The sampler simulation settings. + integrator_settings : IntegratorSettings + Settings defining the behaviour of the integrator. + output_settings : MultiStateOutputSettings + Settings defining the simulation file writing behaviour. + + Raises + ------ + ValueError + * If the + """ steps_per_iteration = settings_validation.convert_steps_per_iteration( simulation_settings=simulation_settings, @@ -671,7 +689,7 @@ def _validate_simulation_settings( if output_settings.positions_write_frequency is not None: _ = settings_validation.divmod_time_and_check( numerator=output_settings.positions_write_frequency, - denominator=sampler_settings.time_per_iteration, + denominator=simulation_settings.time_per_iteration, numerator_name="output settings' position_write_frequency", denominator_name="sampler settings' time_per_iteration", ) @@ -679,13 +697,13 @@ def _validate_simulation_settings( if output_settings.velocities_write_frequency is not None: _ = settings_validation.divmod_time_and_check( numerator=output_settings.velocities_write_frequency, - denominator=sampler_settings.time_per_iteration, + denominator=simulation_settings.time_per_iteration, numerator_name="output settings' velocity_write_frequency", denominator_name="sampler settings' time_per_iteration", ) _, _ = settings_validation.convert_real_time_analysis_iterations( - simulation_settings=sampler_settings, + simulation_settings=simulation_settings, ) def _validate( From 0605d11049934524bafa5d14d4b12362dce7d29b Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 26 Dec 2025 01:03:39 +0000 Subject: [PATCH 10/20] fix a few things --- .../tests/protocols/openmm_rfe/test_hybrid_top_validation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index d0452eb9a..c4a686a01 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -178,7 +178,7 @@ def test_charge_difference_no_corr(benzene_to_benzoic_mapping): def test_charge_difference_no_solvent(benzene_to_benzoic_mapping): errmsg = "Cannot use eplicit charge correction without solvent" - with pytest.raises(ValueError, errmsg): + with pytest.raises(ValueError, match=errmsg): openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference( benzene_to_benzoic_mapping, "pme", @@ -223,7 +223,7 @@ def test_get_charge_difference(mapping_name, result, request, caplog): mapping = request.getfixturevalue(mapping_name) caplog.set_level(logging.INFO) - ion = r"Na\+" if result == -1 else r"Cl\-" + ion = r"Na+" if result == -1 else r"Cl-" msg = ( f"A charge difference of {result} is observed " "between the end states. This will be addressed by " From 48106a297237fffb9e76ca16d2ed99a3d6834bac Mon Sep 17 00:00:00 2001 From: IAlibay Date: Thu, 25 Dec 2025 23:45:07 -0500 Subject: [PATCH 11/20] fix the remaining tests --- .../protocols/openmm_rfe/equil_rfe_methods.py | 10 +- .../openmm_utils/system_validation.py | 33 +- .../openmm_rfe/test_hybrid_top_protocol.py | 94 ----- .../openmm_rfe/test_hybrid_top_validation.py | 363 +++++++++++++++--- 4 files changed, 337 insertions(+), 163 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 356423922..710f4db8e 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -690,7 +690,7 @@ def _validate_simulation_settings( _ = settings_validation.divmod_time_and_check( numerator=output_settings.positions_write_frequency, denominator=simulation_settings.time_per_iteration, - numerator_name="output settings' position_write_frequency", + numerator_name="output settings' positions_write_frequency", denominator_name="sampler settings' time_per_iteration", ) @@ -698,7 +698,7 @@ def _validate_simulation_settings( _ = settings_validation.divmod_time_and_check( numerator=output_settings.velocities_write_frequency, denominator=simulation_settings.time_per_iteration, - numerator_name="output settings' velocity_write_frequency", + numerator_name="output settings' velocities_write_frequency", denominator_name="sampler settings' time_per_iteration", ) @@ -768,8 +768,10 @@ def _validate( # PR #125 temporarily pin lambda schedule spacing to n_replicas if self.settings.simulation_settings.n_replicas != self.settings.lambda_settings.lambda_windows: errmsg = ( - "Number of replicas in simulation_settings must equal " - "number of lambda windows in lambda_settings." + "Number of replicas in ``simulation_settings``: " + f"{self.settings.simulation_settings.n_replicas} must equal " + "the number of lambda windows in lambda_settings: " + f"{self.settings.lambda_settings.lambda_windows}." ) raise ValueError(errmsg) diff --git a/openfe/protocols/openmm_utils/system_validation.py b/openfe/protocols/openmm_utils/system_validation.py index 0fd3c3518..9d67e108f 100644 --- a/openfe/protocols/openmm_utils/system_validation.py +++ b/openfe/protocols/openmm_utils/system_validation.py @@ -95,23 +95,24 @@ def validate_solvent(state: ChemicalSystem, nonbonded_method: str): `nocutoff`. * If the SolventComponent solvent is not water. """ - solv = [comp for comp in state.values() if isinstance(comp, SolventComponent)] + solv_comps = state.get_components_of_type(SolventComponent) - if len(solv) > 0 and nonbonded_method.lower() == "nocutoff": - errmsg = "nocutoff cannot be used for solvent transformations" - raise ValueError(errmsg) + if len(solv_comps) > 0: + if nonbonded_method.lower() == "nocutoff": + errmsg = "nocutoff cannot be used for solvent transformations" + raise ValueError(errmsg) - if len(solv) == 0 and nonbonded_method.lower() == "pme": - errmsg = "PME cannot be used for vacuum transform" - raise ValueError(errmsg) + if len(solv_comps) > 1: + errmsg = "Multiple SolventComponent found, only one is supported" + raise ValueError(errmsg) - if len(solv) > 1: - errmsg = "Multiple SolventComponent found, only one is supported" - raise ValueError(errmsg) - - if len(solv) > 0 and solv[0].smiles != "O": - errmsg = "Non water solvent is not currently supported" - raise ValueError(errmsg) + if solv_comps[0].smiles != "O": + errmsg = "Non water solvent is not currently supported" + raise ValueError(errmsg) + else: + if nonbonded_method.lower() == "pme": + errmsg = "PME cannot be used for vacuum transform" + raise ValueError(errmsg) def validate_protein(state: ChemicalSystem): @@ -129,9 +130,9 @@ def validate_protein(state: ChemicalSystem): ValueError If there are multiple ProteinComponent in the ChemicalSystem. """ - nprot = sum(1 for comp in state.values() if isinstance(comp, ProteinComponent)) + prot_comps = state.get_components_of_type(ProteinComponent) - if nprot > 1: + if len(prot_comps) > 1: errmsg = "Multiple ProteinComponent found, only one is supported" raise ValueError(errmsg) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py index 148f32fe1..d4818f355 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py @@ -948,63 +948,6 @@ def test_lambda_schedule(windows): assert len(lambdas.lambda_schedule) == windows -def test_hightimestep( - benzene_vacuum_system, - toluene_vacuum_system, - benzene_to_toluene_mapping, - vac_settings, - tmpdir, -): - vac_settings.forcefield_settings.hydrogen_mass = 1.0 - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=vac_settings, - ) - - dag = p.create( - stateA=benzene_vacuum_system, - stateB=toluene_vacuum_system, - mapping=benzene_to_toluene_mapping, - ) - dag_unit = list(dag.protocol_units)[0] - - errmsg = "too large for hydrogen mass" - with tmpdir.as_cwd(): - with pytest.raises(ValueError, match=errmsg): - dag_unit.run(dry=True) - - -def test_element_change_warning(atom_mapping_basic_test_files): - # check a mapping with element change gets rejected early - l1 = atom_mapping_basic_test_files["2-methylnaphthalene"] - l2 = atom_mapping_basic_test_files["2-naftanol"] - - # We use the 'old' lomap defaults because the - # basic test files inputs we use aren't fully aligned - mapper = setup.LomapAtomMapper( - time=20, threed=True, max3d=1000.0, element_change=True, seed="", shift=True - ) - - mapping = next(mapper.suggest_mappings(l1, l2)) - - sys1 = openfe.ChemicalSystem( - {"ligand": l1, "solvent": openfe.SolventComponent()}, - ) - sys2 = openfe.ChemicalSystem( - {"ligand": l2, "solvent": openfe.SolventComponent()}, - ) - - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), - ) - with pytest.warns(UserWarning, match="Element change"): - _ = p.create( - stateA=sys1, - stateB=sys2, - mapping=mapping, - ) - - def test_ligand_overlap_warning( benzene_vacuum_system, toluene_vacuum_system, benzene_to_toluene_mapping, vac_settings, tmpdir ): @@ -2023,40 +1966,3 @@ def test_dry_run_vacuum_write_frequency( assert reporter.velocity_interval == velocities_write_frequency.m else: assert reporter.velocity_interval == 0 - - -@pytest.mark.parametrize( - "positions_write_frequency,velocities_write_frequency", - [ - [100.1 * unit.picosecond, 100 * unit.picosecond], - [100 * unit.picosecond, 100.1 * unit.picosecond], - ], -) -def test_pos_write_frequency_not_divisible( - benzene_vacuum_system, - toluene_vacuum_system, - benzene_to_toluene_mapping, - positions_write_frequency, - velocities_write_frequency, - tmpdir, - vac_settings, -): - vac_settings.output_settings.positions_write_frequency = positions_write_frequency - vac_settings.output_settings.velocities_write_frequency = velocities_write_frequency - - protocol = openmm_rfe.RelativeHybridTopologyProtocol( - settings=vac_settings, - ) - - # create DAG from protocol and take first (and only) work unit from within - dag = protocol.create( - stateA=benzene_vacuum_system, - stateB=toluene_vacuum_system, - mapping=benzene_to_toluene_mapping, - ) - dag_unit = list(dag.protocol_units)[0] - - with tmpdir.as_cwd(): - errmsg = "The output settings' " - with pytest.raises(ValueError, match=errmsg): - dag_unit.run(dry=True)["debug"]["sampler"] diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index c4a686a01..984c4d6bb 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -18,7 +18,7 @@ from kartograf.atom_aligner import align_mol_shape from numpy.testing import assert_allclose from openff.toolkit import Molecule -from openff.units import unit +from openff.units import unit as offunit from openff.units.openmm import ensure_quantity, from_openmm, to_openmm from openmm import CustomNonbondedForce, MonteCarloBarostat, NonbondedForce, XmlSerializer, app from openmm import unit as omm_unit @@ -138,6 +138,169 @@ def test_validate_mapping_alchem_not_in(state, benzene_to_toluene_mapping): ) +def test_vaccuum_PME_error( + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + solv_settings +): + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=solv_settings) + + errmsg = "PME cannot be used for vacuum transform" + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + ) + + +def test_solvent_nocutoff_error( + benzene_system, + toluene_system, + benzene_to_toluene_mapping, + vac_settings, +): + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) + + errmsg = "nocutoff cannot be used for solvent transformation" + + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_system, + stateB=toluene_system, + mapping=benzene_to_toluene_mapping, + ) + + +def test_nonwater_solvent_error( + benzene_modifications, + benzene_to_toluene_mapping, + solv_settings, +): + solvent = openfe.SolventComponent(smiles='C') + stateA = openfe.ChemicalSystem( + { + 'ligand': benzene_modifications['benzene'], + 'solvent': solvent, + } + ) + + stateB = openfe.ChemicalSystem( + { + 'ligand': benzene_modifications['toluene'], + 'solvent': solvent + } + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=solv_settings) + + errmsg = "Non water solvent is not currently supported" + + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=stateA, + stateB=stateB, + mapping=benzene_to_toluene_mapping, + ) + + +def test_too_many_solv_comps_error( + benzene_modifications, + benzene_to_toluene_mapping, + solv_settings, +): + stateA = openfe.ChemicalSystem( + { + 'ligand': benzene_modifications['benzene'], + 'solvent!': openfe.SolventComponent(neutralize=True), + 'solvent2': openfe.SolventComponent(neutralize=False), + } + ) + + stateB = openfe.ChemicalSystem( + { + 'ligand': benzene_modifications['toluene'], + 'solvent!': openfe.SolventComponent(neutralize=True), + 'solvent2': openfe.SolventComponent(neutralize=False), + } + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=solv_settings) + + errmsg = "Multiple SolventComponent found, only one is supported" + + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=stateA, + stateB=stateB, + mapping=benzene_to_toluene_mapping, + ) + + +def test_bad_solv_settings( + benzene_system, + toluene_system, + benzene_to_toluene_mapping, + solv_settings, +): + """ + Test a case where the solvent settings would be wrong. + Not doing every cases since those are covered under + ``test_openmmutils.py``. + """ + solv_settings.solvation_settings.solvent_padding = 1.2 * offunit.nanometer + solv_settings.solvation_settings.number_of_solvent_molecules = 20 + + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=solv_settings) + + errmsg = "Only one of solvent_padding, number_of_solvent_molecules," + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_system, + stateB=toluene_system, + mapping=benzene_to_toluene_mapping + ) + + +def test_too_many_prot_comps_error( + benzene_modifications, + benzene_to_toluene_mapping, + T4_protein_component, + eg5_protein, + solv_settings, +): + + stateA = openfe.ChemicalSystem( + { + 'ligand': benzene_modifications['benzene'], + 'solvent': openfe.SolventComponent(), + 'protein1': T4_protein_component, + 'protein2': eg5_protein, + } + ) + + stateB = openfe.ChemicalSystem( + { + 'ligand': benzene_modifications['toluene'], + 'solvent': openfe.SolventComponent(), + 'protein1': T4_protein_component, + 'protein2': eg5_protein, + } + ) + + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=solv_settings) + + errmsg = "Multiple ProteinComponent found, only one is supported" + + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=stateA, + stateB=stateB, + mapping=benzene_to_toluene_mapping, + ) + + def test_element_change_warning(atom_mapping_basic_test_files): # check a mapping with element change gets rejected early l1 = atom_mapping_basic_test_files["2-methylnaphthalene"] @@ -248,81 +411,183 @@ def test_hightimestep( toluene_vacuum_system, benzene_to_toluene_mapping, vac_settings, - tmpdir, ): vac_settings.forcefield_settings.hydrogen_mass = 1.0 - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=vac_settings, - ) - - dag = p.create( - stateA=benzene_vacuum_system, - stateB=toluene_vacuum_system, - mapping=benzene_to_toluene_mapping, - ) - dag_unit = list(dag.protocol_units)[0] + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) errmsg = "too large for hydrogen mass" - with tmpdir.as_cwd(): - with pytest.raises(ValueError, match=errmsg): - dag_unit.run(dry=True) + + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + extends=None + ) -def test_n_replicas_not_n_windows( +def test_time_per_iteration_divmod( benzene_vacuum_system, toluene_vacuum_system, benzene_to_toluene_mapping, vac_settings, - tmpdir, ): - # For PR #125 we pin such that the number of lambda windows - # equals the numbers of replicas used - TODO: remove limitation - # default lambda windows is 11 - vac_settings.simulation_settings.n_replicas = 13 + vac_settings.simulation_settings.time_per_iteration = 10 * offunit.ps + vac_settings.integrator_settings.timestep = 4 * offunit.ps - errmsg = "Number of replicas 13 does not equal the number of lambda windows 11" + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) - with tmpdir.as_cwd(): - with pytest.raises(ValueError, match=errmsg): - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=vac_settings, - ) - dag = p.create( - stateA=benzene_vacuum_system, - stateB=toluene_vacuum_system, - mapping=benzene_to_toluene_mapping, - ) - dag_unit = list(dag.protocol_units)[0] - dag_unit.run(dry=True) + errmsg = "does not evenly divide by the timestep" + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + extends=None + ) -def test_vaccuum_PME_error( - benzene_vacuum_system, benzene_modifications, benzene_to_toluene_mapping + +@pytest.mark.parametrize( + "attribute", ["equilibration_length", "production_length"] +) +def test_simsteps_not_timestep_divisible( + attribute, + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + vac_settings, ): - # state B doesn't have a solvent component (i.e. its vacuum) - stateB = openfe.ChemicalSystem({"ligand": benzene_modifications["toluene"]}) + setattr(vac_settings.simulation_settings, attribute, 102 * offunit.fs) + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) - p = openmm_rfe.RelativeHybridTopologyProtocol( - settings=openmm_rfe.RelativeHybridTopologyProtocol.default_settings(), + errmsg = "Simulation time not divisible by timestep" + + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + extends=None + ) + + +@pytest.mark.parametrize( + "attribute", ["equilibration_length", "production_length"] +) +def test_simsteps_not_mcstep_divisible( + attribute, + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + vac_settings, +): + setattr(vac_settings.simulation_settings, attribute, 102 * offunit.ps) + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) + + errmsg = ( + "should contain a number of steps divisible by the number of " + "integrator timesteps" ) - errmsg = "PME cannot be used for vacuum transform" + with pytest.raises(ValueError, match=errmsg): - _ = p.create( + p.validate( stateA=benzene_vacuum_system, - stateB=stateB, + stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, + extends=None ) -def test_get_alchemical_waters_no_waters( - benzene_solvent_openmm_system, +def test_checkpoint_interval_not_divisible_time_per_iter( + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + vac_settings, ): - system, topology, positions = benzene_solvent_openmm_system + vac_settings.output_settings.checkpoint_interval = 4 * offunit.ps + vac_settings.simulation_settings.time_per_iteration = 2.5 * offunit.ps + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) - errmsg = "There are no waters" + errmsg = "does not evenly divide by the amount of time per state MCMC" with pytest.raises(ValueError, match=errmsg): - topologyhelpers.get_alchemical_waters( - topology, positions, charge_difference=1, distance_cutoff=3.0 * unit.nanometer + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + extends=None + ) + + +@pytest.mark.parametrize( + "attribute", + ["positions_write_frequency", "velocities_write_frequency"] +) +def test_pos_vel_write_frequency_not_divisible( + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + attribute, + vac_settings, +): + setattr(vac_settings.output_settings, attribute, 100.1 * offunit.picosecond) + + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) + + errmsg = f"The output settings' {attribute}" + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + extends=None + ) + + +@pytest.mark.parametrize( + "attribute", + ["real_time_analysis_interval", "real_time_analysis_interval"] +) +def test_pos_vel_write_frequency_not_divisible( + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + attribute, + vac_settings, +): + setattr(vac_settings.simulation_settings, attribute, 100.1 * offunit.picosecond) + + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) + + errmsg = f"The {attribute}" + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + extends=None + ) + +def test_n_replicas_not_n_windows( + benzene_vacuum_system, + toluene_vacuum_system, + benzene_to_toluene_mapping, + vac_settings, + tmpdir, +): + # For PR #125 we pin such that the number of lambda windows + # equals the numbers of replicas used - TODO: remove limitation + vac_settings.simulation_settings.n_replicas = 13 + p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) + + errmsg = "Number of replicas in ``simulation_settings``:" + + with pytest.raises(ValueError, match=errmsg): + p.validate( + stateA=benzene_vacuum_system, + stateB=toluene_vacuum_system, + mapping=benzene_to_toluene_mapping, + extends=None ) From 5af66e81688603d997282176fd4d11c73c50e454 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Thu, 25 Dec 2025 23:50:36 -0500 Subject: [PATCH 12/20] cleanup imports --- .../protocols/openmm_rfe/equil_rfe_methods.py | 1 - .../openmm_rfe/test_hybrid_top_validation.py | 31 +------------------ 2 files changed, 1 insertion(+), 31 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 710f4db8e..189fffe79 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -53,7 +53,6 @@ from openff.units import Quantity, unit from openff.units.openmm import ensure_quantity, from_openmm, to_openmm from openmmtools import multistate -from rdkit import Chem from openfe.due import Doi, due from openfe.protocols.openmm_utils.omm_settings import ( diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 984c4d6bb..0f7db50d6 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -1,42 +1,13 @@ # This code is part of OpenFE and is licensed under the MIT license. # For details, see https://github.com/OpenFreeEnergy/openfe import logging -import copy -import json -import sys -import xml.etree.ElementTree as ET -from importlib import resources -from math import sqrt -from pathlib import Path -from unittest import mock - -import gufe -import mdtraj as mdt -import numpy as np + import pytest -from kartograf import KartografAtomMapper -from kartograf.atom_aligner import align_mol_shape -from numpy.testing import assert_allclose -from openff.toolkit import Molecule from openff.units import unit as offunit -from openff.units.openmm import ensure_quantity, from_openmm, to_openmm -from openmm import CustomNonbondedForce, MonteCarloBarostat, NonbondedForce, XmlSerializer, app -from openmm import unit as omm_unit -from openmmforcefields.generators import SMIRNOFFTemplateGenerator -from openmmtools.multistate.multistatesampler import MultiStateSampler -from rdkit import Chem -from rdkit.Geometry import Point3D import openfe from openfe import setup from openfe.protocols import openmm_rfe -from openfe.protocols.openmm_rfe._rfe_utils import topologyhelpers -from openfe.protocols.openmm_utils import omm_compute, system_creation -from openfe.protocols.openmm_utils.charge_generation import ( - HAS_ESPALOMA_CHARGE, - HAS_NAGL, - HAS_OPENEYE, -) @pytest.fixture() From 792996e9fd5492537f06427a910d818047c3ede4 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 26 Dec 2025 00:30:39 -0500 Subject: [PATCH 13/20] Add news item --- news/validate-rfe.rst | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 news/validate-rfe.rst diff --git a/news/validate-rfe.rst b/news/validate-rfe.rst new file mode 100644 index 000000000..d7036e8d4 --- /dev/null +++ b/news/validate-rfe.rst @@ -0,0 +1,26 @@ +**Added:** + +* The `validate` method for the RelativeHybridTopologyProtocol has been + implemented. This means that settings and system validation can mostly + be done prior to Protocol execution by calling + `RelativeHybridTopologyProtocol.validate(stateA, stateB, mapping)`. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* From 7d179981e86b8d579951f10314ca14107df94f96 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Fri, 26 Dec 2025 20:25:30 -0500 Subject: [PATCH 14/20] fix redefine --- openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 0f7db50d6..349c81d06 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -521,7 +521,7 @@ def test_pos_vel_write_frequency_not_divisible( "attribute", ["real_time_analysis_interval", "real_time_analysis_interval"] ) -def test_pos_vel_write_frequency_not_divisible( +def test_real_time_analysis_not_divisible( benzene_vacuum_system, toluene_vacuum_system, benzene_to_toluene_mapping, From d1bd736414491f41f45ffb1341fca2d0dd36c86f Mon Sep 17 00:00:00 2001 From: IAlibay Date: Sat, 27 Dec 2025 17:23:43 -0500 Subject: [PATCH 15/20] Add charge validation for smcs when dealing with ismorphic molecules --- .../protocols/openmm_rfe/equil_rfe_methods.py | 55 ++++++++++++++++ .../openmm_utils/system_validation.py | 18 +++--- .../openmm_rfe/test_hybrid_top_validation.py | 63 +++++++++++++++++++ 3 files changed, 127 insertions(+), 9 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 189fffe79..6995f82dd 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -560,6 +560,58 @@ def _validate_mapping( logger.warning(wmsg) warnings.warn(wmsg) + @staticmethod + def _validate_smcs( + stateA: ChemicalSystem, + stateB: ChemicalSystem, + ) -> None: + """ + Validates the SmallMoleculeComponents. + + Parameters + ---------- + stateA : ChemicalSystem + The chemical system of end state A. + stateB : ChemicalSystem + The chemical system of end state B. + + Raises + ------ + ValueError + * If there are isomorphic SmallMoleculeComponents with + different charges. + """ + smcs_A = stateA.get_components_of_type(SmallMoleculeComponent) + smcs_B = stateB.get_components_of_type(SmallMoleculeComponent) + smcs_all = list(set(smcs_A).union(set(smcs_B))) + offmols = [m.to_openff() for m in smcs_all] + + def _equal_charges(moli, molj): + # Base case, both molecules don't have charges + if (moli.partial_charges is None) & (molj.partial_charges is None): + return True + # If either is None but not the other + if (moli.partial_charges is None) ^ (molj.partial_charges is None): + return False + # Check if the charges are close to each other + return np.allclose(moli.partial_charges, molj.partial_charges) + + clashes = [] + + for i, moli in enumerate(offmols): + for molj in offmols: + if moli.is_isomorphic_with(molj): + if not _equal_charges(moli, molj): + clashes.append(smcs_all[i]) + + if len(clashes) > 0: + errmsg = ( + "Found SmallMoleculeComponents are are isomorphic " + "but with different charges, this is not currently allowed. " + f"Affected components: {clashes}" + ) + raise ValueError(errmsg) + @staticmethod def _validate_charge_difference( mapping: LigandAtomMapping, @@ -726,6 +778,9 @@ def _validate( alchem_comps = system_validation.get_alchemical_components(stateA, stateB) self._validate_mapping(mapping, alchem_comps) + # Validate the small molecule components + self._validate_smcs(stateA, stateB) + # Validate solvent component nonbond = self.settings.forcefield_settings.nonbonded_method system_validation.validate_solvent(stateA, nonbond) diff --git a/openfe/protocols/openmm_utils/system_validation.py b/openfe/protocols/openmm_utils/system_validation.py index 9d67e108f..750f5f565 100644 --- a/openfe/protocols/openmm_utils/system_validation.py +++ b/openfe/protocols/openmm_utils/system_validation.py @@ -162,24 +162,24 @@ def get_components(state: ChemicalSystem) -> ParseCompRet: small_mols : list[SmallMoleculeComponent] """ - def _get_single_comps(comp_list, comptype): - ret_comps = [comp for comp in comp_list if isinstance(comp, comptype)] - if ret_comps: + def _get_single_comps(state, comptype): + comps = state.get_components_of_type(comptype) + + if len(ret_comps) > 0: return ret_comps[0] else: return None solvent_comp: Optional[SolventComponent] = _get_single_comps( - list(state.values()), SolventComponent + state, + SolventComponent ) protein_comp: Optional[ProteinComponent] = _get_single_comps( - list(state.values()), ProteinComponent + state, + ProteinComponent ) - small_mols = [] - for comp in state.components.values(): - if isinstance(comp, SmallMoleculeComponent): - small_mols.append(comp) + small_mols = state.get_components_of_type(SmallMoleculeComponent) return solvent_comp, protein_comp, small_mols diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 349c81d06..57092ce3c 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -126,6 +126,69 @@ def test_vaccuum_PME_error( ) +@pytest.mark.parametrize('charge', [None, 'gasteiger']) +def test_smcs_same_charge_passes( + charge, + benzene_modifications +): + benzene = benzene_modifications['benzene'] + if charge is None: + smc = benzene + else: + offmol = benzene.to_openff() + offmol.assign_partial_charges(partial_charge_method='gasteiger') + smc = openfe.SmallMoleculeComponent.from_openff(offmol) + + # Just pass the same thing twice + state = openfe.ChemicalSystem({'l': smc}) + openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs(state, state) + + +def test_smcs_different_charges_none_not_none( + benzene_modifications +): + # smcA has no charges + smcA = benzene_modifications['benzene'] + + # smcB has charges + offmol = smcA.to_openff() + offmol.assign_partial_charges(partial_charge_method='gasteiger') + smcB = openfe.SmallMoleculeComponent.from_openff(offmol) + + stateA = openfe.ChemicalSystem({'l': smcA}) + stateB = openfe.ChemicalSystem({'l': smcB}) + + errmsg = "isomorphic but with different charges" + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs( + stateA, stateB + ) + + +def test_smcs_different_charges_all( + benzene_modifications +): + # For this test, we will assign both A and B to both states + # It wouldn't happen in real life, but it tests that within a state + # you can pick up isomorphic molecules with different charges + # create an offmol with gasteiger charges + offmol = benzene_modifications['benzene'].to_openff() + offmol.assign_partial_charges(partial_charge_method='gasteiger') + smcA = openfe.SmallMoleculeComponent.from_openff(offmol) + + # now alter the offmol charges, scaling by 0.1 + offmol.partial_charges *= 0.1 + smcB = openfe.SmallMoleculeComponent.from_openff(offmol) + + state = openfe.ChemicalSystem({'l1': smcA, 'l2': smcB}) + + errmsg = "isomorphic but with different charges" + with pytest.raises(ValueError, match=errmsg): + openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs( + state, state + ) + + def test_solvent_nocutoff_error( benzene_system, toluene_system, From 5270887e58793269d96036281a37d69201c023c6 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Mon, 29 Dec 2025 17:38:05 -0500 Subject: [PATCH 16/20] Fix validation --- openfe/protocols/openmm_utils/system_validation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openfe/protocols/openmm_utils/system_validation.py b/openfe/protocols/openmm_utils/system_validation.py index 750f5f565..7cacaa1f1 100644 --- a/openfe/protocols/openmm_utils/system_validation.py +++ b/openfe/protocols/openmm_utils/system_validation.py @@ -165,8 +165,8 @@ def get_components(state: ChemicalSystem) -> ParseCompRet: def _get_single_comps(state, comptype): comps = state.get_components_of_type(comptype) - if len(ret_comps) > 0: - return ret_comps[0] + if len(comps) > 0: + return comps[0] else: return None From 15c753de60f8536106d738286e1482622d09edcc Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 3 Jan 2026 00:33:24 +0000 Subject: [PATCH 17/20] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../protocols/openmm_rfe/equil_rfe_methods.py | 20 +- .../openmm_utils/system_validation.py | 12 +- .../openmm_rfe/test_hybrid_top_validation.py | 211 ++++++++---------- 3 files changed, 104 insertions(+), 139 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 6995f82dd..346d0aca4 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -536,9 +536,13 @@ def _validate_mapping( # check that the mapping components are in the alchemical components for m in mapping: if m.componentA not in alchemical_components["stateA"]: - raise ValueError(f"Mapping componentA {m.componentA} not in alchemical components of stateA") + raise ValueError( + f"Mapping componentA {m.componentA} not in alchemical components of stateA" + ) if m.componentB not in alchemical_components["stateB"]: - raise ValueError(f"Mapping componentB {m.componentB} not in alchemical components of stateB") + raise ValueError( + f"Mapping componentB {m.componentB} not in alchemical components of stateB" + ) # TODO: remove - this is now the default behaviour? # Check for element changes in mappings @@ -678,10 +682,7 @@ def _validate_charge_difference( ) raise ValueError(errmsg) - ion = { - -1: solvent_component.positive_ion, - 1: solvent_component.negative_ion - }[difference] + ion = {-1: solvent_component.positive_ion, 1: solvent_component.negative_ion}[difference] wmsg = ( f"A charge difference of {difference} is observed " @@ -712,7 +713,7 @@ def _validate_simulation_settings( Raises ------ ValueError - * If the + * If the """ steps_per_iteration = settings_validation.convert_steps_per_iteration( @@ -820,7 +821,10 @@ def _validate( # Validate alchemical settings # PR #125 temporarily pin lambda schedule spacing to n_replicas - if self.settings.simulation_settings.n_replicas != self.settings.lambda_settings.lambda_windows: + if ( + self.settings.simulation_settings.n_replicas + != self.settings.lambda_settings.lambda_windows + ): errmsg = ( "Number of replicas in ``simulation_settings``: " f"{self.settings.simulation_settings.n_replicas} must equal " diff --git a/openfe/protocols/openmm_utils/system_validation.py b/openfe/protocols/openmm_utils/system_validation.py index 7cacaa1f1..3e8ed5c50 100644 --- a/openfe/protocols/openmm_utils/system_validation.py +++ b/openfe/protocols/openmm_utils/system_validation.py @@ -170,15 +170,9 @@ def _get_single_comps(state, comptype): else: return None - solvent_comp: Optional[SolventComponent] = _get_single_comps( - state, - SolventComponent - ) - - protein_comp: Optional[ProteinComponent] = _get_single_comps( - state, - ProteinComponent - ) + solvent_comp: Optional[SolventComponent] = _get_single_comps(state, SolventComponent) + + protein_comp: Optional[ProteinComponent] = _get_single_comps(state, ProteinComponent) small_mols = state.get_components_of_type(SmallMoleculeComponent) diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 57092ce3c..08597df10 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -33,49 +33,54 @@ def test_invalid_protocol_repeats(): settings.protocol_repeats = -1 -@pytest.mark.parametrize('state', ['A', 'B']) +@pytest.mark.parametrize("state", ["A", "B"]) def test_endstate_two_alchemcomp_stateA(state, benzene_modifications): - first_state = openfe.ChemicalSystem({ - 'ligandA': benzene_modifications['benzene'], - 'ligandB': benzene_modifications['toluene'], - 'solvent': openfe.SolventComponent(), - }) - other_state = openfe.ChemicalSystem({ - 'ligandC': benzene_modifications['phenol'], - 'solvent': openfe.SolventComponent(), - }) - - if state == 'A': + first_state = openfe.ChemicalSystem( + { + "ligandA": benzene_modifications["benzene"], + "ligandB": benzene_modifications["toluene"], + "solvent": openfe.SolventComponent(), + } + ) + other_state = openfe.ChemicalSystem( + { + "ligandC": benzene_modifications["phenol"], + "solvent": openfe.SolventComponent(), + } + ) + + if state == "A": args = (first_state, other_state) else: args = (other_state, first_state) with pytest.raises(ValueError, match="Only one alchemical component"): - openmm_rfe.RelativeHybridTopologyProtocol._validate_endstates( - *args - ) + openmm_rfe.RelativeHybridTopologyProtocol._validate_endstates(*args) -@pytest.mark.parametrize('state', ['A', 'B']) + +@pytest.mark.parametrize("state", ["A", "B"]) def test_endstates_not_smc(state, benzene_modifications): - first_state = openfe.ChemicalSystem({ - 'ligand': benzene_modifications['benzene'], - 'foo': openfe.SolventComponent(), - }) - other_state = openfe.ChemicalSystem({ - 'ligand': benzene_modifications['benzene'], - 'foo': benzene_modifications['toluene'], - }) - - if state == 'A': + first_state = openfe.ChemicalSystem( + { + "ligand": benzene_modifications["benzene"], + "foo": openfe.SolventComponent(), + } + ) + other_state = openfe.ChemicalSystem( + { + "ligand": benzene_modifications["benzene"], + "foo": benzene_modifications["toluene"], + } + ) + + if state == "A": args = (first_state, other_state) else: args = (other_state, first_state) errmsg = "only SmallMoleculeComponents transformations" with pytest.raises(ValueError, match=errmsg): - openmm_rfe.RelativeHybridTopologyProtocol._validate_endstates( - *args - ) + openmm_rfe.RelativeHybridTopologyProtocol._validate_endstates(*args) def test_validate_mapping_none_mapping(): @@ -88,12 +93,11 @@ def test_validate_mapping_multi_mapping(benzene_to_toluene_mapping): errmsg = "A single LigandAtomMapping is expected" with pytest.raises(ValueError, match=errmsg): openmm_rfe.RelativeHybridTopologyProtocol._validate_mapping( - [benzene_to_toluene_mapping] * 2, - None + [benzene_to_toluene_mapping] * 2, None ) -@pytest.mark.parametrize('state', ['A', 'B']) +@pytest.mark.parametrize("state", ["A", "B"]) def test_validate_mapping_alchem_not_in(state, benzene_to_toluene_mapping): errmsg = f"not in alchemical components of state{state}" @@ -110,10 +114,7 @@ def test_validate_mapping_alchem_not_in(state, benzene_to_toluene_mapping): def test_vaccuum_PME_error( - benzene_vacuum_system, - toluene_vacuum_system, - benzene_to_toluene_mapping, - solv_settings + benzene_vacuum_system, toluene_vacuum_system, benzene_to_toluene_mapping, solv_settings ): p = openmm_rfe.RelativeHybridTopologyProtocol(settings=solv_settings) @@ -126,67 +127,56 @@ def test_vaccuum_PME_error( ) -@pytest.mark.parametrize('charge', [None, 'gasteiger']) -def test_smcs_same_charge_passes( - charge, - benzene_modifications -): - benzene = benzene_modifications['benzene'] +@pytest.mark.parametrize("charge", [None, "gasteiger"]) +def test_smcs_same_charge_passes(charge, benzene_modifications): + benzene = benzene_modifications["benzene"] if charge is None: smc = benzene else: offmol = benzene.to_openff() - offmol.assign_partial_charges(partial_charge_method='gasteiger') + offmol.assign_partial_charges(partial_charge_method="gasteiger") smc = openfe.SmallMoleculeComponent.from_openff(offmol) # Just pass the same thing twice - state = openfe.ChemicalSystem({'l': smc}) + state = openfe.ChemicalSystem({"l": smc}) openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs(state, state) -def test_smcs_different_charges_none_not_none( - benzene_modifications -): +def test_smcs_different_charges_none_not_none(benzene_modifications): # smcA has no charges - smcA = benzene_modifications['benzene'] + smcA = benzene_modifications["benzene"] # smcB has charges offmol = smcA.to_openff() - offmol.assign_partial_charges(partial_charge_method='gasteiger') + offmol.assign_partial_charges(partial_charge_method="gasteiger") smcB = openfe.SmallMoleculeComponent.from_openff(offmol) - stateA = openfe.ChemicalSystem({'l': smcA}) - stateB = openfe.ChemicalSystem({'l': smcB}) + stateA = openfe.ChemicalSystem({"l": smcA}) + stateB = openfe.ChemicalSystem({"l": smcB}) errmsg = "isomorphic but with different charges" with pytest.raises(ValueError, match=errmsg): - openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs( - stateA, stateB - ) + openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs(stateA, stateB) -def test_smcs_different_charges_all( - benzene_modifications -): +def test_smcs_different_charges_all(benzene_modifications): # For this test, we will assign both A and B to both states # It wouldn't happen in real life, but it tests that within a state # you can pick up isomorphic molecules with different charges # create an offmol with gasteiger charges - offmol = benzene_modifications['benzene'].to_openff() - offmol.assign_partial_charges(partial_charge_method='gasteiger') + offmol = benzene_modifications["benzene"].to_openff() + offmol.assign_partial_charges(partial_charge_method="gasteiger") smcA = openfe.SmallMoleculeComponent.from_openff(offmol) # now alter the offmol charges, scaling by 0.1 offmol.partial_charges *= 0.1 smcB = openfe.SmallMoleculeComponent.from_openff(offmol) - state = openfe.ChemicalSystem({'l1': smcA, 'l2': smcB}) + state = openfe.ChemicalSystem({"l1": smcA, "l2": smcB}) errmsg = "isomorphic but with different charges" with pytest.raises(ValueError, match=errmsg): - openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs( - state, state - ) + openmm_rfe.RelativeHybridTopologyProtocol._validate_smcs(state, state) def test_solvent_nocutoff_error( @@ -212,20 +202,15 @@ def test_nonwater_solvent_error( benzene_to_toluene_mapping, solv_settings, ): - solvent = openfe.SolventComponent(smiles='C') + solvent = openfe.SolventComponent(smiles="C") stateA = openfe.ChemicalSystem( { - 'ligand': benzene_modifications['benzene'], - 'solvent': solvent, + "ligand": benzene_modifications["benzene"], + "solvent": solvent, } ) - stateB = openfe.ChemicalSystem( - { - 'ligand': benzene_modifications['toluene'], - 'solvent': solvent - } - ) + stateB = openfe.ChemicalSystem({"ligand": benzene_modifications["toluene"], "solvent": solvent}) p = openmm_rfe.RelativeHybridTopologyProtocol(settings=solv_settings) @@ -246,17 +231,17 @@ def test_too_many_solv_comps_error( ): stateA = openfe.ChemicalSystem( { - 'ligand': benzene_modifications['benzene'], - 'solvent!': openfe.SolventComponent(neutralize=True), - 'solvent2': openfe.SolventComponent(neutralize=False), + "ligand": benzene_modifications["benzene"], + "solvent!": openfe.SolventComponent(neutralize=True), + "solvent2": openfe.SolventComponent(neutralize=False), } ) stateB = openfe.ChemicalSystem( { - 'ligand': benzene_modifications['toluene'], - 'solvent!': openfe.SolventComponent(neutralize=True), - 'solvent2': openfe.SolventComponent(neutralize=False), + "ligand": benzene_modifications["toluene"], + "solvent!": openfe.SolventComponent(neutralize=True), + "solvent2": openfe.SolventComponent(neutralize=False), } ) @@ -290,11 +275,7 @@ def test_bad_solv_settings( errmsg = "Only one of solvent_padding, number_of_solvent_molecules," with pytest.raises(ValueError, match=errmsg): - p.validate( - stateA=benzene_system, - stateB=toluene_system, - mapping=benzene_to_toluene_mapping - ) + p.validate(stateA=benzene_system, stateB=toluene_system, mapping=benzene_to_toluene_mapping) def test_too_many_prot_comps_error( @@ -304,22 +285,21 @@ def test_too_many_prot_comps_error( eg5_protein, solv_settings, ): - stateA = openfe.ChemicalSystem( { - 'ligand': benzene_modifications['benzene'], - 'solvent': openfe.SolventComponent(), - 'protein1': T4_protein_component, - 'protein2': eg5_protein, + "ligand": benzene_modifications["benzene"], + "solvent": openfe.SolventComponent(), + "protein1": T4_protein_component, + "protein2": eg5_protein, } ) stateB = openfe.ChemicalSystem( { - 'ligand': benzene_modifications['toluene'], - 'solvent': openfe.SolventComponent(), - 'protein1': T4_protein_component, - 'protein2': eg5_protein, + "ligand": benzene_modifications["toluene"], + "solvent": openfe.SolventComponent(), + "protein1": T4_protein_component, + "protein2": eg5_protein, } ) @@ -419,19 +399,16 @@ def test_greater_than_one_charge_difference_error(aniline_to_benzoic_mapping): def test_get_charge_difference(mapping_name, result, request, caplog): mapping = request.getfixturevalue(mapping_name) caplog.set_level(logging.INFO) - + ion = r"Na+" if result == -1 else r"Cl-" msg = ( f"A charge difference of {result} is observed " "between the end states. This will be addressed by " f"transforming a water into a {ion} ion" ) - + openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference( - mapping, - "pme", - True, - openfe.SolventComponent() + mapping, "pme", True, openfe.SolventComponent() ) if result != 0: @@ -457,7 +434,7 @@ def test_hightimestep( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) @@ -479,13 +456,11 @@ def test_time_per_iteration_divmod( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) -@pytest.mark.parametrize( - "attribute", ["equilibration_length", "production_length"] -) +@pytest.mark.parametrize("attribute", ["equilibration_length", "production_length"]) def test_simsteps_not_timestep_divisible( attribute, benzene_vacuum_system, @@ -503,13 +478,11 @@ def test_simsteps_not_timestep_divisible( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) -@pytest.mark.parametrize( - "attribute", ["equilibration_length", "production_length"] -) +@pytest.mark.parametrize("attribute", ["equilibration_length", "production_length"]) def test_simsteps_not_mcstep_divisible( attribute, benzene_vacuum_system, @@ -520,17 +493,14 @@ def test_simsteps_not_mcstep_divisible( setattr(vac_settings.simulation_settings, attribute, 102 * offunit.ps) p = openmm_rfe.RelativeHybridTopologyProtocol(settings=vac_settings) - errmsg = ( - "should contain a number of steps divisible by the number of " - "integrator timesteps" - ) + errmsg = "should contain a number of steps divisible by the number of integrator timesteps" with pytest.raises(ValueError, match=errmsg): p.validate( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) @@ -551,14 +521,11 @@ def test_checkpoint_interval_not_divisible_time_per_iter( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) -@pytest.mark.parametrize( - "attribute", - ["positions_write_frequency", "velocities_write_frequency"] -) +@pytest.mark.parametrize("attribute", ["positions_write_frequency", "velocities_write_frequency"]) def test_pos_vel_write_frequency_not_divisible( benzene_vacuum_system, toluene_vacuum_system, @@ -576,13 +543,12 @@ def test_pos_vel_write_frequency_not_divisible( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) @pytest.mark.parametrize( - "attribute", - ["real_time_analysis_interval", "real_time_analysis_interval"] + "attribute", ["real_time_analysis_interval", "real_time_analysis_interval"] ) def test_real_time_analysis_not_divisible( benzene_vacuum_system, @@ -601,9 +567,10 @@ def test_real_time_analysis_not_divisible( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) + def test_n_replicas_not_n_windows( benzene_vacuum_system, toluene_vacuum_system, @@ -623,5 +590,5 @@ def test_n_replicas_not_n_windows( stateA=benzene_vacuum_system, stateB=toluene_vacuum_system, mapping=benzene_to_toluene_mapping, - extends=None + extends=None, ) From 6a8220ee3d53c5cc88c18421a0d6cb6c5509c72e Mon Sep 17 00:00:00 2001 From: Irfan Alibay Date: Tue, 6 Jan 2026 12:17:03 +0000 Subject: [PATCH 18/20] Apply suggestions from code review Co-authored-by: Hannah Baumann <43765638+hannahbaumann@users.noreply.github.com> Co-authored-by: Josh Horton --- openfe/protocols/openmm_rfe/equil_rfe_methods.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 346d0aca4..bdbe17329 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -610,7 +610,7 @@ def _equal_charges(moli, molj): if len(clashes) > 0: errmsg = ( - "Found SmallMoleculeComponents are are isomorphic " + "Found SmallMoleculeComponents that are isomorphic " "but with different charges, this is not currently allowed. " f"Affected components: {clashes}" ) @@ -644,6 +644,7 @@ def _validate_charge_difference( nonbonded method is not PME. * If the absolute charge difference is greater than one and an explicit charge correction is attempted. + * If an explicit charge correction is attempted and there is no solvent present. UserWarning * If there is any charge difference. """ @@ -664,7 +665,7 @@ def _validate_charge_difference( return if solvent_component is None: - errmsg = "Cannot use eplicit charge correction without solvent" + errmsg = "Cannot use explicit charge correction without solvent" raise ValueError(errmsg) # We implicitly check earlier that we have to have pme for a solvated @@ -775,7 +776,7 @@ def _validate( # Validate the end states self._validate_endstates(stateA, stateB) - # Valildate the mapping + # Validate the mapping alchem_comps = system_validation.get_alchemical_components(stateA, stateB) self._validate_mapping(mapping, alchem_comps) From ccb8c4fb03ec8b139b2878fc41ca39612d464414 Mon Sep 17 00:00:00 2001 From: IAlibay Date: Tue, 6 Jan 2026 07:23:16 -0500 Subject: [PATCH 19/20] address missing docs --- openfe/protocols/openmm_rfe/equil_rfe_methods.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index bdbe17329..080c41258 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -714,7 +714,13 @@ def _validate_simulation_settings( Raises ------ ValueError - * If the + * If any of of the simulation control settings (e.g. + ``equilibration_length`` or ``production_length``) + are not divisible by the timestep or the number of + steps per iteration. + * If the output frequency for position, velocity, or + online analysis are not divisible by the time per + multistate iteration. """ steps_per_iteration = settings_validation.convert_steps_per_iteration( From 03503f4d6bc7b7f9d75f08850a0ef341e73c38df Mon Sep 17 00:00:00 2001 From: IAlibay Date: Tue, 6 Jan 2026 07:34:25 -0500 Subject: [PATCH 20/20] fix test --- openfe/protocols/openmm_rfe/equil_rfe_methods.py | 15 +++++++-------- .../openmm_rfe/test_hybrid_top_validation.py | 2 +- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index 080c41258..524f44ff4 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -535,14 +535,13 @@ def _validate_mapping( # check that the mapping components are in the alchemical components for m in mapping: - if m.componentA not in alchemical_components["stateA"]: - raise ValueError( - f"Mapping componentA {m.componentA} not in alchemical components of stateA" - ) - if m.componentB not in alchemical_components["stateB"]: - raise ValueError( - f"Mapping componentB {m.componentB} not in alchemical components of stateB" - ) + for state in ["A", "B"]: + comp = getattr(m, f"component{state}") + if comp not in alchemical_components[f"state{state}"]: + raise ValueError( + f"Mapping component{state} {comp} not " + f"in alchemical components of state{state}" + ) # TODO: remove - this is now the default behaviour? # Check for element changes in mappings diff --git a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py index 08597df10..35ef57da0 100644 --- a/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py +++ b/openfe/tests/protocols/openmm_rfe/test_hybrid_top_validation.py @@ -353,7 +353,7 @@ def test_charge_difference_no_corr(benzene_to_benzoic_mapping): def test_charge_difference_no_solvent(benzene_to_benzoic_mapping): - errmsg = "Cannot use eplicit charge correction without solvent" + errmsg = "Cannot use explicit charge correction without solvent" with pytest.raises(ValueError, match=errmsg): openmm_rfe.RelativeHybridTopologyProtocol._validate_charge_difference(