AtomMapper fails when scaffolds have radicals

[ijpulidos]

While trying to run the new edges in the https://github.com/openforcefield/protein-ligand-benchmark/tree/new-perses-edges dataset I noticed that at least one edge is failing when creating the mapping with the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 2
      1 atom_mapper = atom_mapping.AtomMapper()
----> 2 atom_mapping = atom_mapper.get_best_mapping(old_molecule, new_molecule)

File ~/workdir/repos/perses/perses/rjmc/atom_mapping.py:993, in AtomMapper.get_best_mapping(self, old_mol, new_mol)
    990 initial_time = time.time()
    992 import numpy as np
--> 993 atom_mappings = self.get_all_mappings(old_mol, new_mol)
    994 if (atom_mappings is None) or len(atom_mappings)==0:
    995     return None

File ~/workdir/repos/perses/perses/rjmc/atom_mapping.py:824, in AtomMapper.get_all_mappings(self, old_mol, new_mol)
    820 else:
    821     # Generate scaffold maps
    822     # TODO: Why are these hard-coded?
    823     from openeye import oechem
--> 824     scaffold_maps = AtomMapper._get_all_maps(old_oescaffold, new_oescaffold,
    825                                              atom_expr=oechem.OEExprOpts_RingMember | oechem.OEExprOpts_IntType,
    826                                              bond_expr=oechem.OEExprOpts_RingMember,
    827                                              external_inttypes=True,
    828                                              unique=False,
    829                                              matching_criterion=self.matching_criterion)
    831     _logger.debug(f'Scaffold mapping generated {len(scaffold_maps)} maps')
    833 if len(scaffold_maps) == 0:
    834     # There are no scaffold maps, so attempt to generate maps between molecules using the factory parameters

File ~/workdir/repos/perses/perses/rjmc/atom_mapping.py:1329, in AtomMapper._get_all_maps(old_oemol, new_oemol, external_inttypes, atom_expr, bond_expr, unique, matching_criterion)
   1327 for match in matches:
   1328     try:
-> 1329         atom_mapping = AtomMapper._create_atom_mapping(old_oemol, new_oemol, match, matching_criterion)
   1330         atom_mappings.add(atom_mapping)
   1331     except InvalidMappingException as e:
   1332         # Mapping is not valid; skip it

File ~/workdir/repos/perses/perses/rjmc/atom_mapping.py:1426, in AtomMapper._create_atom_mapping(old_oemol, new_oemol, match, matching_criterion)
   1422             continue
   1424     new_to_old_atom_map[new_index] = old_index
-> 1426 return AtomMapping(old_oemol, new_oemol, new_to_old_atom_map=new_to_old_atom_map)

File ~/workdir/repos/perses/perses/rjmc/atom_mapping.py:102, in AtomMapping.__init__(self, old_mol, new_mol, new_to_old_atom_map, old_to_new_atom_map)
     99 # Store molecules
    100 # TODO: Should we be allowing unspecified stereochemistry?
    101 from openff.toolkit.topology import Molecule
--> 102 self.old_mol = Molecule(old_mol, allow_undefined_stereo=True)
    103 self.new_mol = Molecule(new_mol, allow_undefined_stereo=True)
    105 # Fix atom name retrieval from OEMol objects
    106 # TODO: When https://github.com/openforcefield/openff-toolkit/issues/1026 is fixed, remove this

File ~/miniconda3/envs/perses-0.10.2/lib/python3.11/site-packages/openff/toolkit/topology/molecule.py:5193, in Molecule.__init__(self, *args, **kwargs)
   5178 def __init__(self, *args, **kwargs):
   5179     """
   5180     See FrozenMolecule.__init__
   5181 
   (...)
   5191 
   5192     """
-> 5193     super(Molecule, self).__init__(*args, **kwargs)

File ~/miniconda3/envs/perses-0.10.2/lib/python3.11/site-packages/openff/toolkit/topology/molecule.py:1045, in FrozenMolecule.__init__(self, other, file_format, toolkit_registry, allow_undefined_stereo)
   1043 for value_error in value_errors:
   1044     msg += str(value_error)
-> 1045 raise ValueError(msg)

ValueError: Cannot construct openff.toolkit.topology.Molecule from <openeye.oechem.OEMol; proxy of <Swig Object of type 'OEMolWrapper *' at 0x7f568e9683c0> >
No registered toolkits can provide the capability "from_object" for args "(<openeye.oechem.OEMol; proxy of <Swig Object of type 'OEMolWrapper *' at 0x7f568e9683c0> >,)" and kwargs "{'allow_undefined_stereo': True, '_cls': <class 'openff.toolkit.topology.molecule.Molecule'>}"
Available toolkits are: [ToolkitWrapper around OpenEye Toolkit version 2023.1.0, ToolkitWrapper around The RDKit version 2023.03.2, ToolkitWrapper around AmberTools version 22.0, ToolkitWrapper around Built-in Toolkit version None]
 ToolkitWrapper around OpenEye Toolkit version 2023.1.0 <class 'openff.toolkit.utils.exceptions.RadicalsNotSupportedError'> : The OpenFF Toolkit does not currently support parsing molecules with radicals. Found 1 radical electrons on molecule [c]1[c][c]c([c][c]1)[N]c2nc3c(c(n2)O[C][C]4[C][C][C][C][C]4)N=[C][N]3.
 ToolkitWrapper around The RDKit version 2023.03.2 <class 'NotImplementedError'> : Cannot create Molecule from <class 'openeye.oechem.OEMol'> object

I wonder if we have to re-evalueate how we are generating the smiles when doing the mapping such that radicals are not generated.

Lastly I just wanted to know that this is not a problem with openff-toolkit 0.10.6. The toolkit stopped supporting radicals from 0.11 forward.

The code to reproduce is as follows

from perses.rjmc import atom_mapping
from openff.toolkit.topology import Molecule

mols = Molecule.from_file("ligands.sdf")
old_molecule = mols[2]
new_molecule = mols[3]
atom_mapper = atom_mapping.AtomMapper()
atom_mapping = atom_mapper.get_best_mapping(old_molecule, new_molecule)

The ligands.sdf file is attached as a tar.gz file here.
ligands.sdf.tar.gz

[ijpulidos]

Related to openforcefield/protein-ligand-benchmark#93 (comment)