UndefinedStereochemistryError for achiral molecule

[kecherry]

Hello, I am trying to load the molecule PFOA from an sdf file but am getting an UndefindedStereochemistryError despite the molecule being achiral. When I load it from a SMILES string, I don't get this error, but I need the coordinates to be correct so that I can later load it back into my protein/solvent system in the correct docked orientation. If you have any insight as to how to fix this or what might be going wrong, that would be much appreciated!

These are the ways I have tried to load it from the sdf file and from SMILES:
ligand = Molecule.from_file(PFOA.sdf)
or
ligand = Molecule.from_smiles("O=C(O)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)F")

Here is the error message I get from loading from the sdf file:

---------------------------------------------------------------------------
UndefinedStereochemistryError             Traceback (most recent call last)
Cell In[51], line 1
----> 1 ligand = Molecule.from_file(LIG_path)
      3 # Print out a SMILES code for the ligand
      4 print(ligand.to_smiles(explicit_hydrogens=False))

File ~/miniforge3/envs/openff_env/lib/python3.12/site-packages/openff/toolkit/topology/molecule.py:3967, in FrozenMolecule.from_file(cls, file_path, file_format, toolkit_registry, allow_undefined_stereo)
   3965     if isinstance(file_path, pathlib.Path):
   3966         file_path = file_path.as_posix()
-> 3967     mols = toolkit.from_file(  # type: ignore[call-arg]
   3968         file_path,
   3969         file_format=file_format,
   3970         allow_undefined_stereo=allow_undefined_stereo,
   3971         _cls=cls,
   3972     )
   3973 elif hasattr(file_path, "read"):
   3974     file_obj = file_path

File ~/miniforge3/envs/openff_env/lib/python3.12/site-packages/openff/toolkit/utils/rdkit_wrapper.py:1075, in RDKitToolkitWrapper.from_file(self, file_path, file_format, allow_undefined_stereo, _cls)
   1071 if (file_format == "MOL") or (file_format == "SDF"):
   1072     sdf_supplier = Chem.ForwardSDMolSupplier(
   1073         file_path, removeHs=False, sanitize=False, strictParsing=True
   1074     )
-> 1075     mols.extend(
   1076         self._process_sdf_supplier(
   1077             sdf_supplier,
   1078             allow_undefined_stereo=allow_undefined_stereo,
   1079             _cls=_cls,
   1080         )
   1081     )
   1083 elif file_format == "SMI":
   1084     # TODO: We have to do some special stuff when we import SMILES (currently
   1085     # just adding H's, but could get fancier in the future). It might be
   1086     # worthwhile to parse the SMILES file ourselves and pass each SMILES
   1087     # through the from_smiles function instead
   1088     for rdmol in Chem.SmilesMolSupplier(file_path, titleLine=False):

File ~/miniforge3/envs/openff_env/lib/python3.12/site-packages/openff/toolkit/utils/rdkit_wrapper.py:1028, in RDKitToolkitWrapper._process_sdf_supplier(self, sdf_supplier, allow_undefined_stereo, _cls)
   1026     continue
   1027 Chem.SetAromaticity(rdmol, Chem.AromaticityModel.AROMATICITY_MDL)
-> 1028 mol = self.from_rdkit(
   1029     rdmol, allow_undefined_stereo=allow_undefined_stereo, _cls=_cls
   1030 )
   1031 yield mol

File ~/miniforge3/envs/openff_env/lib/python3.12/site-packages/openff/toolkit/utils/rdkit_wrapper.py:2307, in RDKitToolkitWrapper.from_rdkit(self, rdmol, allow_undefined_stereo, hydrogens_are_explicit, _cls)
   2305 # Check for undefined stereochemistry.
   2306 if not allow_undefined_stereo:
-> 2307     self._detect_undefined_stereo(
   2308         rdmol,
   2309         err_msg_prefix="Unable to make OFFMol from RDMol: ",
   2310     )
   2312 # Create a new OpenFF Molecule
   2313 offmol = _cls()

