Openff-toolkit: ToolkitWrapper.from_smiles accepts SMILESes with ambiguous stereochemistry

[j-wags]

The SMILES [H]C(F)(Cl)Br shouldn't be accepted when making an OFFMol becsause there is no stereochemistry defined for the C (I think?). However passing it in to different TookitWrappers's from_smiles function doesn't raise an error.
RDKitToolkitWrapper

toolkit_wrapper = RDKitToolkitWrapper()
smiles = '[H]C(F)(Cl)Br'
molecule = toolkit_wrapper.from_smiles(smiles)
print(molecule.to_dict())

The central carbon atom is not assigned a stereochemistry:

OrderedDict([('atomic_number', 6), ('formal_charge', 0), ('is_aromatic', False), ('stereochemistry', None), ('name', '')])

Giving a valid SMILES, [H][C@@](F)(Cl)Br yields

OrderedDict([('atomic_number', 6), ('formal_charge', 0), ('is_aromatic', False), ('stereochemistry', 'R'), ('name', '')])

OpenEyeToolkitWrapper

toolkit_wrapper = OpenEyeToolkitWrapper()
smiles = '[H]C(F)(Cl)Br'
molecule = Molecule.from_smiles(smiles, toolkit_registry=toolkit_wrapper)
print(molecule.to_dict())

The central carbon is not assigned a stereochemistry:

OrderedDict([('atomic_number', 6), ('formal_charge', 0), ('is_aromatic', False), ('stereochemistry', None), ('name', '')])

Giving a valid SMILES, [H][C@@](F)(Cl)Br yields

OrderedDict([('atomic_number', 6), ('formal_charge', 0), ('is_aromatic', False), ('stereochemistry', 'R'), ('name', '')])

[davidlmobley]

Hmm, that's interesting. Probably can be deferred for now if it's an isolated case, but definitely isn't right.

@bannanc did you see this?

[j-wags]

Yeah, this probably isn't a hard fix, just need to track down the appropriate OE and RDK checks and throw an exception. I mostly left this as a reminder to myself to fix while I validate our stereo handling, but I'd love any advice.

Leaving a resource here for when I come back to this: https://docs.eyesopen.com/toolkits/cpp/oechemtk/stereochemistry.html

[bannanc]

@j-wags and @davidlmobley I missed this one. My first question is I thought input molecules were required to have 3D coordinates, if a molecule has coordinates then the OEMol should have coordinates.

If we are no longer requiring 3D coordinates on the input molecule then I think there are options on reading a SMILES where you can require stereochemistry to be specified.

[bannanc]

I'd need to dig into smiles parsing in both packages, I don't know off the top of my head what the API points are for checking on parsing.

[jchodera]

My first question is I thought input molecules were required to have 3D coordinates, if a molecule has coordinates then the OEMol should have coordinates.

I don't think we want to require coordinates, since it should be possible to construct parameters for a system whose coordinates have not yet been specified. Generating charges should also use standard approaches to enumerating conformers so that the charges are not geometry-independent.

Having the ability to associate coordinates with a Topology object will be useful when generating ParmEd Structure objects to convert to popular MD package input formats, but should not be required.

[bannanc]

Ok, in that case, I think I disagree with:

Probably can be deferred for now if it's an isolated case, but definitely isn't right.

Because it won't be an isolated case, as plenty of databases provide SMILES without stereochemistry specified.
@j-wags I can help by looking into the smiles parsing into OEMols and RDMols to see if there is an easy way to raise an error for no stereochemistry other than looping over all atoms and checking if they are a stereocenter.

[davidlmobley]

Ah, I think I was misunderstanding the issue, @j-wags -- I thought you were saying this was losing stereochemistry on a round trip, but it sounds like you're saying it's just not raising an Exception when stereochemistry is undefined. And Caitlin is right -- it very much needs to raise an Exception when stereochemistry is undefined. If people have no idea what molecule they want to be simulating they have no business simulating it. :)

[davidlmobley]

(So I'd bump up the priority relative to what I said earlier.)

[j-wags]

Ah, that's right, @davidlmobley. This is going to be an essential part of our Molecule.from_smiles() functionality.

@bannanc, thanks for the offer. Yes -- it would be really helpful if you found those commands. It's not the end of the world if we have to inefficiently loop over all possible stereocenters, but it would certainly be better if the toolkits have a quick "does-this-smiles-have-ambiguous-stereochemistry" check.

[j-wags]

RDKit unassigned chiral center check:
Chem.FindMolChiralCenters(rdmol, includeUnassigned=True)

OE unassigned chiral center check
bool HasStereoSpecified(unsigned int v=OEAtomStereo::All) const
https://docs.eyesopen.com/toolkits/python/oechemtk/OEChemClasses/OEAtomBase.html

OE unassigned bond stereo check
bool HasStereoSpecified(unsigned int v=OEBondStereo::All) const
https://docs.eyesopen.com/toolkits/python/oechemtk/OEChemClasses/OEBondBase.html#OEChem::OEBondBase::HasStereoSpecified

[bannanc]

I think OE's default behavior is to find the unassigned one, I'm not sure what is already included in smiles parsing in either case, but I've been meaning to look into for my own reasons so I moved it up my list. I'll let you know what I find early next week.

[j-wags]

Fixed in 6a03e2604ffade3bcbb13bb74731db6f4bad0558