Openff-toolkit: RDMol stereochemistry detector too strict

[trevorgokey]

Trying to use:

openforcefield.topology.Molecule.from_rdkit( rdmol, allow_undefined_stereo=False)

where the molecule in question is c1c[cH:1][c:2](cc1)[CH2:3][c:4]2ccccc2:

fails due to the sp3 carbon in the center (index 3 in the smiles above). The function RDKitToolkitWrapper::_find_undefined_stereo_atoms in toolkits.py will fail anything that is CHI_UNSPECIFIED and possibly chiral:

   2304         # Find all atoms with undefined stereo.
   2305         undefined_atom_indices = []
   2306         for atom_idx, atom in enumerate(rdmol.GetAtoms()):
   2307             if (atom.GetChiralTag() == Chem.ChiralType.CHI_UNSPECIFIED and
   2308                     atom.HasProp('_ChiralityPossible')):
   2309                 undefined_atom_indices.append(atom_idx)
   2310         return undefined_atom_indices

However in the case here, the carbon is CHI_UNSPECIFIED, and atom.HasProp('_ChiralityPossible') -> True, causing the failure.

I suggest being a little more relaxed, since chirality "is possible", but it is in fact not chiral. I am not super experienced here, but possibly checking for _CIPCode in the loop might be useful. In this case, for all atoms atom.HasProp('_CIPCode') -> False even after the prior call to Chem.rdmolops.AssignStereochemistry and any other checking function I found, so I assume RDKit is trying to tell us it is not a chiral center after all checks are made. CHI_UNSPECIFIED hints at the fact it's not chiral (i.e. a sufficient condition), but perhaps _CIPCode might be a necessary condition.

[trevorgokey]

Looking just a bit further, the _ChiralityPossible flag is _only_ set immediately prior the above snippet:

2301         # Flag possible chiral centers with the "_ChiralityPossible".
2302         Chem.AssignStereochemistry(rdmol, force=True, flagPossibleStereoCenters=True)
2303

So perhaps a solution can be concocted which is more careful with flagPossibleStereoCenters. If it was False, I'd be inclined to say the above mol would pass through without any other changes.

[j-wags]

Thanks for the report @trevorgokey . This is really strange... That molecule is indeed not chiral. This may be a bug with RDKit. Could you send the file over so I can run some tests?

[trevorgokey]

Which file would you like? To reproduce, I used:

from openforcefield.topology import Molecule
from rdkit import Chem
mol = Chem.MolFromSmiles('[H:14][c:1]1[c:3]([c:7]([c:11]([c:8]([c:4]1[H:17])[H:21])[C:13]([H:24])([H:25])[c:12]2[c:9]([c:5]([c:2]([c:6]([c:10]2[H:23])[H:19])[H:15])[H:18])[H:22])[H:20])[H:16]', sanitize=False)
Chem.SanitizeMol( mol, Chem.SanitizeFlags.SANITIZE_ALL ^ Chem.SanitizeFlags.SANITIZE_ADJUSTHS )
offmol = Molecule.from_rdkit( mol, allow_undefined_stereo=False)

Oddly enough, if I use sanitize=True it works, but then it strips the hydrogen (which I am trying very hard to keep)

[j-wags]

Ahh, I think I see what's happening here. It's treating the labeled atoms as being "different" for the purposes of stereochemistry. That is, an sp2 carbon with index "7" is different from an sp2 carbon with index "2". If I strip the labels off, it loads fine.

from openforcefield.topology import Molecule
from rdkit import Chem
smiles = '[H:14][c:1]1[c:3]([c:7]([c:11]([c:8]([c:4]1[H:17])[H:21])[C:13]([H:24])([H:25])[c:12]2[c:9]([c:5]([c:2]([c:6]([c:10]2[H:23])[H:19])[H:15])[H:18])[H:22])[H:20])[H:16]'
smiles_no_index = '[H][c]1[c]([c]([c]([c]([c]1[H])[H])[C]([H])([H])[c]2[c]([c]([c]([c]([c]2[H])[H])[H])[H])[H])[H])[H]'
mol = Chem.MolFromSmiles(smiles_no_index, sanitize=False)
Chem.SanitizeMol( mol, Chem.SanitizeFlags.SANITIZE_ALL ^ Chem.SanitizeFlags.SANITIZE_ADJUSTHS )
offmol = Molecule.from_rdkit( mol, allow_undefined_stereo=False)

I would think that the "chemically correct" behavior is to not see a difference between atoms which differ only by index, so I'm optimistic that there is a resolution for this.

Secondarily, we should keep our eyes open for a different issue -- your use case probably relies on knowing the indices of the atoms, ~but I'm not at all certain that those indices get transferred into the toolkit in any meaningful way~ (apparently they do, see next comment).

I need to work on some things that are due tomorrow, but I'll try to revisit both of these questions ASAP (maybe early next week)

[j-wags]

Actually, the second part (whether the toolkit does anything with atom indices) may already be handled correctly, so that's a bit of good news. I would test this once or twice to verify first, however.

[trevorgokey]

Hmm, weird. Well the indices are very necessary; RDKit stores them and exposes them via atom.GetAtomMapNum. They are used for example to get the right coordinates from a QCMol if I enumerate the atoms via mol.GetAtoms.

Previous solutions included stripping the hydrogen, and a custom pairwise distance calculation to determine overlap. However, keeping the hydrogen as stated in the QCMol will be important when variable pH is used.

[trevorgokey]

So, in the meantime, I am forced to use allow_undefined_stereo=True, which is uncomfortable, but manageable, since the molecule was already successfully passed through QCArchive, and "should be" problem free.

[trevorgokey]

