6
$\begingroup$

I am having trouble getting the SMILES strings and drawings of the different resonance structutes of a molecule with RDKit. I found an old example online doing this but it doesn't seem to work with the current version of RDKit (2023-03-01).

Here is an example with 1,2-Dichlorobenzene:

from rdkit import Chen
from rdkit.Chem.Draw import rdMolDraw2D

# smiles = ['c1c(Cl)c(Cl)ccc1', 'C1C(Cl)=C(Cl)C=CC=1', 'C1=C(Cl)C(Cl)=CC=C1']
smiles = ['c1c(Cl)c(Cl)ccc1']
smi = smiles[0]
mol = Chem.MolFromSmiles(smi)
suppl = Chem.ResonanceMolSupplier(mol, Chem.KEKULE_ALL)
print(f'# of kekule forms {len(suppl)}\n')

for i, m in enumerate(suppl):
    print(f'Resonance structure #{i}\n')
    smi_kkl = Chem.MolToSmiles(m) # change to turn kekulization on
    print(f'smi_kkl: {smi_kkl}\n') 
    ## newer way of plotting does not work either
    # d = rdMolDraw2D.MolDraw2DCairo(350, 300)
    # d.drawOptions().maxFontSize = 20
    # rdMolDraw2D.PrepareAndDrawMolecule(d, mol, kekulize=False) # change to turn kekulization on
    # d.FinishDrawing()
    # d.WriteDrawingText(os.path.join(parent_path, 'figures', f'{smi}_{i}.png'))

    img = Draw.MolToImage(mol, kekulize=False) # change to turn kekulization on
    img.save(os.path.join(parent_path, 'figures', f'{smi}_{i}.png'))

The output is:

# of kekule forms 2

Resonance structure #0

smi_kkl: Clc1ccccc1Cl # should be C1C(Cl)=C(Cl)C=CC=1

Resonance structure #1

smi_kkl: Clc1c(Cl)cccc1 # should be C1=C(Cl)C(Cl)=CC=C1

And the images both use aromatic rings instead of the kekule forms. If I turn on kekulization where specified in the code, I only get one kekulized structure for both resonance objects (in this case ClC1=CC=CC=C1Cl).

Any ideas how to get the two structures in SMILES and images?

$\endgroup$

3 Answers 3

5
$\begingroup$

I did some digging around and found a hack. The Chem.MolToSmiles function has an argument rootedAtAtom which if non-negative, forces the SMILES to start at a particular atom. It Defaults to -1.

So I tried to set it to the various atom ids present in the molecule like shown below:

from rdkit import Chem
smiles = ['c1c(Cl)c(Cl)ccc1']
smi = smiles[0]
mol = Chem.MolFromSmiles(smi)
suppl = Chem.ResonanceMolSupplier(mol, Chem.KEKULE_ALL)


for i in range(8):
  print (Chem.MolToSmiles(suppl[0], kekuleSmiles=True, rootedAtAtom=i))

This prints out:

C1=C(Cl)C(Cl)=CC=C1
C1(Cl)=CC=CC=C1Cl
ClC1=CC=CC=C1Cl
C1(Cl)=CC=CC=C1Cl
ClC1=CC=CC=C1Cl
C1=C(Cl)C(Cl)=CC=C1
C1=CC=C(Cl)C(Cl)=C1
C1=CC=C(Cl)C(Cl)=C1

Out of which the first one: C1=C(Cl)C(Cl)=CC=C1 is exactly one of the SMILES of the resonance structures which you were expecting.

And regarding drawing the structures, I noticed something hacky, which I will show below:

s1 = Chem.MolToSmiles(suppl[1], kekuleSmiles=True, rootedAtAtom=0)
print(s1)
Chem.MolFromSmiles(s1)

>> C1=C(Cl)C(Cl)=CC=C1

enter image description here

s2 = Chem.MolToSmiles(suppl[1], kekuleSmiles=False, rootedAtAtom=0)
print(s2)
Chem.MolFromSmiles(s2)

>> c1cccc(Cl)c1Cl

enter image description here

$\endgroup$
3
  • $\begingroup$ This is helpful, but not what I am looking for. 'rootedAtAtom' prints a SMILES with the given atom first. This leads to the different kekule forms as byproduct of the atom iterations. I want to automatically print only the valid kekule forms of the given molecules (in image and SMILES). $\endgroup$
    – Njw96
    May 12, 2023 at 14:50
  • $\begingroup$ I would suggest you open an issue on rdkit GitHub repo @Njw96 :) $\endgroup$ May 12, 2023 at 14:54
  • $\begingroup$ yes I did, waiting for one more reply then I can post the answer here. You can follow the thread here if you want: github.com/rdkit/rdkit/discussions/6358 $\endgroup$
    – Njw96
    May 12, 2023 at 15:24
5
$\begingroup$

You have to kekulize your MolFromSmiles and mark all bonds non-aromatic.

from rdkit import Chem
from rdkit.Chem import Draw, rdDepictor
rdDepictor.SetPreferCoordGen(True)
from rdkit.Chem.Draw import IPythonConsole

smiles = 'c1c(Cl)c(Cl)ccc1'

mol = Chem.MolFromSmiles(smiles)
Chem.Kekulize(mol, clearAromaticFlags=True)

Now you get the resonance structures.

suppl = Chem.ResonanceMolSupplier(mol, Chem.KEKULE_ALL)

