ChayaSt
commented
5 years ago Thank you @j-wags for this write-up. I found that if Hydrogens are added to a primary imine bond, the code with RDKit that I was using to flag bonds that are missing stereo information finds an undefined stereobond.
def flag_unspec_stereo(mol):
unspec = False
Chem.FindPotentialStereoBonds(mol)
for bond in mol.GetBonds():
if bond.GetStereo() == Chem.BondStereo.STEREOANY:
print(bond.GetBeginAtom().GetSymbol(), bond.GetSmarts(), bond.GetEndAtom().GetSymbol())
unspec = True
return unspec
mol = Chem.MolFromSmiles('CN(C)C(=N)NC(=N)N')
mol_h = Chem.AddHs(mol)
print(flag_unspec_stereo(mol))
print(flag_unspec_stereo(mol_h))output:
False
C = N
C = N
TrueThe code you provided above does not find unspecified stereochemistry even if I add Hs to the molecule.
def rdkit_unspec_stereo(smiles):
from rdkit import Chem
from rdkit.Chem import EnumerateStereoisomers
rdmol = Chem.MolFromSmiles(smiles)
# Adding H's can hide undefined bond stereochemistry, so we have to test for undefined stereo here
unspec_stereo = False
rdmol_copy = Chem.Mol(rdmol)
rdmol_copy = Chem.AddHs(rdmol)
Chem.FindPotentialStereoBonds(rdmol_copy)
enumsi_opt = EnumerateStereoisomers.StereoEnumerationOptions(maxIsomers=2, onlyUnassigned=True)
stereoisomers = [isomer for isomer in Chem.EnumerateStereoisomers.EnumerateStereoisomers(rdmol_copy, enumsi_opt)]
if len(stereoisomers) != 1:
unspec_stereo = True
if unspec_stereo:
return True
return False
print(rdkit_unspec_stereo('CN(C)C(=N)NC(=N)N'))FalseIt's interesting because if you check the flags in each bond, these bonds are flagged as undefined stereo but when enumerating stereoisomers they are ignored.
jchodera
j-wags
proteneer
coleb
greglandrum
Purpose
Because it's a complicated topic, and to maintain the ability to restart this discussion in the future, I think it's important to record our design decisions with regard to stereochemistry in this repository.
Background
Stereochemical differences are important in the SMIRNOFF spec, as they can lead to different partial charges or partial bond orders on a molecule, and (if the project goes in that direction), different assignment of library parameters using stereochemistry-dependent SMARTS.
However, the definition of meaningful stereochemistry can change depending on experimental conditions and timescale. For example a tertiary amine will invert 10^3 to 10^5 times per second in solution. Therefore, on a wet-lab-experiment timescale, the amine does not have meaningful stereochemistry. However, if one is running a QM calculation on the tertiary amine, there may be meaningful differences in the results depending on its stereochemistry. Thus, whether or not stereochemistry is required in a tertiary amine depends on your use case.
Our definition of stereochemistry will ideally cover the fact that MD simulations can run for micro- and milli-seconds, therefore molecular features that are likely to change stereochemistry on those timescales should not be assigned stereochemistry-specific properties. It is unclear currently whether this decision should be made at the point of SMIRNOFF molecule definition (ie, we enforce that short-lived stereochemical features do not distinguish one OFF molecule from another, thereby preventing our engine from assigning them different simulation parameters), or if the SMARTS that define our parameter libraries should be policed moving forward, to ensure that no two parameters are differentiated only by low-barrier stereochemistry (eg, we will never create a SMARTS that relies on the stereochemistry of a tertiary amine).
Current implementation plans
Read below for details, but I think a good plan is to let the toolkits use their own (slightly different!) standards in choosing how to accept SMILES. While they mostly agree with what "fully specified stereochemistry" is, they differ on some borderline cases.
While this isn't ideal, ensuring the same stereochemistry standard across toolkits will require a lot of engineering, and will restrict us to the union of the most restrictive stereochemistry standards of all of our toolkits (eg. OE omega complains if bond stereo isn't defined for
CC=NH, which is a bit much).3D molecule input
We haven't run into too many problems with sdf and mol2 input. OE and RDK currently use their default stereochemistry detection in these cases and seem to rarely have problems. I haven't dug very deeply into these, so I will update this post as I look further into it.
Potentially relevant to 3D stereochemistry detection: #141, #137
SMILES input
"Common sense" atom and double bond stereochemistry is required, as is implemented by default in OE and RDK. There are, however, a few cases where disagreement occurs, both between OE and RDK, and between a toolkit and our intuition.
This table records the cases of obvious and border-case stereochemistry we've found so far, whether we think molecules should labeled as having undefined stereochemistry (on the MD timescale), and whether the toolkits define the SMILES as having undefined stereochemistry. If applicable, experimental barriers/interconversion rates from literature should be added.
C(Cl)(F)(Br)CC(F)=C(Cl)C\C(F)=C(/F)CC(C)(Cl)BrCC(F)=C(F)C[C@@](C)(Cl)BrC\C(F)=C(/F)C[C@@](C)(Cl)BrCN=C(F)CN=C(F)CN(Cl)(F)CCP(O)(=O)NCNC(C)=OCN=C(C)OC(=[O+][H])CN(C)(CC)(CCC)=OC[C@H]1CN1C, may have tested using different methodC[C@H]1CN1, may have tested using different method[H]c1c(c(c(c(c1[H])[H])P2C(=C3C(=C2c4c(c(c(c(n4)[H])[H])[H])[H])C(C(C(C3([H])[H])([H])[H])([H])[H])([H])[H])c5c(c(c(c(n5)[H])[H])[H])[H])[H])[H]CCS(=O)CCOP(=O)([O-])OCC[1] OE says "Yes" when using method A, but when "No" when method B is given
enum_nitrogens=False. [2] I'm not certain what the rotation timescale for amide bonds is, but this might just be in molecular dynamics range. Letting the enol form of the amide get stereochemistry might give users the choice of whether they want it defined.Code to reproduce the above