File ~/miniforge3/envs/openff_env/lib/python3.12/site-packages/openff/toolkit/utils/rdkit_wrapper.py:3273, in RDKitToolkitWrapper._detect_undefined_stereo(cls, rdmol, err_msg_prefix)
   3267         atom1, atom2 = bond.GetBeginAtom(), bond.GetEndAtom()
   3268         msg += (
   3269             f" - Bond {undefined_bond_idx} (atoms {atom1.GetIdx()}-{atom2.GetIdx()} of element "
   3270             "({atom1.GetSymbol()}-{atom2.GetSymbol()})\n"
   3271         )
-> 3273 raise UndefinedStereochemistryError(err_msg_prefix + msg)

UndefinedStereochemistryError: Unable to make OFFMol from RDMol: RDMol has unspecified stereochemistry. RDMol name: obj01Undefined chiral centers are:
 - Atom C (index 5)

This is the sdf file (sorry, it wouldn't let me upload the file so I copied and pasted it!):

obj01
  PyMOL3.0          3D                             0

 26 25  0  0  0  0  0  0  0  0999 V2000
   38.2630   37.1970   33.2080 C   0  0  0  0  0  0  0  0  0  0  0  0
   36.3070   37.6520   36.2500 O   0  0  0  0  0  0  0  0  0  0  0  0
   38.8130   37.7400   34.5440 C   0  0  0  0  0  0  0  0  0  0  0  0
   37.7950   37.9900   37.9250 O   0  0  0  0  0  0  0  0  0  0  0  0
   38.6110   36.8750   35.8120 C   0  0  0  0  0  0  0  0  0  0  0  0
   37.4470   37.5170   36.6640 C   0  0  0  0  0  0  0  0  0  0  0  0
   38.3530   38.0600   31.9360 C   0  0  0  0  0  0  0  0  0  0  0  0
   38.1980   37.3530   30.5740 C   0  0  0  0  0  0  0  0  0  0  0  0
   37.8750   38.2530   29.3660 C   0  0  0  0  0  0  0  0  0  0  0  0
   38.0710   37.6500   27.9660 C   0  0  0  0  0  0  0  0  0  0  0  0
   38.8740   36.0340   32.9270 F   0  0  0  0  0  0  0  0  0  0  0  0
   37.2310   36.3990   30.6640 F   0  0  0  0  0  0  0  0  0  0  0  0
   36.5940   38.6580   29.5150 F   0  0  0  0  0  0  0  0  0  0  0  0
   38.5940   39.3940   29.4380 F   0  0  0  0  0  0  0  0  0  0  0  0
   37.5170   36.4480   27.8620 F   0  0  0  0  0  0  0  0  0  0  0  0
   37.5500   38.3830   26.9750 F   0  0  0  0  0  0  0  0  0  0  0  0
   39.3680   37.4900   27.6740 F   0  0  0  0  0  0  0  0  0  0  0  0
   36.9770   36.8390   33.3940 F   0  0  0  0  0  0  0  0  0  0  0  0
   38.3060   38.9830   34.7490 F   0  0  0  0  0  0  0  0  0  0  0  0
   40.1440   37.9530   34.3820 F   0  0  0  0  0  0  0  0  0  0  0  0
   39.7600   36.7340   36.5240 F   0  0  0  0  0  0  0  0  0  0  0  0
   38.2460   35.5980   35.4990 F   0  0  0  0  0  0  0  0  0  0  0  0
   37.4360   39.0340   32.0480 F   0  0  0  0  0  0  0  0  0  0  0  0
   39.5290   38.7220   31.8850 F   0  0  0  0  0  0  0  0  0  0  0  0
   39.3460   36.6580   30.3310 F   0  0  0  0  0  0  0  0  0  0  0  0
   38.7210   37.8610   38.1800 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  3  1  0  0  0  0
  1  7  1  0  0  0  0
  1 11  1  0  0  0  0
  1 18  1  0  0  0  0
  3  5  1  0  0  0  0
  3 19  1  0  0  0  0
  3 20  1  0  0  0  0
  4 26  1  0  0  0  0
  5  6  1  0  0  0  0
  5 21  1  0  0  0  0
  5 22  1  0  0  0  0
  2  6  1  0  0  0  0
  4  6  1  0  0  0  0
  7  8  1  0  0  0  0
  7 23  1  0  0  0  0
  7 24  1  0  0  0  0
  8  9  1  0  0  0  0
  8 12  1  0  0  0  0
  8 25  1  0  0  0  0
  9 10  1  0  0  0  0
  9 13  1  0  0  0  0
  9 14  1  0  0  0  0
 10 15  1  0  0  0  0
 10 16  1  0  0  0  0
 10 17  1  0  0  0  0
M  END
$$$$

[mattwthompson]

Looks like you might have found a solution, but just a few notes in passing

• We appreciate the thorough reproduction!
• Passing allow_undefined_stereo=True allows molecules like these to be loaded, losing only stereochemistry
  • You probably don't need stereochemistry, something of a "if you know you need it, you probably do"
• The coordinates should be read precisely even if the stereochemistry is not wholly preserved

[j-wags]

Seconded, thanks for the thorough writeup. I was unable to reproduce the error (both input methods succeeded) but feel free to write in again or reopen this discussion if you run into it again.

[mattwthompson]

<100% sure without re-running but I recall this being a file which RDKit and OEChem (% our wrappers' quirks) disagreed on

[kecherry]

Thanks for the follow up! I am still not sure where the error came from, but I was able to fix it. Appreciate the help!

[j-wags]

Oh, I actually did reproduce it with Matt's clarification that it only reproduces with RDKit. The input SDF is a different molecule than the SMILES.

The issue seems to be that the input SDF has the wrong bond orders (everything is bond order 1, when the intent seems to be to have a double bond to one of the oxygens)

...
  4 26  1  0  0  0  0 <-- Bond between O 4 and H26 (correct)
  5  6  1  0  0  0  0
  5 21  1  0  0  0  0
  5 22  1  0  0  0  0
  2  6  1  0  0  0  0 <-- Bond between C 6 and O 2 (wrong, should be double)
  4  6  1  0  0  0  0 <-- Bond between C 6 and O 4 (correct, should be single)
  7  8  1  0  0  0  0
  7 23  1  0  0  0  0
  7 24  1  0  0  0  0
  8  9  1  0  0  0  0
...

What happens when the SDF with all-1 bond orders gets loaded is that RDKit assumes there's an implicit hydrogen on the C 6 to fill its valence. Here's code to see that in a jupyter notebook:

from openff.toolkit import Molecule
mol = Molecule.from_file('input.sdf', allow_undefined_stereo=True)
mol.visualize(backend='rdkit')

I've circled the newly-added Hs in the output:

So loading that file with allow_undefined_stereo=True kinda masks the real issue, which is that this molecule is chemically different than the one in the SMILES (two Hs are added).

Modifying the SDF to have the correct double bond gets it to load correctly without issue:

...
  5 22  1  0  0  0  0
  2  6  2  0  0  0  0 <-- Changed to double bond
  4  6  1  0  0  0  0
...

mol = Molecule.from_file('input_fixed.sdf')
mol.visualize(backend='rdkit')

It seems like something along the pathway to making that SDF lost track of the bond order to that oxygen, maybe check your workflow or consider other tools that will preserve it.

[kecherry]

Awesome, thanks so much for checking this out! I did make the SDF from a protein-ligand docked structure that was originally a PDB, so I assume that's where the bond order was lost. I'll make sure to edit that bond manually whenever I use this workflow.