Inferring stereochemistry

[orionarcher]

Is there any way in Openff toolkit to infer stereochemistry from molecule coordinates?

If not, are there other tools the devs would recommend that can do this?

Ideally I'd be able to infer stereo labels like E/Z, cis-trans directly from the molecular graph + coordinates, which I am extracting from DFT.

[j-wags]

Not directly in the OpenFF toolkit - We have to be careful not to cross the line into cheminformatics-land, since that's a huge behavior surface area to maintain.

What I'd recommend here is taking a quick round trip with to_rdkit and from_rdkit. I'll put together a code snippet when I have a minute, but roughly something like this, where the important part is telling RDKit to assign stereo from 3D.

[j-wags]

Here you go:

from openff.toolkit import Molecule
from rdkit import Chem

## Make a graph mol with stereo but no conformers

mol = Molecule.from_smiles('[C@H](F)(Cl)\C=C/C')
mol.visualize(backend="rdkit")

# Generate a conformer with the correct stereo and make a copy of it
mol.generate_conformers()
conf = mol.conformers[0]
mol.visualize(backend='nglview')

# Remove the molecule's stereo (sorry, this is hacky but I'm just doing this to get to the starting point of your use case)
mol._invalidate_cached_properties()
for atom in mol.atoms:
    atom.stereochemistry = None
for bond in mol.bonds:
    bond._stereochemistry = None
    
# This removes the molecule's stereo information (and conformers)
mol.visualize(backend="rdkit")

# To verify - Trying to roundtrip this mol leads to an undefinedstereoerror since there's no graph stereo any more
Molecule.from_rdkit(mol.to_rdkit())

UndefinedStereochemistryError: Unable to make OFFMol from RDMol: RDMol has unspecified stereochemistry. RDMol name:
Undefined chiral centers are:

• Atom C (index 0)
  Bonds with undefined stereochemistry are:
• Bond 4 (atoms 4-5 of element (C-C)

# Add the conformer with the correct stereo back, but note that even though the stereo is present in the 3D, OpenFF and RDKit haven't updated the graph representation.
mol.add_conformer(conf)
mol.visualize(backend="rdkit")

# This should be the part you want - We send the molecule to RDKit and tell it to assign stereo from 3D
rdmol = mol.to_rdkit()
Chem.AssignStereochemistryFrom3D(rdmol)
mol2 = Molecule.from_rdkit(rdmol)

# Now the graph stereo has been added back based on the 3d coords :-)
mol2.visualize(backend="rdkit")

[orionarcher]

Awesome! Thank you so much @j-wags!