From e7882f0ca650ea84853436b91a398f362f32ae30 Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Thu, 25 Jun 2026 15:58:46 +0200
Subject: [PATCH 01/19] Initial change of lig and cofactor resname

---
 .../protocols/openmm_rfe/hybridtop_units.py   |  14 ++-
 .../protocols/openmm_utils/molecule_utils.py  | 109 ++++++++++++++++++
 .../openmm_rfe/test_hybrid_top_protocol.py    |  31 +++++
 3 files changed, 153 insertions(+), 1 deletion(-)
 create mode 100644 src/openfe/protocols/openmm_utils/molecule_utils.py

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 5f4c1c9f3..03fbe4e9c 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -49,6 +49,7 @@
 from openfe.protocols.openmm_utils.omm_settings import (
     BasePartialChargeSettings,
 )
+from openfe.protocols.openmm_utils.molecule_utils import _get_offmol_resname, _set_offmol_resname
 
 from ...analysis import plotting
 from ...utils import log_system_probe, without_oechem_backend
@@ -753,6 +754,17 @@ def run(
         alchem_comps = self._inputs["alchemical_components"]
         solvent_comp, protein_comp, small_mols = self._get_components(stateA, stateB)
 
+        alchemical = set(alchem_comps["stateA"]) | set(alchem_comps["stateB"])
+        cofactors = [smc for smc in small_mols if smc not in alchemical]
+
+        for smc, offmol in small_mols.items():
+            if _get_offmol_resname(offmol) is not None:
+                continue
+            if smc in alchemical:
+                _set_offmol_resname(offmol, "LIG")
+            else:
+                _set_offmol_resname(offmol, f"CF{cofactors.index(smc) + 1}")
+
         # Assign partial charges now to avoid any discrepancies later
         self._assign_partial_charges(settings["charge_settings"], small_mols)
 
@@ -1535,7 +1547,7 @@ def _structural_analysis(
         from openfe_analysis import rmsd
 
         try:
-            data = rmsd.gather_rms_data(pdb_file, trj_file)
+            data = rmsd.gather_rms_data(pdb_file, trj_file, ligand_selection="resname LIG")
         # TODO: eventually change this to more specific exception types
         except Exception as e:
             return {"structural_analysis_error": str(e)}
diff --git a/src/openfe/protocols/openmm_utils/molecule_utils.py b/src/openfe/protocols/openmm_utils/molecule_utils.py
new file mode 100644
index 000000000..aec4a21d7
--- /dev/null
+++ b/src/openfe/protocols/openmm_utils/molecule_utils.py
@@ -0,0 +1,109 @@
+# This code is part of OpenFE and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/openfe
+import logging
+from typing import Any
+
+from openff.toolkit import Molecule as OFFMolecule
+
+logger = logging.getLogger(__name__)
+
+
+def _set_offmol_metadata(
+    offmol: OFFMolecule,
+    key: Any,
+    val: Any | None,
+) -> None:
+    """
+    Set a given metadata entry for a whole Molecule.
+
+    Parameters
+    ----------
+    offmol : openff.toolkit.Molecule
+      The Molecule to set the metadata for.
+    key : Any
+      The metadata key.
+    val : Any
+      The value to set the metadata entry to.
+    """
+    if val is None:
+        for a in offmol.atoms:
+            a.metadata.pop(key, None)
+    else:
+        for a in offmol.atoms:
+            a.metadata[key] = val
+
+
+def _get_offmol_metadata(offmol: OFFMolecule, key: Any) -> Any | None:
+    """
+    Get an offmol's given metadata entry and make sure it is
+    consistent across all atoms in the Molecule.
+
+    Parameters
+    ----------
+    offmol : openff.toolkit.Molecule
+      Molecule to get the metadata value from.
+    key: Any
+      The metadata entry key.
+
+    Returns
+    -------
+    value : Any | None
+      Metadata for the given key in the molecule. ``None`` if the
+      Molecule does not have that metadata entry set, or if
+      the value is inconsistent across all the atoms.
+    """
+    value: Any | None = None
+    for a in offmol.atoms:
+        if value is None:
+            try:
+                value = a.metadata[key]
+            except KeyError:
+                return None
+
+        if value != a.metadata[key]:
+            wmsg = f"Inconsistent metadata {key} in OFFMol: {offmol}"
+            logger.warning(wmsg)
+            return None
+
+    return value
+
+
+def _set_offmol_resname(
+    offmol: OFFMolecule,
+    resname: str | None,
+) -> None:
+    """
+    Helper method to set offmol residue names
+
+    Parameters
+    ----------
+    offmol : openff.toolkit.Molecule
+      Molecule to assign a residue name to.
+    resname : str | None
+      Residue name to be set. Set to None to clear it.
+
+    Returns
+    -------
+    None
+    """
+    _set_offmol_metadata(offmol, "residue_name", resname)
+
+
+def _get_offmol_resname(offmol: OFFMolecule) -> str | None:
+    """
+    Helper method to get an offmol's residue name and make sure it is
+    consistent across all atoms in the Molecule.
+
+    Parameters
+    ----------
+    offmol : openff.toolkit.Molecule
+      Molecule to get the residue name from.
+
+    Returns
+    -------
+    resname : Optional[str]
+      Residue name of the molecule. ``None`` if the Molecule
+      does not have a residue name, or if the residue name is
+      inconsistent across all the atoms.
+    """
+    return _get_offmol_metadata(offmol, "residue_name")
\ No newline at end of file
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index a1f9fef2d..bd0ce84e4 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -11,6 +11,7 @@
 from unittest import mock
 
 import gufe
+import MDAnalysis as mda
 import mdtraj as mdt
 import numpy as np
 import openmm
@@ -2401,3 +2402,33 @@ def test_dry_run_vacuum_write_frequency(
         assert reporter.velocity_interval == velocities_write_frequency.m
     else:
         assert reporter.velocity_interval == 0
+
+
+@pytest.fixture
+def eg5_vac_inputs(eg5_ligands, eg5_cofactor):
+    ligA, ligB = eg5_ligands[0], eg5_ligands[1]
+    mapper = openfe.LomapAtomMapper()
+    mapping = next(mapper.suggest_mappings(ligA, ligB))
+    stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": eg5_cofactor})
+    stateB = openfe.ChemicalSystem({"ligand": ligB, "cofactor": eg5_cofactor})
+    return stateA, stateB, mapping
+
+
+def _run_setup_dry(stateA, stateB, mapping, settings, tmp_path):
+    protocol = openmm_rfe.RelativeHybridTopologyProtocol(settings=settings)
+    dag = protocol.create(stateA=stateA, stateB=stateB, mapping=mapping)
+    setup_unit = _get_units(dag.protocol_units, HybridTopologySetupUnit)[0]
+    return setup_unit.run(dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path)
+
+
+def test_ligand_separable_from_cofactor(eg5_vac_inputs, vac_settings, tmp_path):
+    vac_settings.output_settings.output_indices = "all"
+    stateA, stateB, mapping = eg5_vac_inputs
+    out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
+
+    u = mda.Universe(out["pdb_structure"])
+    lig = u.select_atoms("resname LIG")
+    cof = u.select_atoms("resname CF1")
+    assert lig.n_atoms > 0
+    assert cof.n_atoms > 0
+    assert set(lig.indices).isdisjoint(cof.indices)
\ No newline at end of file

From a32cf025c9757c1a10f2fb04e1694a2730a9c10c Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 29 Jun 2026 08:16:13 +0000
Subject: [PATCH 02/19] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/openfe/protocols/openmm_rfe/hybridtop_units.py              | 2 +-
 src/openfe/protocols/openmm_utils/molecule_utils.py             | 2 +-
 .../tests/protocols/openmm_rfe/test_hybrid_top_protocol.py      | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 03fbe4e9c..e92239668 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -46,10 +46,10 @@
 from openmmtools import multistate
 
 import openfe
+from openfe.protocols.openmm_utils.molecule_utils import _get_offmol_resname, _set_offmol_resname
 from openfe.protocols.openmm_utils.omm_settings import (
     BasePartialChargeSettings,
 )
-from openfe.protocols.openmm_utils.molecule_utils import _get_offmol_resname, _set_offmol_resname
 
 from ...analysis import plotting
 from ...utils import log_system_probe, without_oechem_backend
diff --git a/src/openfe/protocols/openmm_utils/molecule_utils.py b/src/openfe/protocols/openmm_utils/molecule_utils.py
index aec4a21d7..a715af4e2 100644
--- a/src/openfe/protocols/openmm_utils/molecule_utils.py
+++ b/src/openfe/protocols/openmm_utils/molecule_utils.py
@@ -106,4 +106,4 @@ def _get_offmol_resname(offmol: OFFMolecule) -> str | None:
       does not have a residue name, or if the residue name is
       inconsistent across all the atoms.
     """
-    return _get_offmol_metadata(offmol, "residue_name")
\ No newline at end of file
+    return _get_offmol_metadata(offmol, "residue_name")
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index bd0ce84e4..0f4259ae6 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2431,4 +2431,4 @@ def test_ligand_separable_from_cofactor(eg5_vac_inputs, vac_settings, tmp_path):
     cof = u.select_atoms("resname CF1")
     assert lig.n_atoms > 0
     assert cof.n_atoms > 0
-    assert set(lig.indices).isdisjoint(cof.indices)
\ No newline at end of file
+    assert set(lig.indices).isdisjoint(cof.indices)

From 369f9d5a07771e4086ad702b7557130b8441d63c Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Mon, 29 Jun 2026 11:42:02 +0200
Subject: [PATCH 03/19] Fix tests

---
 .../hybrid_topology_atoms.csv                 | 32 +++++++++----------
 .../openmm_rfe/test_hybrid_top_protocol.py    |  4 +--
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv b/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv
index 3ed6ee0e8..86ed7f043 100644
--- a/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv
+++ b/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv
@@ -1,17 +1,17 @@
 ,serial,name,element,resSeq,resName,chainID,segmentID,formal_charge
-0,,C1x,C,0,UNK,0,,
-1,,C2x,C,0,UNK,0,,
-2,,C3x,C,0,UNK,0,,
-3,,C4x,C,0,UNK,0,,
-4,,C5x,C,0,UNK,0,,
-5,,C6x,C,0,UNK,0,,
-6,,H1x,H,0,UNK,0,,
-7,,H2x,H,0,UNK,0,,
-8,,H3x,H,0,UNK,0,,
-9,,H4x,H,0,UNK,0,,
-10,,H5x,H,0,UNK,0,,
-11,,H6x,H,0,UNK,0,,
-12,,H1x,H,0,UNK,0,,
-13,,H2x,H,0,UNK,0,,
-14,,C1x,C,0,UNK,0,,
-15,,H3x,H,0,UNK,0,,
+0,,C1x,C,0,LIG,0,,
+1,,C2x,C,0,LIG,0,,
+2,,C3x,C,0,LIG,0,,
+3,,C4x,C,0,LIG,0,,
+4,,C5x,C,0,LIG,0,,
+5,,C6x,C,0,LIG,0,,
+6,,H1x,H,0,LIG,0,,
+7,,H2x,H,0,LIG,0,,
+8,,H3x,H,0,LIG,0,,
+9,,H4x,H,0,LIG,0,,
+10,,H5x,H,0,LIG,0,,
+11,,H6x,H,0,LIG,0,,
+12,,H1x,H,0,LIG,0,,
+13,,H2x,H,0,LIG,0,,
+14,,C1x,C,0,LIG,0,,
+15,,H3x,H,0,LIG,0,,
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 0f4259ae6..6d8e110a4 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -508,7 +508,7 @@ def test_dry_run_ligand(
 ):
     # this might be a bit time consuming
     solv_settings.simulation_settings.sampler_method = method
-    solv_settings.output_settings.output_indices = "resname UNK"
+    solv_settings.output_settings.output_indices = "resname LIG"
 
     protocol = openmm_rfe.RelativeHybridTopologyProtocol(
         settings=solv_settings,
@@ -1092,7 +1092,7 @@ def test_dry_run_complex(
 ):
     # this will be very time consuming
     solv_settings.simulation_settings.sampler_method = method
-    solv_settings.output_settings.output_indices = "protein or resname  UNK"
+    solv_settings.output_settings.output_indices = "protein or resname  LIG"
 
     protocol = openmm_rfe.RelativeHybridTopologyProtocol(
         settings=solv_settings,

From 69a6a0807fe9fbbc5d1cfc023208b2eb73b50cdd Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Mon, 29 Jun 2026 13:48:46 +0200
Subject: [PATCH 04/19] Add news entry

---
 news/fix_unk_resnames.rst | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100644 news/fix_unk_resnames.rst

diff --git a/news/fix_unk_resnames.rst b/news/fix_unk_resnames.rst
new file mode 100644
index 000000000..8f252ccef
--- /dev/null
+++ b/news/fix_unk_resnames.rst
@@ -0,0 +1,23 @@
+**Added:**
+
+* <news item>
+
+**Changed:**
+
+* Small molecules in RelativeHybridTopology topologies (including the output PDB) are now named LIG (alchemical ligand) and CF1, CF2… (cofactors) instead of UNK.
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* Fixed inflated ligand RMSD in the RelativeHybridTopology protocol's structural analysis for systems containing cofactors; the ligand RMSD is now computed for the alchemical ligand alone rather than conflating it with cofactors that shared the UNK residue name.
+
+**Security:**
+
+* <news item>

From 2609bc4fcbe5f117e9de2874b561dcab6e6d9d4d Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Mon, 29 Jun 2026 14:06:11 +0200
Subject: [PATCH 05/19] Add more tests

---
 .../openmm_rfe/test_hybrid_top_protocol.py    | 47 +++++++++++++++++++
 1 file changed, 47 insertions(+)

diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 6d8e110a4..c03cd71c4 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -51,6 +51,12 @@
     HAS_NAGL,
     HAS_OPENEYE,
 )
+from openfe.protocols.openmm_utils.molecule_utils import (
+    _get_offmol_metadata,
+    _set_offmol_metadata,
+    _set_offmol_resname,
+    _get_offmol_resname,
+)
 
 
 def _get_units(protocol_units, unit_type):
@@ -2432,3 +2438,44 @@ def test_ligand_separable_from_cofactor(eg5_vac_inputs, vac_settings, tmp_path):
     assert lig.n_atoms > 0
     assert cof.n_atoms > 0
     assert set(lig.indices).isdisjoint(cof.indices)
+
+
+def test_get_metadata_inconsistent_warns(caplog):
+    mol = Molecule.from_smiles("CC")
+    _set_offmol_metadata(mol, "residue_name", "LIG")
+    mol.atoms[0].metadata["residue_name"] = "COF"
+
+    with caplog.at_level(logging.WARNING):
+        result = _get_offmol_metadata(mol, "residue_name")
+
+    assert result is None
+    assert "Inconsistent metadata" in caplog.text
+
+
+def test_set_metadata_none_clears():
+    mol = Molecule.from_smiles("CC")
+    _set_offmol_metadata(mol, "residue_name", "LIG")
+    _set_offmol_metadata(mol, "residue_name", None)
+    assert all("residue_name" not in a.metadata for a in mol.atoms)
+
+
+def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+    vac_settings.output_settings.output_indices = "all"
+
+    ligA, ligB = eg5_ligands[0], eg5_ligands[1]
+    off = eg5_cofactor.to_openff()
+    _set_offmol_resname(off, "MYC")
+    cof = openfe.SmallMoleculeComponent.from_openff(off)
+
+    mapper = openfe.setup.KartografAtomMapper()
+    mapping = next(mapper.suggest_mappings(ligA, ligB))
+    stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": cof})
+    stateB = openfe.ChemicalSystem({"ligand": ligB, "cofactor": cof})
+
+    out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
+
+    u = mda.Universe(out["pdb_structure"])
+    resnames = {r.resname for r in u.residues}
+    assert "MYC" in resnames
+    assert "CF1" not in resnames
+    assert "LIG" in resnames
\ No newline at end of file

From 0844937cbf1da39371a2fa2274cf39bfd5f67396 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 29 Jun 2026 12:12:56 +0000
Subject: [PATCH 06/19] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 .../tests/protocols/openmm_rfe/test_hybrid_top_protocol.py    | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index c03cd71c4..9a5dec58d 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -53,9 +53,9 @@
 )
 from openfe.protocols.openmm_utils.molecule_utils import (
     _get_offmol_metadata,
+    _get_offmol_resname,
     _set_offmol_metadata,
     _set_offmol_resname,
-    _get_offmol_resname,
 )
 
 
@@ -2478,4 +2478,4 @@ def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp
     resnames = {r.resname for r in u.residues}
     assert "MYC" in resnames
     assert "CF1" not in resnames
-    assert "LIG" in resnames
\ No newline at end of file
+    assert "LIG" in resnames

From 11b9c2c9dcbfb2d213799a2d7713a291f040fac9 Mon Sep 17 00:00:00 2001
From: Hannah Baumann <43765638+hannahbaumann@users.noreply.github.com>
Date: Thu, 2 Jul 2026 11:06:40 +0200
Subject: [PATCH 07/19] Update
 src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py

Co-authored-by: Irfan Alibay <IAlibay@users.noreply.github.com>
---
 .../tests/protocols/openmm_rfe/test_hybrid_top_protocol.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 9a5dec58d..40857e73d 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2475,7 +2475,6 @@ def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
     u = mda.Universe(out["pdb_structure"])
-    resnames = {r.resname for r in u.residues}
-    assert "MYC" in resnames
-    assert "CF1" not in resnames
-    assert "LIG" in resnames
+    assert "MYC" in u.residues.resnames
+    assert "CF1" not in u.residues.resnames
+    assert "LIG" in u.residues.resnames

From 12b9840868838b7cda1e2cf370ae4d5e8280bde6 Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Thu, 2 Jul 2026 11:52:39 +0200
Subject: [PATCH 08/19] Address review comments

---
 news/fix_unk_resnames.rst                     |  2 +-
 .../protocols/openmm_rfe/hybridtop_units.py   | 17 +++-
 ...molecule_utils.py => offmolecule_utils.py} |  0
 .../openmm_rfe/test_hybrid_top_protocol.py    | 80 ++++++++++++++-----
 .../tests/protocols/test_openmmutils.py       | 27 ++++++-
 5 files changed, 101 insertions(+), 25 deletions(-)
 rename src/openfe/protocols/openmm_utils/{molecule_utils.py => offmolecule_utils.py} (100%)

diff --git a/news/fix_unk_resnames.rst b/news/fix_unk_resnames.rst
index 8f252ccef..0dd6b8302 100644
--- a/news/fix_unk_resnames.rst
+++ b/news/fix_unk_resnames.rst
@@ -4,7 +4,7 @@
 
 **Changed:**
 
-* Small molecules in RelativeHybridTopology topologies (including the output PDB) are now named LIG (alchemical ligand) and CF1, CF2… (cofactors) instead of UNK.
+* Small molecules in RelativeHybridTopology topologies (including the output PDB) are now named LIG (alchemical ligand) and CF1, CF2… (cofactors) instead of UNK. If a residue name was already assigned, the assigned one is kept.
 
 **Deprecated:**
 
diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index e92239668..1cefbb3a4 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -46,7 +46,7 @@
 from openmmtools import multistate
 
 import openfe
-from openfe.protocols.openmm_utils.molecule_utils import _get_offmol_resname, _set_offmol_resname
+from openfe.protocols.openmm_utils.offmolecule_utils import _get_offmol_resname, _set_offmol_resname
 from openfe.protocols.openmm_utils.omm_settings import (
     BasePartialChargeSettings,
 )
@@ -765,6 +765,9 @@ def run(
             else:
                 _set_offmol_resname(offmol, f"CF{cofactors.index(smc) + 1}")
 
+        # Record the resname the alchemical ligand actually ended up with
+        alchem_resname = _get_offmol_resname(small_mols[mapping.componentA])
+
         # Assign partial charges now to avoid any discrepancies later
         self._assign_partial_charges(settings["charge_settings"], small_mols)
 
@@ -822,6 +825,7 @@ def run(
             "positions": positions_outfile,
             "pdb_structure": self.shared_basepath / settings["output_settings"].output_structure,
             "selection_indices": selection_indices,
+            "alchemical_resname": alchem_resname,
         }
 
         if dry:
@@ -1517,6 +1521,7 @@ def _structural_analysis(
         trj_file: pathlib.Path,
         output_directory: pathlib.Path,
         dry: bool,
+        ligand_resname: str = "LIG",
     ) -> dict[str, str | pathlib.Path]:
         """
         Run structural analysis using ``openfe-analysis``.
@@ -1532,6 +1537,8 @@ def _structural_analysis(
           will be stored.
         dry : bool
           Whether or not we are running a dry run.
+        ligand_resname: str
+          The residue name of the hybrid topology ligand. Default: LIG
 
         Returns
         -------
@@ -1547,7 +1554,7 @@ def _structural_analysis(
         from openfe_analysis import rmsd
 
         try:
-            data = rmsd.gather_rms_data(pdb_file, trj_file, ligand_selection="resname LIG")
+            data = rmsd.gather_rms_data(pdb_file, trj_file, ligand_selection=f"resname {ligand_resname}")
         # TODO: eventually change this to more specific exception types
         except Exception as e:
             return {"structural_analysis_error": str(e)}
@@ -1585,6 +1592,7 @@ def run(
         pdb_file: pathlib.Path,
         trajectory: pathlib.Path,
         checkpoint: pathlib.Path,
+        ligand_resname: str = "LIG",
         dry: bool = False,
         verbose: bool = True,
         scratch_basepath: pathlib.Path | None = None,
@@ -1600,6 +1608,8 @@ def run(
           Path to the MultiStateReporter generated NetCDF file.
         checkpoint : pathlib.Path
           Path to the checkpoint file generated by MultiStateReporter.
+        ligand_resname: str
+          The residue name of the hybrid topology ligand. Default: LIG
         dry : bool
           Do a dry run of the calculation, creating all necessary hybrid
           system components (topology, system, sampler, etc...) but without
@@ -1653,6 +1663,7 @@ def run(
             trj_file=trajectory,
             output_directory=self.shared_basepath,
             dry=dry,
+            ligand_resname=ligand_resname,
         )
 
         # Return relevant things
@@ -1681,6 +1692,8 @@ def _execute(
             pdb_file=pdb_file,
             trajectory=trajectory,
             checkpoint=checkpoint,
+            ligand_resname=setup_results.outputs.get("alchemical_resname",
+                                                     "LIG"), # Adding this for backward compatibility
             scratch_basepath=ctx.scratch,
             shared_basepath=ctx.shared,
         )
diff --git a/src/openfe/protocols/openmm_utils/molecule_utils.py b/src/openfe/protocols/openmm_utils/offmolecule_utils.py
similarity index 100%
rename from src/openfe/protocols/openmm_utils/molecule_utils.py
rename to src/openfe/protocols/openmm_utils/offmolecule_utils.py
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 40857e73d..516d3e752 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -51,7 +51,7 @@
     HAS_NAGL,
     HAS_OPENEYE,
 )
-from openfe.protocols.openmm_utils.molecule_utils import (
+from openfe.protocols.openmm_utils.offmolecule_utils import (
     _get_offmol_metadata,
     _get_offmol_resname,
     _set_offmol_metadata,
@@ -2440,25 +2440,6 @@ def test_ligand_separable_from_cofactor(eg5_vac_inputs, vac_settings, tmp_path):
     assert set(lig.indices).isdisjoint(cof.indices)
 
 
-def test_get_metadata_inconsistent_warns(caplog):
-    mol = Molecule.from_smiles("CC")
-    _set_offmol_metadata(mol, "residue_name", "LIG")
-    mol.atoms[0].metadata["residue_name"] = "COF"
-
-    with caplog.at_level(logging.WARNING):
-        result = _get_offmol_metadata(mol, "residue_name")
-
-    assert result is None
-    assert "Inconsistent metadata" in caplog.text
-
-
-def test_set_metadata_none_clears():
-    mol = Molecule.from_smiles("CC")
-    _set_offmol_metadata(mol, "residue_name", "LIG")
-    _set_offmol_metadata(mol, "residue_name", None)
-    assert all("residue_name" not in a.metadata for a in mol.atoms)
-
-
 def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
     vac_settings.output_settings.output_indices = "all"
 
@@ -2478,3 +2459,62 @@ def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp
     assert "MYC" in u.residues.resnames
     assert "CF1" not in u.residues.resnames
     assert "LIG" in u.residues.resnames
+
+
+def test_hybrid_ligand_takes_componentA_resname(
+    eg5_ligands, eg5_cofactor, vac_settings, tmp_path
+):
+    vac_settings.output_settings.output_indices = "all"
+
+    ligA, ligB = eg5_ligands[0], eg5_ligands[1]
+    offA, offB = ligA.to_openff(), ligB.to_openff()
+    _set_offmol_resname(offA, "LGA")
+    _set_offmol_resname(offB, "LGB")
+    ligA_named = openfe.SmallMoleculeComponent.from_openff(offA)
+    ligB_named = openfe.SmallMoleculeComponent.from_openff(offB)
+
+    mapper = openfe.setup.KartografAtomMapper()
+    mapping = next(mapper.suggest_mappings(ligA_named, ligB_named))
+    stateA = openfe.ChemicalSystem({"ligand": ligA_named, "cofactor": eg5_cofactor})
+    stateB = openfe.ChemicalSystem({"ligand": ligB_named, "cofactor": eg5_cofactor})
+
+    out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
+
+    u = mda.Universe(out["pdb_structure"])
+    lig = u.select_atoms("resname LGA")
+    cof = u.select_atoms("resname CF1")
+
+    # the setup unit records componentA's resname for the analysis unit
+    assert out["alchemical_resname"] == "LGA"
+
+    # the whole hybrid ligand is under componentA's resname
+    assert lig.n_atoms > 0
+    assert u.select_atoms("resname LGB").n_atoms == 0
+    # user-set names are kept, not overwritten with the LIG default
+    assert u.select_atoms("resname LIG").n_atoms == 0
+    assert "UNK" not in {r.resname for r in u.residues}
+    # cofactor stays separate
+    assert cof.n_atoms > 0
+    assert set(lig.indices).isdisjoint(cof.indices)
+
+
+def test_structural_analysis_uses_ligand_resname():
+    captured = {}
+
+    def fake_gather(pdb, trj, ligand_selection="resname UNK"):
+        captured["ligand_selection"] = ligand_selection
+        return {
+            "protein_RMSD": [], "ligand_RMSD": [], "ligand_wander": [],
+            "protein_2D_RMSD": [], "time(ps)": [],
+        }
+
+    with mock.patch("openfe_analysis.rmsd.gather_rms_data", side_effect=fake_gather):
+        HybridTopologyMultiStateAnalysisUnit._structural_analysis(
+            pdb_file=Path("dummy.pdb"),
+            trj_file=Path("dummy.nc"),
+            output_directory=Path("."),
+            dry=True,
+            ligand_resname="LGA",
+        )
+
+    assert captured["ligand_selection"] == "resname LGA"
\ No newline at end of file
diff --git a/src/openfe/tests/protocols/test_openmmutils.py b/src/openfe/tests/protocols/test_openmmutils.py
index ffacad813..e230afbce 100644
--- a/src/openfe/tests/protocols/test_openmmutils.py
+++ b/src/openfe/tests/protocols/test_openmmutils.py
@@ -2,6 +2,7 @@
 # For details, see https://github.com/OpenFreeEnergy/openfe
 import copy
 import gzip
+import logging
 import os
 import sys
 from importlib import resources
@@ -9,9 +10,7 @@
 from unittest import mock
 
 import numpy as np
-import pooch
 import pytest
-from gufe import BaseSolventComponent
 from gufe.components.errors import ComponentValidationError
 from gufe.settings import OpenMMSystemGeneratorFFSettings, ThermoSettings
 from numpy.testing import assert_allclose, assert_equal
@@ -44,6 +43,11 @@
     HAS_NAGL,
     HAS_OPENEYE,
 )
+from openfe.protocols.openmm_utils.offmolecule_utils import (
+    _get_offmol_metadata,
+    _set_offmol_metadata,
+    _set_offmol_resname,
+)
 
 from ..conftest import HAS_INTERNET
 
@@ -1189,3 +1193,22 @@ def test_forward_backwards_failure(simulation_nc):
         ret = ana.get_forward_and_reverse_analysis()
 
     assert ret is None
+
+
+def test_get_metadata_inconsistent_warns(caplog):
+    mol = OFFMol.from_smiles("CC")
+    _set_offmol_metadata(mol, "residue_name", "LIG")
+    mol.atoms[0].metadata["residue_name"] = "COF"
+
+    with caplog.at_level(logging.WARNING):
+        result = _get_offmol_metadata(mol, "residue_name")
+
+    assert result is None
+    assert "Inconsistent metadata" in caplog.text
+
+
+def test_set_metadata_none_clears():
+    mol = OFFMol.from_smiles("CC")
+    _set_offmol_metadata(mol, "residue_name", "LIG")
+    _set_offmol_metadata(mol, "residue_name", None)
+    assert all("residue_name" not in a.metadata for a in mol.atoms)

From c68b88eca8c8228e26ba4f97b0f3c85d37ebe4c8 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 2 Jul 2026 09:54:04 +0000
Subject: [PATCH 09/19] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/openfe/protocols/openmm_rfe/hybridtop_units.py  |  9 ++++++---
 .../openmm_rfe/test_hybrid_top_protocol.py          | 13 +++++++------
 2 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 1cefbb3a4..2e5436c18 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -1554,7 +1554,9 @@ def _structural_analysis(
         from openfe_analysis import rmsd
 
         try:
-            data = rmsd.gather_rms_data(pdb_file, trj_file, ligand_selection=f"resname {ligand_resname}")
+            data = rmsd.gather_rms_data(
+                pdb_file, trj_file, ligand_selection=f"resname {ligand_resname}"
+            )
         # TODO: eventually change this to more specific exception types
         except Exception as e:
             return {"structural_analysis_error": str(e)}
@@ -1692,8 +1694,9 @@ def _execute(
             pdb_file=pdb_file,
             trajectory=trajectory,
             checkpoint=checkpoint,
-            ligand_resname=setup_results.outputs.get("alchemical_resname",
-                                                     "LIG"), # Adding this for backward compatibility
+            ligand_resname=setup_results.outputs.get(
+                "alchemical_resname", "LIG"
+            ),  # Adding this for backward compatibility
             scratch_basepath=ctx.scratch,
             shared_basepath=ctx.shared,
         )
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 516d3e752..662b46ae1 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2461,9 +2461,7 @@ def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp
     assert "LIG" in u.residues.resnames
 
 
-def test_hybrid_ligand_takes_componentA_resname(
-    eg5_ligands, eg5_cofactor, vac_settings, tmp_path
-):
+def test_hybrid_ligand_takes_componentA_resname(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
     vac_settings.output_settings.output_indices = "all"
 
     ligA, ligB = eg5_ligands[0], eg5_ligands[1]
@@ -2504,8 +2502,11 @@ def test_structural_analysis_uses_ligand_resname():
     def fake_gather(pdb, trj, ligand_selection="resname UNK"):
         captured["ligand_selection"] = ligand_selection
         return {
-            "protein_RMSD": [], "ligand_RMSD": [], "ligand_wander": [],
-            "protein_2D_RMSD": [], "time(ps)": [],
+            "protein_RMSD": [],
+            "ligand_RMSD": [],
+            "ligand_wander": [],
+            "protein_2D_RMSD": [],
+            "time(ps)": [],
         }
 
     with mock.patch("openfe_analysis.rmsd.gather_rms_data", side_effect=fake_gather):
@@ -2517,4 +2518,4 @@ def fake_gather(pdb, trj, ligand_selection="resname UNK"):
             ligand_resname="LGA",
         )
 
-    assert captured["ligand_selection"] == "resname LGA"
\ No newline at end of file
+    assert captured["ligand_selection"] == "resname LGA"

From de84279d938095f675ec3eafadc576b414a62855 Mon Sep 17 00:00:00 2001
From: Hannah Baumann <43765638+hannahbaumann@users.noreply.github.com>
Date: Fri, 3 Jul 2026 15:59:19 +0200
Subject: [PATCH 10/19] Update
 src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py

Co-authored-by: Irfan Alibay <IAlibay@users.noreply.github.com>
---
 .../tests/protocols/openmm_rfe/test_hybrid_top_protocol.py      | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 662b46ae1..9e62ecf16 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -52,9 +52,7 @@
     HAS_OPENEYE,
 )
 from openfe.protocols.openmm_utils.offmolecule_utils import (
-    _get_offmol_metadata,
     _get_offmol_resname,
-    _set_offmol_metadata,
     _set_offmol_resname,
 )
 

From 8fca8353d3309fb0189ea1854ab047367d4bc17b Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Fri, 3 Jul 2026 16:30:24 +0200
Subject: [PATCH 11/19] address review comments

---
 .../protocols/openmm_rfe/hybridtop_units.py   | 62 +++++++++----
 .../openmm_rfe/test_hybrid_top_protocol.py    | 87 ++++++++++++-------
 2 files changed, 103 insertions(+), 46 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 2e5436c18..3307e8af8 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -755,18 +755,46 @@ def run(
         solvent_comp, protein_comp, small_mols = self._get_components(stateA, stateB)
 
         alchemical = set(alchem_comps["stateA"]) | set(alchem_comps["stateB"])
-        cofactors = [smc for smc in small_mols if smc not in alchemical]
 
+        # Store user-provided resnames, so auto-assigned names can't clash with them.
+        used: set[str] = set()
+        for offmol in small_mols.values():
+            name = _get_offmol_resname(offmol)
+            if name is not None:
+                used.add(name)
+
+        def _unique(base: str) -> str:
+            # If the requested name is already in use, append a numeric suffix (e.g.
+            # CF1 -> CF11) until a unique name is found.
+            if base not in used:
+                return base
+            i = 1
+            candidate = f"{base}{i}"
+            while candidate in used:
+                i += 1
+                candidate = f"{base}{i}"
+            return candidate
+
+        cof_counter = 1
         for smc, offmol in small_mols.items():
             if _get_offmol_resname(offmol) is not None:
                 continue
+
             if smc in alchemical:
-                _set_offmol_resname(offmol, "LIG")
+                base = "LIG"
             else:
-                _set_offmol_resname(offmol, f"CF{cofactors.index(smc) + 1}")
+                base = f"CF{cof_counter}"
+                cof_counter += 1
+
+            name = _unique(base)
+            _set_offmol_resname(offmol, name)
+            used.add(name)
 
-        # Record the resname the alchemical ligand actually ended up with
-        alchem_resname = _get_offmol_resname(small_mols[mapping.componentA])
+        # Record the resnames of the ligands
+        alchem_resnames = sorted({
+            _get_offmol_resname(small_mols[comp])
+            for comp in (mapping.componentA, mapping.componentB)
+        })
 
         # Assign partial charges now to avoid any discrepancies later
         self._assign_partial_charges(settings["charge_settings"], small_mols)
@@ -825,7 +853,7 @@ def run(
             "positions": positions_outfile,
             "pdb_structure": self.shared_basepath / settings["output_settings"].output_structure,
             "selection_indices": selection_indices,
-            "alchemical_resname": alchem_resname,
+            "alchemical_resnames": alchem_resnames,
         }
 
         if dry:
@@ -1521,7 +1549,7 @@ def _structural_analysis(
         trj_file: pathlib.Path,
         output_directory: pathlib.Path,
         dry: bool,
-        ligand_resname: str = "LIG",
+        ligand_resnames: list[str],
     ) -> dict[str, str | pathlib.Path]:
         """
         Run structural analysis using ``openfe-analysis``.
@@ -1537,8 +1565,8 @@ def _structural_analysis(
           will be stored.
         dry : bool
           Whether or not we are running a dry run.
-        ligand_resname: str
-          The residue name of the hybrid topology ligand. Default: LIG
+        ligand_resnames: list[str]
+          The residue names of the ligands.
 
         Returns
         -------
@@ -1555,7 +1583,7 @@ def _structural_analysis(
 
         try:
             data = rmsd.gather_rms_data(
-                pdb_file, trj_file, ligand_selection=f"resname {ligand_resname}"
+                pdb_file, trj_file, ligand_selection=f"resname {ligand_resnames}"
             )
         # TODO: eventually change this to more specific exception types
         except Exception as e:
@@ -1594,7 +1622,7 @@ def run(
         pdb_file: pathlib.Path,
         trajectory: pathlib.Path,
         checkpoint: pathlib.Path,
-        ligand_resname: str = "LIG",
+        ligand_resnames: list[str] = ["LIG"],
         dry: bool = False,
         verbose: bool = True,
         scratch_basepath: pathlib.Path | None = None,
@@ -1610,8 +1638,8 @@ def run(
           Path to the MultiStateReporter generated NetCDF file.
         checkpoint : pathlib.Path
           Path to the checkpoint file generated by MultiStateReporter.
-        ligand_resname: str
-          The residue name of the hybrid topology ligand. Default: LIG
+        ligand_resnames: list[str]
+          The residue names of the ligands. Default: ["LIG"]
         dry : bool
           Do a dry run of the calculation, creating all necessary hybrid
           system components (topology, system, sampler, etc...) but without
@@ -1665,7 +1693,7 @@ def run(
             trj_file=trajectory,
             output_directory=self.shared_basepath,
             dry=dry,
-            ligand_resname=ligand_resname,
+            ligand_resnames=ligand_resnames,
         )
 
         # Return relevant things
@@ -1694,9 +1722,9 @@ def _execute(
             pdb_file=pdb_file,
             trajectory=trajectory,
             checkpoint=checkpoint,
-            ligand_resname=setup_results.outputs.get(
-                "alchemical_resname", "LIG"
-            ),  # Adding this for backward compatibility
+            ligand_resnames=setup_results.outputs.get(
+                "alchemical_resnames", ["LIG"]
+            ),  # Adding "LIG" for backward compatibility
             scratch_basepath=ctx.scratch,
             shared_basepath=ctx.shared,
         )
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 9e62ecf16..923da5930 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2425,6 +2425,14 @@ def _run_setup_dry(stateA, stateB, mapping, settings, tmp_path):
     return setup_unit.run(dry=True, scratch_basepath=tmp_path, shared_basepath=tmp_path)
 
 
+def _named_smc(smc, resname):
+    off = smc.to_openff()
+    _set_offmol_resname(off, resname)
+    result = openfe.SmallMoleculeComponent.from_openff(off)
+    assert _get_offmol_resname(result.to_openff()) == resname
+    return result
+
+
 def test_ligand_separable_from_cofactor(eg5_vac_inputs, vac_settings, tmp_path):
     vac_settings.output_settings.output_indices = "all"
     stateA, stateB, mapping = eg5_vac_inputs
@@ -2438,13 +2446,17 @@ def test_ligand_separable_from_cofactor(eg5_vac_inputs, vac_settings, tmp_path):
     assert set(lig.indices).isdisjoint(cof.indices)
 
 
-def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+def test_default_resnames_stored(eg5_vac_inputs, vac_settings, tmp_path):
     vac_settings.output_settings.output_indices = "all"
+    stateA, stateB, mapping = eg5_vac_inputs
+    out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
+    assert out["alchemical_resnames"] == ["LIG"]
+
 
+def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+    vac_settings.output_settings.output_indices = "all"
     ligA, ligB = eg5_ligands[0], eg5_ligands[1]
-    off = eg5_cofactor.to_openff()
-    _set_offmol_resname(off, "MYC")
-    cof = openfe.SmallMoleculeComponent.from_openff(off)
+    cof = _named_smc(eg5_cofactor, "MYC")
 
     mapper = openfe.setup.KartografAtomMapper()
     mapping = next(mapper.suggest_mappings(ligA, ligB))
@@ -2459,15 +2471,10 @@ def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp
     assert "LIG" in u.residues.resnames
 
 
-def test_hybrid_ligand_takes_componentA_resname(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+def test_distinct_endstate_resnames(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
     vac_settings.output_settings.output_indices = "all"
-
-    ligA, ligB = eg5_ligands[0], eg5_ligands[1]
-    offA, offB = ligA.to_openff(), ligB.to_openff()
-    _set_offmol_resname(offA, "LGA")
-    _set_offmol_resname(offB, "LGB")
-    ligA_named = openfe.SmallMoleculeComponent.from_openff(offA)
-    ligB_named = openfe.SmallMoleculeComponent.from_openff(offB)
+    ligA_named = _named_smc(eg5_ligands[0], "LGA")
+    ligB_named = _named_smc(eg5_ligands[1], "LGB")
 
     mapper = openfe.setup.KartografAtomMapper()
     mapping = next(mapper.suggest_mappings(ligA_named, ligB_named))
@@ -2476,35 +2483,57 @@ def test_hybrid_ligand_takes_componentA_resname(eg5_ligands, eg5_cofactor, vac_s
 
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
+    # both endstate names stored
+    assert out["alchemical_resnames"] == ["LGA", "LGB"]
+
     u = mda.Universe(out["pdb_structure"])
     lig = u.select_atoms("resname LGA")
     cof = u.select_atoms("resname CF1")
-
-    # the setup unit records componentA's resname for the analysis unit
-    assert out["alchemical_resname"] == "LGA"
-
-    # the whole hybrid ligand is under componentA's resname
+    # hybrid ligand carries componentA's name; LGB never reaches the hybrid
     assert lig.n_atoms > 0
     assert u.select_atoms("resname LGB").n_atoms == 0
-    # user-set names are kept, not overwritten with the LIG default
     assert u.select_atoms("resname LIG").n_atoms == 0
-    assert "UNK" not in {r.resname for r in u.residues}
-    # cofactor stays separate
+    assert "UNK" not in set(u.residues.resnames)
     assert cof.n_atoms > 0
     assert set(lig.indices).isdisjoint(cof.indices)
+    # the stored union recovers the full ligand
+    sel = u.select_atoms("resname " + " ".join(out["alchemical_resnames"]))
+    assert set(sel.indices) == set(lig.indices)
+
+
+def test_cofactor_ligand_name_clash(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+    """A cofactor pre-named LIG must not collide with the ligand."""
+    vac_settings.output_settings.output_indices = "all"
+    ligA, ligB = eg5_ligands[0], eg5_ligands[1]
+    cof = _named_smc(eg5_cofactor, "LIG")
+
+    mapper = openfe.setup.KartografAtomMapper()
+    mapping = next(mapper.suggest_mappings(ligA, ligB))
+    stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": cof})
+    stateB = openfe.ChemicalSystem({"ligand": ligB, "cofactor": cof})
+
+    out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
+
+    u = mda.Universe(out["pdb_structure"])
+    resnames = list(u.residues.resnames)
+    assert len(resnames) == len(set(resnames))          # all residues distinct
+    # pre-named cofactor keeps LIG; ligand bumps
+    assert u.select_atoms("resname LIG").n_atoms > 0     # the cofactor
+    ligand_names = out["alchemical_resnames"]
+    assert "LIG" not in ligand_names
+    lig = u.select_atoms("resname " + " ".join(ligand_names))
+    cof = u.select_atoms("resname LIG")
+    assert set(lig.indices).isdisjoint(cof.indices)
 
 
-def test_structural_analysis_uses_ligand_resname():
+def test_structural_analysis_uses_ligand_resnames():
     captured = {}
 
-    def fake_gather(pdb, trj, ligand_selection="resname UNK"):
+    def fake_gather(*args, ligand_selection="resname UNK"):
         captured["ligand_selection"] = ligand_selection
         return {
-            "protein_RMSD": [],
-            "ligand_RMSD": [],
-            "ligand_wander": [],
-            "protein_2D_RMSD": [],
-            "time(ps)": [],
+            "protein_RMSD": [], "ligand_RMSD": [], "ligand_wander": [],
+            "protein_2D_RMSD": [], "time(ps)": [],
         }
 
     with mock.patch("openfe_analysis.rmsd.gather_rms_data", side_effect=fake_gather):
@@ -2513,7 +2542,7 @@ def fake_gather(pdb, trj, ligand_selection="resname UNK"):
             trj_file=Path("dummy.nc"),
             output_directory=Path("."),
             dry=True,
-            ligand_resname="LGA",
+            ligand_resnames=["LGA", "LGB"],
         )
 
-    assert captured["ligand_selection"] == "resname LGA"
+    assert captured["ligand_selection"] == "resname LGA LGB"

From 9b21e8057116cbd6322df748675d8dc9a680d9e2 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 3 Jul 2026 14:31:56 +0000
Subject: [PATCH 12/19] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/openfe/protocols/openmm_rfe/hybridtop_units.py    | 10 ++++++----
 .../protocols/openmm_rfe/test_hybrid_top_protocol.py  | 11 +++++++----
 2 files changed, 13 insertions(+), 8 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 3307e8af8..22fe22f27 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -791,10 +791,12 @@ def _unique(base: str) -> str:
             used.add(name)
 
         # Record the resnames of the ligands
-        alchem_resnames = sorted({
-            _get_offmol_resname(small_mols[comp])
-            for comp in (mapping.componentA, mapping.componentB)
-        })
+        alchem_resnames = sorted(
+            {
+                _get_offmol_resname(small_mols[comp])
+                for comp in (mapping.componentA, mapping.componentB)
+            }
+        )
 
         # Assign partial charges now to avoid any discrepancies later
         self._assign_partial_charges(settings["charge_settings"], small_mols)
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 923da5930..172a4a71f 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2516,9 +2516,9 @@ def test_cofactor_ligand_name_clash(eg5_ligands, eg5_cofactor, vac_settings, tmp
 
     u = mda.Universe(out["pdb_structure"])
     resnames = list(u.residues.resnames)
-    assert len(resnames) == len(set(resnames))          # all residues distinct
+    assert len(resnames) == len(set(resnames))  # all residues distinct
     # pre-named cofactor keeps LIG; ligand bumps
-    assert u.select_atoms("resname LIG").n_atoms > 0     # the cofactor
+    assert u.select_atoms("resname LIG").n_atoms > 0  # the cofactor
     ligand_names = out["alchemical_resnames"]
     assert "LIG" not in ligand_names
     lig = u.select_atoms("resname " + " ".join(ligand_names))
@@ -2532,8 +2532,11 @@ def test_structural_analysis_uses_ligand_resnames():
     def fake_gather(*args, ligand_selection="resname UNK"):
         captured["ligand_selection"] = ligand_selection
         return {
-            "protein_RMSD": [], "ligand_RMSD": [], "ligand_wander": [],
-            "protein_2D_RMSD": [], "time(ps)": [],
+            "protein_RMSD": [],
+            "ligand_RMSD": [],
+            "ligand_wander": [],
+            "protein_2D_RMSD": [],
+            "time(ps)": [],
         }
 
     with mock.patch("openfe_analysis.rmsd.gather_rms_data", side_effect=fake_gather):

From 6ae227d8ffa3818bc66ecb3a4538ffbf96e5a989 Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Fri, 3 Jul 2026 16:49:11 +0200
Subject: [PATCH 13/19] Make mypy happy

---
 src/openfe/protocols/openmm_rfe/hybridtop_units.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 3307e8af8..1c5cf6feb 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -791,10 +791,12 @@ def _unique(base: str) -> str:
             used.add(name)
 
         # Record the resnames of the ligands
-        alchem_resnames = sorted({
-            _get_offmol_resname(small_mols[comp])
-            for comp in (mapping.componentA, mapping.componentB)
-        })
+        names = set()
+        for comp in (mapping.componentA, mapping.componentB):
+            name = _get_offmol_resname(small_mols[comp])
+            assert name is not None, "alchemical ligand has no resname"
+            names.add(name)
+        alchem_resnames = sorted(names)
 
         # Assign partial charges now to avoid any discrepancies later
         self._assign_partial_charges(settings["charge_settings"], small_mols)

From 3bcbca833c3be832563a5e4a0a3488a642bf0282 Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Mon, 6 Jul 2026 15:21:06 +0200
Subject: [PATCH 14/19] Address review comments

---
 .../protocols/openmm_rfe/hybridtop_units.py   |  62 ++++-----
 .../hybrid_topology_atoms.csv                 |  32 ++---
 .../openmm_rfe/test_hybrid_top_protocol.py    | 119 ++++++++----------
 3 files changed, 102 insertions(+), 111 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 1c5cf6feb..3dad38476 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -693,9 +693,14 @@ def _subsample_topology(
         bfactors[np.isin(selection_indices, list(atom_classes["unique_new_atoms"]))] = 0.75
 
         if len(selection_indices) > 0:
+            sub_top = hybrid_topology.subset(selection_indices)
+            # Renumber sequentially so same-named residues (e.g. two COF cofactors,
+            # both resSeq=0) don't merge into one residue on PDB write/read.
+            for i, res in enumerate(sub_top.residues, start=1):
+                res.resSeq = i
             traj = mdt.Trajectory(
                 hybrid_positions[selection_indices, :],
-                hybrid_topology.subset(selection_indices),
+                sub_top,
             ).save_pdb(
                 self.shared_basepath / output_filename,
                 bfactors=bfactors,
@@ -755,46 +760,43 @@ def run(
         solvent_comp, protein_comp, small_mols = self._get_components(stateA, stateB)
 
         alchemical = set(alchem_comps["stateA"]) | set(alchem_comps["stateB"])
+        stateA_comps = set(alchem_comps["stateA"])
+
+        def _unique(used: set[str]) -> str:
+            """Next free 'LG{n}' name (n = 1..9) for the alchemical ligands."""
+            for i in range(1, 10):
+                candidate = f"LG{i}"
+                if candidate not in used:
+                    return candidate
+            raise ValueError(
+                "Could not assign a unique ligand residue name; too many "
+                "colliding 'LG#' names."
+            )
 
-        # Store user-provided resnames, so auto-assigned names can't clash with them.
+        # Seed with user-provided resnames so auto-assigned names avoid them.
         used: set[str] = set()
         for offmol in small_mols.values():
             name = _get_offmol_resname(offmol)
             if name is not None:
                 used.add(name)
 
-        def _unique(base: str) -> str:
-            # If the requested name is already in use, append a numeric suffix (e.g.
-            # CF1 -> CF11) until a unique name is found.
-            if base not in used:
-                return base
-            i = 1
-            candidate = f"{base}{i}"
-            while candidate in used:
-                i += 1
-                candidate = f"{base}{i}"
-            return candidate
-
-        cof_counter = 1
+        # ligands take LG1 (endstate A) and LG2 (endstate B) by default, only
+        # bumping to LG3+ if a user pre-named something LG1/LG2. All cofactors
+        # share "COF" and are distinguished by residue index, not name.
         for smc, offmol in small_mols.items():
             if _get_offmol_resname(offmol) is not None:
                 continue
-
             if smc in alchemical:
-                base = "LIG"
+                base = "LG1" if smc in stateA_comps else "LG2"
+                name = base if base not in used else _unique(used)
+                used.add(name)
             else:
-                base = f"CF{cof_counter}"
-                cof_counter += 1
-
-            name = _unique(base)
+                name = "COF"
             _set_offmol_resname(offmol, name)
-            used.add(name)
 
-        # Record the resnames of the ligands
         names = set()
         for comp in (mapping.componentA, mapping.componentB):
             name = _get_offmol_resname(small_mols[comp])
-            assert name is not None, "alchemical ligand has no resname"
             names.add(name)
         alchem_resnames = sorted(names)
 
@@ -1551,7 +1553,7 @@ def _structural_analysis(
         trj_file: pathlib.Path,
         output_directory: pathlib.Path,
         dry: bool,
-        ligand_resnames: list[str],
+        ligand_resnames: list[str] | None = None,
     ) -> dict[str, str | pathlib.Path]:
         """
         Run structural analysis using ``openfe-analysis``.
@@ -1585,7 +1587,7 @@ def _structural_analysis(
 
         try:
             data = rmsd.gather_rms_data(
-                pdb_file, trj_file, ligand_selection=f"resname {ligand_resnames}"
+                pdb_file, trj_file, ligand_selection="resname " + " ".join(ligand_resnames)
             )
         # TODO: eventually change this to more specific exception types
         except Exception as e:
@@ -1624,7 +1626,7 @@ def run(
         pdb_file: pathlib.Path,
         trajectory: pathlib.Path,
         checkpoint: pathlib.Path,
-        ligand_resnames: list[str] = ["LIG"],
+        ligand_resnames: list[str] = ["LG1", "LG2"],
         dry: bool = False,
         verbose: bool = True,
         scratch_basepath: pathlib.Path | None = None,
@@ -1641,7 +1643,7 @@ def run(
         checkpoint : pathlib.Path
           Path to the checkpoint file generated by MultiStateReporter.
         ligand_resnames: list[str]
-          The residue names of the ligands. Default: ["LIG"]
+          The residue names of the ligands. Default: ["LG1", "LG2"]
         dry : bool
           Do a dry run of the calculation, creating all necessary hybrid
           system components (topology, system, sampler, etc...) but without
@@ -1725,8 +1727,8 @@ def _execute(
             trajectory=trajectory,
             checkpoint=checkpoint,
             ligand_resnames=setup_results.outputs.get(
-                "alchemical_resnames", ["LIG"]
-            ),  # Adding "LIG" for backward compatibility
+                "alchemical_resnames", ["LG1", "LG2"]
+            ),
             scratch_basepath=ctx.scratch,
             shared_basepath=ctx.shared,
         )
diff --git a/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv b/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv
index 86ed7f043..9c6e9c7ce 100644
--- a/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv
+++ b/src/openfe/tests/data/openmm_rfe/benzene_toluene_hybrid_top/hybrid_topology_atoms.csv
@@ -1,17 +1,17 @@
 ,serial,name,element,resSeq,resName,chainID,segmentID,formal_charge
-0,,C1x,C,0,LIG,0,,
-1,,C2x,C,0,LIG,0,,
-2,,C3x,C,0,LIG,0,,
-3,,C4x,C,0,LIG,0,,
-4,,C5x,C,0,LIG,0,,
-5,,C6x,C,0,LIG,0,,
-6,,H1x,H,0,LIG,0,,
-7,,H2x,H,0,LIG,0,,
-8,,H3x,H,0,LIG,0,,
-9,,H4x,H,0,LIG,0,,
-10,,H5x,H,0,LIG,0,,
-11,,H6x,H,0,LIG,0,,
-12,,H1x,H,0,LIG,0,,
-13,,H2x,H,0,LIG,0,,
-14,,C1x,C,0,LIG,0,,
-15,,H3x,H,0,LIG,0,,
+0,,C1x,C,0,LG1,0,,
+1,,C2x,C,0,LG1,0,,
+2,,C3x,C,0,LG1,0,,
+3,,C4x,C,0,LG1,0,,
+4,,C5x,C,0,LG1,0,,
+5,,C6x,C,0,LG1,0,,
+6,,H1x,H,0,LG1,0,,
+7,,H2x,H,0,LG1,0,,
+8,,H3x,H,0,LG1,0,,
+9,,H4x,H,0,LG1,0,,
+10,,H5x,H,0,LG1,0,,
+11,,H6x,H,0,LG1,0,,
+12,,H1x,H,0,LG1,0,,
+13,,H2x,H,0,LG1,0,,
+14,,C1x,C,0,LG1,0,,
+15,,H3x,H,0,LG1,0,,
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 172a4a71f..c9cae68a4 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -512,7 +512,7 @@ def test_dry_run_ligand(
 ):
     # this might be a bit time consuming
     solv_settings.simulation_settings.sampler_method = method
-    solv_settings.output_settings.output_indices = "resname LIG"
+    solv_settings.output_settings.output_indices = "resname LG1 LG2"
 
     protocol = openmm_rfe.RelativeHybridTopologyProtocol(
         settings=solv_settings,
@@ -1096,7 +1096,7 @@ def test_dry_run_complex(
 ):
     # this will be very time consuming
     solv_settings.simulation_settings.sampler_method = method
-    solv_settings.output_settings.output_indices = "protein or resname  LIG"
+    solv_settings.output_settings.output_indices = "protein or resname  LG1 or resname LG2"
 
     protocol = openmm_rfe.RelativeHybridTopologyProtocol(
         settings=solv_settings,
@@ -2426,87 +2426,84 @@ def _run_setup_dry(stateA, stateB, mapping, settings, tmp_path):
 
 
 def _named_smc(smc, resname):
-    off = smc.to_openff()
+    off = copy.deepcopy(smc.to_openff())
     _set_offmol_resname(off, resname)
     result = openfe.SmallMoleculeComponent.from_openff(off)
     assert _get_offmol_resname(result.to_openff()) == resname
     return result
 
 
-def test_ligand_separable_from_cofactor(eg5_vac_inputs, vac_settings, tmp_path):
+def test_default_resnames(eg5_vac_inputs, vac_settings, tmp_path):
+    """Unnamed ligands -> LG1/LG2 stored; cofactor -> COF; all separable."""
     vac_settings.output_settings.output_indices = "all"
     stateA, stateB, mapping = eg5_vac_inputs
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
+    assert out["alchemical_resnames"] == ["LG1", "LG2"]
+
     u = mda.Universe(out["pdb_structure"])
-    lig = u.select_atoms("resname LIG")
-    cof = u.select_atoms("resname CF1")
+    lig = u.select_atoms("resname " + " ".join(out["alchemical_resnames"]))
+    cof = u.select_atoms("resname COF")
     assert lig.n_atoms > 0
     assert cof.n_atoms > 0
     assert set(lig.indices).isdisjoint(cof.indices)
 
 
-def test_default_resnames_stored(eg5_vac_inputs, vac_settings, tmp_path):
-    vac_settings.output_settings.output_indices = "all"
-    stateA, stateB, mapping = eg5_vac_inputs
-    out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
-    assert out["alchemical_resnames"] == ["LIG"]
-
-
-def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+def test_cofactors_share_cof_distinct_resindices(
+    eg5_ligands, eg5_cofactor, benzene_modifications, vac_settings, tmp_path
+):
+    """Two cofactors share the name COF but are separate residues."""
     vac_settings.output_settings.output_indices = "all"
     ligA, ligB = eg5_ligands[0], eg5_ligands[1]
-    cof = _named_smc(eg5_cofactor, "MYC")
-
+    cof2 = benzene_modifications["toluene"]  # stand-in second cofactor
     mapper = openfe.setup.KartografAtomMapper()
     mapping = next(mapper.suggest_mappings(ligA, ligB))
-    stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": cof})
-    stateB = openfe.ChemicalSystem({"ligand": ligB, "cofactor": cof})
+    stateA = openfe.ChemicalSystem(
+        {"ligand": ligA, "cofactor": eg5_cofactor, "cofactor2": cof2}
+    )
+    stateB = openfe.ChemicalSystem(
+        {"ligand": ligB, "cofactor": eg5_cofactor, "cofactor2": cof2}
+    )
 
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
     u = mda.Universe(out["pdb_structure"])
-    assert "MYC" in u.residues.resnames
-    assert "CF1" not in u.residues.resnames
-    assert "LIG" in u.residues.resnames
+    cof = u.select_atoms("resname COF")
+    assert cof.n_atoms > 0
+    assert cof.residues.n_residues == 2
+    assert len(set(cof.residues.resindices)) == 2
 
 
-def test_distinct_endstate_resnames(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+def test_cofactor_clash_worked_around(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+    """A cofactor pre-named LG1 doesn't give an error; the ligand names work around it."""
     vac_settings.output_settings.output_indices = "all"
-    ligA_named = _named_smc(eg5_ligands[0], "LGA")
-    ligB_named = _named_smc(eg5_ligands[1], "LGB")
-
+    ligA, ligB = eg5_ligands[0], eg5_ligands[1]
+    cof = _named_smc(eg5_cofactor, "LG1")   # deliberately clashes with a ligand default
     mapper = openfe.setup.KartografAtomMapper()
-    mapping = next(mapper.suggest_mappings(ligA_named, ligB_named))
-    stateA = openfe.ChemicalSystem({"ligand": ligA_named, "cofactor": eg5_cofactor})
-    stateB = openfe.ChemicalSystem({"ligand": ligB_named, "cofactor": eg5_cofactor})
+    mapping = next(mapper.suggest_mappings(ligA, ligB))
+    stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": cof})
+    stateB = openfe.ChemicalSystem({"ligand": ligB, "cofactor": cof})
 
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
-    # both endstate names stored
-    assert out["alchemical_resnames"] == ["LGA", "LGB"]
+    # ligands avoid LG1 (taken by the cofactor)
+    assert out["alchemical_resnames"] == ["LG2", "LG3"]
 
     u = mda.Universe(out["pdb_structure"])
-    lig = u.select_atoms("resname LGA")
-    cof = u.select_atoms("resname CF1")
-    # hybrid ligand carries componentA's name; LGB never reaches the hybrid
+    resnames = list(u.residues.resnames)
+    assert len(resnames) == len(set(resnames))   # all residues distinct-named
+    assert "LG1" in resnames                      # the cofactor kept its name
+    lig = u.select_atoms("resname " + " ".join(out["alchemical_resnames"]))
+    cof_sel = u.select_atoms("resname LG1")
     assert lig.n_atoms > 0
-    assert u.select_atoms("resname LGB").n_atoms == 0
-    assert u.select_atoms("resname LIG").n_atoms == 0
-    assert "UNK" not in set(u.residues.resnames)
-    assert cof.n_atoms > 0
-    assert set(lig.indices).isdisjoint(cof.indices)
-    # the stored union recovers the full ligand
-    sel = u.select_atoms("resname " + " ".join(out["alchemical_resnames"]))
-    assert set(sel.indices) == set(lig.indices)
+    assert set(lig.indices).isdisjoint(cof_sel.indices)
 
-
-def test_cofactor_ligand_name_clash(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
-    """A cofactor pre-named LIG must not collide with the ligand."""
+    
+def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
+    """A pre-set cofactor name survives, not overwritten with COF."""
     vac_settings.output_settings.output_indices = "all"
     ligA, ligB = eg5_ligands[0], eg5_ligands[1]
-    cof = _named_smc(eg5_cofactor, "LIG")
-
+    cof = _named_smc(eg5_cofactor, "MYC")
     mapper = openfe.setup.KartografAtomMapper()
     mapping = next(mapper.suggest_mappings(ligA, ligB))
     stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": cof})
@@ -2514,38 +2511,30 @@ def test_cofactor_ligand_name_clash(eg5_ligands, eg5_cofactor, vac_settings, tmp
 
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
-    u = mda.Universe(out["pdb_structure"])
-    resnames = list(u.residues.resnames)
-    assert len(resnames) == len(set(resnames))  # all residues distinct
-    # pre-named cofactor keeps LIG; ligand bumps
-    assert u.select_atoms("resname LIG").n_atoms > 0  # the cofactor
-    ligand_names = out["alchemical_resnames"]
-    assert "LIG" not in ligand_names
-    lig = u.select_atoms("resname " + " ".join(ligand_names))
-    cof = u.select_atoms("resname LIG")
-    assert set(lig.indices).isdisjoint(cof.indices)
+    resnames = set(mda.Universe(out["pdb_structure"]).residues.resnames)
+    assert "MYC" in resnames        # preserved
+    assert "COF" not in resnames    # not auto-assigned over the user's name
+    assert out["alchemical_resnames"] == ["LG1", "LG2"]
 
 
-def test_structural_analysis_uses_ligand_resnames():
+def test_structural_analysis_uses_ligand_resnames(tmp_path):
+    """The analysis builds a union selection from the stored resname list."""
     captured = {}
 
     def fake_gather(*args, ligand_selection="resname UNK"):
         captured["ligand_selection"] = ligand_selection
         return {
-            "protein_RMSD": [],
-            "ligand_RMSD": [],
-            "ligand_wander": [],
-            "protein_2D_RMSD": [],
-            "time(ps)": [],
+            "protein_RMSD": [], "ligand_RMSD": [], "ligand_wander": [],
+            "protein_2D_RMSD": [], "time(ps)": [],
         }
 
     with mock.patch("openfe_analysis.rmsd.gather_rms_data", side_effect=fake_gather):
         HybridTopologyMultiStateAnalysisUnit._structural_analysis(
             pdb_file=Path("dummy.pdb"),
             trj_file=Path("dummy.nc"),
-            output_directory=Path("."),
+            output_directory=tmp_path,
             dry=True,
-            ligand_resnames=["LGA", "LGB"],
+            ligand_resnames=["LG1", "LG2"],
         )
 
-    assert captured["ligand_selection"] == "resname LGA LGB"
+    assert captured["ligand_selection"] == "resname LG1 LG2"

From 9deae8a13e22d1d2be254c875abec0b4cd697009 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 6 Jul 2026 13:23:00 +0000
Subject: [PATCH 15/19] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 .../protocols/openmm_rfe/hybridtop_units.py   |  7 ++---
 .../openmm_rfe/test_hybrid_top_protocol.py    | 27 +++++++++----------
 2 files changed, 15 insertions(+), 19 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 3dad38476..067a0453e 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -769,8 +769,7 @@ def _unique(used: set[str]) -> str:
                 if candidate not in used:
                     return candidate
             raise ValueError(
-                "Could not assign a unique ligand residue name; too many "
-                "colliding 'LG#' names."
+                "Could not assign a unique ligand residue name; too many colliding 'LG#' names."
             )
 
         # Seed with user-provided resnames so auto-assigned names avoid them.
@@ -1726,9 +1725,7 @@ def _execute(
             pdb_file=pdb_file,
             trajectory=trajectory,
             checkpoint=checkpoint,
-            ligand_resnames=setup_results.outputs.get(
-                "alchemical_resnames", ["LG1", "LG2"]
-            ),
+            ligand_resnames=setup_results.outputs.get("alchemical_resnames", ["LG1", "LG2"]),
             scratch_basepath=ctx.scratch,
             shared_basepath=ctx.shared,
         )
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index c9cae68a4..170f4e94b 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2458,12 +2458,8 @@ def test_cofactors_share_cof_distinct_resindices(
     cof2 = benzene_modifications["toluene"]  # stand-in second cofactor
     mapper = openfe.setup.KartografAtomMapper()
     mapping = next(mapper.suggest_mappings(ligA, ligB))
-    stateA = openfe.ChemicalSystem(
-        {"ligand": ligA, "cofactor": eg5_cofactor, "cofactor2": cof2}
-    )
-    stateB = openfe.ChemicalSystem(
-        {"ligand": ligB, "cofactor": eg5_cofactor, "cofactor2": cof2}
-    )
+    stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": eg5_cofactor, "cofactor2": cof2})
+    stateB = openfe.ChemicalSystem({"ligand": ligB, "cofactor": eg5_cofactor, "cofactor2": cof2})
 
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
@@ -2478,7 +2474,7 @@ def test_cofactor_clash_worked_around(eg5_ligands, eg5_cofactor, vac_settings, t
     """A cofactor pre-named LG1 doesn't give an error; the ligand names work around it."""
     vac_settings.output_settings.output_indices = "all"
     ligA, ligB = eg5_ligands[0], eg5_ligands[1]
-    cof = _named_smc(eg5_cofactor, "LG1")   # deliberately clashes with a ligand default
+    cof = _named_smc(eg5_cofactor, "LG1")  # deliberately clashes with a ligand default
     mapper = openfe.setup.KartografAtomMapper()
     mapping = next(mapper.suggest_mappings(ligA, ligB))
     stateA = openfe.ChemicalSystem({"ligand": ligA, "cofactor": cof})
@@ -2491,14 +2487,14 @@ def test_cofactor_clash_worked_around(eg5_ligands, eg5_cofactor, vac_settings, t
 
     u = mda.Universe(out["pdb_structure"])
     resnames = list(u.residues.resnames)
-    assert len(resnames) == len(set(resnames))   # all residues distinct-named
-    assert "LG1" in resnames                      # the cofactor kept its name
+    assert len(resnames) == len(set(resnames))  # all residues distinct-named
+    assert "LG1" in resnames  # the cofactor kept its name
     lig = u.select_atoms("resname " + " ".join(out["alchemical_resnames"]))
     cof_sel = u.select_atoms("resname LG1")
     assert lig.n_atoms > 0
     assert set(lig.indices).isdisjoint(cof_sel.indices)
 
-    
+
 def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
     """A pre-set cofactor name survives, not overwritten with COF."""
     vac_settings.output_settings.output_indices = "all"
@@ -2512,8 +2508,8 @@ def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp
     out = _run_setup_dry(stateA, stateB, mapping, vac_settings, tmp_path)
 
     resnames = set(mda.Universe(out["pdb_structure"]).residues.resnames)
-    assert "MYC" in resnames        # preserved
-    assert "COF" not in resnames    # not auto-assigned over the user's name
+    assert "MYC" in resnames  # preserved
+    assert "COF" not in resnames  # not auto-assigned over the user's name
     assert out["alchemical_resnames"] == ["LG1", "LG2"]
 
 
@@ -2524,8 +2520,11 @@ def test_structural_analysis_uses_ligand_resnames(tmp_path):
     def fake_gather(*args, ligand_selection="resname UNK"):
         captured["ligand_selection"] = ligand_selection
         return {
-            "protein_RMSD": [], "ligand_RMSD": [], "ligand_wander": [],
-            "protein_2D_RMSD": [], "time(ps)": [],
+            "protein_RMSD": [],
+            "ligand_RMSD": [],
+            "ligand_wander": [],
+            "protein_2D_RMSD": [],
+            "time(ps)": [],
         }
 
     with mock.patch("openfe_analysis.rmsd.gather_rms_data", side_effect=fake_gather):

From 6fdcb749affc3c6d8f48da9059a0cce75fa5d349 Mon Sep 17 00:00:00 2001
From: Hannah Baumann <43765638+hannahbaumann@users.noreply.github.com>
Date: Mon, 6 Jul 2026 15:24:27 +0200
Subject: [PATCH 16/19] Apply suggestion from @hannahbaumann

---
 news/fix_unk_resnames.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/news/fix_unk_resnames.rst b/news/fix_unk_resnames.rst
index 0dd6b8302..ac068b046 100644
--- a/news/fix_unk_resnames.rst
+++ b/news/fix_unk_resnames.rst
@@ -4,7 +4,7 @@
 
 **Changed:**
 
-* Small molecules in RelativeHybridTopology topologies (including the output PDB) are now named LIG (alchemical ligand) and CF1, CF2… (cofactors) instead of UNK. If a residue name was already assigned, the assigned one is kept.
+* Small molecules in RelativeHybridTopology topologies (including the output PDB) are now named LG1 amd LG2 (alchemical ligand) and COF (cofactors) instead of UNK. If a residue name was already assigned, the assigned one is kept.
 
 **Deprecated:**
 

From 465045e1accd9d93d046f01b7262c26d60277c96 Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Mon, 6 Jul 2026 15:35:09 +0200
Subject: [PATCH 17/19] Make mypy happy

---
 src/openfe/protocols/openmm_rfe/hybridtop_units.py           | 5 +++--
 .../tests/protocols/openmm_rfe/test_hybrid_top_protocol.py   | 3 ++-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index 3dad38476..96c71d05b 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -794,9 +794,10 @@ def _unique(used: set[str]) -> str:
                 name = "COF"
             _set_offmol_resname(offmol, name)
 
-        names = set()
+        names: set[str] = set()
         for comp in (mapping.componentA, mapping.componentB):
             name = _get_offmol_resname(small_mols[comp])
+            assert name is not None # for typing reasons, cannot be None
             names.add(name)
         alchem_resnames = sorted(names)
 
@@ -1553,7 +1554,7 @@ def _structural_analysis(
         trj_file: pathlib.Path,
         output_directory: pathlib.Path,
         dry: bool,
-        ligand_resnames: list[str] | None = None,
+        ligand_resnames: list[str],
     ) -> dict[str, str | pathlib.Path]:
         """
         Run structural analysis using ``openfe-analysis``.
diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index c9cae68a4..a4e3da554 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2340,6 +2340,7 @@ def test_structural_analysis_error(tmp_path):
         Path(tmp_path),
         Path(tmp_path),
         True,
+        ["LG1", "LG2"],
     )
 
     assert "structural_analysis_error" in ret
@@ -2498,7 +2499,7 @@ def test_cofactor_clash_worked_around(eg5_ligands, eg5_cofactor, vac_settings, t
     assert lig.n_atoms > 0
     assert set(lig.indices).isdisjoint(cof_sel.indices)
 
-    
+
 def test_existing_resname_preserved(eg5_ligands, eg5_cofactor, vac_settings, tmp_path):
     """A pre-set cofactor name survives, not overwritten with COF."""
     vac_settings.output_settings.output_indices = "all"

From e88e61b45294fe36b20913f24eb198962df4b21a Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 6 Jul 2026 13:41:34 +0000
Subject: [PATCH 18/19] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/openfe/protocols/openmm_rfe/hybridtop_units.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/openfe/protocols/openmm_rfe/hybridtop_units.py b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
index cc503ece6..d753516e4 100644
--- a/src/openfe/protocols/openmm_rfe/hybridtop_units.py
+++ b/src/openfe/protocols/openmm_rfe/hybridtop_units.py
@@ -796,7 +796,7 @@ def _unique(used: set[str]) -> str:
         names: set[str] = set()
         for comp in (mapping.componentA, mapping.componentB):
             name = _get_offmol_resname(small_mols[comp])
-            assert name is not None # for typing reasons, cannot be None
+            assert name is not None  # for typing reasons, cannot be None
             names.add(name)
         alchem_resnames = sorted(names)
 

From de270f0549ffd8767eb642da02ca750d2c805e59 Mon Sep 17 00:00:00 2001
From: hannahbaumann <hbaumann@uci.edu>
Date: Mon, 6 Jul 2026 16:49:04 +0200
Subject: [PATCH 19/19] Small fix

---
 .../tests/protocols/openmm_rfe/test_hybrid_top_protocol.py  | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
index 8518ed2be..32f14da00 100644
--- a/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
+++ b/src/openfe/tests/protocols/openmm_rfe/test_hybrid_top_protocol.py
@@ -2427,11 +2427,9 @@ def _run_setup_dry(stateA, stateB, mapping, settings, tmp_path):
 
 
 def _named_smc(smc, resname):
-    off = copy.deepcopy(smc.to_openff())
+    off = Molecule.from_rdkit(copy.deepcopy(smc.to_rdkit()))
     _set_offmol_resname(off, resname)
-    result = openfe.SmallMoleculeComponent.from_openff(off)
-    assert _get_offmol_resname(result.to_openff()) == resname
-    return result
+    return openfe.SmallMoleculeComponent.from_openff(off)
 
 
 def test_default_resnames(eg5_vac_inputs, vac_settings, tmp_path):