KeyError when removing constraints from mapping

[jthorton]

Describe the bug

While trying to run a RBFE calculation on some MALT1 ligands I get the following error when trying to remove parts of the mapping that involve a change in constraints:

Traceback (most recent call last):
File "/Users/joshua/Documents/Projects/openfe/dev/femto_benchmarking/malt1/reproduce_error.py", line 32, in
vac_unit.run(dry=True, verbose=True)
File "/Users/joshua/Documents/Software/openfe/openfe/protocols/openmm_rfe/equil_rfe_methods.py", line 844, in run
ligand_mappings = _rfe_utils.topologyhelpers.get_system_mappings(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/joshua/Documents/Software/openfe/openfe/protocols/openmm_rfe/_rfe_utils/topologyhelpers.py", line 599, in get_system_mappings
adjusted_old_to_new_map = _remove_constraints(
^^^^^^^^^^^^^^^^^^^^
File "/Users/joshua/Documents/Software/openfe/openfe/protocols/openmm_rfe/_rfe_utils/topologyhelpers.py", line 513, in _remove_constraints
del no_const_old_to_new_atom_map[idx]
~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^
KeyError: 34

I have narrowed the cause of this down to the index 34 being added to the to_del list twice. The overlapping ligands are shown here (the pymol index is +1 on the mapping index):

and the atom mapping made by Kartograf is

{26: 28, 27: 30, 28: 29, 29: 31, 30: 27, 31: 32, 32: 38, 34: 13, 39: 40, 40: 41, 0: 0, 1: 2, 2: 1, 3: 3, 4: 6, 5: 4, 6: 5, 7: 7, 8: 8, 9: 9, 10: 10, 11: 11, 12: 23, 14: 14, 15: 15, 16: 16, 17: 17, 18: 18, 19: 19, 20: 20, 21: 21, 22: 22}

I think this mapping may have an issue (so maybe a kartograf issue?) as 32(H) -> 38(H) means that a H will be left on the ether oxygen at lambda=1 which is not correct, and this issue goes away if you set atom_map_hydrogens=False in the script below, which we intend to make the default going forward and should resolve this.

To Reproduce
Steps to reproduce the behavior (ideally a minimally reproducible example):
See the attached ligand sdf files and this script:

from gufe import Transformation, SmallMoleculeComponent

from kartograf import KartografAtomMapper
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol
from openfe import ChemicalSystem

lig1 = SmallMoleculeComponent.from_sdf_file("shapefit_Pfizer-01-01.sdf")
lig2 = SmallMoleculeComponent.from_sdf_file("shapefit_1832577-09-9.sdf")

settings = RelativeHybridTopologyProtocol._default_settings()
settings.forcefield_settings.nonbonded_method = "nocutoff"
protocol = RelativeHybridTopologyProtocol(settings)

atom_mapper = KartografAtomMapper(atom_map_hydrogens=False)

mapping = next(atom_mapper.suggest_mappings(lig1, lig2))

print(mapping.componentA_to_componentB)
print(mapping.componentA.name, mapping.componentB.name)
sys_a_dict = {"ligand": mapping.componentA}
sys_b_dict = {"ligand": mapping.componentB}
system_a = ChemicalSystem(sys_a_dict, name=f"{mapping.componentA.name}_vac")
system_b = ChemicalSystem(sys_b_dict, name=f"{mapping.componentB.name}_vac")
transform = Transformation(
    stateA=system_a,
    stateB=system_b,
    mapping={"ligand": mapping},
    protocol=protocol
)
vac_leg_dag = transform.create()
vac_unit = list(vac_leg_dag.protocol_units)[0]
vac_unit.run(dry=True, verbose=True)

ligands.zip

Software versions

• Which operating system and version did you use? (e.g. ubuntu 22.04.5) mac
• Which method did you use to install this package? (e.g. conda-forge) conda-forge
• Copy/paste the output of conda list (or the equivalent for your package manager):

Output

Expected behavior

Maybe we should raise a more descriptive error if we hit this case where an index is in the to delete list more than once, I think this generally indicates that an atom in a constraint is being transformed into another atom involved in a different constraint.

Additional context