diff --git a/news/abfe_host_smc.rst b/news/abfe_host_smc.rst
new file mode 100644
index 000000000..b7c31a871
--- /dev/null
+++ b/news/abfe_host_smc.rst
@@ -0,0 +1,25 @@
+**Added:**
+
+* Add support for SmallMoleculeComponents to act
+  as host in the ``AbsoluteBindingProtocol``
+  (`PR 2020 <https://github.com/OpenFreeEnergy/openfe/pull/2020>`_).
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>
diff --git a/news/abfe_user_dfn_restraints.rst b/news/abfe_user_dfn_restraints.rst
new file mode 100644
index 000000000..65c1d3b75
--- /dev/null
+++ b/news/abfe_user_dfn_restraints.rst
@@ -0,0 +1,24 @@
+**Added:**
+
+* Added support for user-defined Boresch restraints in the
+  ABFE Protocol (`PR #2019 <https://github.com/OpenFreeEnergy/openfe/pull/2019>`_).
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>
diff --git a/news/net_charge_abfe.rst b/news/net_charge_abfe.rst
new file mode 100644
index 000000000..7adbffaaa
--- /dev/null
+++ b/news/net_charge_abfe.rst
@@ -0,0 +1,25 @@
+**Added:**
+
+* Add support for transformations involving net
+  charges in the ABFE Protocol
+  (`PR #2020 <https://github.com/OpenFreeEnergy/openfe/pull/2020>`_).
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>
diff --git a/src/openfe/protocols/openmm_afe/abfe_units.py b/src/openfe/protocols/openmm_afe/abfe_units.py
index 343cbd2e8..d3a475e2e 100644
--- a/src/openfe/protocols/openmm_afe/abfe_units.py
+++ b/src/openfe/protocols/openmm_afe/abfe_units.py
@@ -9,6 +9,7 @@
 import logging
 import pathlib
 from collections.abc import Iterable
+from copy import deepcopy
 
 import MDAnalysis as mda
 import numpy as np
@@ -17,22 +18,30 @@
     SolventComponent,
 )
 from gufe.components import Component, SolvatedPDBComponent
+from MDAnalysis.lib.distances import calc_bonds
 from openff.units import Quantity
+from openff.units import unit as offunit
 from openff.units.openmm import to_openmm
-from openmm import System
+from openmm import (
+    HarmonicBondForce,
+    NonbondedForce,
+    System,
+)
 from openmm import unit as ommunit
 from openmm.app import Topology as omm_topology
 from openmmtools.states import ThermodynamicState
 from rdkit import Chem
 
 from openfe.protocols.openmm_afe.equil_afe_settings import (
-    BoreschRestraintSettings,
+    ABFEBoreschRestraintSettings,
     SettingsBaseModel,
 )
 from openfe.protocols.openmm_utils import system_validation
 from openfe.protocols.restraint_utils import geometry
 from openfe.protocols.restraint_utils.geometry.boresch import BoreschRestraintGeometry
+from openfe.protocols.restraint_utils.geometry.utils import FindHostAtoms, get_central_atom_idx
 from openfe.protocols.restraint_utils.openmm import omm_restraints
