Boresch restraint search crashes on ligands with three or more fused aromatic rings (`get_aromatic_rings`: `list.remove(x): x not in list`)

Start with openfe/protocols/restraint_utils/geometry/utils.py and the get_aromatic_rings function, then run the minimal RDKit reproducer from the issue. Add the suggested regression cases and verify that fused aromatic rings produce the expected systems without crashing, while separate rings remain separate.

Summary

openfe.protocols.restraint_utils.geometry.utils.get_aromatic_rings raises
ValueError: list.remove(x): x not in list for any molecule with three or more fused
aromatic rings, such as anthracene, phenanthrene or carbazole.

It is called during automatic Boresch guest-atom selection, so for these ligands the
AbsoluteBindingProtocol complex-leg setup unit fails. The failure only appears once
execution has started: the protocol validates, the solvent leg runs, and the complex leg
then fails at setup. SepTopProtocol uses the same restraint search, so it is likely
affected as well; I have not run it to confirm.

Environment

• openfe 1.12.0, gufe 1.12.0, RDKit 2025.09.3, Python 3.13, Linux
• The same code is still present on main as of 2026-09-18

Minimal reproducer

No input files are needed:

from rdkit import Chem
from openfe.protocols.restraint_utils.geometry.utils import get_aromatic_rings

for name, smiles in {
    "benzene": "c1ccccc1",
    "naphthalene": "c1ccc2ccccc2c1",
    "anthracene": "c1ccc2cc3ccccc3cc2c1",
    "phenanthrene": "c1ccc2c(c1)ccc1ccccc12",
}.items():
    try:
        rings = get_aromatic_rings(Chem.MolFromSmiles(smiles))
        print(f"OK     {name}: {len(rings)} ring system(s)")
    except ValueError as error:
        print(f"CRASH  {name}: {error}")

Output:

OK     benzene: 1 ring system(s)
OK     naphthalene: 1 ring system(s)
CRASH  anthracene: list.remove(x): x not in list
CRASH  phenanthrene: list.remove(x): x not in list

Traceback from an ABFE run

From openfe quickrun on an AbsoluteBindingProtocol transformation whose ligand has a
fused tricyclic aromatic core:

Error: The protocol unit 'ABFE Setup: <ligand> complex leg: repeat 0 generation 0' failed with the error message:
ValueError: list.remove(x): x not in list

  File ".../openfe/protocols/openmm_afe/abfe_units.py", line 252, in _get_boresch_restraint
    geom = geometry.boresch.find_boresch_restraint(
  File ".../openfe/protocols/restraint_utils/geometry/boresch/geometry.py", line 246, in find_boresch_restraint
    guest_anchors = find_guest_atom_candidates(
  File ".../openfe/protocols/restraint_utils/geometry/boresch/guest.py", line 217, in find_guest_atom_candidates
    atom_pool, rings_only = _get_guest_atom_pool(rdmol, rmsf, rmsf_cutoff)
  File ".../openfe/protocols/restraint_utils/geometry/boresch/guest.py", line 157, in _get_guest_atom_pool
    for ring in get_aromatic_rings(rdmol):
  File ".../openfe/protocols/restraint_utils/geometry/utils.py", line 105, in get_aromatic_rings
    aromatic_rings.remove(y)
ValueError: list.remove(x): x not in list

Expected behaviour

Fused aromatic rings are merged into one ring system — anthracene gives a single
14-atom system, as naphthalene already gives a single 10-atom one — and restraint
selection continues.

Cause

The merge loop iterates over combinations(aromatic_rings, 2) while removing from
aromatic_rings:

for x, y in combinations(aromatic_rings, 2):
    if not x.isdisjoint(y):
        x.update(y)
        aromatic_rings.remove(y)

itertools.combinations takes a snapshot of the list when it is created, so it keeps
yielding pairs containing rings that have already been merged and removed. For three
linearly fused rings A–B–C:

1. (A, B) overlap: B is merged into A and removed.
2. (A, C) overlap, since A now contains B's atoms: C is merged into A and removed.
3. (B, C) still overlap, but C has already been removed, so remove(C) raises.

Any molecule where three or more aromatic rings form one connected fused system reaches
step 3. Two fused rings never do, which is why naphthalene works.

Proposed fix

Merge each ring into every existing group it touches, building connected components
without mutating the sequence being iterated:

def get_aromatic_rings(rdmol: Chem.Mol) -> list[set[int]]:
    ringinfo = rdmol.GetRingInfo()
    arom_idxs = get_aromatic_atom_idxs(rdmol)

    aromatic_rings = [
        set(ring) for ring in ringinfo.AtomRings() if all(a in arom_idxs for a in ring)
    ]

    # Merge rings that share atoms into fused ring systems
    merged: list[set[int]] = []
    for ring in aromatic_rings:
        overlapping = [group for group in merged if not group.isdisjoint(ring)]
        for group in overlapping:
            merged.remove(group)
            ring |= group
        merged.append(ring)

    return merged

Checked against these cases, all of which give the expected grouping:

A suggested regression test:

@pytest.mark.parametrize(
    "smiles, sizes",
    [
        ("c1ccc2cc3ccccc3cc2c1", [14]),       # anthracene
        ("c1ccc2c(c1)ccc1ccccc12", [14]),     # phenanthrene
        ("c1cc2ccc3cccc4ccc(c1)c2c34", [16]), # pyrene
        ("c1ccc(cc1)-c1ccccc1", [6, 6]),      # biphenyl stays two systems
    ],
)
def test_get_aromatic_rings_fused(smiles, sizes):
    rings = get_aromatic_rings(Chem.MolFromSmiles(smiles))
    assert sorted(len(r) for r in rings) == sizes

Impact

Fused tricyclic and larger aromatic cores are common in drug-like molecules — carbazoles,
acridines, phenanthridines, fused quinoxalines — and none of these can currently be run
with automatic restraint selection.

I have not checked whether the unreleased user-defined Boresch restraint support (#2019)
bypasses this code path, which would give a workaround once released.