On Tue, Apr 15, 2008 at 12:41 PM, Noel O'Boyle <[email protected]> wrote: > I am running into trouble building a molecule from scratch. Let's say > the molecule represented by "F/C=C/F" (taken from the Daylight SMILES > tutorial). Here's the code to build it (according to me): > > import Chem > > # Set some lookups > _bondtypes = {1: Chem.BondType.SINGLE, > 2: Chem.BondType.DOUBLE, > 3: Chem.BondType.TRIPLE} > _bondstereo = {0: Chem.rdchem.BondStereo.STEREONONE, > 1: Chem.rdchem.BondStereo.STEREOE, > 2: Chem.rdchem.BondStereo.STEREOZ} > > # Build F/C=C/F > rdmol = Chem.Mol() > rdedmol = Chem.EditableMol(rdmol)
ah... my friend EditableMol. As you are discovering, this is a very good way to shoot yourself in the foot. (Particularly if there's no documentation telling you how to aim the gun) :-) > for atomnum in [9, 6, 6, 9]: > rdatom = Chem.Atom(atomnum) > rdedmol.AddAtom(rdatom) > > for bond in [(0,1,1), (1,2,2), (2,3,1)]: > rdedmol.AddBond(bond[0], > bond[1], > _bondtypes[bond[2]]) > > rdmol = rdedmol.GetMol() > for stereoID, newbond in zip([0, 1, 0], rdmol.GetBonds()): > newbond.SetStereo(_bondstereo[stereoID]) The problem is here. I shouldn't even have exposed SetStereo to Python. To indicate bond stereochemistry, you need to set the directionalities of the neighboring single bonds. So you'd do: [13] >>> rdmol.GetBondWithIdx(0).SetBondDir(Chem.BondDir.ENDUPRIGHT) [15] >>> rdmol.GetBondWithIdx(2).SetBondDir(Chem.BondDir.ENDUPRIGHT) Notice that this does not set the stereochemistry of the double bond: [16] >>> rdmol.GetBondWithIdx(1).GetStereo() Out[16]: Chem.rdchem.BondStereo.STEREONONE The RDKit is, at the moment, lazy about stereochemistry : it doesn't do the perception until it's required. So if you ask for SMILES: [17] >>> Chem.MolToSmiles(rdmol,True) Out[17]: 'F/C=C/F' you get stereochemistry calculated: [18] >>> rdmol.GetBondWithIdx(1).GetStereo() Out[18]: Chem.rdchem.BondStereo.STEREOE If you want the stereochem information earlier, you can call: Chem.AssignBondStereoCodes(rdmol) > Chem.SanitizeMol(rdmol) > # Gets here fine > > # Print the result > print Chem.MolToSmiles(rdmol, isomericSmiles=True) > # RuntimeError: Pre-condition Violation There actually is an error message associated with this, but it doesn't make it through the C++ / Python wrapper. That's something for me to look into. > What am I missing? It works fine if I use [0, 0, 0] for stereoID, so > it's something to do with the way I handle the stereoisomer. If you're interested in the actual reason behind the error, I'm happy to explain, but I'd understand if you just want to know how to make it go away. :-) -greg
