Merge pull request #123 from openforcefield/stereo_fix

Fix #118

openff/fragmenter/fragment.py

@@ -61,7 +61,7 @@ class Fragment(BaseModel):
     @property
     def molecule(self) -> Molecule:
         """The fragment represented as an OpenFF molecule object."""
-        return Molecule.from_smiles(self.smiles)
+        return Molecule.from_smiles(self.smiles, allow_undefined_stereo=True)
 
 
 class FragmentationResult(BaseModel):
@@ -217,7 +217,7 @@ def _check_stereo(
     @classmethod
     def _fix_stereo(
         cls, fragment: Molecule, parent_stereo: Stereochemistries
-    ) -> Molecule:
+    ) -> Optional[Molecule]:
         """Flip all stereocenters and find the stereoisomer that matches the parent
 
         Parameters
@@ -229,7 +229,7 @@ def _fix_stereo(
 
         Returns
         -------
-            A fragment with the same stereochemistry as parent molecule
+            A fragment with the same stereochemistry as parent molecule if possible else ``None``.
         """
 
         for stereoisomer in _enumerate_stereoisomers(fragment, 200, True):
@@ -239,9 +239,7 @@ def _fix_stereo(
 
             return stereoisomer
 
-        raise RuntimeError(
-            f"The stereochemistry of {fragment.to_smiles()} could not be fixed."
-        )
+        return None
 
     @classmethod
     def _find_functional_groups(
@@ -370,7 +368,7 @@ def _atom_bond_set_to_mol(
         parent_stereo: Stereochemistries,
         atoms: Set[int],
         bonds: Set[BondTuple],
-    ) -> Molecule:
+    ) -> Tuple[Molecule, bool]:
         """Extracts a subset of a molecule based on a set of atom and bond indices.
 
         Parameters
@@ -387,15 +385,19 @@ def _atom_bond_set_to_mol(
 
         Returns
         -------
-            The subset molecule.
+            The subset molecule and a flag to indicate whether any new stereocenters are present in the fragment.
         """
 
         fragment = extract_fragment(parent, atoms, bonds)
 
         if not cls._check_stereo(fragment, parent_stereo):
-            fragment = cls._fix_stereo(fragment, parent_stereo)
+            fixed_fragment = cls._fix_stereo(fragment, parent_stereo)
+            if fixed_fragment is None:
+                return fragment, True
+            else:
+                return fixed_fragment, False
 
-        return fragment
+        return fragment, False
 
     @classmethod
     def _get_torsion_quartet(
@@ -1186,13 +1188,20 @@ def _build_fragment(
         if cap:
             atoms, bonds = cls._cap_open_valence(parent, parent_groups, atoms, bonds)
 
-        fragment = cls._atom_bond_set_to_mol(parent, parent_stereo, atoms, bonds)
+        fragment, has_new_stereocenter = cls._atom_bond_set_to_mol(
+            parent, parent_stereo, atoms, bonds
+        )
 
-        wbo_difference = cls._compare_wbo(fragment, bond_tuple, parent_wbo, **kwargs)
+        if has_new_stereocenter:
+            wbo_difference = threshold + 1.0
+        else:
+            wbo_difference = cls._compare_wbo(
+                fragment, bond_tuple, parent_wbo, **kwargs
+            )
 
         while fragment is not None and wbo_difference > threshold:
 
-            fragment = cls._add_next_substituent(
+            fragment, has_new_stereocenter = cls._add_next_substituent(
                 parent,
                 parent_stereo,
                 parent_groups,
@@ -1206,16 +1215,21 @@ def _build_fragment(
             if fragment is None:
                 break
 
-            wbo_difference = cls._compare_wbo(
-                fragment, bond_tuple, parent_wbo, **kwargs
-            )
+            if has_new_stereocenter:
+                # For now keep growing the fragment as it is not yet clear how to handle such cases robustly.
+                wbo_difference = threshold + 1.0
+
+            else:
+                wbo_difference = cls._compare_wbo(
+                    fragment, bond_tuple, parent_wbo, **kwargs
+                )
 
         # A work around for a known bug where if stereochemistry changes or gets removed,
         # the WBOs can change more than the threshold (this will sometimes happen if a
         # very small threshold is chosen) and even the parent will have a WBO difference
         # greater than the threshold. In this case, return the molecule
         if fragment is None:
-            fragment = cls._atom_bond_set_to_mol(parent, parent_stereo, atoms, bonds)
+            fragment, _ = cls._atom_bond_set_to_mol(parent, parent_stereo, atoms, bonds)
 
         return fragment
 
@@ -1357,7 +1371,7 @@ def _add_next_substituent(
         bonds: Set[BondTuple],
         target_bond: BondTuple,
         heuristic: Heuristic = "path_length",
-    ) -> Optional[Molecule]:
+    ) -> Tuple[Optional[Molecule], bool]:
         """Expand the fragment to include the next set of substituents / ring systems.
 
         Parameters
@@ -1401,7 +1415,7 @@ def _add_next_substituent(
             )
 
         if neighbour_atom_and_bond is None:
-            return None
+            return None, False
 
         neighbour_atom, neighbour_bond = neighbour_atom_and_bond
 
@@ -1463,7 +1477,7 @@ def _fragment(
 
             atoms, bonds = self._cap_open_valence(parent, parent_groups, atoms, bonds)
 
-            fragments[bond] = self._atom_bond_set_to_mol(
+            fragments[bond], _ = self._atom_bond_set_to_mol(
                 parent, parent_stereo, atoms, bonds
             )
         return FragmentationResult(

openff/fragmenter/tests/test_fragmenter.py

@@ -207,7 +207,9 @@ def test_atom_bond_set_to_mol(abemaciclib):
         and get_map_index(molecule, bond.atom2_index) in atoms
     }
 
-    fragment = Fragmenter._atom_bond_set_to_mol(molecule, {}, atoms=atoms, bonds=bonds)
+    fragment, _ = Fragmenter._atom_bond_set_to_mol(
+        molecule, {}, atoms=atoms, bonds=bonds
+    )
 
     for bond in fragment.bonds:
 
@@ -600,7 +602,7 @@ def test_add_substituent():
         if bond.atom1.atomic_number != 1 and bond.atom2.atomic_number != 1
     )
 
-    fragment = WBOFragmenter._add_next_substituent(
+    fragment, _ = WBOFragmenter._add_next_substituent(
         result.parent_molecule, {}, {}, {}, atoms, bonds, target_bond=(3, 5)
     )
 

openff/fragmenter/tests/test_regression.py

@@ -43,3 +43,15 @@ def test_wbo_regression(parent, fragments):
     for fragment in fragments:
         fragment_mol = Molecule.from_mapped_smiles(fragment)
         assert fragment_mol in result_fragments
+
+
+def test_stereo_error():
+    """
+    Make sure the molecule can be fragmented, before PR123 this molecule would result in a runtime error as new
+    stereocenters would be created.
+    """
+    parent = Molecule.from_smiles(
+        "[H][O][C]([H])([H])[C]([H])([H])[C]([H])([N]([H])[C]1=[N][C]([H])=[C]2[C](=[N]1)[N]([C]([H])([H])[H])[C](=[O])[C]([O][C]1=[C]([F])[C]([H])=[C]([F])[C]([H])=[C]1[H])=[C]2[H])[C]([H])([H])[C]([H])([H])[O][H]"
+    )
+    engine = WBOFragmenter()
+    _ = engine.fragment(parent)