+from openfe.protocols.restraint_utils.openmm.omm_forces import add_force_in_separate_group
 from openfe.protocols.restraint_utils.openmm.omm_restraints import BoreschRestraint
 
 from .base_afe_units import (
@@ -44,6 +53,290 @@
 logger = logging.getLogger(__name__)
 
 
+def _get_mda_universe(
+    topology: omm_topology,
+    positions: ommunit.Quantity | None,
+    trajectory: pathlib.Path | None,
+) -> mda.Universe:
+    """
+    Helper method to get a Universe from an openmm Topology,
+    and either an input trajectory or a set of positions.
+
+    Parameters
+    ----------
+    topology : openmm.app.Topology
+      An OpenMM Topology that defines the System.
+    positions: openmm.unit.Quantity | None
+      The System's current positions.
+      Used if a trajectory file is None or is not a file.
+    trajectory: pathlib.Path | None
+      A Path to a trajectory file to read positions from.
+
+    Returns
+    -------
+    mda.Universe
+      An MDAnalysis Universe of the System.
+    """
+    from MDAnalysis.coordinates.memory import MemoryReader
+
+    # If the trajectory file doesn't exist, then we use positions
+    if trajectory is not None and trajectory.is_file():
+        return mda.Universe(
+            topology,
+            trajectory,
+            topology_format="OPENMMTOPOLOGY",
+        )
+    else:
+        if positions is None:
+            raise ValueError("No positions to create the Universe with")
+
+        # Positions is an openmm Quantity in nm we need
+        # to convert to angstroms
+        return mda.Universe(
+            topology,
+            np.array(positions._value) * 10,
+            topology_format="OPENMMTOPOLOGY",
+            trajectory_format=MemoryReader,
+        )
+
+
+def _get_idxs_from_residxs(
+    topology: omm_topology,
+    residxs: Iterable[int],
+) -> list[int]:
+    """
+    Helper method to get the a list of atom indices which belong to a list
+    of residues.
+
+    Parameters
+    ----------
+    topology : openmm.app.Topology
+      An OpenMM Topology that defines the System.
+    residxs : Iterable[int]
+      A list of residue numbers who's atoms we should get atom indices.
+
+    Returns
+    -------
+    atom_ids : list[int]
+      A list of atom indices.
+
+    TODO
+    ----
+    * Check how this works when we deal with virtual sites.
+    """
+    atom_ids = []
+
+    for r in topology.residues():
+        if r.index in residxs:
+            atom_ids.extend([at.index for at in r.atoms()])
+
+    return atom_ids
+
+
+def _get_minimum_image_distance(box_dimensions: npt.NDArray) -> Quantity:
+    """
+    Get the minimum image distance using the minimum perpendicular width
+    of the triclinic vectors.
+
+    Parameters
+    ----------
+    box_dimensions : npt.NDArray
+      The box dimensions as obtained from an MDAnalysis Universe.
+
+    Returns
+    -------
+    openff.units.Quantity
+      The minimum perpendicular width in units of Angstrom.
+
+    Acknowledgements
+    ----------------
+    Originally contributed by Bendict Tan (@aqemia-benedict-tan).
+    """
+    from MDAnalysis.lib import mdamath
+
+    box_vectors = mdamath.triclinic_vectors(box_dimensions)
+
+    # Calculate the volume based on the scalar triple product
+    volume = mdamath.stp(box_vectors[0], box_vectors[1], box_vectors[2])
+
+    # Now calculate the perpendicular widths using perp_width_i = Volume / Area_of_face_i
+    # Where Area_of_face_i is |box_vectors_{i+1} × box_vectors_{i+2}|.
+    areas = np.cross(box_vectors[[1, 2, 0]], box_vectors[[2, 0, 1]])
+    perp_widths = volume / np.linalg.norm(areas, axis=1)
+
+    return perp_widths.min() * offunit.angstrom
+
+
+def _find_most_common_ions(
+    openmm_topology: omm_topology,
+    openmm_system: System,
+    target_charge: int,
+) -> list[int] | None:
+    """
+    Get the most common ions of a given net charge in a system.
+
+    Parameters
+    ----------
+    openmm_topology : openmm.app.Topology
+      The Topology of the OpenMM System.
+    openmm_system : openmm.System
+      The OpenMM System.
+    target_charge : int
+      The charge the ion should have.
+
+    Returns
+    -------
+    list[int] | None
+      If present, the list of indices matching the most common ion.
+
+    Notes
+    -----
+    This is similar to what is done in ``_get_ion_parameters`` in
+    :mod:`openfe.protocols.openmm_rfe._rfe_utils.topologyhelpers`.
+    """
+    from collections import Counter, defaultdict
+
+    nbf = [i for i in openmm_system.getForces() if isinstance(i, NonbondedForce)][0]
+
+    ion_counts: Counter = Counter()
+    ion_atom_indices: dict[str, list[int]] = defaultdict(list)
+
+    for residue in openmm_topology.residues():
+        atoms = list(residue.atoms())
+
+        # We are only interested in single atom counterions
+        if len(atoms) != 1:
+            continue
+
+        charge, _, _ = nbf.getParticleParameters(atoms[0].index)
+        charge_val = charge.value_in_unit(ommunit.elementary_charge)
+
+        if np.isclose(charge_val, target_charge, atol=0.01):
+            ion_counts[residue.name] += 1
+            ion_atom_indices[residue.name].append(atoms[0].index)
+
+    if ion_counts:
+        best_resname = ion_counts.most_common(1)[0][0]
+        return ion_atom_indices[best_resname]
+    else:
+        return None
+
+
+class ABFESetupUnitMixin:
+    """
+    Mixin for common class methods between Units
+    """
+
+    def _get_alchemical_ions(
+        self,
+        alchemical_components: dict[str, list[Component]],
+        comp_resids: dict[Component, npt.NDArray],
+        openmm_topology: omm_topology,
+        openmm_system: System,
+        positions: ommunit.Quantity,
+        settings: dict[str, SettingsBaseModel],
+        dry: bool,
+    ) -> list[int] | None:
+        """
+        Find a suitable alchemical ion for a net charge transformation.
+
+        Parameters
+        ----------
+        alchemical_components: dict[str, list[Component]]
+          A dictionary with a list of alchemical components
+          in both state A and B.
+        comp_resids: dict[Component, npt.NDArray]
+          A dictionary keyed by each Component in the System
+          which contains arrays with the residue indices that is contained
+          by that Component.
+        openmm_topology : openmm.app.Topology
+          The OpenMM Topology of the system.
+        openmm_system : openmm.System
+          The OpenMM System to work on.
+        positions : openmm.unit.Quantity
+          The positions of the system.
+        settings : dict[str, SettingsBaseModel]
+          A dictionary of settings that defines how to find and set
+          the restraint.
+        dry: bool
+          ``True`` if we are dry-running.
+
+        Returns
+        -------
+        list[int]
+          The indices of the alchemical ions.
+
+        """
+        total_charge = alchemical_components["stateA"][0].total_charge
+
+        # Don't add an alchemical ion if we have zero net charge
+        # or we didn't request it.
+        if total_charge == 0 or not settings["alchemical_settings"].explicit_charge_correction:
+            return None
+
+        # TODO: For now, let's stick with a single -1/+1 case, but we should expand to more
+        if abs(total_charge) > 1:
+            errmsg = "Cannot handle net charge correction on charges greater than one"
+            raise ValueError(errmsg)
+
+        # Get the indices of the most common ion type that can act as a counterion
+        ion_indices = _find_most_common_ions(openmm_topology, openmm_system, -total_charge)
+
+        if ion_indices is None:
+            errmsg = "No suitable ions could be found to act as counterion in the system"
+            raise ValueError(errmsg)
+
+        univ = _get_mda_universe(
+            openmm_topology,
+            positions,
+            self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename,  # type: ignore[attr-defined]
+        )
+
+        # get an atomgroup of the possible alchemical ions
+        ions_atomgroup = univ.atoms[ion_indices]
+
+        # get the alchemical atoms
+        residxs = np.concatenate([comp_resids[key] for key in alchemical_components["stateA"]])
+        alchem_idxs = _get_idxs_from_residxs(topology=openmm_topology, residxs=residxs)
+        alchem_atomgroup = univ.atoms[alchem_idxs]
+
+        # Get the maximum distance we can use to find ions
+        univ.trajectory[-1]  # use the box dimensions from the last frame
+        box = univ.dimensions
+
+        if box is None or np.all(np.isinf(box)) or np.any(box[:3] <= 0.0):
+            # If it's not a dry simulation then error out
+            if not dry:
+                errmsg = f"Invalid box for co-alchemical ion search: {box}"
+                raise ValueError(errmsg)
+
+            # For a dry execution, just assign a reasonably high value
+            max_search_distance = 1.5 * offunit.nanometer
+        else:
+            # Set the max search distance to half the smallest perpendicular width
+            # with a 1 Angstrom padding
+            max_search_distance = (_get_minimum_image_distance(box) * 0.5) - 1 * offunit.angstrom  # type: ignore[operator]
+
+        # Re-using a utility from the restraints utilities
+        # TODO: rename this class!
+        atom_finder = FindHostAtoms(
+            host_atoms=ions_atomgroup,
+            guest_atoms=alchem_atomgroup,
+            min_search_distance=settings["alchemical_settings"].alchemical_ion_min_distance,
+            max_search_distance=max_search_distance,
+        )
+
+        # only run on the final frame
+        atom_finder.run(frames=[-1])
+
+        if len(atom_finder.results.host_idxs) == 0:
+            errmsg = "No suitable alchemical ion was found"
+            raise ValueError(errmsg)
+
+        # Just use the first one that comes back ok
+        return [atom_finder.results.host_idxs[0]]
+
+
 class ComplexComponentsMixin:
     def _get_components(self):
         """
@@ -123,7 +416,9 @@ def _get_settings(self) -> dict[str, SettingsBaseModel]:
         return settings
 
 
-class ABFEComplexSetupUnit(ComplexComponentsMixin, ComplexSettingsMixin, BaseAbsoluteSetupUnit):
+class ABFEComplexSetupUnit(
+    ABFESetupUnitMixin, ComplexComponentsMixin, ComplexSettingsMixin, BaseAbsoluteSetupUnit
+):
     """
     Setup unit for the complex phase of absolute binding free energy
     transformations.
@@ -131,86 +426,6 @@ class ABFEComplexSetupUnit(ComplexComponentsMixin, ComplexSettingsMixin, BaseAbs
 
     simtype = "complex"
 
-    @staticmethod
-    def _get_mda_universe(
-        topology: omm_topology,
-        positions: ommunit.Quantity | None,
-        trajectory: pathlib.Path | None,
-    ) -> mda.Universe:
-        """
-        Helper method to get a Universe from an openmm Topology,
-        and either an input trajectory or a set of positions.
-
-        Parameters
-        ----------
-        topology : openmm.app.Topology
-          An OpenMM Topology that defines the System.
-        positions: openmm.unit.Quantity | None
-          The System's current positions.
-          Used if a trajectory file is None or is not a file.
-        trajectory: pathlib.Path | None
-          A Path to a trajectory file to read positions from.
-
-        Returns
-        -------
-        mda.Universe
-          An MDAnalysis Universe of the System.
-        """
-        from MDAnalysis.coordinates.memory import MemoryReader
-
-        # If the trajectory file doesn't exist, then we use positions
-        if trajectory is not None and trajectory.is_file():
-            return mda.Universe(
-                topology,
-                trajectory,
-                topology_format="OPENMMTOPOLOGY",
-            )
-        else:
-            if positions is None:
-                raise ValueError("No positions to create the Universe with")
-
-            # Positions is an openmm Quantity in nm we need
-            # to convert to angstroms
-            return mda.Universe(
-                topology,
-                np.array(positions._value) * 10,
-                topology_format="OPENMMTOPOLOGY",
-                trajectory_format=MemoryReader,
-            )
-
-    @staticmethod
-    def _get_idxs_from_residxs(
-        topology: omm_topology,
-        residxs: Iterable[int],
-    ) -> list[int]:
-        """
-        Helper method to get the a list of atom indices which belong to a list
-        of residues.
-
-        Parameters
-        ----------
-        topology : openmm.app.Topology
-          An OpenMM Topology that defines the System.
-        residxs : Iterable[int]
-          A list of residue numbers who's atoms we should get atom indices.
-
-        Returns
-        -------
-        atom_ids : list[int]
-          A list of atom indices.
-
-        TODO
-        ----
-        * Check how this works when we deal with virtual sites.
-        """
-        atom_ids = []
-
-        for r in topology.residues():
-            if r.index in residxs:
-                atom_ids.extend([at.index for at in r.atoms()])
-
-        return atom_ids
-
     @staticmethod
     def _get_boresch_restraint(
         universe: mda.Universe,
@@ -218,7 +433,7 @@ def _get_boresch_restraint(
         guest_atom_ids: list[int],
         host_atom_ids: list[int],
         temperature: Quantity,
-        settings: BoreschRestraintSettings,
+        settings: ABFEBoreschRestraintSettings,
     ) -> tuple[BoreschRestraintGeometry, BoreschRestraint]:
         """
         Get a Boresch-like restraint Geometry and OpenMM restraint force
@@ -236,7 +451,7 @@ def _get_boresch_restraint(
           A list of atom indices defining the host molecules in the universe.
         temperature : openff.units.Quantity
           The temperature of the simulation where the restraint will be added.
-        settings : BoreschRestraintSettings
+        settings : ABFEBoreschRestraintSettings
           Settings on how the Boresch-like restraint should be defined.
 
         Returns
@@ -254,6 +469,12 @@ def _get_boresch_restraint(
             guest_rdmol=guest_rdmol,
             guest_idxs=guest_atom_ids,
             host_idxs=host_atom_ids,
+            guest_restraint_atoms_idxs=list(settings.guest_restraint_ids)
+            if settings.guest_restraint_ids is not None
+            else None,
+            host_restraint_atoms_idxs=list(settings.host_restraint_ids)
+            if settings.host_restraint_ids is not None
+            else None,
             host_selection=settings.host_selection,
             anchor_finding_strategy=settings.anchor_finding_strategy,
             dssp_filter=settings.dssp_filter,
@@ -275,6 +496,7 @@ def _add_restraints(
         alchem_comps: dict[str, list[Component]],
         comp_resids: dict[Component, npt.NDArray],
         settings: dict[str, SettingsBaseModel],
+        alchemical_ions: list[int] | None,
     ) -> tuple[
         Quantity,
         System,
@@ -306,6 +528,8 @@ def _add_restraints(
         settings : dict[str, SettingsBaseModel]
           A dictionary of settings that defines how to find and set
           the restraint.
+        alchemical_ions : list[int] | None
+          The alchemical ion indices, if they exist.
 
         Returns
         -------
@@ -315,6 +539,10 @@ def _add_restraints(
           A copy of the System with the restraint added.
         rest_geom : geometry.HostGuestRestraintGeometry
           The restraint Geometry object.
+
+        TODO
+        ----
+        Add a restraint between the alchemical ion and the guest molecule?
         """
         if self.verbose:
             self.logger.info("Generating restraints")
@@ -334,7 +562,7 @@ def _add_restraints(
         residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]])
 
         # get the alchemicical atom ids
-        guest_atom_ids = self._get_idxs_from_residxs(topology, residxs)
+        guest_atom_ids = _get_idxs_from_residxs(topology, residxs)
 
         # Now get the host idxs
         # We assume this is everything but the alchemical component
@@ -343,19 +571,19 @@ def _add_restraints(
         exclude_comps = [alchem_comps["stateA"]] + solv_comps
         residxs = np.concatenate([v for i, v in comp_resids.items() if i not in exclude_comps])
 
-        host_atom_ids = self._get_idxs_from_residxs(topology, residxs)
+        host_atom_ids = _get_idxs_from_residxs(topology, residxs)
 
         # Finally create an MDAnalysis Universe
         # We try to pass the equilibration production file path through
         # In some cases (debugging / dry runs) this won't be available
         # so we'll default to using input positions.
-        univ = self._get_mda_universe(
+        univ = _get_mda_universe(
             topology,
             positions,
             self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename,
         )
 
-        if isinstance(settings["restraint_settings"], BoreschRestraintSettings):
+        if isinstance(settings["restraint_settings"], ABFEBoreschRestraintSettings):
             rest_geom, restraint = self._get_boresch_restraint(
                 univ,
                 guest_rdmol,
@@ -383,7 +611,7 @@ def _add_restraints(
         restraint.add_force(
             thermodynamic_state,
             rest_geom,
-            controlling_parameter_name="lambda_restraints",
+            controlling_parameter_name="lambda_restraints_A",
         )
         # Get the standard state correction as a unit.Quantity
         correction = restraint.get_standard_state_correction(
@@ -491,7 +719,9 @@ def _get_settings(self) -> dict[str, SettingsBaseModel]:
         return settings
 
 
-class ABFESolventSetupUnit(SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteSetupUnit):
+class ABFESolventSetupUnit(
+    ABFESetupUnitMixin, SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteSetupUnit
+):
     """
     Setup unit for the solvent phase of absolute binding free energy
     transformations.
@@ -499,6 +729,120 @@ class ABFESolventSetupUnit(SolventComponentsMixin, SolventSettingsMixin, BaseAbs
 
     simtype = "solvent"
 
+    def _add_restraints(
+        self,
+        system: System,
+        topology: omm_topology,
+        positions: ommunit.Quantity,
+        alchem_comps: dict[str, list[Component]],
+        comp_resids: dict[Component, npt.NDArray],
+        settings: dict[str, SettingsBaseModel],
+        alchemical_ions: list[int] | None,
+    ) -> tuple[
+        Quantity | None,
+        System,
+        geometry.HostGuestRestraintGeometry | None,
+    ]:
+        """
+        Find and add restraints to the OpenMM System.
+
+        Notes
+        -----
+        Currently, only Boresch-like restraints are supported.
+
+        Parameters
+        ----------
+        system : openmm.System
+          The System to add the restraint to.
+        topology : openmm.app.Topology
+          An OpenMM Topology that defines the System.
+        positions: openmm.unit.Quantity
+          The System's current positions.
+          Used if a trajectory file isn't found.
+        alchem_comps: dict[str, list[Component]]
+          A dictionary with a list of alchemical components
+          in both state A and B.
+        comp_resids: dict[Component, npt.NDArray]
+          A dictionary keyed by each Component in the System
+          which contains arrays with the residue indices that is contained
+          by that Component.
+        settings : dict[str, SettingsBaseModel]
+          A dictionary of settings that defines how to find and set
+          the restraint.
+        alchemical_ions : list[int] | None
+          The alchemical ion indices, if they exist.
+
+        Returns
+        -------
+        correction : openff.units.Quantity | None
+          The standard state correction for the restraint.
+        system : openmm.System
+          A copy of the System with the restraint added.
+        rest_geom : geometry.HostGuestRestraintGeometry | None
+          The restraint Geometry object.
+
+        TODO
+        ----
+        Expand to support restraining multiple ions.
+        """
+        if alchemical_ions is None:
+            return None, system, None
+
+        if len(alchemical_ions) > 1:
+            errmsg = "Currently cannot handle more than one alchemical ion"
+            raise ValueError(errmsg)
+
+        restrained_system = deepcopy(system)
+
+        if self.verbose:
+            self.logger.info("Generating restraints for alchemical ions")
+
+        universe = _get_mda_universe(
+            topology,
+            positions,
+            self.shared_basepath / settings["equil_output_settings"].production_trajectory_filename,
+        )
+
+        # alchemical ion atom atomgroup
+        alchem_ion_ag = universe.atoms[alchemical_ions]
+
+        # get the alchemical ligand atoms
+        ligand_rdmol = alchem_comps["stateA"][0].to_rdkit()
+        residxs = np.concatenate([comp_resids[key] for key in alchem_comps["stateA"]])
+        ligand_alchem_idxs = _get_idxs_from_residxs(topology=topology, residxs=residxs)
+        ligand_central_atom = ligand_alchem_idxs[get_central_atom_idx(ligand_rdmol)]
+        ligand_central_atom_ag = universe.atoms[ligand_central_atom]
+
+        # Get the ligand-ion distance based on the final frame
+        universe.trajectory[-1]
+
+        distance = float(
+            calc_bonds(
+                alchem_ion_ag.atoms[0].position,
+                ligand_central_atom_ag.position,
+                box=universe.dimensions,
+            )
+        )
+
+        spring_constant = to_openmm(
+            settings["alchemical_settings"].alchemical_ion_solvent_spring_constant
+        )
+
+        force = HarmonicBondForce()
+
+        force.addBond(
+            ligand_central_atom,
+            alchemical_ions[0],
+            distance * ommunit.angstrom,
+            spring_constant,
+        )
+
+        force.setName("ion_restraint")
+        # TODO: Temporarily disabled this for testing
+        # add_force_in_separate_group(restrained_system, force)
+
+        return None, restrained_system, None
+
 
 class ABFESolventSimUnit(
     SolventComponentsMixin, SolventSettingsMixin, BaseAbsoluteMultiStateSimulationUnit
diff --git a/src/openfe/protocols/openmm_afe/base_afe_units.py b/src/openfe/protocols/openmm_afe/base_afe_units.py
index 676ba4c51..5a5379d47 100644
--- a/src/openfe/protocols/openmm_afe/base_afe_units.py
+++ b/src/openfe/protocols/openmm_afe/base_afe_units.py
@@ -46,7 +46,6 @@
 from openmmtools.alchemy import (
     AbsoluteAlchemicalFactory,
     AlchemicalRegion,
-    AlchemicalState,
 )
 from openmmtools.states import (
     GlobalParameterState,
@@ -85,6 +84,10 @@
     make_vec3_box,
     serialize,
 )
+from openfe.protocols.openmm_utils.states import (
+    DualRegionAlchemicalState,
+    SingleRegionAlchemicalState,
+)
 from openfe.protocols.restraint_utils import geometry
 from openfe.protocols.restraint_utils.openmm import omm_restraints
 from openfe.utils import log_system_probe, without_oechem_backend
@@ -534,6 +537,25 @@ def _get_omm_objects(
 
         return topology, system, positions, comp_resids
 
+    def _get_alchemical_ions(
+        self,
+        alchemical_components: dict[str, list[Component]],
+        comp_resids: dict[Component, npt.NDArray],
+        openmm_topology: app.Topology,
+        openmm_system: openmm.System,
+        positions: openmm.unit.Quantity,
+        settings: dict[str, SettingsBaseModel],
+        dry: bool,
+    ) -> list[int] | None:
+        """
+        Placeholder method for finding alchemical ions.
+
+        Returns
+        -------
+        None
+        """
+        return None
+
     def _add_restraints(
         self,
         system: openmm.System,
@@ -542,9 +564,10 @@ def _add_restraints(
         alchem_comps: dict[str, list[Component]],
         comp_resids: dict[Component, npt.NDArray],
         settings: dict[str, SettingsBaseModel],
+        alchemical_ions: list[int] | None,
     ) -> tuple[
         Quantity | None,
-        openmm.System | None,
+        openmm.System,
         geometry.BaseRestraintGeometry | None,
     ]:
         """
@@ -557,9 +580,10 @@ def _get_alchemical_system(
         topology: app.Topology,
         system: openmm.System,
         comp_resids: dict[Component, npt.NDArray],
-        alchem_comps: dict[str, list[Component]],
+        alchemical_components: dict[str, list[Component]],
+        alchemical_ions: list[int] | None,
         alchemical_settings: AlchemicalSettings,
-    ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, list[int]]:
+    ) -> tuple[AbsoluteAlchemicalFactory, openmm.System, dict[str, list[int]]]:
         """
         Get an alchemically modified system and its associated factory
 
@@ -571,8 +595,10 @@ def _get_alchemical_system(
           System to alchemically modify.
         comp_resids : dict[str, npt.NDArray]
           A dictionary of residues for each component in the System.
-        alchem_comps : dict[str, list[Component]]
+        alchemical_components : dict[str, list[Component]]
           A dictionary of alchemical components for each end state.
+        alchemical_ions : list[int] | None
+          List of indices for any alchemical ions, there are any.
         alchemical_settings : AlchemicalSettings
           Settings controlling how the alchemical system is built.
 
@@ -582,29 +608,41 @@ def _get_alchemical_system(
           Factory for creating an alchemically modified system.
         alchemical_system : openmm.System
           Alchemically modified system
-        alchemical_indices : list[int]
-          A list of atom indices for the alchemically modified
-          species in the system.
+        alchemical_indices : dict[str, list[int]]
+          A dictionary containing a list of atom indices
+          for each independent alchemically modified species in the system.
 
         TODO
         ----
         * Add support for all alchemical factory options
         """
-        alchemical_indices = self._get_alchemical_indices(topology, comp_resids, alchem_comps)
-
-        alchemical_region = AlchemicalRegion(
-            alchemical_atoms=alchemical_indices,
-            softcore_alpha=alchemical_settings.softcore_alpha,
-            annihilate_electrostatics=True,
-            annihilate_sterics=alchemical_settings.annihilate_sterics,
-            softcore_a=alchemical_settings.softcore_a,
-            softcore_b=alchemical_settings.softcore_b,
-            softcore_c=alchemical_settings.softcore_c,
-            softcore_beta=0.0,
-            softcore_d=1.0,
-            softcore_e=1.0,
-            softcore_f=2.0,
-        )
+        # first region is the alchemically changing species
+        alchemical_indices = {
+            "A": self._get_alchemical_indices(topology, comp_resids, alchemical_components)
+        }
+        # second region is any alchemical ions
+        if alchemical_ions is not None:
+            alchemical_indices["B"] = alchemical_ions
+
+        alchemical_regions = []
+
+        for region in alchemical_indices:
+            alchemical_regions.append(
+                AlchemicalRegion(
+                    alchemical_atoms=alchemical_indices[region],
+                    name=region,
+                    softcore_alpha=alchemical_settings.softcore_alpha,
+                    annihilate_electrostatics=True,
+                    annihilate_sterics=alchemical_settings.annihilate_sterics,
+                    softcore_a=alchemical_settings.softcore_a,
+                    softcore_b=alchemical_settings.softcore_b,
+                    softcore_c=alchemical_settings.softcore_c,
+                    softcore_beta=0.0,
+                    softcore_d=1.0,
+                    softcore_e=1.0,
+                    softcore_f=2.0,
+                )
+            )
 
         alchemical_factory = AbsoluteAlchemicalFactory(
             consistent_exceptions=False,
@@ -614,7 +652,10 @@ def _get_alchemical_system(
             disable_alchemical_dispersion_correction=alchemical_settings.disable_alchemical_dispersion_correction,
             split_alchemical_forces=True,
         )
-        alchemical_system = alchemical_factory.create_alchemical_system(system, alchemical_region)
+        alchemical_system = alchemical_factory.create_alchemical_system(
+            reference_system=system,
+            alchemical_regions=alchemical_regions,
+        )
 
         return alchemical_factory, alchemical_system, alchemical_indices
 
@@ -715,6 +756,17 @@ def run(
             omm_system, omm_topology, positions, settings, dry
         )
 
+        # Get alchemical ions, if needed / allowed
+        alchem_ions = self._get_alchemical_ions(
+            alchemical_components=alchem_comps,
+            comp_resids=comp_resids,
+            openmm_topology=omm_topology,
+            openmm_system=omm_system,
+            positions=positions,
+            settings=settings,
+            dry=dry,
+        )
+
         # Add restraints
         # Note: when no restraint is applied, restrained_omm_system == omm_system
         (
@@ -728,6 +780,7 @@ def run(
             alchem_comps,
             comp_resids,
             settings,
+            alchem_ions,
         )
 
         # Get alchemical system
@@ -735,7 +788,8 @@ def run(
             topology=omm_topology,
             system=restrained_omm_system,
             comp_resids=comp_resids,
-            alchem_comps=alchem_comps,
+            alchemical_components=alchem_comps,
+            alchemical_ions=alchem_ions,
             alchemical_settings=settings["alchemical_settings"],
         )
 
@@ -767,6 +821,7 @@ def run(
             "pdb_structure": pdb_structure,
             "selection_indices": selection_indices,
             "box_vectors": from_openmm(box_vectors),
+            "alchemical_indices": alchem_indices,
         }
 
         if standard_state_corr is not None:
@@ -786,7 +841,6 @@ def run(
                 "standard_system": omm_system,
                 "restrained_system": restrained_omm_system,
                 "alchem_system": alchem_system,
-                "alchem_indices": alchem_indices,
                 "alchem_factory": alchem_factory,
                 "debug_positions": positions,
             }
@@ -919,6 +973,7 @@ def _get_states(
         lambdas: dict[str, list[float]],
         solvent_component: BaseSolventComponent | None,
         alchemically_restrained: bool,
+        alchemical_indices: dict[str, list[int]],
     ) -> tuple[list[SamplerState], list[ThermodynamicState]]:
         """
         Get a list of sampler and thermodynmic states from an
@@ -941,6 +996,9 @@ def _get_states(
         alchemically_restrained : bool
           Whether or not the system requires a control parameter
           for any alchemical restraints.
+        alchemical_indices : dict[str, list[int]]
+          Dictionary of alchemical indices for each alchemical
+          region in the system.
 
         Returns
         -------
@@ -950,7 +1008,13 @@ def _get_states(
           A list of ThermodynamicState for each replica in the system.
         """
         # Fetch an alchemical state
-        alchemical_state = AlchemicalState.from_system(alchemical_system)
+        if len(alchemical_indices) == 1:
+            alchemical_state = SingleRegionAlchemicalState.from_system(alchemical_system)
+        elif len(alchemical_indices) == 2:
+            alchemical_state = DualRegionAlchemicalState.from_system(alchemical_system)
+        else:
+            errmsg = "More than two alchemical regions are not supported"
+            raise ValueError(errmsg)
 
         # Set up the system constants
         temperature = thermodynamic_settings.temperature
@@ -962,24 +1026,29 @@ def _get_states(
             constants["pressure"] = ensure_quantity(pressure, "openmm")
 
         # Get the thermodynamic parameter protocol
-        param_protocol = copy.deepcopy(lambdas)
+        # We populate the protocol from lambdas based on the number of regions
+        # we have.
+        def _add_lambdas_to_protocol(protocol, lambdas, region_name):
+            for key in ["lambda_electrostatics", "lambda_sterics"]:
+                protocol[f"{key}_{region_name}"] = copy.deepcopy(lambdas[key])
 
-        # Get the composable states
-        if alchemically_restrained:
-            restraint_state = omm_restraints.RestraintParameterState(lambda_restraints=1.0)
-            composable_states = [alchemical_state, restraint_state]
-        else:
-            composable_states = [alchemical_state]
+        param_protocol: dict[str, list[float]] = {}
+        # Always add for the first region
+        _add_lambdas_to_protocol(param_protocol, lambdas, "A")
+
+        # If we have two regions (i.e. an alchemical ion) add
+        if len(alchemical_indices) == 2:
+            _add_lambdas_to_protocol(param_protocol, lambdas, "B")
 
-            # In this case we also don't have a restraint being controlled
-            # so we drop it from the protocol
-            param_protocol.pop("lambda_restraints", None)
+        # Only the first region (ligand) can be restrained
+        if alchemically_restrained:
+            param_protocol["lambda_restraints_A"] = copy.deepcopy(lambdas["lambda_restraints"])
 
         cmp_states = create_thermodynamic_state_protocol(
             alchemical_system,
             protocol=param_protocol,
             constants=constants,
-            composable_states=composable_states,
+            composable_states=[alchemical_state],
         )
 
         sampler_state = SamplerState(positions=positions)
@@ -1355,6 +1424,7 @@ def run(
         positions: openmm.unit.Quantity,
         box_vectors: Quantity,
         selection_indices: npt.NDArray,
+        alchemical_indices: dict[str, list[int]],
         alchemical_restraints: bool,
         dry: bool = False,
         verbose: bool = True,
@@ -1374,6 +1444,8 @@ def run(
           The box vectors of the System.
         selection_indices : npt.NDArray
           Indices of the System particles to write to file.
+        alchemical_indices : dict[str, list[int]]
+          Dictionary of alchemical indices for each alchemical region.
         alchemical_restraints: bool,
           Whether or not the system has alchemical restraints.
         dry: bool
@@ -1434,6 +1506,7 @@ def run(
             lambdas=lambdas,
             solvent_component=solv_comp,
             alchemically_restrained=alchemical_restraints,
+            alchemical_indices=alchemical_indices,
         )
 
         # Get the integrator
@@ -1528,6 +1601,7 @@ def _execute(
         system = deserialize(setup_results.outputs["system"])
         positions = to_openmm(np.load(setup_results.outputs["positions"]) * offunit.nanometer)
         selection_indices = setup_results.outputs["selection_indices"]
+        alchemical_indices = setup_results.outputs["alchemical_indices"]
         box_vectors = setup_results.outputs["box_vectors"]
 
         if setup_results.outputs["restraint_geometry"] is not None:
@@ -1540,6 +1614,7 @@ def _execute(
             positions=positions,
             box_vectors=box_vectors,
             selection_indices=selection_indices,
+            alchemical_indices=alchemical_indices,
             alchemical_restraints=alchemical_restraints,
             scratch_basepath=ctx.scratch,
             shared_basepath=ctx.shared,
diff --git a/src/openfe/protocols/openmm_afe/equil_afe_settings.py b/src/openfe/protocols/openmm_afe/equil_afe_settings.py
index 0941b6784..dcdbdc002 100644
--- a/src/openfe/protocols/openmm_afe/equil_afe_settings.py
+++ b/src/openfe/protocols/openmm_afe/equil_afe_settings.py
@@ -22,6 +22,10 @@
     SettingsBaseModel,
     ThermoSettings,
 )
+from gufe.settings.typing import (
+    NanometerQuantity,
+)
+from openff.units import unit as offunit
 from pydantic import field_validator
 
 from openfe.protocols.openmm_utils.omm_settings import (
@@ -38,6 +42,7 @@
 from openfe.protocols.restraint_utils.settings import (
     BaseRestraintSettings,
     BoreschRestraintSettings,
+    SpringConstantLinearQuantity,
 )
 
 
@@ -87,6 +92,27 @@ class AlchemicalSettings(SettingsBaseModel):
     """
 
 
+class ABFEAlchemicalSettings(AlchemicalSettings):
+    # Dev note: we make a separate class for ABFEs so that SepTop and AHFE can
+    # use the parent class without having ABFE-specific net charge settings
+    explicit_charge_correction: bool = True
+    """
+    Whether or not to use explicit charge correction using
+    a co-alchemical ion.
+    """
+    alchemical_ion_min_distance: NanometerQuantity = 1.0 * offunit.nanometer
+    """
+    The minimum distance to search for a co-alchemical ion.
+    """
+    alchemical_ion_solvent_spring_constant: SpringConstantLinearQuantity = (
+        1000.0 * offunit.kilojoule_per_mole / offunit.nm**2
+    )
+    """
+    The spring constant holding the ion away from the alchemical solute
+    in the solvent leg.
+    """
+
+
 class LambdaSettings(SettingsBaseModel):
     """Lambda schedule settings.
 
@@ -213,6 +239,29 @@ def must_be_all(cls, v):
         return v
 
 
+class ABFEBoreschRestraintSettings(BoreschRestraintSettings):
+    host_restraint_ids: tuple[int, int, int] | None = None
+    """
+    The indices of the host component atoms to restrain.
+    The entries define the H0, H1, and H2 atoms in order.
+    If defined, these will override any automatic selection.
+    """
+    guest_restraint_ids: tuple[int, int, int] | None = None
+    """
+    The indices of the guest component atoms to restraint.
+    The entries define the G0, G1, and G2 atoms in order.
+    If defined, these will override any automatic selection.
+    """
+
+    @field_validator("guest_restraint_ids", "host_restraint_ids")
+    def positive_idxs_three_tuple(cls, v):
+        if v is not None:
+            if any([i < 0 for i in v]):
+                errmsg = "``guest_atoms`` and ``host_atoms`` cannot have negative indices."
+                raise ValueError(errmsg)
+        return v
+
+
 # This subclasses from SettingsBaseModel as it has vacuum_forcefield and
 # solvent_forcefield fields, not just a single forcefield_settings field
 class AbsoluteSolvationSettings(SettingsBaseModel):
@@ -360,10 +409,21 @@ def must_be_positive(cls, v):
     """Settings for solvating the system in the complex."""
 
     # Alchemical settings
-    alchemical_settings: AlchemicalSettings
+    alchemical_settings: ABFEAlchemicalSettings
     """
     Alchemical protocol settings.
     """
+
+    @field_validator("alchemical_settings", mode="before")
+    @classmethod
+    def coerce_alchemical_settings(cls, v):
+        """
+        Migration from previous default to the new one
+        """
+        if isinstance(v, AlchemicalSettings) and not isinstance(v, ABFEAlchemicalSettings):
+            return ABFEAlchemicalSettings(**v.model_dump())
+        return v
+
     complex_lambda_settings: LambdaSettings
     """
     Settings for controlling the complex transformation leg
diff --git a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py
index 06dcbaed0..e68dd2699 100644
--- a/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py
+++ b/src/openfe/protocols/openmm_afe/equil_binding_afe_method.py
@@ -43,10 +43,10 @@
 
 from openfe.due import Doi, due
 from openfe.protocols.openmm_afe.equil_afe_settings import (
+    ABFEAlchemicalSettings,
+    ABFEBoreschRestraintSettings,
     ABFEPreEquilOutputSettings,
     AbsoluteBindingSettings,
-    AlchemicalSettings,
-    BoreschRestraintSettings,
     IntegratorSettings,
     LambdaSettings,
     MDSimulationSettings,
@@ -134,7 +134,7 @@ def _default_settings(cls):
                 temperature=298.15 * offunit.kelvin,
                 pressure=1 * offunit.bar,
             ),
-            alchemical_settings=AlchemicalSettings(),
+            alchemical_settings=ABFEAlchemicalSettings(),
             solvent_lambda_settings=LambdaSettings(
                 lambda_elec=[
                     0.0, 0.25, 0.5, 0.75, 1.0,
@@ -174,7 +174,7 @@ def _default_settings(cls):
             engine_settings=OpenMMEngineSettings(),
             solvent_integrator_settings=IntegratorSettings(),
             complex_integrator_settings=IntegratorSettings(),
-            restraint_settings=BoreschRestraintSettings(),
+            restraint_settings=ABFEBoreschRestraintSettings(),
             solvent_equil_simulation_settings=MDSimulationSettings(
                 equilibration_length_nvt=0.1 * offunit.nanosecond,
                 equilibration_length=0.2 * offunit.nanosecond,
@@ -254,6 +254,7 @@ def _adaptive_settings(
     def _validate_endstates(
         stateA: ChemicalSystem,
         stateB: ChemicalSystem,
+        charge_correction: bool,
     ) -> None:
         """
         A binding transformation is defined (in terms of gufe components)
@@ -266,6 +267,8 @@ def _validate_endstates(
           The chemical system of end state A
         stateB : ChemicalSystem
           The chemical system of end state B
+        charge_correction : bool
+          If we are using a charge correction scheme.
 
         Raises
         ------
@@ -275,11 +278,18 @@ def _validate_endstates(
           If stateA has more than one unique Component.
           If the stateA unique Component is not a SmallMoleculeComponent.
           If stateB contains any unique Components.
-          If the alchemical species is charged.
+          If the alchemical species is has greater than one unit charge.
+        UserWarning
+          If the alchemical species has a net charge but ``charge_correction``
+          is ``False``.
         """
         if not (stateA.contains(ProteinComponent) and stateB.contains(ProteinComponent)):
-            errmsg = "No ProteinComponent found"
-            raise ValueError(errmsg)
+            # Check if there is a suitable SmallMoleculeComponent that could
+            # be the host molecule instead.
+            smcs = stateA.get_components_of_type(SmallMoleculeComponent)
+            if all(smc in stateA.component_diff(stateB)[0] for smc in smcs):
+                errmsg = "No suitable host molecule found"
+                raise ValueError(errmsg)
 
         if not (stateA.contains(SolventComponent) and stateB.contains(SolventComponent)):
             errmsg = "No SolventComponent found"
@@ -302,14 +312,25 @@ def _validate_endstates(
             )
             raise ValueError(errmsg)
 
-        # Check that the state A unique isn't charged
+        # Check that the state A unique total charge
         if diff[0][0].total_charge != 0:
-            errmsg = (
-                "Charged alchemical molecules are not currently "
-                "supported for solvation free energies. "
-                f"Molecule total charge: {diff[0][0].total_charge}."
-            )
-            raise ValueError(errmsg)
+            # Error if the total charge is greater than 1
+            if abs(diff[0][0].total_charge) > 1:
+                errmsg = (
+                    "Alchemical molecules with a formal charge of greater than 1 "
+                    "are not currently supported in binding free energy calculations. "
+                    f"Molecule total charge: {diff[0][0].total_charge}."
+                )
+                raise ValueError(errmsg)
+
+            # Warn if we aren't using a charge correction
+            if not charge_correction:
+                wmsg = (
+                    "Ligand has a net charge but no charge correction scheme has been requested. "
+                    "Please note that you will need to apply your own correction scheme "
+                    "to account for finite-size effects."
+                )
+                warnings.warn(wmsg)
 
         # If there are any alchemical Components in state B
         if len(diff[1]) > 0:
@@ -401,7 +422,11 @@ def _validate(
         # Validate the end states & alchemical components
         system_validation.validate_chemical_system(stateA)
         system_validation.validate_chemical_system(stateB)
-        self._validate_endstates(stateA, stateB)
+        self._validate_endstates(
+            stateA,
+            stateB,
+            self.settings.alchemical_settings.explicit_charge_correction,
+        )
 
         # Validate the complex lambda schedule
         self._validate_lambda_schedule(
diff --git a/src/openfe/protocols/openmm_septop/base_units.py b/src/openfe/protocols/openmm_septop/base_units.py
index bf19f8dd5..173499ba0 100644
--- a/src/openfe/protocols/openmm_septop/base_units.py
+++ b/src/openfe/protocols/openmm_septop/base_units.py
@@ -76,6 +76,7 @@
 from openfe.protocols.openmm_utils import omm_compute
 from openfe.protocols.openmm_utils.omm_settings import MultiStateAnalysisSettings, SettingsBaseModel
 from openfe.protocols.openmm_utils.serialization import deserialize
+from openfe.protocols.openmm_utils.states import DualRegionAlchemicalState
 from openfe.utils import log_system_probe, without_oechem_backend
 
 from ..openmm_utils import (
@@ -86,7 +87,6 @@
     system_validation,
 )
 from ..openmm_utils.mdtraj_utils import mdtraj_from_openmm
-from .utils import SepTopParameterState
 
 logger = logging.getLogger(__name__)
 
@@ -879,7 +879,7 @@ def _get_states(
         cmp_states : list[ThermodynamicState]
           A list of ThermodynamicState for each replica in the system.
         """
-        alchemical_state = SepTopParameterState.from_system(alchemical_system)
+        alchemical_state = DualRegionAlchemicalState.from_system(alchemical_system)
 
         # Set up the system constants
         temperature = settings["thermo_settings"].temperature
diff --git a/src/openfe/protocols/openmm_septop/utils.py b/src/openfe/protocols/openmm_septop/utils.py
deleted file mode 100644
index 38d7f1d34..000000000
--- a/src/openfe/protocols/openmm_septop/utils.py
+++ /dev/null
@@ -1,85 +0,0 @@
-from openmmtools import states
-from openmmtools.states import GlobalParameterState
-
-
-class SepTopParameterState(GlobalParameterState):
-    """
-    Composable state to control lambda parameters for two ligands.
-    See :class:`openmmtools.states.GlobalParameterState` for more details.
-    Parameters
-    ----------
-    parameters_name_suffix : Optional[str]
-      If specified, the state will control a modified version of the parameter
-      ``lambda_restraints_{parameters_name_suffix}` instead of just
-      ``lambda_restraints``.
-    lambda_sterics_A : Optional[float]
-      The value for the vdW interactions for ligand A.
-      If defined, must be between 0 and 1.
-    lambda_electrosterics_A : Optional[float]
-      The value for the electrostatics interactions for ligand A.
-      If defined, must be between 0 and 1.
-    lambda_restraints_A : Optional[float]
-      The strength of the restraint for ligand A.
-      If defined, must be between 0 and 1.
-    lambda_bonds_A : Optional[float]
-      The value for modifying bonds for ligand A.
-      If defined, must be between 0 and 1.
-    lambda_angles_A : Optional[float]
-      The value for modifying angles for ligand A.
-      If defined, must be between 0 and 1.
-    lambda_dihedrals_A : Optional[float]
-      The value for modifying dihedrals for ligand A.
-      If defined, must be between 0 and 1.
-    lambda_sterics_B : Optional[float]
-      The value for the vdW interactions for ligand B.
-      If defined, must be between 0 and 1.
-    lambda_electrosterics_B : Optional[float]
-      The value for the electrostatics interactions for ligand B.
-      If defined, must be between 0 and 1.
-    lambda_restraints_B : Optional[float]
-      The strength of the restraint for ligand B.
-      If defined, must be between 0 and 1.
-    lambda_bonds_B : Optional[float]
-      The value for modifying bonds for ligand B.
-      If defined, must be between 0 and 1.
-    lambda_angles_B : Optional[float]
-      The value for modifying angles for ligand B.
-      If defined, must be between 0 and 1.
-    lambda_dihedrals_B : Optional[float]
-      The value for modifying dihedrals for ligand B.
-      If defined, must be between 0 and 1.
-    """
-
-    class _LambdaParameter(states.GlobalParameterState.GlobalParameter):
-        """A global parameter in the interval [0, 1] with standard
-        value 1."""
-
-        def __init__(self, parameter_name):
-            super().__init__(parameter_name, standard_value=1.0, validator=self.lambda_validator)
-
-        @staticmethod
-        def lambda_validator(self, instance, parameter_value):
-            if parameter_value is None:
-                return parameter_value
-            if not (0.0 <= parameter_value <= 1.0):
-                raise ValueError("{} must be between 0 and 1.".format(self.parameter_name))
-            return float(parameter_value)
-
-    # Lambda parameters for ligand A
-    lambda_sterics_A = _LambdaParameter("lambda_sterics_A")
-    lambda_electrostatics_A = _LambdaParameter("lambda_electrostatics_A")
-    lambda_restraints_A = _LambdaParameter("lambda_restraints_A")
-    lambda_bonds_A = _LambdaParameter("lambda_bonds_A")
-    lambda_angles_A = _LambdaParameter("lambda_angles_A")
-    lambda_torsions_A = _LambdaParameter("lambda_torsions_A")
-
-    # Lambda parameters for ligand B
-    lambda_sterics_B = _LambdaParameter("lambda_sterics_B")
-    lambda_electrostatics_B = _LambdaParameter("lambda_electrostatics_B")
-    lambda_restraints_B = _LambdaParameter("lambda_restraints_B")
-    lambda_bonds_B = _LambdaParameter("lambda_bonds_B")
-    lambda_angles_B = _LambdaParameter("lambda_angles_B")
-    lambda_torsions_B = _LambdaParameter("lambda_torsions_B")
-
-    # # Restraints solvent
-    # lambda_restraints = _LambdaParameter('lambda_restraints')
diff --git a/src/openfe/protocols/openmm_utils/states.py b/src/openfe/protocols/openmm_utils/states.py
new file mode 100644
index 000000000..c79bfacd9
--- /dev/null
+++ b/src/openfe/protocols/openmm_utils/states.py
@@ -0,0 +1,128 @@
+# This code is part of OpenFE and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/openfe
+#
+# Acknowledgements:
+# This module derives heavily from the AlchemicalState class
+# in openmmtools.alchemy (https://github.com/choderalab/openmmtools).
+
+from openmmtools.states import GlobalParameterState
+
+
+class _LambdaParameter(GlobalParameterState.GlobalParameter):
+    """
+    A global parameter in the interval [0, 1] with standard value 1.
+    """
+
+    def __init__(self, parameter_name):
+        super().__init__(parameter_name, standard_value=1.0, validator=self.lambda_validator)
+
+    @staticmethod
+    def lambda_validator(self, instance, parameter_value):
+        if parameter_value is None:
+            return parameter_value
+        if not (0.0 <= parameter_value <= 1.0):
+            raise ValueError("{} must be between 0 and 1.".format(self.parameter_name))
+        return float(parameter_value)
+
+
+class SingleRegionAlchemicalState(GlobalParameterState):
+    """
+    Composable state to control lambda parameters for a single
+    alchemical molecule / region (``ligand_A``).
+
+    Parameters
+    ----------
+    parameters_name_suffix : str | None
+      If specified, the state will control a modified version of the parameter
+      ``lambda_restraints_{parameters_name_suffix}` instead of just
+      ``lambda_restraints``.
+    lambda_sterics_A : float | None
+      Control parameter for the vdW interactions for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_electrostatics_A : float | None
+      Control parameter for the electrostatics interactions for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_bonds_A : float | None
+      Control parameter for alchemically modified bonds for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_angles_A : float | None
+      Control parameter for alchemically modified angles for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_torsions_A : float | None
+      Control parameter for alchemically modified dihedrals for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_restraints_A : float | None
+      Control parameter for alchemically modified restraints for ligand A.
+
+    See Also
+    --------
+    :class:`openmmtools.states.GlobalParameterState`
+    :class:`openfe.protocols.openmm_utils.states.DualRegionAlchemicalState`
+    """
+
+    lambda_sterics_A = _LambdaParameter("lambda_sterics_A")
+    lambda_electrostatics_A = _LambdaParameter("lambda_electrostatics_A")
+    lambda_bonds_A = _LambdaParameter("lambda_bonds_A")
+    lambda_angles_A = _LambdaParameter("lambda_angles_A")
+    lambda_torsions_A = _LambdaParameter("lambda_torsions_A")
+    lambda_restraints_A = _LambdaParameter("lambda_restraints_A")
+
+
+class DualRegionAlchemicalState(SingleRegionAlchemicalState):
+    """
+    Composable state to control lambda parameters for a system
+    with two alchemical molecules / regions (``ligand_A`` and ``ligand_B``).
+
+    Parameters
+    ----------
+    parameters_name_suffix : str | None
+      If specified, the state will control a modified version of the parameter
+      ``lambda_restraints_A_{parameters_name_suffix}` instead of just
+      ``lambda_restraints_A``.
+    lambda_sterics_A : float | None
+      Control parameter for the vdW interactions for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_electrostatics_A : float | None
+      Control parameter for the electrostatics interactions for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_bonds_A : float | None
+      Control parameter for alchemically modified bonds for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_angles_A : float | None
+      Control parameter for alchemically modified angles for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_torsions_A : float | None
+      Control parameter for alchemically modified dihedrals for ligand A.
+      If defined, must be between 0 and 1.
+    lambda_restraints_A : float | None
+      Control parameter for alchemically modified restraints for ligand A.
+    lambda_sterics_B : float | None
+      Control parameter for the vdW interactions for ligand B.
+      If defined, must be between 0 and 1.
+    lambda_electrostatics_B : float | None
+      Control parameter for the electrostatics interactions for ligand B.
+      If defined, must be between 0 and 1.
+    lambda_bonds_B : float | None
+      Control parameter for alchemically modified bonds for ligand B.
+      If defined, must be between 0 and 1.
+    lambda_angles_B : float | None
+      Control parameter for alchemically modified angles for ligand B.
+      If defined, must be between 0 and 1.
+    lambda_torsions_B : float | None
+      Control parameter for alchemically modified dihedrals for ligand B.
+      If defined, must be between 0 and 1.
+    lambda_restraints_B : float | None
+      Control parameter for alchemically modified restraints for ligand B.
+
+    See Also
+    --------
+    :class:`openmmtools.states.GlobalParameterState`
+    :class:`openfe.protocols.openmm_utils.states.SingleRegionAlchemicalState`
+    """
+
+    lambda_sterics_B = _LambdaParameter("lambda_sterics_B")
+    lambda_electrostatics_B = _LambdaParameter("lambda_electrostatics_B")
+    lambda_bonds_B = _LambdaParameter("lambda_bonds_B")
+    lambda_angles_B = _LambdaParameter("lambda_angles_B")
+    lambda_torsions_B = _LambdaParameter("lambda_torsions_B")
+    lambda_restraints_B = _LambdaParameter("lambda_restraints_B")
diff --git a/src/openfe/protocols/restraint_utils/settings.py b/src/openfe/protocols/restraint_utils/settings.py
index 231b6c141..8ab8dcd86 100644
--- a/src/openfe/protocols/restraint_utils/settings.py
+++ b/src/openfe/protocols/restraint_utils/settings.py
@@ -174,17 +174,6 @@ class BoreschRestraintSettings(BaseRestraintSettings):
     Boresch-like restraint search parameter.
     The maximum distance between any host atom and the guest G0 atom. Must be in units compatible with nanometer.
     """
-    # TODO: re-enable this (Issue #1555)
-    # host_atoms: Optional[list[int]] = None
-    # """
-    # The indices of the host component atoms to restrain.
-    # If defined, these will override any automatic selection.
-    # """
-    # guest_atoms: Optional[list[int]] = None
-    # """
-    # The indices of the guest component atoms to restraint.
-    # If defined, these will override any automatic selection.
-    # """
     anchor_finding_strategy: Literal["multi-residue", "bonded"] = "bonded"
     """
     The Boresch atom picking strategy to use.
@@ -193,11 +182,3 @@ class BoreschRestraintSettings(BaseRestraintSettings):
       * `bonded`: pick host atoms that are bonded to each other.
       * `multi-residue`: pick host atoms which can span multiple residues.
     """
-
-
-#     @field_validator("guest_atoms", "host_atoms")
-#     def positive_idxs_list(cls, v):
-#         if v is not None and any([i < 0 for i in v]):
-#             errmsg = "negative indices passed"
-#             raise ValueError(errmsg)
-#         return v
diff --git a/src/openfe/tests/conftest.py b/src/openfe/tests/conftest.py
index fcd41412f..10e7047b0 100644
--- a/src/openfe/tests/conftest.py
+++ b/src/openfe/tests/conftest.py
@@ -233,6 +233,17 @@ def benzene_modifications():
     return files
 
 
+@pytest.fixture(scope="session")
+def benzene_modifications_am1bcc():
+    files = {}
+    with resources.as_file(resources.files("openfe.tests.data")) as d:
+        fn = str(d / "benzene_modifications_am1bcc.sdf")
+        supp = Chem.SDMolSupplier(str(fn), removeHs=False)
+        for rdmol in supp:
+            files[rdmol.GetProp("_Name")] = SmallMoleculeComponent(rdmol)
+    return files
+
+
 @pytest.fixture(scope="session")
 def charged_benzene_modifications():
     files = {}
@@ -244,6 +255,25 @@ def charged_benzene_modifications():
     return files
 
 
+@pytest.fixture(scope="session")
+def OA_guests_charged() -> dict[str, SmallMoleculeComponent]:
+    files = {}
+    with resources.as_file(resources.files("openfe.tests.data.host_guest")) as d:
+        fn = str(d / "OA_guests_nagl_charged.sdf")
+        supp = Chem.SDMolSupplier(str(fn), removeHs=False)
+        for rdmol in supp:
+            files[rdmol.GetProp("_Name")] = SmallMoleculeComponent(rdmol)
+    return files
+
+
+@pytest.fixture(scope="session")
+def OA_host_charged() -> SmallMoleculeComponent:
+    with resources.as_file(resources.files("openfe.tests.data.host_guest")) as d:
+        fn = str(d / "OA_host_nagl_charges.sdf")
+        mol = [m for m in Chem.SDMolSupplier(str(fn), removeHs=False)][0]
+        yield SmallMoleculeComponent(mol)
+
+
 @pytest.fixture(scope="session")
 def T4L_septop_reference_xml():
     with resources.as_file(resources.files("openfe.tests.data.openmm_septop")) as d:
diff --git a/src/openfe/tests/data/benzene_modifications_am1bcc.sdf b/src/openfe/tests/data/benzene_modifications_am1bcc.sdf
new file mode 100644
index 000000000..206c82d61
--- /dev/null
+++ b/src/openfe/tests/data/benzene_modifications_am1bcc.sdf
@@ -0,0 +1,295 @@
+benzene
+     RDKit          3D
+
+ 12 12  0  0  0  0  0  0  0  0999 V2000
+   25.9780    5.3270    4.7790 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.3950    5.0740    3.4990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.3400    5.8600    2.9020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8370    6.9210    3.5690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.4200    7.1960    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.4980    6.3790    5.4690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.2298    4.6859    5.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9676    4.2351    2.9497 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.6890    5.6311    1.8951 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.5730    7.5660    3.0889 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8209    8.0598    5.3863 H   0  0  0  0  0  0  0  0  0  0  0  0
+   26.1874    6.5720    6.4958 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  2  0
+  1  6  1  0
+  1  7  1  0
+  2  3  1  0
+  2  8  1  0
+  3  4  2  0
+  3  9  1  0
+  4  5  1  0
+  4 10  1  0
+  5  6  2  0
+  5 11  1  0
+  6 12  1  0
+M  END
+
+>  <ofe-name>
+benzene
+
+>  <atom.dprop.PartialCharge>
+-0.13 -0.13 -0.13 -0.13 -0.13 -0.13 0.13 0.13 0.13 0.13 0.13 0.13
+
+$$$$
+toluene
+     RDKit          3D
+
+ 15 15  0  0  0  0  0  0  0  0999 V2000
+   28.9072    8.7434    5.1220 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.1966    8.1433    6.6393 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.9864    8.4164    5.6052 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.2579    9.2269    5.5838 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9780    5.3270    4.7790 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.3950    5.0740    3.4990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.3400    5.8600    2.9020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8370    6.9210    3.5690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.4200    7.1960    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.4980    6.3790    5.4690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.2298    4.6859    5.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9676    4.2351    2.9497 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.6890    5.6311    1.8951 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.5730    7.5660    3.0889 H   0  0  0  0  0  0  0  0  0  0  0  0
+   26.1874    6.5720    6.4958 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  3  1  0
+  2  3  1  0
+  3  4  1  0
+  5  6  2  0
+  5 10  1  0
+  5 11  1  0
+  6  7  1  0
+  6 12  1  0
+  7  8  2  0
+  7 13  1  0
+  8  9  1  0
+  8 14  1  0
+  3  9  1  0
+  9 10  2  0
+ 10 15  1  0
+M  END
+
+>  <ofe-name>
+toluene
+
+>  <atom.dprop.PartialCharge>
+0.044033066666666669 0.044033066666666669 -0.053799933333333334 0.044033066666666669 -0.12699993333333334 -0.13499993333333335 -0.12699993333333334 -0.13099993333333335
+-0.077299933333333321 -0.13099993333333335 0.13000006666666666 0.13000006666666666 0.13000006666666666 0.13000006666666666 0.13000006666666666
+
+$$$$
+phenol
+     RDKit          3D
+
+ 13 13  0  0  0  0  0  0  0  0999 V2000
+   25.9780    5.3270    4.7790 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.3950    5.0740    3.4990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.3400    5.8600    2.9020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8370    6.9210    3.5690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.4200    7.1960    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.4980    6.3790    5.4690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.2298    4.6859    5.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9676    4.2351    2.9497 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.6890    5.6311    1.8951 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.5730    7.5660    3.0889 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.1311    8.0887    6.4624 H   0  0  0  0  0  0  0  0  0  0  0  0
+   26.1874    6.5720    6.4958 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.9460    8.3293    5.5517 O   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  2  0
+  1  6  1  0
+  1  7  1  0
+  2  3  1  0
+  2  8  1  0
+  3  4  2  0
+  3  9  1  0
+  4  5  1  0
+  4 10  1  0
+  5  6  2  0
+  5 13  1  0
+  6 12  1  0
+ 11 13  1  0
+M  END
+
+>  <ofe-name>
+phenol
+
+>  <atom.dprop.PartialCharge>
+-0.094423076923076915 -0.16592307692307692 -0.094423076923076915 -0.18492307692307691 0.12317692307692309 -0.18492307692307691 0.13307692307692309 0.13307692307692309 0.13307692307692309
+0.14157692307692307 0.41807692307692307 0.14157692307692307 -0.4990230769230769
+
+$$$$
+benzonitrile
+     RDKit          3D
+
+ 13 13  0  0  0  0  0  0  0  0999 V2000
+   28.5559    9.5700    6.2831 N   0  0  0  0  0  0  0  0  0  0  0  0
+   27.9981    8.4043    5.5824 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9780    5.3270    4.7790 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.3950    5.0740    3.4990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.3400    5.8600    2.9020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8370    6.9210    3.5690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.4200    7.1960    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.4980    6.3790    5.4690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.2298    4.6859    5.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9676    4.2351    2.9497 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.6890    5.6311    1.8951 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.5730    7.5660    3.0889 H   0  0  0  0  0  0  0  0  0  0  0  0
+   26.1874    6.5720    6.4958 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  3  0
+  3  4  2  0
+  3  8  1  0
+  3  9  1  0
+  4  5  1  0
+  4 10  1  0
+  5  6  2  0
+  5 11  1  0
+  6  7  1  0
+  6 12  1  0
+  2  7  1  0
+  7  8  2  0
+  8 13  1  0
+M  END
+
+>  <ofe-name>
+benzonitrile
+
+>  <atom.dprop.PartialCharge>
+-0.36380000000000001 0.23380000000000001 -0.13500000000000001 -0.107 -0.13500000000000001 -0.090999999999999998 -0.019000000000000003 -0.090999999999999998 0.14000000000000001
+0.13800000000000001 0.14000000000000001 0.14499999999999999 0.14499999999999999
+
+$$$$
+anisole
+     RDKit          3D
+
+ 16 16  0  0  0  0  0  0  0  0999 V2000
+   29.2873    8.8784    4.9226 C   0  0  0  0  0  0  0  0  0  0  0  0
+   29.5502    9.7990    5.4437 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.0548    8.3459    5.4720 O   0  0  0  0  0  0  0  0  0  0  0  0
+   30.0866    8.1484    5.0502 H   0  0  0  0  0  0  0  0  0  0  0  0
+   29.1525    9.0868    3.8612 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9780    5.3270    4.7790 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.3950    5.0740    3.4990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.3400    5.8600    2.9020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8370    6.9210    3.5690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.4200    7.1960    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.4980    6.3790    5.4690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.2298    4.6859    5.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9676    4.2351    2.9497 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.6890    5.6311    1.8951 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.5730    7.5660    3.0889 H   0  0  0  0  0  0  0  0  0  0  0  0
+   26.1874    6.5720    6.4958 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  1  0
+  1  4  1  0
+  1  5  1  0
+  1  3  1  0
+  6  7  2  0
+  6 11  1  0
+  6 12  1  0
+  7  8  1  0
+  7 13  1  0
+  8  9  2  0
+  8 14  1  0
+  9 10  1  0
+  9 15  1  0
+  3 10  1  0
+ 10 11  2  0
+ 11 16  1  0
+M  END
+
+>  <ofe-name>
+anisole
+
+>  <atom.dprop.PartialCharge>
+0.113825 0.043825000000000003 -0.32877500000000004 0.043825000000000003 0.043825000000000003 -0.097875000000000004 -0.16487500000000002 -0.097875000000000004 -0.17987500000000001 0.123225
+-0.17987500000000001 0.13212499999999999 0.13212499999999999 0.13212499999999999 0.14212499999999997 0.14212499999999997
+
+$$$$
+benzaldehyde
+     RDKit          3D
+
+ 14 14  0  0  0  0  0  0  0  0999 V2000
+   29.2079    8.8492    4.9632 O   0  0  0  0  0  0  0  0  0  0  0  0
+   27.5482    8.8691    6.4597 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.9981    8.4043    5.5824 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9780    5.3270    4.7790 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.3950    5.0740    3.4990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.3400    5.8600    2.9020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8370    6.9210    3.5690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.4200    7.1960    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.4980    6.3790    5.4690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.2298    4.6859    5.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9676    4.2351    2.9497 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.6890    5.6311    1.8951 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.5730    7.5660    3.0889 H   0  0  0  0  0  0  0  0  0  0  0  0
+   26.1874    6.5720    6.4958 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  3  2  0
+  2  3  1  0
+  4  5  2  0
+  4  9  1  0
+  4 10  1  0
+  5  6  1  0
+  5 11  1  0
+  6  7  2  0
+  6 12  1  0
+  7  8  1  0
+  7 13  1  0
+  3  8  1  0
+  8  9  2  0
+  9 14  1  0
+M  END
+
+>  <ofe-name>
+benzaldehyde
+
+>  <atom.dprop.PartialCharge>
+-0.52817142857142862 -0.0028714285714285795 0.5754285714285714 -0.14507142857142857 -0.098071428571428587 -0.14507142857142857 -0.078071428571428583 -0.19767142857142858
+-0.078071428571428583 0.13742857142857143 0.13492857142857143 0.13742857142857143 0.14392857142857141 0.14392857142857141
+
+$$$$
+styrene
+     RDKit          3D
+
+ 16 16  0  0  0  0  0  0  0  0999 V2000
+   29.2873    8.8784    4.9226 C   0  0  0  0  0  0  0  0  0  0  0  0
+   29.6609    8.3486    4.0463 H   0  0  0  0  0  0  0  0  0  0  0  0
+   29.8344    9.7353    5.3157 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.5365    8.8812    6.4825 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.9864    8.4164    5.6052 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9780    5.3270    4.7790 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.3950    5.0740    3.4990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.3400    5.8600    2.9020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.8370    6.9210    3.5690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   27.4200    7.1960    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   26.4980    6.3790    5.4690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   25.2298    4.6859    5.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
+   25.9676    4.2351    2.9497 H   0  0  0  0  0  0  0  0  0  0  0  0
+   27.6890    5.6311    1.8951 H   0  0  0  0  0  0  0  0  0  0  0  0
+   28.5730    7.5660    3.0889 H   0  0  0  0  0  0  0  0  0  0  0  0
+   26.1874    6.5720    6.4958 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  1  0
+  1  3  1  0
+  1  5  2  0
+  4  5  1  0
+  6  7  2  0
+  6 11  1  0
+  6 12  1  0
+  7  8  1  0
+  7 13  1  0
+  8  9  2  0
+  8 14  1  0
+  9 10  1  0
+  9 15  1  0
+  5 10  1  0
+ 10 11  2  0
+ 11 16  1  0
+M  END
+
+>  <ofe-name>
+styrene
+
+>  <atom.dprop.PartialCharge>
+-0.20899999999999999 0.11349999999999999 0.11349999999999999 0.123 -0.1132 -0.13100000000000001 -0.127 -0.13100000000000001 -0.11849999999999999 -0.057800000000000004 -0.11849999999999999
+0.13100000000000001 0.13100000000000001 0.13100000000000001 0.13150000000000001 0.13150000000000001
+
+$$$$
diff --git a/src/openfe/tests/data/host_guest/OA_guests_nagl_charged.sdf b/src/openfe/tests/data/host_guest/OA_guests_nagl_charged.sdf
new file mode 100644
index 000000000..b1f356172
--- /dev/null
+++ b/src/openfe/tests/data/host_guest/OA_guests_nagl_charged.sdf
@@ -0,0 +1,414 @@
+OA-G0
+     RDKit          3D
+
+ 20 20  0  0  0  0  0  0  0  0999 V2000
+   -0.7955    1.2297   -3.7854 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3091   -0.3847   -0.7486 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5096    0.2717    0.3768 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2838   -0.8860   -1.7613 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9597    0.1420   -0.0138 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9546   -0.0222   -1.5321 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.8697    1.3180   -2.2595 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2049    0.6341   -4.2709 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.7693    1.7767   -4.3800 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0136    0.3154   -1.2099 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.8934   -1.2223   -0.3497 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7004   -0.2524    1.3210 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7925    1.3196    0.5218 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.6733   -0.7466   -2.7812 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0732   -1.9566   -1.6356 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4456   -0.7024    0.4948 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5023    1.0547    0.2785 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8449   -0.5532   -1.8532 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0170    1.8572   -1.9078 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.7461    1.9209   -1.9970 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  7  1  0
+  1  8  1  0
+  1  9  2  0
+  2  3  1  0
+  2  4  1  0
+  3  5  1  0
+  4  6  1  0
+  5  6  1  0
+  6  7  1  0
+  2 10  1  0
+  2 11  1  0
+  3 12  1  0
+  3 13  1  0
+  4 14  1  0
+  4 15  1  0
+  5 16  1  0
+  5 17  1  0
+  6 18  1  0
+  7 19  1  0
+  7 20  1  0
+M  CHG  1   8  -1
+M  END
+
+>  <Chemgauss4>
+-7.9565659999999996
+
+>  <ofe-name>
+OA-G0
+
+>  <atom.dprop.PartialCharge>
+0.90913933776319023 -0.080904421582818034 -0.080904421582818034 -0.08442279435694218 -0.08442279435694218 -0.062293932959437373 -0.19467459358274936 -0.846026298776269 -0.846026298776269
+0.025389452278614045 0.025389452278614045 0.025389452278614045 0.025389452278614045 0.043140637502074239 0.043140637502074239 0.043140637502074239 0.043140637502074239 0.053682352229952809
+0.021366753429174424 0.021366753429174424
+
+$$$$
+OA-G1
+     RDKit          3D
+
+ 17 16  0  0  0  0  0  0  0  0999 V2000
+   -1.0006   -0.7469   -3.8833 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.0084   -0.5611   -3.0020 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.8433   -0.3922   -5.3688 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.0381    0.0211    0.8063 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1259   -0.8917   -1.5409 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.1069    0.3450   -0.6725 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2681    0.0857   -5.7303 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8736   -0.6292   -6.0659 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.9447   -1.1634   -3.5343 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.9430   -0.1434   -3.3206 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0312   -0.3847    1.0235 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0983    0.9255    1.4073 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7104   -0.7126    1.1224 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1300   -1.2920   -1.3561 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5854   -1.6783   -1.2604 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1106    0.7483   -0.8545 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6092    1.1326   -0.9379 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  2  0
+  1  3  1  0
+  2  5  1  0
+  3  7  1  0
+  3  8  2  0
+  4  6  1  0
+  5  6  1  0
+  1  9  1  0
+  2 10  1  0
+  4 11  1  0
+  4 12  1  0
+  4 13  1  0
+  5 14  1  0
+  5 15  1  0
+  6 16  1  0
+  6 17  1  0
+M  CHG  1   7  -1
+M  END
+
+>  <Chemgauss4>
+-6.8312020000000002
+
+>  <ofe-name>
+OA-G1
+
+>  <atom.dprop.PartialCharge>
+-0.17128288493875196 -0.23447102056268385 0.90111100925680465 -0.10154167310718228 -0.040073477937018168 -0.082682838006054651 -0.84050458417657548 -0.84050458417657548 0.07679569973226856
+0.11721654488321613 0.027137529762352213 0.027137529762352213 0.027137529762352213 0.031102622585261568 0.031102622585261568 0.036159987287486303 0.036159987287486303
+
+$$$$
+OA-G2
+     RDKit          3D
+
+ 25 25  0  0  0  0  0  0  0  0999 V2000
+    0.2241   -1.1357   -4.2804 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4652    0.1472   -4.5856 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.7894   -0.5542   -1.5717 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1825    0.4275   -5.8319 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5566   -0.2851   -1.1171 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4939   -1.5520   -3.0259 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0275    1.3042   -3.7403 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0051    0.9021   -2.6793 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6591   -0.4247   -1.9881 C   0  0  2  0  0  0  0  0  0  0  0  0
+    0.3129    0.1732    0.2902 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.6132    0.0157   -6.8813 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.2758    1.0443   -5.6932 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5476   -1.9343   -4.9434 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.6613   -0.4502   -0.9346 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.9554   -0.8848   -2.5917 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4888   -1.9348   -3.3003 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0542   -2.4021   -2.5903 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8919    1.7520   -3.2369 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4008    2.0735   -4.3938 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0802    1.7062   -1.9245 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.0002    0.8451   -3.1539 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5025   -0.6880   -1.3241 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.3537    1.2659    0.3463 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6741   -0.1560    0.6311 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0569   -0.2247    0.9886 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  2  0
+  1  6  1  0
+  2  4  1  0
+  2  7  1  0
+  3  5  2  0
+  4 11  1  0
+  4 12  2  0
+  5  9  1  0
+  5 10  1  0
+  6  9  1  0
+  7  8  1  0
+  8  9  1  0
+  1 13  1  0
+  3 14  1  0
+  3 15  1  0
+  6 16  1  0
+  6 17  1  0
+  7 18  1  0
+  7 19  1  0
+  8 20  1  0
+  8 21  1  0
+  9 22  1  1
+ 10 23  1  0
+ 10 24  1  0
+ 10 25  1  0
+M  CHG  1  11  -1
+M  END
+
+>  <Chemgauss4>
+-9.8733339999999998
+
+>  <ofe-name>
+OA-G2
+
+>  <atom.dprop.PartialCharge>
+-0.24234159156680107 -0.13431934878230095 -0.236714898198843 0.91454404726624494 -0.1058359132707119 -0.0076737576723098751 -0.029078564196825026 -0.066287456601858141
+-0.033445252627134325 -0.06594990059733391 -0.83899009093642229 -0.83899009093642229 0.11854273959994316 0.10526807740330696 0.10526807740330696 0.021749406531453134 0.021749406531453134
+0.047736430019140241 0.047736430019140241 0.03281161695718765 0.03281161695718765 0.046814215779304502 0.034864944815635679 0.034864944815635679 0.034864911288022993
+
+$$$$
+OA-G3
+     RDKit          3D
+
+ 17 16  0  0  0  0  0  0  0  0999 V2000
+   -0.2421   -0.0718    0.3818 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8393   -0.2077   -0.3937 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3051    0.2551   -3.3541 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1850    0.7201   -1.5211 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1018   -0.6816   -3.2022 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2330   -0.0262   -2.8578 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5879    0.9108   -2.3102 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8592    0.2770   -4.4873 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4542   -0.7664    1.1869 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.9353    0.7462    0.2188 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.5207   -1.0316   -0.1977 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.1638    1.1681   -1.3121 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4650    1.5450   -1.5737 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.3602   -1.3853   -2.4008 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0074   -1.2732   -4.1197 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4879    0.6810   -3.6576 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0316   -0.7782   -2.8399 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  2  0
+  2  4  1  0
+  3  5  1  0
+  3  7  1  0
+  3  8  2  0
+  4  6  1  0
+  5  6  1  0
+  1  9  1  0
+  1 10  1  0
+  2 11  1  0
+  4 12  1  0
+  4 13  1  0
+  5 14  1  0
+  5 15  1  0
+  6 16  1  0
+  6 17  1  0
+M  CHG  1   7  -1
+M  END
+
+>  <Chemgauss4>
+-7.7345790000000001
+
+>  <ofe-name>
+OA-G3
+
+>  <atom.dprop.PartialCharge>
+-0.25688242101494002 -0.14571022176567247 0.90786136007484264 -0.056502479402457964 -0.21427863025489977 -0.055603332160150301 -0.84741937303367787 -0.84741937303367787 0.10201528734144043
+0.10201528734144043 0.10772500997957062 0.038299766095245588 0.038299766095245588 0.016157646389568552 0.016157646389568552 0.047642030479276884 0.047642030479276884
+
+$$$$
+OA-G4
+     RDKit          3D
+
+ 29 28  0  0  0  0  0  0  0  0999 V2000
+   -0.3023   -0.8943   -1.9056 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4373   -0.5293   -0.8413 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0556    1.4485   -5.1381 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1104   -0.2112    0.5212 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.9350   -0.4263   -0.9512 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.9787    1.0573   -5.6003 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.1827   -1.2245   -3.2898 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.0475    1.5230   -4.0753 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.8663   -0.8574   -4.3441 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2784    0.6305   -4.3059 C   0  0  2  0  0  0  0  0  0  0  0  0
+    2.0098    2.2581   -4.9567 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8803    0.6125   -6.0669 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3802   -0.9841   -1.7762 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1488   -1.1197    1.1320 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1238    0.1959    0.4987 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5344    0.5154    1.0259 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.2639    0.5732   -1.2544 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.3111   -1.1387   -1.6928 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.4012   -0.6500    0.0139 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3241    0.9347   -6.4681 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.2613    2.1149   -5.5589 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.8879    0.4712   -5.7716 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.3887   -2.3007   -3.3170 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1330   -0.7307   -3.5143 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.3900    1.2984   -3.0936 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.3600    2.5753   -4.0472 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4864   -1.1087   -5.3434 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.7539   -1.4826   -4.1850 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.9857    0.7778   -3.4799 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  2  0
+  1  7  1  0
+  2  4  1  0
+  2  5  1  0
+  3  8  1  0
+  3 11  1  0
+  3 12  2  0
+  6 10  1  0
+  7  9  1  0
+  8 10  1  0
+  9 10  1  0
+  1 13  1  0
+  4 14  1  0
+  4 15  1  0
+  4 16  1  0
+  5 17  1  0
+  5 18  1  0
+  5 19  1  0
+  6 20  1  0
+  6 21  1  0
+  6 22  1  0
+  7 23  1  0
+  7 24  1  0
+  8 25  1  0
+  8 26  1  0
+  9 27  1  0
+  9 28  1  0
+ 10 29  1  1
+M  CHG  1  11  -1
+M  END
+
+>  <Chemgauss4>
+-9.1737099999999998
+
+>  <ofe-name>
+OA-G4
+
+>  <atom.dprop.PartialCharge>
+-0.16766302782142983 -0.12061256664837229 0.90602863830482139 -0.070390922306426643 -0.070390922306426643 -0.078514670310863136 -0.042349432884105323 -0.19354121881569253
+-0.064572726665385841 -0.05895461264098513 -0.84640865522468911 -0.84640865522468911 0.11097088706647527 0.034675578297726037 0.034675578297726037 0.034675578297726037 0.034675578297726037
+0.034675578297726037 0.034675578297726037 0.024976786868325596 0.024976786868325596 0.024976786868325596 0.042367330463281991 0.042367330463281991 0.018510187687031155 0.018510187687031155
+0.037587809087387444 0.037587809087387444 0.062893400611034753
+
+$$$$
+OA-G5
+     RDKit          3D
+
+ 17 16  0  0  0  0  0  0  0  0999 V2000
+    0.4721   -0.3596    0.1049 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6140    0.0086   -0.5838 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.3138   -0.7087   -3.7837 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.6475    1.3143   -4.0310 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.5825    0.9091   -1.7860 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2344    0.2559   -3.0098 C   0  0  2  0  0  0  0  0  0  0  0  0
+    0.5382   -1.3387   -3.0968 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.5250   -0.7473   -5.0305 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4022   -1.0124    0.9678 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4547   -0.0073   -0.1892 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5854   -0.3514   -0.2552 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7758    1.8832   -4.3717 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.3607    2.0256   -3.5996 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.1200    0.8627   -4.9107 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4463    1.2159   -2.0129 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1278    1.8263   -1.5296 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.1320   -0.2920   -2.6927 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  2  0
+  2  5  1  0
+  3  6  1  0
+  3  7  1  0
+  3  8  2  0
+  4  6  1  0
+  5  6  1  0
+  1  9  1  0
+  1 10  1  0
+  2 11  1  0
+  4 12  1  0
+  4 13  1  0
+  4 14  1  0
+  5 15  1  0
+  5 16  1  0
+  6 17  1  1
+M  CHG  1   7  -1
+M  END
+
+>  <Chemgauss4>
+-7.6388220000000002
+
+>  <ofe-name>
+OA-G5
+
+>  <atom.dprop.PartialCharge>
+-0.2979469792369534 -0.10777685506378903 0.91644908027613869 -0.067875248325221682 -0.053967733812682772 -0.20121777517830625 -0.84354089660679588 -0.84354089660679588 0.094638693201191282
+0.094638693201191282 0.11095144884551272 0.018116577583200792 0.018116577583200792 0.018116577583200792 0.05721881898010478 0.05721881898010478 0.030401098596699098
+
+$$$$
+OA-G6
+     RDKit          3D
+
+ 19 18  0  0  0  0  0  0  0  0999 V2000
+   -1.3428    0.3977   -4.1555 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3494    0.8562   -1.6206 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1465    0.0031    0.1950 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4127   -0.7040   -3.6398 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4163   -0.9172   -2.1271 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.0498    0.3297   -1.2978 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4251    0.5122   -5.4098 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.9062    1.0746   -3.2485 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.1072    0.0897   -1.4250 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.6012    1.7234   -1.0027 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4304    1.1499   -2.6726 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5555   -0.7928    0.4655 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1538   -0.3402    0.4542 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0812    0.8805    0.8099 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.6080   -0.5200   -3.9934 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7392   -1.6346   -4.1217 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2862   -1.7273   -1.8917 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4072   -1.2739   -1.8180 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7712    1.1284   -1.5083 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  4  1  0
+  1  7  1  0
+  1  8  2  0
+  2  6  1  0
+  3  6  1  0
+  4  5  1  0
+  5  6  1  0
+  2  9  1  0
+  2 10  1  0
+  2 11  1  0
+  3 12  1  0
+  3 13  1  0
+  3 14  1  0
+  4 15  1  0
+  4 16  1  0
+  5 17  1  0
+  5 18  1  0
+  6 19  1  0
+M  CHG  1   7  -1
+M  END
+
+>  <Chemgauss4>
+-7.5684560000000003
+
+>  <ofe-name>
+OA-G6
+
+>  <atom.dprop.PartialCharge>
+0.90856356252180903 -0.083272322228080342 -0.083272322228080342 -0.21121091317189367 -0.054996575375920849 -0.083306490590697835 -0.84895735155594976 -0.84895735155594976
+0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022056599000566883 0.022712499687546177 0.022712499687546177 0.034720391819351597
+0.034720391819351597 0.0582043871675667
+
+$$$$
diff --git a/src/openfe/tests/data/host_guest/OA_host_nagl_charges.sdf b/src/openfe/tests/data/host_guest/OA_host_nagl_charges.sdf
new file mode 100644
index 000000000..c49fe096f
--- /dev/null
+++ b/src/openfe/tests/data/host_guest/OA_host_nagl_charges.sdf
@@ -0,0 +1,422 @@
+
+     RDKit          3D
+
+184204  0  0  0  0  0  0  0  0999 V2000
+    1.9780   -2.9130    2.5560 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.7890   -0.8910    1.7980 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0950   -1.7370    3.3190 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.8090   -3.1130    1.4400 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.7500   -2.1260    1.0940 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.8730   -0.6510    2.8620 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.5480    0.2040    1.3130 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.9840    0.5840    3.7170 C   0  0  1  0  0  0  0  0  0  0  0  0
+    3.9320    1.3630    0.6550 C   0  0  1  0  0  0  0  0  0  0  0  0
+    3.6580    2.3950    1.6720 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.3640    2.5270    2.1880 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0150    1.6640    3.2190 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7110    1.8220    3.7590 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.5080    3.4420    1.6300 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.1860    2.7330    3.1850 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2210    3.5300    2.1100 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6470    4.3030    1.4010 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.6530    2.7550    3.6810 C   0  0  1  0  0  0  0  0  0  0  0  0
+   -2.5740    1.8700    2.8770 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.1320    2.4120    1.6520 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5380    3.6600    0.4250 C   0  0  1  0  0  0  0  0  0  0  0  0
+   -2.7830    3.6460    1.1110 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.5440   -4.1280    0.4970 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8190   -3.8400    2.8300 C   0  0  2  0  0  0  0  0  0  0  0  0
+   -0.4320   -3.4440    2.0490 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2820   -4.7530   -0.0200 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6540   -3.9910    0.7700 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.9750    1.6730    0.8700 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.2640    0.3670    1.1370 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8020   -3.7270    0.0330 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.7750   -2.9180    0.5540 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.8390   -0.4100    0.1310 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.8440   -2.5740   -0.2240 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.9300   -1.1640   -0.7330 C   0  0  1  0  0  0  0  0  0  0  0  0
+   -1.4680   -2.7260    2.6070 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.6720   -2.4640    1.9030 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.9270    0.5100    3.1980 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.9200   -1.7400    2.6170 C   0  0  2  0  0  0  0  0  0  0  0  0
+   -3.7670   -0.2670    2.3280 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.4480   -3.9480   -0.4130 C   0  0  2  0  0  0  0  0  0  0  0  0
+   -1.7230    2.5680    5.1830 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.1160    2.7320    5.7950 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.0940    2.4680    7.3220 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.1470    1.9470    7.8190 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.9710    0.3330    5.2440 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.4740    1.5040    6.0680 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.9650    1.5110    5.9710 C   0  0  0  0  0  0  0  0  0  0  0  0
+    5.6550    0.8410    6.7490 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.6190   -4.2620    4.3570 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.4420   -5.3440    4.4880 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2540   -6.6130    4.0330 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0480   -7.1760    4.7900 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.9730   -2.2770    4.1800 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.1780   -1.5800    4.8560 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.1500   -1.9220    6.3420 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.2720   -1.6220    7.0980 O   0  0  0  0  0  0  0  0  0  0  0  0
+    4.8530    1.8910   -0.4090 C   0  0  0  0  0  0  0  0  0  0  0  0
+    6.4190    2.7220   -2.5260 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.8060    3.1880   -0.8380 C   0  0  0  0  0  0  0  0  0  0  0  0
+    5.7010    1.0210   -1.0690 C   0  0  0  0  0  0  0  0  0  0  0  0
+    6.5160    1.4180   -2.1280 C   0  0  0  0  0  0  0  0  0  0  0  0
+    5.5390    3.5760   -1.8830 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.9240   -4.1690   -1.8430 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.5740   -4.1720   -4.5420 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0450   -4.6830   -2.7060 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.1640   -3.6710   -2.3470 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.5120   -3.7000   -3.6730 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3590   -4.6800   -4.0870 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.3640   -1.2630   -2.1710 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.1030   -1.2900   -4.8080 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.1920   -0.2470   -2.6580 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.9340   -2.3460   -3.0260 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.2790   -2.3050   -4.3980 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.5090   -0.2170   -4.0010 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5850    4.4050   -0.9270 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5140    5.2850   -3.5240 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.7380    4.3570   -1.7190 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.3570    4.8830   -1.4230 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.3650    5.2850   -2.7420 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.6800    4.7550   -2.9980 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.6560   -3.1320   -4.2310 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.9210   -3.3160   -5.2690 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -6.1530    0.7740   -4.5970 O   0  0  0  0  0  0  0  0  0  0  0  0
+    7.3710    0.5580   -2.7490 O   0  0  0  0  0  0  0  0  0  0  0  0
+    5.5640    4.8870   -2.3690 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.8000    4.8130   -3.8840 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.8470    5.5530   -3.3740 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.1250   -4.3820   -4.8430 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3880   -6.3310   -3.9960 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.7660   -4.2420   -5.0990 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.6200   -5.4880   -4.0900 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.7470   -6.4710   -3.6320 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.8840   -5.2010   -4.6980 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4360   -5.0770   -4.9660 O   0  0  0  0  0  0  0  0  0  0  0  0
+    7.0140   -0.7940   -2.8120 C   0  0  0  0  0  0  0  0  0  0  0  0
+    6.6030   -3.4800   -2.7760 C   0  0  0  0  0  0  0  0  0  0  0  0
+    7.8680   -1.6190   -2.0560 C   0  0  0  0  0  0  0  0  0  0  0  0
+    5.8570   -1.2990   -3.4370 C   0  0  0  0  0  0  0  0  0  0  0  0
+    5.6910   -2.6750   -3.4480 C   0  0  0  0  0  0  0  0  0  0  0  0
+    7.6440   -2.9510   -2.0130 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -6.1080    2.0390   -4.1700 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -6.0420    4.5660   -3.0460 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -7.2560    2.5070   -3.5050 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.9680    2.7640   -4.2470 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.9020    4.0350   -3.6650 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -7.2160    3.7920   -2.9240 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -8.3420    4.2640   -2.0140 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -9.2730    3.5770   -1.6770 O   0  0  0  0  0  0  0  0  0  0  0  0
+    3.5430    9.3010   -1.4370 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.6040    9.6750   -1.0440 O   0  0  0  0  0  0  0  0  0  0  0  0
+    8.5220   -3.8700   -1.0650 C   0  0  0  0  0  0  0  0  0  0  0  0
+    8.5050   -5.0240   -1.2770 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.0100    2.7470    7.9930 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -6.1860   -2.4440    6.7280 O   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0950   -6.9830    2.9000 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.1550   -7.5970   -2.7690 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.3410   -8.3000   -2.1910 O   0  0  0  0  0  0  0  0  0  0  0  0
+    2.5090    9.9980   -1.3550 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -8.3240    5.4290   -1.6400 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.3750   -7.7630   -2.5820 O   0  0  0  0  0  0  0  0  0  0  0  0
+    9.0540   -3.4610   -0.0610 O   0  0  0  0  0  0  0  0  0  0  0  0
+    5.5220    2.2790    5.1950 O   0  0  0  0  0  0  0  0  0  0  0  0
+    4.4110    5.6700   -2.4190 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.1750    7.2800   -2.2470 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.1540    5.1300   -2.8380 C   0  0  0  0  0  0  0  0  0  0  0  0
+    4.5390    7.0130   -1.9880 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.3890    7.8860   -1.8580 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.0370    5.9900   -2.7670 C   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3140   -1.5380    4.0470 H   0  0  0  0  0  0  0  0  0  0  0  0
+    4.3950   -2.2610    0.2290 H   0  0  0  0  0  0  0  0  0  0  0  0
+    3.9410    1.0910    3.5770 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.9520    1.0690    0.2500 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.4670    1.0980    4.5330 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.8430    4.0220    0.7740 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.9970    3.7840    3.5520 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2860    2.6080    0.2230 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1220   -4.7900    2.3850 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.3700    2.1250   -0.0360 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.8910   -4.1970   -0.9440 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.8730   -0.8670   -0.6830 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5000   -2.4000    3.6440 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.5320    0.0400    4.0950 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.8830   -1.9290    2.1380 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.0410   -2.9290   -0.3430 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2640    1.6140    5.4480 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0850    3.3470    5.6060 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.4140    3.7700    5.6300 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.7740    1.9840    5.3470 H   0  0  0  0  0  0  0  0  0  0  0  0
+    3.4060   -0.6560    5.4020 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.9960    0.0970    5.6760 H   0  0  0  0  0  0  0  0  0  0  0  0
+    3.0110    2.4500    5.7790 H   0  0  0  0  0  0  0  0  0  0  0  0
+    3.2700    1.1450    7.0800 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.5540   -4.5780    4.8230 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.2700   -3.3680    4.8770 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6790   -5.3400    5.5540 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.3680   -5.2640    3.9150 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.1990   -3.3400    4.2770 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.1350   -2.0140    4.8300 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.0920   -0.4960    4.7620 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -6.1050   -1.8690    4.3570 H   0  0  0  0  0  0  0  0  0  0  0  0
+    6.9260    3.1600   -3.3810 H   0  0  0  0  0  0  0  0  0  0  0  0
+    4.0220    3.8410   -0.4650 H   0  0  0  0  0  0  0  0  0  0  0  0
+    5.8180    0.0340   -0.6280 H   0  0  0  0  0  0  0  0  0  0  0  0
+    2.7470   -4.1360   -5.6140 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.0270   -4.9510   -2.4310 H   0  0  0  0  0  0  0  0  0  0  0  0
+    3.8240   -3.2270   -1.6060 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.3290   -1.2620   -5.8720 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -5.4360    0.6030   -2.0260 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.2850   -3.1120   -2.6070 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.5800    5.8250   -4.4650 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -3.6250    3.8900   -1.2980 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.5490    4.8370   -0.8240 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.6830   -7.0780   -3.6390 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4480   -3.3280   -5.5950 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.6980   -5.5820   -3.9820 H   0  0  0  0  0  0  0  0  0  0  0  0
+    6.2860   -4.5190   -2.7310 H   0  0  0  0  0  0  0  0  0  0  0  0
+    8.6710   -1.1640   -1.4820 H   0  0  0  0  0  0  0  0  0  0  0  0
+    5.0750   -0.6430   -3.8120 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -6.1300    5.5770   -2.6560 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -8.1780    1.9330   -3.5290 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -4.1080    2.2910   -4.7150 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2730    7.8750   -2.3700 H   0  0  0  0  0  0  0  0  0  0  0  0
+    3.0630    4.0900   -3.1410 H   0  0  0  0  0  0  0  0  0  0  0  0
+    5.4990    7.4650   -1.7500 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  3  2  0
+  1  4  1  0
+  4  5  2  0
+  2  5  1  0
+  3  6  1  0
+  2  6  2  0
+  2  7  1  0
+  6  8  1  0
+  7  9  1  0
+  9 10  1  0
+ 10 11  1  0
+ 11 12  2  0
+  8 12  1  0
+ 12 13  1  0
+ 11 14  1  0
+ 13 15  2  0
+ 15 16  1  0
+ 14 16  2  0
+ 16 17  1  0
+ 15 18  1  0
+ 18 19  1  0
+ 19 20  2  0
+ 17 21  1  0
+ 21 22  1  0
+ 20 22  1  0
+  4 23  1  0
+  1 24  1  0
+ 24 25  1  0
+ 26 27  1  0
+ 25 27  2  0
+ 20 28  1  0
+ 28 29  2  0
+ 27 30  1  0
+ 30 31  2  0
+ 29 32  1  0
+ 31 33  1  0
+ 32 34  1  0
+ 33 34  1  0
+ 25 35  1  0
+ 35 36  2  0
+ 31 36  1  0
+ 19 37  1  0
+ 36 38  1  0
+ 38 39  1  0
+ 37 39  2  0
+ 29 39  1  0
+ 26 40  1  0
+ 23 40  1  0
+ 18 41  1  0
+ 41 42  1  0
+ 42 43  1  0
+ 43 44  1  0
+  8 45  1  0
+ 45 46  1  0
+ 46 47  1  0
+ 47 48  1  0
+ 24 49  1  0
+ 49 50  1  0
+ 50 51  1  0
+ 51 52  1  0
+ 38 53  1  0
+ 53 54  1  0
+ 54 55  1  0
+ 55 56  2  0
+  9 57  1  0
+ 57 59  2  0
+ 57 60  1  0
+ 60 61  2  0
+ 58 61  1  0
+ 59 62  1  0
+ 58 62  2  0
+ 40 63  1  0
+ 63 65  2  0
+ 63 66  1  0
+ 64 67  1  0
+ 66 67  2  0
+ 65 68  1  0
+ 64 68  2  0
+ 34 69  1  0
+ 69 71  2  0
+ 69 72  1  0
+ 70 73  1  0
+ 72 73  2  0
+ 70 74  2  0
+ 71 74  1  0
+ 21 75  1  0
+ 75 77  2  0
+ 75 78  1  0
+ 78 79  2  0
+ 76 79  1  0
+ 77 80  1  0
+ 76 80  2  0
+ 67 81  1  0
+ 73 82  1  0
+ 74 83  1  0
+ 61 84  1  0
+ 62 85  1  0
+ 80 86  1  0
+ 79 87  1  0
+ 82 88  1  0
+ 88 90  2  0
+ 88 91  1  0
+ 91 92  2  0
+ 89 92  1  0
+ 89 93  2  0
+ 90 93  1  0
+ 93 94  1  0
+ 68 94  1  0
+ 84 95  1  0
+ 95 97  2  0
+ 95 98  1  0
+ 98 99  2  0
+ 96 99  1  0
+ 81 99  1  0
+ 97100  1  0
+ 96100  2  0
+ 83101  1  0
+101103  2  0
+101104  1  0
+104105  2  0
+102105  1  0
+ 86105  1  0
+103106  1  0
+102106  2  0
+106107  1  0
+107108  2  0
+109110  2  0
+100111  1  0
+111112  2  0
+ 43113  2  0
+ 55114  1  0
+ 51115  2  0
+ 92116  1  0
+116117  2  0
+109118  1  0
+107119  1  0
+116120  1  0
+111121  1  0
+ 47122  2  0
+ 85123  1  0
+123125  2  0
+123126  1  0
+124127  1  0
+109127  1  0
+126127  2  0
+125128  1  0
+124128  2  0
+ 87128  1  0
+  3129  1  0
+  5130  1  0
+  8131  1  1
+  9132  1  1
+ 13133  1  0
+ 14134  1  0
+ 18135  1  6
+ 21136  1  1
+ 24137  1  6
+ 28138  1  0
+ 30139  1  0
+ 34140  1  1
+ 35141  1  0
+ 37142  1  0
+ 38143  1  6
+ 40144  1  1
+ 41145  1  0
+ 41146  1  0
+ 42147  1  0
+ 42148  1  0
+ 45149  1  0
+ 45150  1  0
+ 46151  1  0
+ 46152  1  0
+ 49153  1  0
+ 49154  1  0
+ 50155  1  0
+ 50156  1  0
+ 53157  1  0
+ 53158  1  0
+ 54159  1  0
+ 54160  1  0
+ 58161  1  0
+ 59162  1  0
+ 60163  1  0
+ 64164  1  0
+ 65165  1  0
+ 66166  1  0
+ 70167  1  0
+ 71168  1  0
+ 72169  1  0
+ 76170  1  0
+ 77171  1  0
+ 78172  1  0
+ 89173  1  0
+ 90174  1  0
+ 91175  1  0
+ 96176  1  0
+ 97177  1  0
+ 98178  1  0
+102179  1  0
+103180  1  0
+104181  1  0
+124182  1  0
+125183  1  0
+126184  1  0
+M  CHG  8  44  -1  48  -1  52  -1 114  -1 118  -1 119  -1 120  -1 121  -1
+M  END
+
+>  <ofe-name>
+
+
+>  <atom.dprop.PartialCharge>
+-0.091408472103269203 0.10657664056381454 -0.067896451396138771 0.10657664056381454 -0.16556027118602526 -0.091408472103269203 -0.3436575490657402 -0.02315850887933503 0.31904197152218094
+-0.3436575490657402 0.10657664056381454 -0.091408531707913979 -0.067896451396138771 -0.16556027118602526 -0.091408472103269203 0.10657664056381454 -0.3436575490657402 -0.02315850887933503
+-0.091408531707913979 0.10657664056381454 0.31904197152218094 -0.3436575490657402 -0.3436575490657402 -0.02315850887933503 -0.091408531707913979 -0.3436575490657402 0.10657664056381454
+-0.16556027118602526 0.10657664056381454 -0.16556027118602526 0.10657664056381454 -0.3436575490657402 -0.3436575490657402 0.31904197152218094 -0.067896451396138771 -0.091408472103269203
+-0.067896451396138771 -0.02315850887933503 -0.091408472103269203 0.31904197152218094 -0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777
+-0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777 -0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777
+-0.037677105111272438 -0.19577824656406176 0.91143131669124833 -0.84737205092349777 -0.091567691010625465 -0.16220025305190813 -0.14929413382449877 -0.14929413382449877 0.1104655456121849
+0.1104655456121849 -0.091567691010625465 -0.16220025305190813 -0.14929410402217638 -0.14929413382449877 0.1104655456121849 0.1104655456121849 -0.091567691010625465 -0.16220025305190813
+-0.14929413382449877 -0.14929413382449877 0.1104655456121849 0.1104655456121849 -0.091567691010625465 -0.16220025305190813 -0.14929413382449877 -0.14929413382449877 0.1104655456121849
+0.1104655456121849 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 -0.26024511043468246 0.085777629571764366
+-0.11220825795570145 -0.18432724062839281 -0.11220825795570145 -0.082448165458829506 0.085777629571764366 -0.26024511043468246 0.085777629571764366 -0.11220825795570145
+-0.11220825795570145 -0.18432724062839281 0.085777629571764366 -0.082448165458829506 0.085777629571764366 -0.11220825795570145 -0.11220825795570145 -0.18432724062839281
+0.085777629571764366 -0.082448165458829506 0.89783174212536088 -0.81537824456134567 0.89783174212536088 -0.81537824456134567 0.89783174212536088 -0.81537824456134567 -0.84737205092349777
+-0.84737205092349777 -0.84737205092349777 0.89783174212536088 -0.81537824456134567 -0.81537824456134567 -0.81537824456134567 -0.81537824456134567 -0.81537824456134567 -0.84737205092349777
+0.085777629571764366 -0.11220825795570145 -0.1843272704307152 -0.11220825795570145 -0.082448165458829506 0.085777629571764366 0.12293871160110702 0.1574754160221504 0.077845622258989708
+0.048985436963646309 0.12293871160110702 0.1574754160221504 0.077845622258989708 0.048985436963646309 0.077845622258989708 0.1574754160221504 0.1574754160221504 0.048985436963646309
+0.12293871160110702 0.12293871160110702 0.077845622258989708 0.048985436963646309 0.066841882127134697 0.066841882127134697 0.022772260415165321 0.022772260415165321 0.066841882127134697
+0.066841882127134697 0.022772260415165321 0.022772260415165321 0.066841882127134697 0.066841882127134697 0.022772260415165321 0.022772260415165321 0.066841882127134697 0.066841882127134697
+0.022772260415165321 0.022772260415165321 0.15592389221748579 0.14880136068424452 0.14880136068424452 0.15592389221748579 0.14880136068424452 0.14880136068424452 0.15592389221748579
+0.14880136068424452 0.14880136068424452 0.15592389221748579 0.14880136068424452 0.14880136068424452 0.1685743671234535 0.13864818628391493 0.16857433732113111 0.16857433732113111
+0.16857433732113111 0.13864818628391493 0.16857433732113111 0.16857433732113111 0.13864818628391493 0.16857433732113111 0.13864812667927015 0.16857439692577589
+
+$$$$
diff --git a/src/openfe/tests/data/host_guest/__init__.py b/src/openfe/tests/data/host_guest/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf b/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf
index ea2735835..18a65a49c 100644
--- a/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf
+++ b/src/openfe/tests/data/openmm_rfe/charged_benzenes.sdf
@@ -98,3 +98,41 @@ benzoic_acid
   9 14  1  0  0  0  0
 M  END
 $$$$
+phthalic_acid
+  PyMOL3.1          3D                             0
+
+ 16 16  0  0  0  0  0  0  0  0999 V2000
+    1.4045   -0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.7445   -0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.2607   -1.0833    0.0000 O   0  5  0  0  0  0  0  0  0  0  0  0
+    0.7022    1.2164    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7023    1.2164    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    3.4265    1.0159    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.4045   -0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7023   -1.2164    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7023   -1.2164    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2540    2.1720    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -2.5079   -0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.2540   -2.1719    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2540   -2.1720    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.3723    2.3768    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
+    2.5685    2.2821    0.0000 O   0  5  0  0  0  0  0  0  0  0  0  0
+    0.8335    3.4754    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  1  0  0  0  0
+  1  4  2  0  0  0  0
+  1  9  1  0  0  0  0
+  2  3  1  0  0  0  0
+  2  6  2  0  0  0  0
+  4  5  1  0  0  0  0
+  5  7  2  0  0  0  0
+  5 10  1  0  0  0  0
+  7  8  1  0  0  0  0
+  7 11  1  0  0  0  0
+  8  9  2  0  0  0  0
+  8 12  1  0  0  0  0
+  9 13  1  0  0  0  0
+ 14 15  1  0  0  0  0
+ 14 16  2  0  0  0  0
+  4 14  1  0  0  0  0
+M  END
+$$$$
diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py
index 14539d968..22e487d09 100644
--- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py
+++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_energies.py
@@ -13,7 +13,6 @@
 from openmmtools.alchemy import (
     AbsoluteAlchemicalFactory,
     AlchemicalRegion,
-    AlchemicalState,
 )
 
 from openfe.protocols import openmm_afe
@@ -22,17 +21,7 @@
 )
 from openfe.protocols.openmm_utils.omm_settings import OpenMMSolvationSettings
 from openfe.protocols.openmm_utils.serialization import deserialize
-
-
-class AlchemStateRest(AlchemicalState):
-    """
-    A modified AlchemicalState for testing.
-
-    Note: we don't need this in the main protocol since we use composable
-    thermodynamic states.
-    """
-
-    lambda_restraints = AlchemicalState._LambdaParameter("lambda_restraints")
+from openfe.protocols.openmm_utils.states import SingleRegionAlchemicalState
 
 
 def get_alchemical_energy_components(alchemical_system, alchemical_state, positions, platform):
@@ -166,13 +155,11 @@ def get_energy_components(
         platform = Platform.getPlatformByName("Reference")
         alchemical_region = AlchemicalRegion(alchemical_atoms=indices)
 
-        alchemical_state = AlchemStateRest.from_system(
-            system, parameters_name_suffix=alchemical_region.name
-        )
+        alchemical_state = SingleRegionAlchemicalState.from_system(system)
 
-        alchemical_state.lambda_sterics = lambda_sterics
-        alchemical_state.lambda_electrostatics = lambda_electrostatics
-        alchemical_state.lambda_restraints = lambda_restraints
+        alchemical_state.lambda_sterics_A = lambda_sterics
+        alchemical_state.lambda_electrostatics_A = lambda_electrostatics
+        alchemical_state.lambda_restraints_A = lambda_restraints
 
         return get_alchemical_energy_components(
             system, alchemical_state, positions, platform=platform
@@ -182,7 +169,7 @@ def get_energy_components(
     def test_energies_regression(self, lambda_val, t4_reference_system, t4_validation_data):
         energies_ref = self.get_energy_components(
             t4_reference_system,
-            t4_validation_data["alchem_indices"],
+            t4_validation_data["alchemical_indices"],
             t4_validation_data["debug_positions"],
             lambda_val,
             lambda_val,
@@ -191,7 +178,7 @@ def test_energies_regression(self, lambda_val, t4_reference_system, t4_validatio
 
         energies_val = self.get_energy_components(
             t4_validation_data["alchem_system"],
-            t4_validation_data["alchem_indices"],
+            t4_validation_data["alchemical_indices"],
             t4_validation_data["debug_positions"],
             lambda_val,
             lambda_val,
@@ -223,7 +210,7 @@ def assert_energies(actual, expected, nonbonded_lower: bool):
         # lambda 1 on all
         energies = self.get_energy_components(
             t4_validation_data["alchem_system"],
-            t4_validation_data["alchem_indices"],
+            t4_validation_data["alchemical_indices"],
             t4_validation_data["debug_positions"],
             lambda_sterics=1.0,
             lambda_electrostatics=1.0,
@@ -236,7 +223,7 @@ def assert_energies(actual, expected, nonbonded_lower: bool):
 
         energies = self.get_energy_components(
             t4_validation_data["alchem_system"],
-            t4_validation_data["alchem_indices"],
+            t4_validation_data["alchemical_indices"],
             t4_validation_data["debug_positions"],
             lambda_sterics=1.0,
             lambda_electrostatics=1.0,
@@ -249,7 +236,7 @@ def assert_energies(actual, expected, nonbonded_lower: bool):
 
         energies = self.get_energy_components(
             t4_validation_data["alchem_system"],
-            t4_validation_data["alchem_indices"],
+            t4_validation_data["alchemical_indices"],
             t4_validation_data["debug_positions"],
             lambda_sterics=1.0,
             lambda_electrostatics=0.0,
@@ -261,16 +248,16 @@ def assert_energies(actual, expected, nonbonded_lower: bool):
 
         # turn off sterics
         expected = copy.deepcopy(energies)
-        expected["alchemically modified NonbondedForce for non-alchemical/alchemical sterics"] = (
-            0 * ommunit.kilojoule_per_mole
-        )
         expected[
-            "alchemically modified BondForce for non-alchemical/alchemical sterics exceptions"
+            "alchemically modified NonbondedForce for non-alchemical/alchemical sterics for region A"
+        ] = 0 * ommunit.kilojoule_per_mole
+        expected[
+            "alchemically modified BondForce for non-alchemical/alchemical sterics exceptions for region A"
         ] = 0 * ommunit.kilojoule_per_mole
 
         energies = self.get_energy_components(
             t4_validation_data["alchem_system"],
-            t4_validation_data["alchem_indices"],
+            t4_validation_data["alchemical_indices"],
             t4_validation_data["debug_positions"],
             lambda_sterics=0.0,
             lambda_electrostatics=0.0,
diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py
index 6ba27ca2b..ceb421e1c 100644
--- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py
+++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol.py
@@ -30,8 +30,8 @@
 from openmmtools.multistate.multistatesampler import MultiStateSampler
 from openmmtools.tests import test_alchemy
 from openmmtools.tests.test_alchemy import (
-    check_interacting_energy_components,
-    check_noninteracting_energy_components,
+    check_multi_interacting_energy_components,
+    check_multi_noninteracting_energy_components,
     compare_system_energies,
 )
 
@@ -41,6 +41,8 @@
 from openfe.protocols.openmm_afe import (
     AbsoluteBindingProtocol,
 )
+from openfe.protocols.openmm_afe.abfe_units import _get_mda_universe
+from openfe.protocols.restraint_utils.geometry import BoreschRestraintGeometry
 
 from .utils import UNIT_TYPES, _get_units
 
@@ -168,15 +170,13 @@ def test_mda_universe_error():
     when calling the mda Universe getter.
     """
     with pytest.raises(ValueError, match="No positions to create"):
-        _ = openmm_afe.ABFEComplexSetupUnit._get_mda_universe(
-            topology="foo", positions=None, trajectory=None
-        )
+        _ = _get_mda_universe(topology="foo", positions=None, trajectory=None)
 
 
 class TestT4LysozymeDryRun:
     solvent = SolventComponent(ion_concentration=0 * offunit.molar)
     num_all_not_water = 2634  # 9 counterions to neutralize
-    num_protein_component_atoms = 2614
+    num_host_component_atoms = 2614
     # No ions
     num_ligand_atoms = 12
     expected_complex_particles = 32607
@@ -436,14 +436,14 @@ def _test_energies(reference_system, alchemical_system, alchemical_regions, posi
             positions=positions,
         )
 
-        check_noninteracting_energy_components(
+        check_multi_noninteracting_energy_components(
             reference_system=reference_system,
             alchemical_system=alchemical_system,
             alchemical_regions=alchemical_regions,
             positions=positions,
         )
 
-        check_interacting_energy_components(
+        check_multi_interacting_energy_components(
             reference_system=reference_system,
             alchemical_system=alchemical_system,
             alchemical_regions=alchemical_regions,
@@ -461,6 +461,7 @@ def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings,
             system=setup_results["alchem_system"],
             positions=setup_results["debug_positions"],
             selection_indices=setup_results["selection_indices"],
+            alchemical_indices=setup_results["alchemical_indices"],
             box_vectors=setup_results["box_vectors"],
             alchemical_restraints=True,
             dry=True,
@@ -480,9 +481,9 @@ def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings,
 
         # Check the alchemical indices
         expected_indices = [
-            i + self.num_protein_component_atoms - 1 for i in range(self.num_ligand_atoms)
+            i + self.num_host_component_atoms - 1 for i in range(self.num_ligand_atoms)
         ]
-        assert expected_indices == setup_results["alchem_indices"]
+        assert expected_indices == setup_results["alchemical_indices"]["A"]
 
         # Check the non-alchemical system
         assert setup_results["standard_system"].getNumParticles() == self.expected_complex_particles
@@ -503,11 +504,19 @@ def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings,
         assert pdb.n_atoms == self.num_all_not_water
 
         # Check energies
-        alchem_region = AlchemicalRegion(alchemical_atoms=setup_results["alchem_indices"])
+        alchemical_regions = []
+        for key in setup_results["alchemical_indices"]:
+            alchemical_regions.append(
+                AlchemicalRegion(
+                    alchemical_atoms=setup_results["alchemical_indices"][key],
+                    name=key,
+                )
+            )
+
         self._test_energies(
             reference_system=setup_results["standard_system"],
             alchemical_system=setup_results["alchem_system"],
-            alchemical_regions=alchem_region,
+            alchemical_regions=alchemical_regions,
             positions=setup_results["debug_positions"],
         )
 
@@ -522,6 +531,7 @@ def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings,
             system=setup_results["alchem_system"],
             positions=setup_results["debug_positions"],
             selection_indices=setup_results["selection_indices"],
+            alchemical_indices=setup_results["alchemical_indices"],
             box_vectors=setup_results["box_vectors"],
             alchemical_restraints=False,
             dry=True,
@@ -544,7 +554,7 @@ def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings,
 
         # Check the alchemical indices
         expected_indices = [i for i in range(self.num_ligand_atoms)]
-        assert expected_indices == setup_results["alchem_indices"]
+        assert expected_indices == setup_results["alchemical_indices"]["A"]
 
         # Check the non-alchemical system
         assert (
@@ -567,12 +577,19 @@ def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings,
         assert pdb.n_atoms == self.num_ligand_atoms
 
         # Check energies
-        alchem_region = AlchemicalRegion(alchemical_atoms=setup_results["alchem_indices"])
+        alchemical_regions = []
+        for key in setup_results["alchemical_indices"]:
+            alchemical_regions.append(
+                AlchemicalRegion(
+                    alchemical_atoms=setup_results["alchemical_indices"][key],
+                    name=key,
+                )
+            )
 
         self._test_energies(
             reference_system=setup_results["standard_system"],
             alchemical_system=setup_results["alchem_system"],
-            alchemical_regions=alchem_region,
+            alchemical_regions=alchemical_regions,
             positions=setup_results["debug_positions"],
         )
 
@@ -614,6 +631,299 @@ def settings(self):
         return s
 
 
+class TestHostGuestCharged(TestT4LysozymeDryRun):
+    """
+    A host guest system with user defined Boresch restraints
+    and a net charge.
+    """
+    solvent = SolventComponent(ion_concentration=0.15 * offunit.molar)
+    num_all_not_water = 223
+    num_host_component_atoms = 192
+    # No ions
+    num_ligand_atoms = 75
+    expected_complex_particles = 6166
+    expected_ligand_solvent_particles = 29910
+
+    barostat_by_phase = {
+        "complex": MonteCarloBarostat,
+        "solvent": MonteCarloBarostat,
+    }
+
+    @pytest.fixture(scope="class")
+    def settings(self):
+        s = openmm_afe.AbsoluteBindingProtocol.default_settings()
+        s.protocol_repeats = 1
+        s.engine_settings.compute_platform = "cpu"
+        s.forcefield_settings.nonbonded_cutoff = 1.2 * offunit.nanometer
+        s.complex_output_settings.output_indices = "not water"
+        s.complex_solvation_settings.solvent_padding = None
+        s.complex_solvation_settings.number_of_solvent_molecules = 2000
+        s.solvent_solvation_settings.box_shape = "cube"
+        s.solvent_solvation_settings.solvent_padding = None
+        s.solvent_solvation_settings.number_of_solvent_molecules = 10000
+        s.restraint_settings.guest_restraint_ids = [1, 2, 4]
+        s.restraint_settings.host_restraint_ids = [1, 4, 3]
+        s.alchemical_settings.alchemical_ion_min_distance = 1.0 * offunit.nanometer
+        return s
+
+
+    @pytest.fixture(scope="class")
+    def dag(self, protocol, OA_guests_charged, OA_host_charged):
+        stateA = ChemicalSystem(
+            {
+                "ligand": OA_guests_charged["OA-G0"],
+                "host": OA_host_charged,
+                "solvent": self.solvent,
+            }
+        )
+
+        stateB = ChemicalSystem(
+            {
+                "host": OA_host_charged,
+                "solvent": self.solvent,
+            }
+        )
+
+        return protocol.create(
+            stateA=stateA,
+            stateB=stateB,
+            mapping=None,
+        )
+
+    def _assert_expected_alchemical_forces(self, system, phase: str, settings):
+        """
+        Assert the forces expected in the alchemical system.
+        """
+        barostat_type = self.barostat_by_phase[phase]
+        self._assert_force_num(system, NonbondedForce, 1)
+        self._assert_force_num(system, CustomNonbondedForce, 4)
+        self._assert_force_num(system, CustomBondForce, 4)
+        self._assert_force_num(system, HarmonicAngleForce, 1)
+        self._assert_force_num(system, PeriodicTorsionForce, 1)
+        self._assert_force_num(system, barostat_type, 1)
+
+        if phase == "complex":
+            self._assert_force_num(system, HarmonicBondForce, 1)
+            self._assert_force_num(system, CustomCompoundBondForce, 1)
+        else:
+            self._assert_force_num(system, HarmonicBondForce, 2)
+        
+        assert len(system.getForces()) == 14
+
+        # Check the nonbonded force has the right contents
+        nonbond = [f for f in system.getForces() if isinstance(f, NonbondedForce)]
+        assert len(nonbond) == 1
+        assert nonbond[0].getNonbondedMethod() == NonbondedForce.PME
+        assert (
+            from_openmm(nonbond[0].getCutoffDistance())
+            == settings.forcefield_settings.nonbonded_cutoff
+        )
+
+        # Check the barostat made it all the way through
+        barostat = [f for f in system.getForces() if isinstance(f, barostat_type)]
+        assert len(barostat) == 1
+        expected_frequency = int(
+            (
+                settings.complex_integrator_settings
+                if phase == "complex"
+                else settings.solvent_integrator_settings
+            ).barostat_frequency.m
+        )
+        assert barostat[0].getFrequency() == expected_frequency
+        assert barostat[0].getDefaultPressure() == to_openmm(settings.thermo_settings.pressure)
+        assert barostat[0].getDefaultTemperature() == to_openmm(
+            settings.thermo_settings.temperature
+        )
+
+    def test_complex_dry_run(self, complex_setup_units, complex_sim_units, settings, tmp_path):
+        setup_results = complex_setup_units[0].run(
+            dry=True,
+            verbose=True,
+            scratch_basepath=tmp_path,
+            shared_basepath=tmp_path,
+        )
+        sim_results = complex_sim_units[0].run(
+            system=setup_results["alchem_system"],
+            positions=setup_results["debug_positions"],
+            selection_indices=setup_results["selection_indices"],
+            alchemical_indices=setup_results["alchemical_indices"],
+            box_vectors=setup_results["box_vectors"],
+            alchemical_restraints=True,
+            dry=True,
+            scratch_basepath=tmp_path,
+            shared_basepath=tmp_path,
+        )
+
+        # Check the sampler
+        self._verify_sampler(sim_results["sampler"], phase="complex", settings=settings)
+
+        # Check the alchemical system
+        assert setup_results["alchem_system"].getNumParticles() == self.expected_complex_particles
+        self._assert_expected_alchemical_forces(
+            setup_results["alchem_system"], phase="complex", settings=settings
+        )
+        self._check_box_vectors(setup_results["alchem_system"])
+
+        # # Check the alchemical indices
+        # expected_indices = [
+        #     i + self.num_host_component_atoms - 1 for i in range(self.num_ligand_atoms)
+        # ]
+        # assert expected_indices == setup_results["alchemical_indices"]["A"]
+
+        # Check the non-alchemical system
+        assert setup_results["standard_system"].getNumParticles() == self.expected_complex_particles
+        self._assert_expected_nonalchemical_forces(
+            setup_results["standard_system"],
+            "complex",
+            settings=settings,
+        )
+        self._check_box_vectors(setup_results["standard_system"])
+        # Check the box vectors haven't changed (they shouldn't have because we didn't do MD)
+        assert_allclose(
+            from_openmm(setup_results["alchem_system"].getDefaultPeriodicBoxVectors()),
+            from_openmm(setup_results["standard_system"].getDefaultPeriodicBoxVectors()),
+        )
+
+        # Check the PDB
+        pdb = mdt.load_pdb(setup_results["pdb_structure"])
+        assert pdb.n_atoms == self.num_all_not_water
+
+        # Check energies
+        alchemical_regions = []
+        for key in setup_results["alchemical_indices"]:
+            alchemical_regions.append(
+                AlchemicalRegion(
+                    alchemical_atoms=setup_results["alchemical_indices"][key],
+                    name=key,
+                )
+            )
+
+        self._test_energies(
+            reference_system=setup_results["standard_system"],
+            alchemical_system=setup_results["alchem_system"],
+            alchemical_regions=alchemical_regions,
+            positions=setup_results["debug_positions"],
+        )
+
+    def test_solvent_dry_run(self, solvent_setup_units, solvent_sim_units, settings, tmp_path):
+        setup_results = solvent_setup_units[0].run(
+            dry=True,
+            verbose=True,
+            scratch_basepath=tmp_path,
+            shared_basepath=tmp_path,
+        )
+        sim_results = solvent_sim_units[0].run(
+            system=setup_results["alchem_system"],
+            positions=setup_results["debug_positions"],
+            selection_indices=setup_results["selection_indices"],
+            alchemical_indices=setup_results["alchemical_indices"],
+            box_vectors=setup_results["box_vectors"],
+            alchemical_restraints=False,
+            dry=True,
+            scratch_basepath=tmp_path,
+            shared_basepath=tmp_path,
+        )
+
+        # Check the sampler
+        self._verify_sampler(sim_results["sampler"], phase="solvent", settings=settings)
+
+        # Check the alchemical system
+        assert (
+            setup_results["alchem_system"].getNumParticles()
+            == self.expected_ligand_solvent_particles
+        )
+        self._assert_expected_alchemical_forces(
+            setup_results["alchem_system"], phase="solvent", settings=settings
+        )
+        self._test_cubic_vectors(setup_results["alchem_system"])
+
+        # # Check the alchemical indices
+        # expected_indices = [i for i in range(self.num_ligand_atoms)]
+        # assert expected_indices == setup_results["alchemical_indices"]["A"]
+
+        # Check the non-alchemical system
+        assert (
+            setup_results["standard_system"].getNumParticles()
+            == self.expected_ligand_solvent_particles
+        )
+        self._assert_expected_nonalchemical_forces(
+            setup_results["standard_system"], phase="solvent", settings=settings
+        )
+        self._test_cubic_vectors(setup_results["standard_system"])
+
+        # Check the box vectors haven't changed (they shouldn't have because we didn't do MD)
+        assert_allclose(
+            from_openmm(setup_results["alchem_system"].getDefaultPeriodicBoxVectors()),
+            from_openmm(setup_results["standard_system"].getDefaultPeriodicBoxVectors()),
+        )
+
+        # Check the PDB
+        pdb = mdt.load_pdb(setup_results["pdb_structure"])
+        assert pdb.n_atoms == self.num_ligand_atoms
+
+        # Check energies
+        alchemical_regions = []
+        for key in setup_results["alchemical_indices"]:
+            alchemical_regions.append(
+                AlchemicalRegion(
+                    alchemical_atoms=setup_results["alchemical_indices"][key],
+                    name=key,
+                )
+            )
+
+        self._test_energies(
+            reference_system=setup_results["standard_system"],
+            alchemical_system=setup_results["alchem_system"],
+            alchemical_regions=alchemical_regions,
+            positions=setup_results["debug_positions"],
+        )
+
+
+def test_user_restraint(benzene_modifications_am1bcc, T4_protein_component, tmp_path):
+    s = openmm_afe.AbsoluteBindingProtocol.default_settings()
+    s.protocol_repeats = 1
+    s.engine_settings.compute_platform = "cpu"
+    s.restraint_settings.guest_restraint_ids = [0, 1, 2]
+    # Ca and C from VAL 87, and N from TYR 88
+    s.restraint_settings.host_restraint_ids = [1383, 1384, 1398]
+
+    protocol = openmm_afe.AbsoluteBindingProtocol(settings=s)
+
+    stateA = gufe.ChemicalSystem(
+        {
+            "protein": T4_protein_component,
+            "benzene": benzene_modifications_am1bcc["benzene"],
+            "solvent": gufe.SolventComponent(),
+        }
+    )
+
+    stateB = gufe.ChemicalSystem(
+        {
+            "protein": T4_protein_component,
+            "solvent": gufe.SolventComponent(),
+        }
+    )
+
+    dag = protocol.create(stateA=stateA, stateB=stateB, mapping=None)
+
+    complex_setup_units = _get_units(dag.protocol_units, UNIT_TYPES["complex"]["setup"])
+
+    results = complex_setup_units[0].run(
+        dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path
+    )
+
+    geom = BoreschRestraintGeometry.model_validate(results["restraint_geometry"])
+    # This should be C1, C2, and C3 on the benzene
+    assert geom.guest_atoms == [2613, 2614, 2615]
+    assert geom.host_atoms == [1383, 1384, 1398]
+    assert pytest.approx(geom.r_aA0.to("nanometer").m, rel=1e-4) == 0.510798
+    assert pytest.approx(geom.theta_A0.to("radians").m, rel=1e-4) == 1.20278
+    assert pytest.approx(geom.theta_B0.to("radians").m, rel=1e-4) == 1.25705
+    assert pytest.approx(geom.phi_A0.to("radians").m, rel=1e-4) == 0.86035
+    assert pytest.approx(geom.phi_B0.to("radians").m, rel=1e-4) == 1.59444
+    assert pytest.approx(geom.phi_C0.to("radians").m, rel=1e-4) == 2.92365
+
+
 def test_user_charges(benzene_modifications, T4_protein_component, tmp_path):
     s = openmm_afe.AbsoluteBindingProtocol.default_settings()
     s.protocol_repeats = 1
@@ -682,6 +992,12 @@ def assign_fictitious_charges(offmol):
 
         c, s, e = system_nbf.getParticleParameters(index)
         assert pytest.approx(prop_chgs[i]) == c.value_in_unit(ommunit.elementary_charge)
+
+        # Alchemical system should be 0 charge in standard parameters
+        # and all charge in the offset
+        c, s, e = alchem_system_nbf.getParticleParameters(index)
+        assert pytest.approx(0.0) == c.value_in_unit(ommunit.elementary_charge)
+
         offsets = alchem_system_nbf.getParticleParameterOffset(i)
         assert pytest.approx(prop_chgs[i]) == offsets[2]
 
@@ -699,7 +1015,7 @@ class TestA2AMembraneDryRun(TestT4LysozymeDryRun):
 
     solvent = SolventComponent(ion_concentration=0 * offunit.molar)
     num_all_not_water = 16080
-    num_protein_component_atoms = 39391
+    num_host_component_atoms = 39391
     expected_complex_particles = 39426
     expected_ligand_solvent_particles = 3036
 
diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py
index 2c63ea819..9586cc8c7 100644
--- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py
+++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_protocol_results.py
@@ -29,6 +29,7 @@ def patcher():
                 "positions": Path("positions.npy"),
                 "pdb_structure": Path("hybrid_system.pdb"),
                 "selection_indices": np.zeros(100),
+                "alchemical_indices": {"A": list(range(10))},
                 "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm,
                 "standard_state_correction": 0 * offunit.kilocalorie_per_mole,
                 "restraint_geometry": None,
@@ -44,6 +45,7 @@ def patcher():
                 "positions": Path("positions.npy"),
                 "pdb_structure": Path("hybrid_system.pdb"),
                 "selection_indices": np.zeros(100),
+                "alchemical_indices": {"A": list(range(10))},
                 "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm,
                 "standard_state_correction": 0 * offunit.kilocalorie_per_mole,
                 "restraint_geometry": True,
diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py
index fad97e6ad..013c775ff 100644
--- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py
+++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_settings.py
@@ -7,6 +7,7 @@
     AbsoluteBindingProtocol,
     AbsoluteBindingSettings,
 )
+from openfe.protocols.openmm_afe.equil_afe_settings import ABFEBoreschRestraintSettings
 
 
 @pytest.fixture()
@@ -70,6 +71,15 @@ def test_monotonic_lambda_windows(val, default_settings):
         lambda_settings.lambda_restraints = val["restraints"]
 
 
+@pytest.mark.parametrize("parameter", ["host_restraint_ids", "guest_restraint_ids"])
+def test_boresch_restraints_restraint_negative_ids(parameter):
+    setting = ABFEBoreschRestraintSettings()
+
+    errmsg = "``guest_atoms`` and ``host_atoms`` cannot have negative indices."
+    with pytest.raises(ValueError, match=errmsg):
+        setattr(setting, parameter, [1, 2, -3])
+
+
 def test_equil_not_all_complex(default_settings):
     with pytest.raises(ValueError, match="output_indices must be all"):
         default_settings.complex_equil_output_settings.output_indices = "not water"
diff --git a/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py
index 378702434..eb2a4d79d 100644
--- a/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py
+++ b/src/openfe/tests/protocols/openmm_abfe/test_abfe_validation.py
@@ -121,14 +121,37 @@ def test_validate_no_protcomp(
 
     stateB = ChemicalSystem(
         {
-            "benzene": benzene_modifications["benzene"],
             "solvent": SolventComponent(),
         }
     )
 
-    errmsg = "No ProteinComponent found"
+    errmsg = "No suitable host molecule found"
     with pytest.raises(ValueError, match=errmsg):
-        AbsoluteBindingProtocol._validate_endstates(stateA, stateB)
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
+
+
+def test_validate_smc_host(
+    benzene_modifications,
+):
+    """
+    Should pass if there's another smc present that could be a host.
+    """
+    stateA = ChemicalSystem(
+        {
+            "benzene": benzene_modifications["benzene"],
+            "toluene": benzene_modifications["toluene"],
+            "solvent": SolventComponent(),
+        }
+    )
+
+    stateB = ChemicalSystem(
+        {
+            "toluene": benzene_modifications["toluene"],
+            "solvent": SolventComponent(),
+        }
+    )
+
+    AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
 
 
 def test_validate_endstates_nosolvcomp_stateA(benzene_modifications, T4_protein_component):
@@ -148,7 +171,7 @@ def test_validate_endstates_nosolvcomp_stateA(benzene_modifications, T4_protein_
     )
 
     with pytest.raises(ValueError, match="No SolventComponent found"):
-        AbsoluteBindingProtocol._validate_endstates(stateA, stateB)
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
 
 
 def test_validate_endstates_nosolvcomp_stateB(benzene_modifications, T4_protein_component):
@@ -168,7 +191,7 @@ def test_validate_endstates_nosolvcomp_stateB(benzene_modifications, T4_protein_
     )
 
     with pytest.raises(ValueError, match="No SolventComponent found"):
-        AbsoluteBindingProtocol._validate_endstates(stateA, stateB)
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
 
 
 def test_validate_endstates_multiple_uniqueA(benzene_modifications, T4_protein_component):
@@ -189,7 +212,7 @@ def test_validate_endstates_multiple_uniqueA(benzene_modifications, T4_protein_c
     )
 
     with pytest.raises(ValueError, match="Only one alchemical species"):
-        AbsoluteBindingProtocol._validate_endstates(stateA, stateB)
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
 
 
 def test_validate_solvent_endstates_solvent_dissapearing(
@@ -214,7 +237,7 @@ def test_validate_solvent_endstates_solvent_dissapearing(
 
     errmsg = "Only disappearing small molecule components"
     with pytest.raises(ValueError, match=errmsg):
-        AbsoluteBindingProtocol._validate_endstates(stateA, stateB)
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
 
 
 def test_validate_endstates_unique_stateB(benzene_modifications, T4_protein_component):
@@ -235,13 +258,13 @@ def test_validate_endstates_unique_stateB(benzene_modifications, T4_protein_comp
     )
 
     with pytest.raises(ValueError, match="appearing in state B"):
-        AbsoluteBindingProtocol._validate_endstates(stateA, stateB)
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
 
 
-def test_charged_endstate(charged_benzene_modifications, T4_protein_component):
+def test_highly_charged_endstate_error(charged_benzene_modifications, T4_protein_component):
     stateA = ChemicalSystem(
         {
-            "benzene": charged_benzene_modifications["benzoic_acid"],
+            "benzene": charged_benzene_modifications["phthalic_acid"],
             "protein": T4_protein_component,
             "solvent": SolventComponent(),
         }
@@ -254,9 +277,30 @@ def test_charged_endstate(charged_benzene_modifications, T4_protein_component):
         }
     )
 
-    errmsg = "Charged alchemical molecules are not currently supported"
+    errmsg = "Alchemical molecules with a formal charge of greater than 1"
     with pytest.raises(ValueError, match=errmsg):
-        AbsoluteBindingProtocol._validate_endstates(stateA, stateB)
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, True)
+
+
+def test_charge_endstate_warning(charged_benzene_modifications, T4_protein_component):
+    stateA = ChemicalSystem(
+        {
+            "benzene": charged_benzene_modifications["benzoic_acid"],
+            "protein": T4_protein_component,
+            "solvent": SolventComponent(),
+        }
+    )
+
+    stateB = ChemicalSystem(
+        {
+            "protein": T4_protein_component,
+            "solvent": SolventComponent(),
+        }
+    )
+
+    msg = "Ligand has a net charge but no charge correction scheme has been requested."
+    with pytest.warns(UserWarning, match=msg):
+        AbsoluteBindingProtocol._validate_endstates(stateA, stateB, False)
 
 
 def test_validate_fail_extends(benzene_modifications, T4_protein_component, default_settings):
diff --git a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py
index a5193a530..2a047ead0 100644
--- a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py
+++ b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol.py
@@ -230,6 +230,7 @@ def test_setup_dry_sim_vac_benzene(benzene_system, method, protocol_dry_settings
         selection_indices=setup_results["selection_indices"],
         box_vectors=setup_results["box_vectors"],
         alchemical_restraints=False,
+        alchemical_indices=setup_results["alchemical_indices"],
         dry=True,
     )
 
@@ -398,6 +399,7 @@ def test_setup_solv_benzene(benzene_system, protocol_dry_settings, tmp_path):
         selection_indices=setup_results["selection_indices"],
         box_vectors=setup_results["box_vectors"],
         alchemical_restraints=False,
+        alchemical_indices=setup_results["alchemical_indices"],
         dry=True,
         scratch_basepath=tmp_path,
         shared_basepath=tmp_path,
@@ -451,6 +453,7 @@ def test_dry_run_vsite_fail(benzene_system, tmp_path, protocol_dry_settings):
             selection_indices=setup_results["selection_indices"],
             box_vectors=setup_results["box_vectors"],
             alchemical_restraints=False,
+            alchemical_indices=setup_results["alchemical_indices"],
             dry=True,
             scratch_basepath=tmp_path,
             shared_basepath=tmp_path,
@@ -498,6 +501,7 @@ def test_setup_dry_sim_solv_benzene_tip4p(benzene_system, protocol_dry_settings,
         selection_indices=setup_results["selection_indices"],
         box_vectors=setup_results["box_vectors"],
         alchemical_restraints=False,
+        alchemical_indices=setup_results["alchemical_indices"],
         dry=True,
         scratch_basepath=tmp_path,
         shared_basepath=tmp_path,
@@ -608,7 +612,9 @@ def assign_fictitious_charges(offmol):
     nonbond = [f for f in system.getForces() if isinstance(f, CustomNonbondedForce)]
     assert len(nonbond) == 4
 
-    custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics"][0]
+    custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][
+        0
+    ]
 
     # loop through the 12 benzene atoms
     for i in range(12):
@@ -675,7 +681,9 @@ def test_dry_run_charge_backends(
     nonbond = [f for f in system.getForces() if isinstance(f, CustomNonbondedForce)]
     assert len(nonbond) == 4
 
-    custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics"][0]
+    custom_elec = [n for n in nonbond if n.getGlobalParameterName(0) == "lambda_electrostatics_A"][
+        0
+    ]
 
     charges = []
     for i in range(system.getNumParticles()):
@@ -756,6 +764,7 @@ def test_dry_run_vacuum_write_frequency(
             selection_indices=setup_results["selection_indices"],
             box_vectors=setup_results["box_vectors"],
             alchemical_restraints=False,
+            alchemical_indices=setup_results["alchemical_indices"],
             dry=True,
             scratch_basepath=tmp_path,
             shared_basepath=tmp_path,
diff --git a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py
index c8ca84020..87cd90aa1 100644
--- a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py
+++ b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_protocol_results.py
@@ -52,6 +52,7 @@ def patcher():
                 "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm,
                 "standard_state_correction": 0 * offunit.kilocalorie_per_mole,
                 "restraint_geometry": None,
+                "alchemical_indices": [1, 2, 3, 4, 5, 6],
                 "gufe_version": gufe.__version__,
                 "openfe_version": openfe.__version__,
                 "openmm_version": openmm.__version__,
@@ -64,6 +65,7 @@ def patcher():
                 "positions": Path("positions.npy"),
                 "pdb_structure": Path("hybrid_system.pdb"),
                 "selection_indices": np.zeros(100),
+                "alchemical_indices": [1, 2, 3, 4, 5, 6],
                 "box_vectors": [np.zeros(3), np.zeros(3), np.zeros(3)] * offunit.nm,
                 "standard_state_correction": 0 * offunit.kilocalorie_per_mole,
                 "restraint_geometry": None,
diff --git a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py
index e62767f3d..2ebd8cbbe 100644
--- a/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py
+++ b/src/openfe/tests/protocols/openmm_ahfe/test_ahfe_resume.py
@@ -228,6 +228,7 @@ def test_resume(
             selection_indices=setup_results["selection_indices"],
             box_vectors=setup_results["box_vectors"],
             alchemical_restraints=False,
+            alchemical_indices=setup_results["alchemical_indices"],
             scratch_basepath=tmp_path,
             shared_basepath=tmp_path,
         )
@@ -303,6 +304,7 @@ def test_resume_fail_particles(
                 selection_indices=setup_results["selection_indices"],
                 box_vectors=setup_results["box_vectors"],
                 alchemical_restraints=False,
+                alchemical_indices=setup_results["alchemical_indices"],
                 scratch_basepath=tmp_path,
                 shared_basepath=tmp_path,
             )
@@ -345,6 +347,7 @@ def test_resume_fail_constraints(
                 selection_indices=setup_results["selection_indices"],
                 box_vectors=setup_results["box_vectors"],
                 alchemical_restraints=False,
+                alchemical_indices=setup_results["alchemical_indices"],
                 scratch_basepath=tmp_path,
                 shared_basepath=tmp_path,
             )
@@ -385,6 +388,7 @@ def test_resume_fail_forces(
                 selection_indices=setup_results["selection_indices"],
                 box_vectors=setup_results["box_vectors"],
                 alchemical_restraints=False,
+                alchemical_indices=setup_results["alchemical_indices"],
                 scratch_basepath=tmp_path,
                 shared_basepath=tmp_path,
             )
@@ -441,6 +445,7 @@ def test_resume_differ_barostat(
                 selection_indices=setup_results["selection_indices"],
                 box_vectors=setup_results["box_vectors"],
                 alchemical_restraints=False,
+                alchemical_indices=setup_results["alchemical_indices"],
                 scratch_basepath=tmp_path,
                 shared_basepath=tmp_path,
             )
@@ -486,7 +491,7 @@ def test_resume_differ_forces(
         # Now add a fake force
         new_force = openmm.NonbondedForce()
         new_force.setNonbondedMethod(openmm.NonbondedForce.PME)
-        new_force.addGlobalParameter("lambda_electrostatics", 1.0)
+        new_force.addGlobalParameter("lambda_electrostatics_A", 1.0)
 
         fake_system.addForce(new_force)
 
@@ -500,6 +505,7 @@ def test_resume_differ_forces(
             selection_indices=setup_results["selection_indices"],
             box_vectors=setup_results["box_vectors"],
             alchemical_restraints=False,
+            alchemical_indices=setup_results["alchemical_indices"],
             scratch_basepath=tmp_path,
             shared_basepath=tmp_path,
             dry=True,
@@ -546,6 +552,7 @@ def test_resume_bad_files(
                 selection_indices=setup_results["selection_indices"],
                 box_vectors=setup_results["box_vectors"],
                 alchemical_restraints=False,
+                alchemical_indices=setup_results["alchemical_indices"],
                 scratch_basepath=tmp_path,
                 shared_basepath=tmp_path,
             )
diff --git a/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py b/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py
index c003432a0..022b57b6d 100644
--- a/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py
+++ b/src/openfe/tests/protocols/openmm_septop/test_septop_slow.py
@@ -18,7 +18,7 @@
     SepTopProtocol,
     SepTopSolventSetupUnit,
 )
-from openfe.protocols.openmm_septop.utils import SepTopParameterState
+from openfe.protocols.openmm_utils.states import DualRegionAlchemicalState
 from openfe.protocols.openmm_utils.serialization import deserialize
 
 
@@ -30,7 +30,7 @@ def default_settings():
 
 
 def compare_energies(alchemical_system, positions):
-    alchemical_state = SepTopParameterState.from_system(alchemical_system)
+    alchemical_state = DualRegionAlchemicalState.from_system(alchemical_system)
 
     from openmmtools.alchemy import AbsoluteAlchemicalFactory