To further the quality of the cleanup of this issue, can we do something with the warning message? Right now I get, for every occurrence:

2409         if msg is not None:                                                      
   1             if raise_warning:                                                    
   2                 msg = 'Warning (not error because allow_undefined_stereo=True): '
   3                 logger.warning(msg)                                              
   4             else:                                                                
   5                 msg = 'Unable to make OFFMol from RDMol: ' + msg                 
   6                 raise UndefinedStereochemistryError(msg)

With the colon at the end, it makes me think more stuff will print (the intention was probably msg = 'Warning (not error because allow_undefined_stereo=True): ' + msg). Regardless, can we suppress this warning, given the fact I explicitly set allow_undefined_stereo=True?

[j-wags]

We talked about this today. Basically, what we EVENTUALLY want is a native way to turn QCA molecules into OFF molecules. There are several intermediate goals that we could snag along the way. But in the interest of time, I'm defining the minimum thing we could do to unblock @trevorgokey.

So, this issue will be closed when the following code block does not return an undefined stereochemistry error.

from openforcefield.topology import Molecule
from rdkit import Chem
smiles = '[H:14][c:1]1[c:3]([c:7]([c:11]([c:8]([c:4]1[H:17])[H:21])[C:13]([H:24])([H:25])[c:12]2[c:9]([c:5]([c:2]([c:6]([c:10]2[H:23])[H:19])[H:15])[H:18])[H:22])[H:20])[H:16]'
offmol=Molecule.from_smiles(smiles)

[jchodera]

Just so I understand what is being requested: Do you mean you want to integrate the cmiles functionality into Molecule.from_smiles() so that it will automatically handled index-tagged SMILES properly?

@chayast: Which method(s) from cmiles could we use to ingest an index-tagged SMILES and return a property-ordered RDMol without getting an undefined-stereochemistry error? Or if that doesn't exist, what parts of cmiles could we use to help with this?

[j-wags]

Just so I understand what is being requested: Do you mean you want to integrate the cmiles functionality into Molecule.from_smiles() so that it will automatically handled index-tagged SMILES properly?

Eventually, yes. But not now. There are a few questions about what behavior we would want out of that -- Like, what if there are gaps in the indices? What if only some atoms have indices? What if the indices don't start at 1 (or zero)? What if two atoms have the same index?

One solution that we could implement quickly is to make a copy of the rdmol, strip all of the copy's map indices, run AssignStereochemistry on the copy, and then copy back the _CHI_UNSPECIFIED values from the copy to the original.

[jchodera]

what if there are gaps in the indices? What if only some atoms have indices? What if the indices don't start at 1 (or zero)? What if two atoms have the same index?

It's pretty clear that the answer to all of those questions should be "throw an exception". (Same is true of an implicit hydrogen SMILES with tagged atom indices.)

Being able to use a CMILES-generated index-tagged explicit-hydrogen canonical isomeric SMILES string is _very_ useful for creating molecules and manipulating topologies, but the utility is basically almost zero once you leave the CMILES-generated land.

[j-wags]

Actually, yeah. I agree with all of that, but we should call the function with that behavior Molecule.from_cmiles.

I don't want to encode that as the behavior of Molecule.from_smiles since there are all sorts of reasons/methods people use to assign map indices outside of our ecosystem, and we shouldn't play favorites. Actually, even IN our ecosystem, the tagged SMILES we assign to a TorsionDrive only has four atoms with indices -- The atoms defining the driven torsion.

[jchodera]

Actually, yeah. I agree with all of that, but we should call the function with that behavior Molecule.from_cmiles.

That could create more confusion, not less, since cmiles outputs a bunch of different forms, including InChI, and we are just talking about one of them. Would from_cmiles accept each output from cmiles? What if cmiles is extended to output more types?

How about Molecule.from_indexed_smiles()? Molecule.from_smiles() could just strip off any indices and ignore them.

[trevorgokey]

Can we do Molecule.from_smiles( string, ignore_index=False)? This keeps current behavior for non-indexed smiles whether ignore_index is True or False, and adds new functionality when there are indices. This way, the default signature will give me expected results regardless if the string is indexed or not. Setting ignore_index=True can just be an override for the user in case they want to explicitly ignore their indices (e.g. similar to allow_undefined_stereo).

[ChayaSt]

cmiles does not respect tags in a SMILES when it ingests it because it will generate tagged SMILES with the map indices in a canonical order (function to canonicalize order https://github.com/openforcefield/cmiles/blob/master/cmiles/_cmiles_oe.py#L68-L109). However, if you give it a qcschema molecule, it will create an oemol or rdkit mol with atoms in the order of qcschema xyz coordinates. See https://github.com/openforcefield/cmiles/blob/master/cmiles/_cmiles_oe.py#L11-L65

It can be pretty challenging to create an oemol or rdkit mol with the atoms ordered according to the map indices if you don't have the qcschema molecule. Because then you rely on rdkit's internal SMILES parser and it will order the atoms according to its algorithm. However, you can get the mapping from the map indices to the atom indices pretty easily by iterating over the atoms and mapping the atom index to the map index.

@trevorgokey, when would you have a situation where you would want to create a molecule just from a mapped SMILES without the associated coordinates?

[ChayaSt]

Another point relating to rdkit's steroechmistry detection. If the atom has map indices, it will include that information in the weight of the atom and then complain that it should be chiral. This is cmiles handles this feature/bug
https://github.com/openforcefield/cmiles/blob/master/cmiles/_cmiles_rd.py#L321-L327

[j-wags]

I just happened upon the explanation for this -- RDKit considers atom map indices to be equivalent to isotopic labels.
https://github.com/rdkit/rdkit/issues/553