for resMol in suppl:
    res_smi = Chem.MolToSmiles(resMol, kekuleSmiles=True)
    print(res_smi)

ClC1=CC=CC=C1Cl
ClC1=C(Cl)C=CC=C1

suppl[0]

enter image description here

suppl[1]

enter image description here

$\endgroup$
3
  • $\begingroup$ Thanks @rapelpy, this works for the smiles but not the images... If you just add this plot to the loop you get only one kekule form: img = Draw.MolToImage(mol, kekulize=True) # change to turn kekulization on img.save(os.path.join(parent_path, 'figures', f'{smi}_{i}.png')) Any tips for this? $\endgroup$
    – Njw96
    May 14, 2023 at 13:44
  • $\begingroup$ @Njw96 You are drawing the images from the starting mol and not from the resonance mols. $\endgroup$
    – rapelpy
    May 14, 2023 at 15:24
  • $\begingroup$ no actually in my code mol is your resMol because I have the plotting logic in a separate function. What was missing is to unset the aromatic flag, I will post the answer now! $\endgroup$
    – Njw96
    May 15, 2023 at 7:29
5
$\begingroup$

One of the contributors of RDKit gave a solution here. I am paraphrasing it below.

The problem occurs because RDKit by default assigns aromatic labels to atom and bonds as part of its sanitization process, and therefore loses the location of the double bonds in the original input. If you request a kekule form (in SMILES or image), RDKit uses its own kekulizer to generate a valid representation.

Below I present work-arounds for two scenarios. First a plotting helper function:

def generate_images(mols, should_kekulize_aromatic_smiles=False):
    imgs = []
    
    for i, m in enumerate(mols):
        smi = Chem.MolToSmiles(m)
        print(f'smi {smi}\n')
        ## only kekulize molecule if it has aromatic atoms
        ## used to choose whether to plot kekulized or aromatic smiles for 'c1c(Cl)c(Cl)ccc1'
        has_aromatic_atoms = any(a.GetIsAromatic() for a in m.GetAtoms())
        if not has_aromatic_atoms: should_kekulize_aromatic_smiles = False
        drawer = rdMolDraw2D.MolDraw2DCairo(200, 100)
        drawer.drawOptions().prepareMolsBeforeDrawing = False
        m = rdMolDraw2D.PrepareMolForDrawing(m, kekulize=should_kekulize_aromatic_smiles)
        drawer.DrawMolecule(m)
        drawer.FinishDrawing()
        img = drawer.GetDrawingText()
        imgs.append(img)
        with open(os.path.join(parent_path, 'figures', f'{smi}_{i}.png'), "wb") as hnd:
            hnd.write(img)
    return imgs

Plotting with given SMILES

Let's look at the example of plotting from defined kekule and aromatic SMILES, e.g. ['c1c(Cl)c(Cl)ccc1', 'C1C(Cl)=C(Cl)C=CC=1', 'C1=C(Cl)C(Cl)=CC=C1']. In this case, we need to sanitize the molecules without kekulization/aromaticity.

def sanitize_without_aromatic(smiles):
    mols = []
    smiles_parser_params = Chem.SmilesParserParams()
    smiles_parser_params.sanitize = False   

    for i, smi in enumerate(smiles):
        m = Chem.MolFromSmiles(smi, smiles_parser_params)
        sanitize_flags = Chem.SanitizeFlags.SANITIZE_ALL
        ## turn off kekulization and aromaticity for molecules with no aromatic atoms
        has_aromatic_atoms = any(a.GetIsAromatic() for a in m.GetAtoms())
        if not has_aromatic_atoms:
            sanitize_flags ^= (Chem.SanitizeFlags.SANITIZE_KEKULIZE | Chem.SanitizeFlags.SANITIZE_SETAROMATICITY)
        Chem.SanitizeMol(m, sanitize_flags)

        mols.append(m)

    return mols

putting the two functions together:

smiles = ['c1c(Cl)c(Cl)ccc1', 'C1C(Cl)=C(Cl)C=CC=1', 'C1=C(Cl)C(Cl)=CC=C1']
mols = sanitize_without_aromatic(smiles)
generate_images(mols=mols, should_kekulize_aromatic_smiles=False)

gives the following output:

smi Clc1ccccc1Cl

imgs[0]

aromatic_smi

smi ClC1=C(Cl)C=CC=C1

enter image description here

smi ClC1=CC=CC=C1Cl

enter image description here

Plotting the output of ResonanceMolSupplier

ResonanceMolSupplier returns sanitized molecules, and we need to unset the aromatic flags manually. Here is a function to do that.

def unset_aromatic_flags(m):
    for b in m.GetBonds():
        if b.GetBondType() != Chem.BondType.AROMATIC and b.GetIsAromatic():
            b.SetIsAromatic(False)
            b.GetBeginAtom().SetIsAromatic(False)
            b.GetEndAtom().SetIsAromatic(False)
    return m

You can test it with the following example:

smi = 'c1c(Cl)c(Cl)ccc1'
mol = Chem.MolFromSmiles(smi)
suppl = Chem.ResonanceMolSupplier(mol, Chem.KEKULE_ALL)
mols = [mol]+[unset_aromatic_flags(m) for m in suppl]
imgs = generate_images(mols=mols, should_kekulize_aromatic_smiles=False)

The output is the same as above.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .