# RDKit Cookbook

## Introduction

### What is this?

This document provides example recipes of how to carry out particular tasks using the RDKit functionality from Python. The contents have been contributed by the RDKit community, tested with the latest RDKit release, and then compiled into this document. The RDKit Cookbook is written in reStructuredText, which supports Sphinx doctests, allowing for easier validation and maintenance of the RDKit Cookbook code examples, where appropriate.

### What gets included?

The examples included come from various online sources such as blogs, shared gists, and the RDKit mailing lists. Generally, only minimal editing is added to the example code/notes for formatting consistency and to incorporate the doctests. We have made a conscious effort to appropriately credit the original source and authors. One of the first priorities of this document is to compile useful short examples shared on the RDKit mailing lists, as these can be difficult to discover. It will take some time, but we hope to expand this document into 100s of examples. As the document grows, it may make sense to prioritize examples included in the RDKit Cookbook based on community demand.

### Feedback and Contributing

If you have suggestions for how to improve the Cookbook and/or examples you would like included, please contribute directly in the source document (the .rst file). Alternatively, you can also send Cookbook revisions and addition requests to the mailing list: <rdkit-discuss@lists.sourceforge.net> (you will need to subscribe first).

Note

The Index ID# (e.g., RDKitCB_##) is simply a way to track Cookbook entries and image file names. New Cookbook additions are sequentially index numbered, regardless of where they are placed within the document. As such, for reference, the next Cookbook entry is RDKitCB_40.

## Drawing Molecules (Jupyter)

### Include an Atom Index

Author: Takayuki Serizawa
Original Source: https://iwatobipen.wordpress.com/2017/02/25/draw-molecule-with-atom-index-in-rdkit/
Index ID#: RDKitCB_0
Summary: Draw a molecule with atom index numbers.

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    from rdkit.Chem import Draw
    IPythonConsole.ipython_useSVG=True  #< set this to False if you want PNGs instead of SVGs

    def mol_with_atom_index(mol):
        for atom in mol.GetAtoms():
            atom.SetAtomMapNum(atom.GetIdx())
        return mol

    # Test in a kinase inhibitor
    mol = Chem.MolFromSmiles("C1CC2=C3C(=CC=C2)C(=CN3C1)[C@H]4[C@@H](C(=O)NC4=O)C5=CNC6=CC=CC=C65")
    # Default
    mol

    # With atom index
    mol_with_atom_index(mol)

A simpler way to add atom indices is to adjust the IPythonConsole properties. This produces a similar image to the example above, the difference being that the atom indices are now near the atom, rather than at the atom position.

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    from rdkit.Chem import Draw
    IPythonConsole.drawOptions.addAtomIndices = True
    IPythonConsole.molSize = 300,300

    mol = Chem.MolFromSmiles("C1CC2=C3C(=CC=C2)C(=CN3C1)[C@H]4[C@@H](C(=O)NC4=O)C5=CNC6=CC=CC=C65")
    mol

### Include a Calculation

Author: Greg Landrum
Original Source: https://sourceforge.net/p/rdkit/mailman/message/36457619/
Index ID#: RDKitCB_23
Summary: Draw a molecule with a calculation value displayed (e.g., Gasteiger Charge)

    from rdkit import Chem
    from rdkit.Chem import AllChem
    from rdkit.Chem.Draw import IPythonConsole
    IPythonConsole.molSize = 250,250

    m = Chem.MolFromSmiles('c1ncncc1C(=O)[O-]')
    AllChem.ComputeGasteigerCharges(m)
    m

    m2 = Chem.Mol(m)
    for at in m2.GetAtoms():
        lbl = '%.2f'%(at.GetDoubleProp("_GasteigerCharge"))
        at.SetProp('atomNote',lbl)
    m2

### Include Stereo Annotations

Author: Greg Landrum
Source: https://github.com/rdkit/UGM_2020/blob/master/Notebooks/Landrum_WhatsNew.ipynb
Index ID#: RDKitCB_32
Summary: Draw a molecule with stereochemistry annotations displayed.

    from rdkit import Chem
    from rdkit.Chem import Draw
    from rdkit.Chem.Draw import IPythonConsole
    IPythonConsole.drawOptions.addAtomIndices = False
    IPythonConsole.drawOptions.addStereoAnnotation = True

    # Default Representation uses legacy FindMolChiralCenters() code
    m1 = Chem.MolFromSmiles('C1CC1[C@H](F)C1CCC1')
    m2 = Chem.MolFromSmiles('F[C@H]1CC[C@H](O)CC1')
    Draw.MolsToGridImage((m1,m2), subImgSize=(250,250))

    # new stereochemistry code with more accurate CIP labels, 2020.09 release
    from rdkit.Chem import rdCIPLabeler
    rdCIPLabeler.AssignCIPLabels(m1)
    rdCIPLabeler.AssignCIPLabels(m2)
    Draw.MolsToGridImage((m1,m2), subImgSize=(250,250))

### Black and White Molecules

Author: Greg Landrum and Vincent Scalfani
Source: https://gist.github.com/greglandrum/d85d5693e57c306e30057ec4d4d11342 and https://github.com/rdkit/rdkit/discussions/5885
Index ID#: RDKitCB_1
Summary: Draw a molecule in black and white.

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    from rdkit.Chem import Draw

    ms = [Chem.MolFromSmiles(x) for x in ('Cc1onc(-c2ccccc2)c1C(=O)N[C@@H]1C(=O)N2[C@@H](C(=O)O)C(C)(C)S[C@H]12','CC1(C)SC2C(NC(=O)Cc3ccccc3)C(=O)N2C1C(=O)O.[Na]')]
    Draw.MolsToGridImage(ms)

    IPythonConsole.drawOptions.useBWAtomPalette()
    Draw.MolsToGridImage(ms)

    # Alternatively, use the rdMolDraw2D package
    from rdkit.Chem.Draw import rdMolDraw2D
    import io
    from PIL import Image

    drawer = rdMolDraw2D.MolDraw2DCairo(500,180,200,180)
    drawer.drawOptions().useBWAtomPalette()
    drawer.DrawMolecules(ms)
    drawer.FinishDrawing()
    bio = io.BytesIO(drawer.GetDrawingText())
    Image.open(bio)

    # works for reactions too:
    # rxn is from https://github.com/rdkit/UGM_2020/blob/master/Notebooks/Landrum_WhatsNew.ipynb
    from rdkit.Chem import rdChemReactions
    rxn = rdChemReactions.ReactionFromSmarts("[cH:1]:1:[cH:2]:[cH:3]:[cH:4]:[cH:5](-[C:6]#[N:7]):[c:8]:1-[Cl].\
    [cH:10]:1:[cH:11]:[cH:12](-[Cl:16]):[cH:13]:[cH:14]:[cH:15]:1-B(-O)-O>>\
    [cH:1]:1:[cH:2]:[cH:3]:[cH:4]:[cH:5](-[C:6]#[N:7]):[c:8]:1-[cH:15]:1[cH:10]:[cH:11]:[cH:12](-[Cl:16]):[cH:13]:[cH:14]:1")
    drawer = rdMolDraw2D.MolDraw2DCairo(700,300)
    drawer.drawOptions().useBWAtomPalette()
    drawer.DrawReaction(rxn)
    drawer.FinishDrawing()
    bio = io.BytesIO(drawer.GetDrawingText())
    Image.open(bio)

### Highlight a Substructure in a Molecule

Author: Greg Landrum
Source: https://gist.github.com/greglandrum/5d45b56afe75603b955103cdd0d8e038
Index ID#: RDKitCB_2
Summary: Draw a molecule with a substructure highlight in Jupyter.

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole

    m = Chem.MolFromSmiles('c1cc(C(=O)O)c(OC(=O)C)cc1')
    substructure = Chem.MolFromSmarts('C(=O)O')
    print(m.GetSubstructMatches(substructure))

    ((3, 4, 5), (8, 9, 7))

    m

    # you can also manually set the atoms that should be highlighted:
    m.__sssAtoms = [0,1,2,6,11,12]
    m

### Highlight Molecule Differences

Author: Takayuki Serizawa
Original Source: https://gist.github.com/iwatobipen/6d8708d8c77c615cfffbb89409be730d
Index ID#: RDKitCB_36
Summary: Highlight molecule differences based on maximum common substructure

    from rdkit import Chem
    from rdkit.Chem import Draw
    from rdkit.Chem.Draw import IPythonConsole
    from rdkit.Chem import rdFMCS
    from rdkit.Chem.Draw import rdDepictor
    rdDepictor.SetPreferCoordGen(True)
    IPythonConsole.drawOptions.minFontSize=20

    mol1 = Chem.MolFromSmiles('FC1=CC=C2C(=C1)C=NN2')
    mol2 = Chem.MolFromSmiles('CCC1=C2NN=CC2=CC(Cl)=C1')

    Draw.MolsToGridImage([mol1, mol2])

    def view_difference(mol1, mol2):
        mcs = rdFMCS.FindMCS([mol1,mol2])
        mcs_mol = Chem.MolFromSmarts(mcs.smartsString)
        match1 = mol1.GetSubstructMatch(mcs_mol)
        target_atm1 = []
        for atom in mol1.GetAtoms():
            if atom.GetIdx() not in match1:
                target_atm1.append(atom.GetIdx())
        match2 = mol2.GetSubstructMatch(mcs_mol)
        target_atm2 = []
        for atom in mol2.GetAtoms():
            if atom.GetIdx() not in match2:
                target_atm2.append(atom.GetIdx())
        return Draw.MolsToGridImage([mol1, mol2],highlightAtomLists=[target_atm1, target_atm2])

    view_difference(mol1,mol2)

### Highlight Entire Molecule

Author: Vincent Scalfani
Original Source: https://github.com/vfscalfani/CSN_tutorial
Index ID#: RDKitCB_38
Summary: Highlight all atoms and bonds

    from rdkit import Chem
    from rdkit.Chem.Draw import rdMolDraw2D
    import io
    from PIL import Image

    mol = Chem.MolFromSmiles('CC(C)CN1C(=O)COC2=C1C=CC(=C2)NC(=O)/C=C/C3=CC=CC=C3')
    rgba_color = (0.0, 0.0, 1.0, 0.1) # transparent blue

    atoms = []
    for a in mol.GetAtoms():
        atoms.append(a.GetIdx())

    bonds = []
    for bond in mol.GetBonds():
        aid1 = atoms[bond.GetBeginAtomIdx()]
        aid2 = atoms[bond.GetEndAtomIdx()]
        bonds.append(mol.GetBondBetweenAtoms(aid1,aid2).GetIdx())

    drawer = rdMolDraw2D.MolDraw2DCairo(350,300)
    drawer.drawOptions().fillHighlights=True
    drawer.drawOptions().setHighlightColour((rgba_color))
    drawer.drawOptions().highlightBondWidthMultiplier=20
    drawer.drawOptions().clearBackground = False
    rdMolDraw2D.PrepareAndDrawMolecule(drawer, mol, highlightAtoms=atoms, highlightBonds=bonds)
    bio = io.BytesIO(drawer.GetDrawingText())
    Image.open(bio)

### Highlight Molecule with Multiple Colors

Author: Vincent Scalfani
Original Source: Adapted from http://rdkit.blogspot.com/2020/04/new-drawing-options-in-202003-release.html
Index ID#: RDKitCB_39
Summary: Highlight a molecule with different colors based on if the atom/bond is aromatic.

    from rdkit import Chem
    from rdkit.Chem.Draw import rdMolDraw2D
    import io
    from PIL import Image
    from collections import defaultdict

    mol = Chem.MolFromSmiles('CC1=CC(=CC=C1)NC(=O)CCC2=CC=CC=C2')
    colors = [(0.0, 0.0, 1.0, 0.1), (1.0, 0.0, 0.0, 0.2)]

    athighlights = defaultdict(list)
    arads = {}
    for a in mol.GetAtoms():
        if a.GetIsAromatic():
            aid = a.GetIdx()
            athighlights[aid].append(colors[0])
            arads[aid] = 0.3
        else:
            aid = a.GetIdx()
            athighlights[aid].append(colors[1])
            arads[aid] = 0.3

    bndhighlights = defaultdict(list)
    for bond in mol.GetBonds():
        aid1 = bond.GetBeginAtomIdx()
        aid2 = bond.GetEndAtomIdx()

        if bond.GetIsAromatic():
            bid = mol.GetBondBetweenAtoms(aid1,aid2).GetIdx()
            bndhighlights[bid].append(colors[0])
        else:
            bid = mol.GetBondBetweenAtoms(aid1,aid2).GetIdx()
            bndhighlights[bid].append(colors[1])

    d2d = rdMolDraw2D.MolDraw2DCairo(350,400)
    d2d.DrawMoleculeWithHighlights(mol,"",dict(athighlights),dict(bndhighlights),arads,{})
    d2d.FinishDrawing()
    bio = io.BytesIO(d2d.GetDrawingText())
    Image.open(bio)

### Without Implicit Hydrogens

Author: Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/36746387/ and https://gist.github.com/greglandrum/9283aeadfb66d0fe8a2900e63fb10f3e
Index ID#: RDKitCB_17
Summary: Draw a molecule without implicit hydrogens

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    m = Chem.MolFromSmiles('[Pt](Cl)(Cl)(N)N')
    m

    for atom in m.GetAtoms():
        atom.SetProp("atomLabel", atom.GetSymbol())
    m

### With Abbreviations

Author: Greg Landrum
Source: https://github.com/rdkit/UGM_2020/blob/master/Notebooks/Landrum_WhatsNew.ipynb
Index ID#: RDKitCB_34
Summary: Draw a molecule with functional group abbreviations

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    from rdkit.Chem import Draw
    from rdkit.Chem import rdAbbreviations

    m = Chem.MolFromSmiles('COc1ccc(C(=O)[O-])cc1')
    m

    abbrevs = rdAbbreviations.GetDefaultAbbreviations()
    nm = rdAbbreviations.CondenseMolAbbreviations(m,abbrevs)
    nm

    # abbreviations that cover more than 40% of the molecule won't be applied by default
    m = Chem.MolFromSmiles('c1c[nH]cc1C(F)(F)F')
    nm1 = rdAbbreviations.CondenseMolAbbreviations(m,abbrevs)
    nm2 = rdAbbreviations.CondenseMolAbbreviations(m,abbrevs,maxCoverage=0.8)
    Draw.MolsToGridImage((m,nm1,nm2),legends=('','default','maxCoverage=0.8'))

    # See available abbreviations and their SMILES
    # where * is the dummy atom that the group would attach to
    abbrevs = rdAbbreviations.GetDefaultAbbreviations()
    labels = ["Abbrev", "SMILES"]
    line = '--------'

    print(f"{labels[0]:<10} {labels[1]}")
    print(f"{line:<10} {line}")
    for a in abbrevs:
       print(f"{a.label:<10} {Chem.MolToSmiles(a.mol)}")

    Abbrev     SMILES
    --------   --------
    CO2Et      *C(=O)OCC
    COOEt      *C(=O)OCC
    OiBu       *OCC(C)C
    nDec       *CCCCCCCCCC
    nNon       *CCCCCCCCC
    nOct       *CCCCCCCC
    nHept      *CCCCCCC
    nHex       *CCCCCC
    nPent      *CCCCC
    iPent      *C(C)CCC
    tBu        *C(C)(C)C
    iBu        *C(C)CC
    nBu        *CCCC
    iPr        *C(C)C
    nPr        *CCC
    Et         *CC
    NCF3       *NC(F)(F)F
    CF3        *C(F)(F)F
    CCl3       *C(Cl)(Cl)Cl
    CN         *C#N
    NC         *[N+]#[C-]
    N(OH)CH3   *N(C)[OH]
    NO2        *[N+](=O)[O-]
    NO         *N=O
    SO3H       *S(=O)(=O)[OH]
    CO2H       *C(=O)[OH]
    COOH       *C(=O)[OH]
    OEt        *OCC
    OAc        *OC(C)=O
    NHAc       *NC(C)=O
    Ac         *C(C)=O
    CHO        *C=O
    NMe        *NC
    SMe        *SC
    OMe        *OC
    CO2-       *C(=O)[O-]
    COO-       *C(=O)[O-]

### Using CoordGen Library

Author: Greg Landrum
Source: https://github.com/rdkit/UGM_2020/blob/master/Notebooks/Landrum_WhatsNew.ipynb
Index ID#: RDKitCB_37
Summary: Draw a molecule using CoordGen Library

Some molecules like macrocycles are not represented well using the default RDKit drawing code. As a result, it may be preferable to use the CoordGen integration.

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    IPythonConsole.molSize = 350,300
    from rdkit.Chem import Draw

    # default drawing
    mol = Chem.MolFromSmiles("C/C=C/CC(C)C(O)C1C(=O)NC(CC)C(=O)N(C)CC(=O)N(C)C(CC(C)C)C(=O)NC(C(C)C)C(=O)N(C)C(CC(C)C)C(=O)NC(C)C(=O)NC(C)C(=O)N(C)C(CC(C)C)C(=O)N(C)C(CC(C)C)C(=O)N(C)C(C(C)C)C(=O)N1C")
    mol

    # with CoordGen
    from rdkit.Chem import rdCoordGen
    rdCoordGen.AddCoords(mol)
    mol

It is also possible to use CoordGen with the MolDraw2D class. Here is one way to do that:

    from rdkit import Chem
    from rdkit.Chem import Draw
    from rdkit.Chem.Draw import rdMolDraw2D
    from rdkit.Chem import rdDepictor
    rdDepictor.SetPreferCoordGen(True)
    from rdkit.Chem.Draw import IPythonConsole
    from IPython.display import SVG

    mol = Chem.MolFromSmiles("C/C=C/CC(C)C(O)C1C(=O)NC(CC)C(=O)N(C)CC(=O)N(C)C(CC(C)C)C(=O)NC(C(C)C)C(=O)N(C)C(CC(C)C)C(=O)NC(C)C(=O)NC(C)C(=O)N(C)C(CC(C)C)C(=O)N(C)C(CC(C)C)C(=O)N(C)C(C(C)C)C(=O)N1C")
    drawer = rdMolDraw2D.MolDraw2DSVG(300,300)
    drawer.drawOptions().addStereoAnnotation = False
    drawer.DrawMolecule(mol)
    drawer.FinishDrawing()
    SVG(drawer.GetDrawingText())

### On a Plot

Author: Takayuki Serizawa
Original Source: https://gist.github.com/iwatobipen/1b384d145024663151b3252bf16d2aa8
Index ID#: RDKitCB_35
Summary: Draw a molecule on a matplotlib plot.

    import matplotlib.pyplot as plt
    import numpy as np
    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole

    x = np.arange(0, 180, 1)
    y = np.sin(x)

    mol = Chem.MolFromSmiles('C1CNCCC1C(=O)C')
    im = Chem.Draw.MolToImage(mol)

    fig = plt.figure(figsize=(10,5))
    plt.plot(x, y)
    plt.ylim(-1, 5)
    ax = plt.axes([0.6, 0.47, 0.38, 0.38], frameon=True)
    ax.imshow(im)
    ax.axis('off')
    # plt.show() # commented out to avoid creating plot with doctest

## Bonds and Bonding

### Hybridization Type and Count

Author: Jean-Marc Nuzillard and Andrew Dalke
Source: https://sourceforge.net/p/rdkit/mailman/message/37021222/ and https://sourceforge.net/p/rdkit/mailman/message/37025886/
Index ID#: RDKitCB_26
Summary: Get hybridization type and count

    from rdkit import Chem
    m = Chem.MolFromSmiles("CN1C=NC2=C1C(=O)N(C(=O)N2C)C")
    for x in m.GetAtoms():
        print(x.GetIdx(), x.GetHybridization())

    0 SP3
    1 SP2
    2 SP2
    3 SP2
    4 SP2
    5 SP2
    6 SP2
    7 SP2
    8 SP2
    9 SP2
    10 SP2
    11 SP2
    12 SP3
    13 SP3

    # if you want to count hybridization type (e.g., SP3):
    from rdkit import Chem
    m = Chem.MolFromSmiles("CN1C=NC2=C1C(=O)N(C(=O)N2C)C")
    print(sum((x.GetHybridization() == Chem.HybridizationType.SP3) for x in m.GetAtoms()))

    3

## Rings, Aromaticity, and Kekulization

### Count Ring Systems

Author: Greg Landrum
Source: https://gist.github.com/greglandrum/de1751a42b3cae54011041dd67ae7415
Index ID#: RDKitCB_3
Summary: Count ring systems in a molecule

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole

    def GetRingSystems(mol, includeSpiro=False):
        ri = mol.GetRingInfo()
        systems = []
        for ring in ri.AtomRings():
            ringAts = set(ring)
            nSystems = []
            for system in systems:
                nInCommon = len(ringAts.intersection(system))
                if nInCommon and (includeSpiro or nInCommon>1):
                    ringAts = ringAts.union(system)
                else:
                    nSystems.append(system)
            nSystems.append(ringAts)
            systems = nSystems
        return systems
    mol = Chem.MolFromSmiles('CN1C(=O)CN=C(C2=C1C=CC(=C2)Cl)C3=CC=CC=C3')
    print(GetRingSystems(mol))

    [{1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 12}, {14, 15, 16, 17, 18, 19}]

    # Draw molecule with atom index (see RDKitCB_0)
    def mol_with_atom_index(mol):
        for atom in mol.GetAtoms():
            atom.SetAtomMapNum(atom.GetIdx())
        return mol
    mol_with_atom_index(mol)

### Identify Aromatic Rings

Author: Benjamin Datko and Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/36860045/ and https://sourceforge.net/p/rdkit/mailman/message/23801106/
Index ID#: RDKitCB_8
Summary: Identify which rings are aromatic in a molecule

    from rdkit import Chem
    m = Chem.MolFromSmiles('c1cccc2c1CCCC2')
    m

    ri = m.GetRingInfo()
    # You can interrogate the RingInfo object to tell you the atoms that make up each ring:
    print(ri.AtomRings())

    ((0, 5, 4, 3, 2, 1), (6, 7, 8, 9, 4, 5))

    # or the bonds that make up each ring:
    print(ri.BondRings())

    ((9, 4, 3, 2, 1, 0), (6, 7, 8, 10, 4, 5))

    # To detect aromatic rings, I would loop over the bonds in each ring and
    # flag the ring as aromatic if all bonds are aromatic:
    def isRingAromatic(mol, bondRing):
            for id in bondRing:
                if not mol.GetBondWithIdx(id).GetIsAromatic():
                    return False
            return True

    print(isRingAromatic(m, ri.BondRings()[0]))

    True

    print(isRingAromatic(m, ri.BondRings()[1]))

    False

### Identify Aromatic Atoms

Author: Paolo Tosco
Original Source: https://sourceforge.net/p/rdkit/mailman/message/36862879/
Index ID#: RDKitCB_9
Summary: Differentiate aromatic carbon from olefinic carbon with SMARTS

    from rdkit import Chem
    mol = Chem.MolFromSmiles("c1ccccc1C=CCC")
    aromatic_carbon = Chem.MolFromSmarts("c")
    print(mol.GetSubstructMatches(aromatic_carbon))

    ((0,), (1,), (2,), (3,), (4,), (5,))

    # The RDKit includes a SMARTS extension that allows hybridization queries,
    # here we query for SP2 aliphatic carbons:
    olefinic_carbon = Chem.MolFromSmarts("[C^2]")
    print(mol.GetSubstructMatches(olefinic_carbon))

    ((6,), (7,))

There is also an alternative, more efficient approach, using the rdqueries module:

    from rdkit import Chem
    from rdkit.Chem import rdqueries

    mol = Chem.MolFromSmiles("c1ccccc1C=CCC")
    q = rdqueries.IsAromaticQueryAtom()
    print([x.GetIdx() for x in mol.GetAtomsMatchingQuery(q)])

    [0, 1, 2, 3, 4, 5]

    q = rdqueries.HybridizationEqualsQueryAtom(Chem.HybridizationType.SP2)
    print([x.GetIdx() for x in mol.GetAtomsMatchingQuery(q)])

    [0, 1, 2, 3, 4, 5, 6, 7]

    qcombined = rdqueries.IsAliphaticQueryAtom()
    qcombined.ExpandQuery(q)
    print([x.GetIdx() for x in mol.GetAtomsMatchingQuery(qcombined)])

    [6, 7]

## Stereochemistry

### Identifying Stereochemistry

Author: Vincent Scalfani
Source: https://github.com/rdkit/UGM_2020/blob/master/Notebooks/Landrum_WhatsNew.ipynb
Index ID#: RDKitCB_30
Summary: Find chiral centers and double bond stereochemistry.

    from rdkit import Chem
    from rdkit.Chem import Draw
    from rdkit.Chem.Draw import IPythonConsole
    IPythonConsole.drawOptions.addAtomIndices = True
    IPythonConsole.drawOptions.addStereoAnnotation = False
    IPythonConsole.molSize = 200,200

    m = Chem.MolFromSmiles("C[C@H]1CCC[C@@H](C)[C@@H]1Cl")
    m

    # legacy FindMolChiralCenters()
    print(Chem.FindMolChiralCenters(m,force=True,includeUnassigned=True,useLegacyImplementation=True))

    [(1, 'S'), (5, 'R'), (7, 'R')]

    # new stereochemistry code
    print(Chem.FindMolChiralCenters(m,force=True,includeUnassigned=True,useLegacyImplementation=False))

    [(1, 'S'), (5, 'R'), (7, 'r')]

    # Identifying Double Bond Stereochemistry
    IPythonConsole.molSize = 250,250
    mol = Chem.MolFromSmiles("C\C=C(/F)\C(=C\F)\C=C")
    mol

    # Using GetStereo()
    for b in mol.GetBonds():
        print(b.GetBeginAtomIdx(),b.GetEndAtomIdx(),
              b.GetBondType(),b.GetStereo())

    0 1 SINGLE STEREONONE
    1 2 DOUBLE STEREOZ
    2 3 SINGLE STEREONONE
    2 4 SINGLE STEREONONE
    4 5 DOUBLE STEREOE
    5 6 SINGLE STEREONONE
    4 7 SINGLE STEREONONE
    7 8 DOUBLE STEREONONE

    # Double bond configuration can also be identified with new
    # stereochemistry code using Chem.FindPotentialStereo()
    si = Chem.FindPotentialStereo(mol)
    for element in si:
        print(f'  Type: {element.type}, Which: {element.centeredOn}, Specified: {element.specified}, Descriptor: {element.descriptor} ')

    Type: Bond_Double, Which: 1, Specified: Specified, Descriptor: Bond_Cis
    Type: Bond_Double, Which: 4, Specified: Specified, Descriptor: Bond_Trans

## Manipulating Molecules

### Create Fragments

Author: Paulo Tosco
Source: https://sourceforge.net/p/rdkit/mailman/message/36895168/ and https://gist.github.com/ptosco/3fb93b7c09dac15b6d355eb0ad29f532
Index ID#: RDKitCB_7
Summary: Create fragments of molecules on bonds

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole, MolsToGridImage
    # I have put explicit bonds in the SMILES definition to facilitate comprehension:
    mol = Chem.MolFromSmiles("O-C-C-C-C-N")
    mol1 = Chem.Mol(mol)
    mol2 = Chem.Mol(mol)
    mol1

    # Chem.FragmentOnBonds() will fragment all specified bond indices at once, and return a single molecule
    # with all specified cuts applied. By default, addDummies=True, so empty valences are filled with dummy atoms:
    mol1_f = Chem.FragmentOnBonds(mol1, (0, 2, 4))
    mol1_f

    # This molecule can be split into individual fragments using Chem.GetMolFrags():
    MolsToGridImage(Chem.GetMolFrags(mol1_f, asMols=True))

    # Chem.FragmentOnSomeBonds() will fragment according to all permutations of numToBreak bonds at a time
    # (numToBreak defaults to 1), and return tuple of molecules with numToBreak cuts applied. By default,
    # addDummies=True, so empty valences are filled with dummy atoms:
    mol2_f_tuple = Chem.FragmentOnSomeBonds(mol2, (0, 2, 4))

    mol2_f_tuple[0]

    mol2_f_tuple[1]

    mol2_f_tuple[2]

    # Finally, you can manually cut bonds using Chem.RWMol.RemoveBonds:
    with Chem.RWMol(mol) as rwmol:
      for b_idx in [0, 2, 4]:
        b = rwmol.GetBondWithIdx(b_idx)
        rwmol.RemoveBond(b.GetBeginAtomIdx(), b.GetEndAtomIdx())
    # And then call Chem.GetMolFrags() to get sanitized fragments where empty valences were filled with implicit hydrogens:
    MolsToGridImage(Chem.GetMolFrags(rwmol, asMols=True))

### Largest Fragment

Author: Andrew Dalke and Susan Leung
Source: https://sourceforge.net/p/rdkit/mailman/message/36355644/ and https://github.com/susanhleung/rdkit/blob/dev/GSOC2018_MolVS_Integration/rdkit/Chem/MolStandardize/tutorial/MolStandardize.ipynb
Index ID#: RDKitCB_31
Summary: Select largest fragment from a molecule

    from rdkit import Chem
    from rdkit.Chem import rdmolops
    mol = Chem.MolFromSmiles('CCOC(=O)C(C)(C)OC1=CC=C(C=C1)Cl.CO.C1=CC(=CC=C1C(=O)N[C@@H](CCC(=O)O)C(=O)O)NCC2=CN=C3C(=N2)C(=O)NC(=N3)N')

    mol_frags = rdmolops.GetMolFrags(mol, asMols = True)
    largest_mol = max(mol_frags, default=mol, key=lambda m: m.GetNumAtoms())
    print(Chem.MolToSmiles(largest_mol))

    Nc1nc2ncc(CNc3ccc(C(=O)N[C@@H](CCC(=O)O)C(=O)O)cc3)nc2c(=O)[nH]1

The same result can also be achieved with MolStandardize:

    from rdkit import Chem
    from rdkit.Chem.MolStandardize import rdMolStandardize
    mol = Chem.MolFromSmiles('CCOC(=O)C(C)(C)OC1=CC=C(C=C1)Cl.CO.C1=CC(=CC=C1C(=O)N[C@@H](CCC(=O)O)C(=O)O)NCC2=CN=C3C(=N2)C(=O)NC(=N3)N')

    # setup standardization module
    largest_Fragment = rdMolStandardize.LargestFragmentChooser()
    largest_mol = largest_Fragment.choose(mol)
    print(Chem.MolToSmiles(largest_mol))

    Nc1nc2ncc(CNc3ccc(C(=O)N[C@@H](CCC(=O)O)C(=O)O)cc3)nc2c(=O)[nH]1

### Sidechain-Core Enumeration

Author: Chris Earnshaw, Stephen Roughley, Greg Landrum (Vincent Scalfani added loop example)
Source: https://sourceforge.net/p/rdkit/mailman/message/35730514/ and https://www.rdkit.org/docs/source/rdkit.Chem.rdChemReactions.html
Index ID#: RDKitCB_29
Summary: Replace sidechains on a core and enumerate the combinations.

    from rdkit import Chem
    from rdkit.Chem import Draw
    from rdkit.Chem import AllChem

    # core is '*c1c(C)cccc1(O)'
    # chain is 'CN*'

    rxn = AllChem.ReactionFromSmarts('[c:1][#0].[#0][*:2]>>[c:1]-[*:2]')
    reacts = (Chem.MolFromSmiles('*c1c(C)cccc1(O)'),Chem.MolFromSmiles('CN*'))
    products = rxn.RunReactants(reacts) # tuple
    print(len(products))

    1

    print(len(products[0]))

    1

    print(Chem.MolToSmiles(products[0][0])) # [0][0] to index out the rdchem mol object

    CNc1c(C)cccc1O

    # The above reaction-based approach is flexible, however if you can generate your
    # sidechains in such a way that the atom you want to attach to the core
    # is the first one (atom zero), there's a somewhat easier way to do this
    # kind of simple replacement:

    core = Chem.MolFromSmiles('*c1c(C)cccc1(O)')
    chain = Chem.MolFromSmiles('NC')
    products = Chem.ReplaceSubstructs(core,Chem.MolFromSmarts('[#0]'),chain) # tuple
    print(Chem.MolToSmiles(products[0]))

    CNc1c(C)cccc1O

    # Here is an example in a loop for an imidazolium core with alkyl chains

    core = Chem.MolFromSmiles('*[n+]1cc[nH]c1')
    chains = ['C','CC','CCC','CCCC','CCCCC','CCCCCC']
    chainMols = [Chem.MolFromSmiles(chain) for chain in chains]

    product_smi = []
    for chainMol in chainMols:
        product_mol = Chem.ReplaceSubstructs(core,Chem.MolFromSmarts('[#0]'),chainMol)
        product_smi.append(Chem.MolToSmiles(product_mol[0]))
    print(product_smi)

    ['C[n+]1cc[nH]c1', 'CC[n+]1cc[nH]c1', 'CCC[n+]1cc[nH]c1', 'CCCC[n+]1cc[nH]c1', 'CCCCC[n+]1cc[nH]c1', 'CCCCCC[n+]1cc[nH]c1']

    # View the enumerated molecules:
    Draw.MolsToGridImage([Chem.MolFromSmiles(smi) for smi in product_smi])

### Neutralizing Molecules

Author: Noel O’Boyle (Vincent Scalfani adapted code for RDKit)
Source: https://baoilleach.blogspot.com/2019/12/no-charge-simple-approach-to.html
Index ID#: RDKitCB_33
Summary: Neutralize charged molecules by atom.

This `neutralize_atoms()` algorithm is adapted from Noel O’Boyle’s nocharge code. It is a neutralization by atom approach and neutralizes atoms with a +1 or -1 charge by removing or adding hydrogen where possible. The SMARTS pattern checks for a hydrogen in +1 charged atoms and checks for no neighbors with a negative charge (for +1 atoms) and no neighbors with a positive charge (for -1 atoms), this is to avoid altering molecules with charge separation (e.g., nitro groups).

The `neutralize_atoms()` function differs from the `rdMolStandardize.Uncharger` behavior. See the MolVS documentation for Uncharger:

https://molvs.readthedocs.io/en/latest/api.html#molvs-charge

“This class uncharges molecules by adding and/or removing hydrogens. In cases where there is a positive charge that is not neutralizable, any corresponding negative charge is also preserved.”

As an example, `rdMolStandardize.Uncharger` will not change charges on `C[N+](C)(C)CCC([O-])=O`, as there is a positive charge that is not neutralizable. In contrast, the `neutralize_atoms()` function will attempt to neutralize any atoms it can (in this case to `C[N+](C)(C)CCC(=O)O`). That is, `neutralize_atoms()` ignores the overall charge on the molecule, and attempts to neutralize charges even if the neutralization introduces an overall formal charge on the molecule. See below for a comparison.

    from rdkit import Chem
    from rdkit.Chem import AllChem
    from rdkit.Chem import Draw

    # list of SMILES
    smiList = ['CC(CNC[O-])[N+]([O-])=O',
           'C[N+](C)(C)CCC([O-])=O',
           '[O-]C1=CC=[N+]([O-])C=C1',
           '[O-]CCCN=[N+]=[N-]',
           'C[NH+](C)CC[S-]',
           'CP([O-])(=O)OC[NH3+]']

    # Create RDKit molecular objects
    mols = [Chem.MolFromSmiles(m) for m in smiList]

    # display
    Draw.MolsToGridImage(mols,molsPerRow=3,subImgSize=(200,200))

    def neutralize_atoms(mol):
        pattern = Chem.MolFromSmarts("[+1!h0!$([*]~[-1,-2,-3,-4]),-1!$([*]~[+1,+2,+3,+4])]")
        at_matches = mol.GetSubstructMatches(pattern)
        at_matches_list = [y[0] for y in at_matches]
        if len(at_matches_list) > 0:
            for at_idx in at_matches_list:
                atom = mol.GetAtomWithIdx(at_idx)
                chg = atom.GetFormalCharge()
                hcount = atom.GetTotalNumHs()
                atom.SetFormalCharge(0)
                atom.SetNumExplicitHs(hcount - chg)
                atom.UpdatePropertyCache()
        return mol

    # Neutralize molecules by atom
    for mol in mols:
        neutralize_atoms(mol)
        print(Chem.MolToSmiles(mol))

    CC(CNCO)[N+](=O)[O-]
    C[N+](C)(C)CCC(=O)O
    [O-][n+]1ccc(O)cc1
    [N-]=[N+]=NCCCO
    CN(C)CCS
    CP(=O)(O)OCN

    Draw.MolsToGridImage(mols,molsPerRow=3, subImgSize=(200,200))

Compare to `rdMolStandardize.Uncharger` results:

    from rdkit.Chem.MolStandardize import rdMolStandardize
    un = rdMolStandardize.Uncharger()
    mols2 = [Chem.MolFromSmiles(m) for m in smiList]

    for mol2 in mols2:
        un.uncharge(mol2)
        print(Chem.MolToSmiles(mol2))

    CC(CNC[O-])[N+](=O)[O-]
    C[N+](C)(C)CCC(=O)[O-]
    [O-]c1cc[n+]([O-])cc1
    [N-]=[N+]=NCCC[O-]
    C[NH+](C)CC[S-]
    CP(=O)([O-])OC[NH3+]

    Draw.MolsToGridImage(mols2,molsPerRow=3,subImgSize=(200,200))

## Substructure Matching

### Functional Group with SMARTS queries

Author: Paulo Tosco
Source: https://sourceforge.net/p/rdkit/mailman/message/36810326/
Index ID#: RDKitCB_10
Summary: Match a functional group (e.g., alcohol) with a SMARTS query

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    sucrose = "C([C@@H]1[C@H]([C@@H]([C@H]([C@H](O1)O[C@]2([C@H]([C@@H]([C@H](O2)CO)O)O)CO)O)O)O)O"
    sucrose_mol = Chem.MolFromSmiles(sucrose)
    primary_alcohol = Chem.MolFromSmarts("[CH2][OH1]")
    print(sucrose_mol.GetSubstructMatches(primary_alcohol))

    ((0, 22), (13, 14), (17, 18))

    secondary_alcohol = Chem.MolFromSmarts("[CH1][OH1]")
    print(sucrose_mol.GetSubstructMatches(secondary_alcohol))

    ((2, 21), (3, 20), (4, 19), (9, 16), (10, 15))

### Macrocycles with SMARTS queries

Author: Ivan Tubert-Brohman and David Cosgrove (Vincent Scalfani added example)
Source: https://sourceforge.net/p/rdkit/mailman/message/36781480/
Index ID#: RDKitCB_13
Summary: Match a macrocycle ring with a SMARTS query

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole
    from rdkit.Chem import Draw
    erythromycin = Chem.MolFromSmiles("CC[C@@H]1[C@@]([C@@H]([C@H](C(=O)[C@@H](C[C@@]([C@@H]([C@H]([C@@H]([C@H](C(=O)O1)C)O[C@H]2C[C@@]([C@H]([C@@H](O2)C)O)(C)OC)C)O[C@H]3[C@@H]([C@H](C[C@H](O3)C)N(C)C)O)(C)O)C)C)O)(C)O")
    erythromycin

    # Define SMARTS pattern with ring size > 12
    # This is an RDKit SMARTS extension
    macro = Chem.MolFromSmarts("[r{12-}]")
    print(erythromycin.GetSubstructMatches(macro))

    ((2,), (3,), (4,), (5,), (6,), (8,), (9,), (10,), (11,), (12,), (13,), (14,), (15,), (17,))

    erythromycin

### Returning Substructure Matches as SMILES

Author: Andrew Dalke
Source: https://sourceforge.net/p/rdkit/mailman/message/36735316/
Index ID#: RDKitCB_18
Summary: Obtain SMILES of the matched substructure.

    from rdkit import Chem
    pat = Chem.MolFromSmarts("[NX1]#[CX2]") #matches nitrile
    mol = Chem.MolFromSmiles("CCCC#N") # Butyronitrile
    atom_indices = mol.GetSubstructMatch(pat)
    print(atom_indices)

    (4, 3)

    print(Chem.MolFragmentToSmiles(mol, atom_indices)) # returns the nitrile

    C#N

    # Note however that if only the atom indices are given then Chem.MolFragmentToSmiles() will include all bonds
    # which connect those atoms, even if the original SMARTS does not match those bonds. For example:
    pat = Chem.MolFromSmarts("*~*~*~*") # match 4 linear atoms
    mol = Chem.MolFromSmiles("C1CCC1") # ring of size 4
    atom_indices = mol.GetSubstructMatch(pat)
    print(atom_indices)

    (0, 1, 2, 3)

    print(Chem.MolFragmentToSmiles(mol, atom_indices))  # returns the ring

    C1CCC1

    # If this is important, then you need to pass the correct bond indices to MolFragmentToSmiles().
    # This can be done by using the bonds in the query graph to get the bond indices in the molecule graph.
    def get_match_bond_indices(query, mol, match_atom_indices):
        bond_indices = []
        for query_bond in query.GetBonds():
            atom_index1 = match_atom_indices[query_bond.GetBeginAtomIdx()]
            atom_index2 = match_atom_indices[query_bond.GetEndAtomIdx()]
            bond_indices.append(mol.GetBondBetweenAtoms(
                 atom_index1, atom_index2).GetIdx())
        return bond_indices

    bond_indices = get_match_bond_indices(pat, mol, atom_indices)
    print(bond_indices)

    [0, 1, 2]

    print(Chem.MolFragmentToSmiles(mol, atom_indices, bond_indices))

    CCCC

### Within the Same Fragment

Author: Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/36942946/
Index ID#: RDKitCB_20
Summary: Match a pattern only within the same fragment.

    from rdkit import Chem

    p = Chem.MolFromSmarts('O.N')

    # define a function where matches are contained in a single fragment
    def fragsearch(m,p):
        matches = [set(x) for x in m.GetSubstructMatches(p)]
        frags = [set(y) for y in Chem.GetMolFrags(m)] # had to add this line for code to work
        for frag in frags:
            for match in matches:
                if match.issubset(frag):
                    return match
        return False

    m1 = Chem.MolFromSmiles('OCCCN.CCC')
    m1

    m2 = Chem.MolFromSmiles('OCCC.CCCN')
    m2

    print(m1.HasSubstructMatch(p))

    True

    print(m2.HasSubstructMatch(p))

    True

    print(fragsearch(m1,p))

    {0, 4}

    print(fragsearch(m2,p))

    False

## Descriptor Calculations

### Molecule Hash Strings

Author: Vincent Scalfani and Takayuki Serizawa
Source: https://gist.github.com/vfscalfani/f77d90f9f27e0f820b966882cdadccd0 and https://iwatobipen.wordpress.com/2019/10/27/a-new-function-of-rdkit201909-rdkit-chemoinformatics/
Index ID#: RDKitCB_21
Summary: Calculate hash strings for molecules with the NextMove MolHash functionality within RDKit.
Reference Note: Examples from O’Boyle and Sayle [2]

    from rdkit import Chem
    from rdkit.Chem import rdMolHash
    import rdkit

    s = Chem.MolFromSmiles('CC(C(C1=CC(=C(C=C1)O)O)O)N(C)C(=O)OCC2=CC=CC=C2')
    s

    #  View all of the MolHash hashing functions types with the names method.
    molhashf = rdMolHash.HashFunction.names
    print(molhashf)

    {'AnonymousGraph': rdkit.Chem.rdMolHash.HashFunction.AnonymousGraph,
     'ElementGraph': rdkit.Chem.rdMolHash.HashFunction.ElementGraph,
     'CanonicalSmiles': rdkit.Chem.rdMolHash.HashFunction.CanonicalSmiles,
     'MurckoScaffold': rdkit.Chem.rdMolHash.HashFunction.MurckoScaffold,
     'ExtendedMurcko': rdkit.Chem.rdMolHash.HashFunction.ExtendedMurcko,
     'MolFormula': rdkit.Chem.rdMolHash.HashFunction.MolFormula,
     'AtomBondCounts': rdkit.Chem.rdMolHash.HashFunction.AtomBondCounts,
     'DegreeVector': rdkit.Chem.rdMolHash.HashFunction.DegreeVector,
     'Mesomer': rdkit.Chem.rdMolHash.HashFunction.Mesomer,
     'HetAtomTautomer': rdkit.Chem.rdMolHash.HashFunction.HetAtomTautomer,
     'HetAtomProtomer': rdkit.Chem.rdMolHash.HashFunction.HetAtomProtomer,
     'RedoxPair': rdkit.Chem.rdMolHash.HashFunction.RedoxPair,
     'Regioisomer': rdkit.Chem.rdMolHash.HashFunction.Regioisomer,
     'NetCharge': rdkit.Chem.rdMolHash.HashFunction.NetCharge,
     'SmallWorldIndexBR': rdkit.Chem.rdMolHash.HashFunction.SmallWorldIndexBR,
     'SmallWorldIndexBRL': rdkit.Chem.rdMolHash.HashFunction.SmallWorldIndexBRL,
     'ArthorSubstructureOrder': rdkit.Chem.rdMolHash.HashFunction.ArthorSubstructureOrder,
     'HetAtomTautomerv2': rdkit.Chem.rdMolHash.HashFunction.HetAtomTautomerv2}

    # Generate MolHashes for molecule 's' with all defined hash functions.
    for i, j in molhashf.items():
        print(i, rdMolHash.MolHash(s, j))

    AnonymousGraph **(***1*****1)*(*)*(*)*(*)*1***(*)*(*)*1
    ElementGraph CC(C(O)C1CCC(O)C(O)C1)N(C)C(O)OCC1CCCCC1
    CanonicalSmiles CC(C(O)c1ccc(O)c(O)c1)N(C)C(=O)OCc1ccccc1
    MurckoScaffold c1ccc(CCNCOCc2ccccc2)cc1
    ExtendedMurcko *c1ccc(C(*)C(*)N(*)C(=*)OCc2ccccc2)cc1*
    MolFormula C18H21NO5
    AtomBondCounts 24,25
    DegreeVector 0,8,10,6
    Mesomer CC(C(O)[C]1[CH][CH][C](O)[C](O)[CH]1)N(C)[C]([O])OC[C]1[CH][CH][CH][CH][CH]1_0
    HetAtomTautomer CC(C([O])[C]1[CH][CH][C]([O])[C]([O])[CH]1)N(C)[C]([O])OC[C]1[CH][CH][CH][CH][CH]1_3_0
    HetAtomProtomer CC(C([O])[C]1[CH][CH][C]([O])[C]([O])[CH]1)N(C)[C]([O])OC[C]1[CH][CH][CH][CH][CH]1_3
    RedoxPair CC(C(O)[C]1[CH][CH][C](O)[C](O)[CH]1)N(C)[C]([O])OC[C]1[CH][CH][CH][CH][CH]1
    Regioisomer *C.*CCC.*O.*O.*O.*OC(=O)N(*)*.C.c1ccccc1.c1ccccc1
    NetCharge 0
    SmallWorldIndexBR B25R2
    SmallWorldIndexBRL B25R2L10
    ArthorSubstructureOrder 00180019010012000600009b000000
    HetAtomTautomerv2 [CH3]-[CH](-[CH](-[OH])-[C]1:[C]:[C]:[C](:[O]):[C](:[O]):[C]:1)-[N](-[CH3])-[C](=[O])-[O]-[CH2]-[c]1:[cH]:[cH]:[cH]:[cH]:[cH]:1_5_0

    # Murcko Scaffold Hashes (from slide 16 in [ref2])
    # Create a list of SMILES
    mList = ['CCC1CC(CCC1=O)C(=O)C1=CC=CC(C)=C1','CCC1CC(CCC1=O)C(=O)C1=CC=CC=C1',\
             'CC(=C)C(C1=CC=CC=C1)S(=O)CC(N)=O','CC1=CC(=CC=C1)C(C1CCC(N)CC1)C(F)(F)F',\
             'CNC1CCC(C2=CC(Cl)=C(Cl)C=C2)C2=CC=CC=C12','CCCOC(C1CCCCC1)C1=CC=C(Cl)C=C1']

    # Loop through the SMILES mList and create RDKit molecular objects
    mMols = [Chem.MolFromSmiles(m) for m in mList]
    # Calculate Murcko Scaffold Hashes
    murckoHashList = [rdMolHash.MolHash(mMol, rdkit.Chem.rdMolHash.HashFunction.MurckoScaffold) for mMol in mMols]
    print(murckoHashList)

    ['c1ccc(CC2CCCCC2)cc1',
     'c1ccc(CC2CCCCC2)cc1',
     'c1ccccc1',
     'c1ccc(CC2CCCCC2)cc1',
     'c1ccc(C2CCCc3ccccc32)cc1',
     'c1ccc(CC2CCCCC2)cc1']

    # Get the most frequent Murcko Scaffold Hash
    def mostFreq(list):
        return max(set(list), key=list.count)
    mostFreq_murckoHash = mostFreq(murckoHashList)
    print(mostFreq_murckoHash)

    c1ccc(CC2CCCCC2)cc1

    mostFreq_murckoHash_mol = Chem.MolFromSmiles('c1ccc(CC2CCCCC2)cc1')
    mostFreq_murckoHash_mol

    # Display molecules with MurkoHash as legends and highlight the mostFreq_murckoHash
    highlight_mostFreq_murckoHash = [mMol.GetSubstructMatch(mostFreq_murckoHash_mol) for mMol in mMols]
    Draw.MolsToGridImage(mMols,legends=[murckoHash for murckoHash in murckoHashList],
                         highlightAtomLists = highlight_mostFreq_murckoHash,
                         subImgSize=(250,250), useSVG=False)

    # Regioisomer Hashes (from slide 17 in [ref2])
    # Find Regioisomer matches for this molecule
    r0 = Chem.MolFromSmiles('CC1=CC2=C(C=C1)C(=CN2CCN1CCOCC1)C(=O)C1=CC=CC2=C1C=CC=C2')
    r0

    # Calculate the regioisomer hash for r0
    r0_regioHash = rdMolHash.MolHash(r0,rdkit.Chem.rdMolHash.HashFunction.Regioisomer)
    print(r0_regioHash)

    *C.*C(*)=O.*CC*.C1COCCN1.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1

    r0_regioHash_mol = Chem.MolFromSmiles('*C.*C(*)=O.*CC*.C1COCCN1.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1')
    r0_regioHash_mol

    # Create a list of SMILES
    rList = ['CC1=CC2=C(C=C1)C(=CN2CCN1CCOCC1)C(=O)C1=CC=CC2=C1C=CC=C2',\
            'CCCCCN1C=C(C2=CC=CC=C21)C(=O)C3=CC=CC4=CC=CC=C43',\
            'CC1COCCN1CCN1C=C(C(=O)C2=CC=CC3=C2C=CC=C3)C2=C1C=CC=C2',\
             'CC1=CC=C(C(=O)C2=CN(CCN3CCOCC3)C3=C2C=CC=C3)C2=C1C=CC=C2',\
            'CC1=C(CCN2CCOCC2)C2=C(C=CC=C2)N1C(=O)C1=CC=CC2=CC=CC=C12',\
            'CN1CCN(C(C1)CN2C=C(C3=CC=CC=C32)C(=O)C4=CC=CC5=CC=CC=C54)C']
    # Loop through the SMILES rList and create RDKit molecular objects
    rMols = [Chem.MolFromSmiles(r) for r in rList]

    # Calculate Regioisomer Hashes
    regioHashList = [rdMolHash.MolHash(rMol, rdkit.Chem.rdMolHash.HashFunction.Regioisomer) for rMol in rMols]
    print(regioHashList)

    ['*C.*C(*)=O.*CC*.C1COCCN1.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1',
     '*C(*)=O.*CCCCC.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1',
     '*C.*C(*)=O.*CC*.C1COCCN1.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1',
     '*C.*C(*)=O.*CC*.C1COCCN1.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1',
     '*C.*C(*)=O.*CC*.C1COCCN1.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1',
     '*C.*C.*C(*)=O.*C*.C1CNCCN1.c1ccc2[nH]ccc2c1.c1ccc2ccccc2c1']

    rmatches =[]
    for regioHash in regioHashList:
        if regioHash == r0_regioHash:
            print('Regioisomer: True')
            rmatches.append('Regioisomer: True')
        else:
            print('Regioisomer: False')
            rmatches.append('Regioisomer: False')

    Regioisomer: True
    Regioisomer: False
    Regioisomer: True
    Regioisomer: True
    Regioisomer: True
    Regioisomer: False

    # Create some labels
    index = ['r0: ','r1: ','r2: ','r3: ','r4: ','r5: ']
    labelList = [rmatches + index for rmatches,index in zip(index,rmatches)]
    # Display molecules with labels
    Draw.MolsToGridImage(rMols,legends=[label for label in labelList],
                        subImgSize=(250,250), useSVG=False)
    # note, that r0 is the initial molecule we were interested in.

### Contiguous Rotatable Bonds

Author: Paulo Tosco
Source: https://sourceforge.net/p/rdkit/mailman/message/36405144/
Index ID#: RDKitCB_22
Summary: Calculate the largest number of contiguous rotable bonds.

    from rdkit import Chem
    from rdkit.Chem.Lipinski import RotatableBondSmarts

    mol = Chem.MolFromSmiles('CCC(CC(C)CC1CCC1)C(CC(=O)O)N')
    mol

    def find_bond_groups(mol):
        """Find groups of contiguous rotatable bonds and return them sorted by decreasing size"""
        rot_atom_pairs = mol.GetSubstructMatches(RotatableBondSmarts)
        rot_bond_set = set([mol.GetBondBetweenAtoms(*ap).GetIdx() for ap in rot_atom_pairs])
        rot_bond_groups = []
        while (rot_bond_set):
            i = rot_bond_set.pop()
            connected_bond_set = set([i])
            stack = [i]
            while (stack):
                i = stack.pop()
                b = mol.GetBondWithIdx(i)
                bonds = []
                for a in (b.GetBeginAtom(), b.GetEndAtom()):
                    bonds.extend([b.GetIdx() for b in a.GetBonds() if (
                        (b.GetIdx() in rot_bond_set) and (not (b.GetIdx() in connected_bond_set)))])
                connected_bond_set.update(bonds)
                stack.extend(bonds)
            rot_bond_set.difference_update(connected_bond_set)
            rot_bond_groups.append(tuple(connected_bond_set))
        return tuple(sorted(rot_bond_groups, reverse = True, key = lambda x: len(x)))

    # Find groups of contiguous rotatable bonds in mol
    bond_groups = find_bond_groups(mol)
    # As bond groups are sorted by decreasing size, the size of the first group (if any)
    # is the largest number of contiguous rotatable bonds in mol
    largest_n_cont_rot_bonds = len(bond_groups[0]) if bond_groups else 0

    print(largest_n_cont_rot_bonds)

    8

    print(bond_groups)

    ((1, 2, 3, 5, 6, 10, 11, 12),)

    mol

## Writing Molecules

### Kekule SMILES

Author: Paulo Tosco
Source: https://sourceforge.net/p/rdkit/mailman/message/36893087/
Index ID#: RDKitCB_4
Summary: Kekulize a molecule and write Kekule SMILES

    from rdkit import Chem
    smi = "CN1C(NC2=NC=CC=C2)=CC=C1"
    mol = Chem.MolFromSmiles(smi)
    print(Chem.MolToSmiles(mol))

    Cn1cccc1Nc1ccccn1

    Chem.Kekulize(mol)
    print(Chem.MolToSmiles(mol, kekuleSmiles=True))

    CN1C=CC=C1NC1=NC=CC=C1

### Isomeric SMILES without isotopes

Author: Andrew Dalke
Source: https://sourceforge.net/p/rdkit/mailman/message/36877847/
Index ID#: RDKitCB_5
Summary: Write Isomeric SMILES without isotope information (i.e., only stereochemistry)

    from rdkit import Chem
    def MolWithoutIsotopesToSmiles(mol):
       atom_data = [(atom, atom.GetIsotope()) for atom in mol.GetAtoms()]
       for atom, isotope in atom_data:
       # restore original isotope values
           if isotope:
               atom.SetIsotope(0)
       smiles = Chem.MolToSmiles(mol)
       for atom, isotope in atom_data:
           if isotope:
              atom.SetIsotope(isotope)
       return smiles

    mol = Chem.MolFromSmiles("[19F][13C@H]([16OH])[35Cl]")
    print(MolWithoutIsotopesToSmiles(mol))

    O[C@@H](F)Cl

N.B. There are two limitations noted with this Isomeric SMILES without isotopes method including with isotopic hydrogens, and a requirement to recalculate stereochemistry. See the source discussion linked above for further explanation and examples.

## Reactions

### Reversing Reactions

Author: Greg Landrum
Source: https://gist.github.com/greglandrum/5ca4eebbe78f4d6d9b8cb03f401ad9cd and https://sourceforge.net/p/rdkit/mailman/message/36867857/
Index ID#: RDKitCB_6
Summary: Decompose a reaction product into its reactants
Reference Note: Example reaction from Hartenfeller [1]

    from rdkit import Chem
    from rdkit.Chem import AllChem
    from rdkit.Chem import Draw

    # Pictet-Spengler rxn
    rxn = AllChem.ReactionFromSmarts('[cH1:1]1:[c:2](-[CH2:7]-[CH2:8]-[NH2:9]):[c:3]:[c:4]:[c:5]:[c:6]:1.[#6:11]-[CH1;R0:10]=[OD1]>>[c:1]12:[c:2](-[CH2:7]-[CH2:8]-[NH1:9]-[C:10]-2(-[#6:11])):[c:3]:[c:4]:[c:5]:[c:6]:1')
    rxn

    rxn2 = AllChem.ChemicalReaction()
    for i in range(rxn.GetNumReactantTemplates()):
        rxn2.AddProductTemplate(rxn.GetReactantTemplate(i))
    for i in range(rxn.GetNumProductTemplates()):
        rxn2.AddReactantTemplate(rxn.GetProductTemplate(i))
    rxn2.Initialize()

    reacts = [Chem.MolFromSmiles(x) for x in ('NCCc1ccccc1','C1CC1C(=O)')]
    ps = rxn.RunReactants(reacts)
    ps0 = ps[0]
    for p in ps0:
        Chem.SanitizeMol(p)
    Draw.MolsToGridImage(ps0)

    reacts = ps0
    rps = rxn2.RunReactants(reacts)
    rps0 = rps[0]
    for rp in rps0:
        Chem.SanitizeMol(rp)
    Draw.MolsToGridImage(rps0)

N.B. This approach isn’t perfect and won’t work for every reaction. Reactions that include extensive query information in the original reactants are very likely to be problematic.

### Reaction Fingerprints and Similarity

Author: Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/37034626/
Index ID#: RDKitCB_27
Summary: Construct a reaction fingerprint and compute similarity
Reference Note: Reaction fingerprinting algorithm [6]

    from rdkit import Chem
    from rdkit.Chem import rdChemReactions
    from rdkit.Chem import DataStructs

    # construct the chemical reactions
    rxn1 = rdChemReactions.ReactionFromSmarts('CCCO>>CCC=O')
    rxn2 = rdChemReactions.ReactionFromSmarts('CC(O)C>>CC(=O)C')
    rxn3 = rdChemReactions.ReactionFromSmarts('NCCO>>NCC=O')

    # construct difference fingerprint (subtracts reactant fingerprint from product)
    fp1 = rdChemReactions.CreateDifferenceFingerprintForReaction(rxn1)
    fp2 = rdChemReactions.CreateDifferenceFingerprintForReaction(rxn2)
    fp3 = rdChemReactions.CreateDifferenceFingerprintForReaction(rxn3)

    print(DataStructs.TanimotoSimilarity(fp1,fp2))

    0.0

    # The similarity between fp1 and fp2 is zero because as far as the reaction
    # fingerprint is concerned, the parts which change within the reactions have
    # nothing in common with each other.
    # In contrast, fp1 and fp3 have some common parts
    print(DataStructs.TanimotoSimilarity(fp1,fp3))

    0.42857142857142855

## Error Messages

### Explicit Valence Error - Partial Sanitization

Author: Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/32599798/
Index ID#: RDKitCB_15
Summary: Create a mol object with skipping valence check, followed by a partial sanitization. N.B. Use caution, and make sure your molecules actually make sense before doing this!

    from rdkit import Chem
    # default RDKit behavior is to reject hypervalent P, so you need to set sanitize=False
    m = Chem.MolFromSmiles('F[P-](F)(F)(F)(F)F.CN(C)C(F)=[N+](C)C',sanitize=False)

    # next, you probably want to at least do a partial sanitization so that the molecule is actually useful:
    m.UpdatePropertyCache(strict=False)
    Chem.SanitizeMol(m,Chem.SanitizeFlags.SANITIZE_FINDRADICALS|Chem.SanitizeFlags.SANITIZE_KEKULIZE|Chem.SanitizeFlags.SANITIZE_SETAROMATICITY|Chem.SanitizeFlags.SANITIZE_SETCONJUGATION|Chem.SanitizeFlags.SANITIZE_SETHYBRIDIZATION|Chem.SanitizeFlags.SANITIZE_SYMMRINGS,catchErrors=True)

### Detect Chemistry Problems

Author: Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/36779572/
Index ID#: RDKitCB_14
Summary: Identify and capture error messages when creating mol objects.

    from rdkit import Chem
    m = Chem.MolFromSmiles('CN(C)(C)C', sanitize=False)
    problems = Chem.DetectChemistryProblems(m)
    print(len(problems))

    1

    print(problems[0].GetType())
    print(problems[0].GetAtomIdx())
    print(problems[0].Message())

    AtomValenceException
    1
    Explicit valence for atom # 1 N, 4, is greater than permitted

    m2 = Chem.MolFromSmiles('c1cncc1',sanitize=False)
    problems = Chem.DetectChemistryProblems(m2)
    print(len(problems))

    1

    print(problems[0].GetType())
    print(problems[0].GetAtomIndices())
    print(problems[0].Message())

    KekulizeException
    (0, 1, 2, 3, 4)
    Can't kekulize mol.  Unkekulized atoms: 0 1 2 3 4

## Miscellaneous Topics

### Explicit Valence and Number of Hydrogens

Author: Michael Palmer and Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/29679834/
Index ID#: RDKitCB_11
Summary: Calculate the explicit valence, number of explicit and implicit hydrogens, and total number of hydrogens on an atom. See the link for an important explanation about terminology and implementation of these methods in RDKit. Highlights are presented below.

Most of the time (exception is explained below), explicit refers to atoms that are in the graph and implicit refers to atoms that are not in the graph (i.e., Hydrogens). So given that the ring is aromatic (e.g.,in pyrrole), the explicit valence of each of the atoms (ignoring the Hs that are not present in the graph) in pyrrole is 3. If you want the Hydrogen count, use GetTotalNumHs(); the total number of Hs for each atom is one:

    from rdkit import Chem
    pyrrole = Chem.MolFromSmiles('C1=CNC=C1')
    for atom in pyrrole.GetAtoms():
       print(atom.GetSymbol(), atom.GetExplicitValence(), atom.GetTotalNumHs())

    C 3 1
    C 3 1
    N 3 1
    C 3 1
    C 3 1

In RDKit, there is overlapping nomenclature around the use of the words “explicit” and “implicit” when it comes to Hydrogens. When you specify the Hydrogens for an atom inside of square brackets in the SMILES, it becomes an “explicit” hydrogen as far as atom.GetNumExplicitHs() is concerned. Here is an example:

    pyrrole = Chem.MolFromSmiles('C1=CNC=C1')
    mol1 = Chem.MolFromSmiles('C1=CNCC1')
    mol2 = Chem.MolFromSmiles('C1=C[NH]CC1')

    for atom in pyrrole.GetAtoms():
        print(atom.GetSymbol(), atom.GetExplicitValence(), atom.GetNumImplicitHs(), atom.GetNumExplicitHs(), atom.GetTotalNumHs())

    C 3 1 0 1
    C 3 1 0 1
    N 3 0 1 1
    C 3 1 0 1
    C 3 1 0 1

    for atom in mol1.GetAtoms():
       print(atom.GetSymbol(), atom.GetExplicitValence(), atom.GetNumImplicitHs(), atom.GetNumExplicitHs(), atom.GetTotalNumHs())

    C 3 1 0 1
    C 3 1 0 1
    N 2 1 0 1
    C 2 2 0 2
    C 2 2 0 2

    for atom in mol2.GetAtoms():
       print(atom.GetSymbol(), atom.GetExplicitValence(), atom.GetNumImplicitHs(), atom.GetNumExplicitHs(), atom.GetTotalNumHs())

    C 3 1 0 1
    C 3 1 0 1
    N 3 0 1 1
    C 2 2 0 2
    C 2 2 0 2

### Wiener Index

Author: Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/36802142/
Index ID#: RDKitCB_12
Summary: Calculate the Wiener index (a topological index of a molecule)

    from rdkit import Chem
    def wiener_index(m):
        res = 0
        amat = Chem.GetDistanceMatrix(m)
        num_atoms = m.GetNumAtoms()
        for i in range(num_atoms):
            for j in range(i+1,num_atoms):
                res += amat[i][j]
        return res

    butane = Chem.MolFromSmiles('CCCC')
    print(wiener_index(butane))

    10.0

    isobutane = Chem.MolFromSmiles('CC(C)C')
    print(wiener_index(isobutane))

    9.0

### Organometallics with Dative Bonds

Author: Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/36727044/ and https://gist.github.com/greglandrum/6cd7aadcdedb1ebcafa9537e8a47e3a4
Index ID#: RDKitCB_19
Summary: Process organometallic SMILES by detecting single bonds between metals and replacing with dative bonds.

    from rdkit import Chem
    from rdkit.Chem.Draw import IPythonConsole

    def is_transition_metal(at):
        n = at.GetAtomicNum()
        return (n>=22 and n<=29) or (n>=40 and n<=47) or (n>=72 and n<=79)
    def set_dative_bonds(mol, fromAtoms=(7,8)):
        """ convert some bonds to dative

        Replaces some single bonds between metals and atoms with atomic numbers in fomAtoms
        with dative bonds. The replacement is only done if the atom has "too many" bonds.

        Returns the modified molecule.

        """
        pt = Chem.GetPeriodicTable()
        rwmol = Chem.RWMol(mol)
        rwmol.UpdatePropertyCache(strict=False)
        metals = [at for at in rwmol.GetAtoms() if is_transition_metal(at)]
        for metal in metals:
            for nbr in metal.GetNeighbors():
                if nbr.GetAtomicNum() in fromAtoms and \
                   nbr.GetExplicitValence()>pt.GetDefaultValence(nbr.GetAtomicNum()) and \
                   rwmol.GetBondBetweenAtoms(nbr.GetIdx(),metal.GetIdx()).GetBondType() == Chem.BondType.SINGLE:
                    rwmol.RemoveBond(nbr.GetIdx(),metal.GetIdx())
                    rwmol.AddBond(nbr.GetIdx(),metal.GetIdx(),Chem.BondType.DATIVE)
        return rwmol

    m = Chem.MolFromSmiles('CN(C)(C)[Pt]', sanitize=False)
    m2 = set_dative_bonds(m)
    m2

    # we can check the bond between nitrogen and platinum
    print(m2.GetBondBetweenAtoms(1,4).GetBondType())

    DATIVE

    # It also shows up in the output SMILES
    # This is an RDKit extension to SMILES
    print(Chem.MolToSmiles(m2))

    CN(C)(C)->[Pt]

### Enumerate SMILES

Author: Guillaume Godin and Greg Landrum
Source: https://sourceforge.net/p/rdkit/mailman/message/36591616/ and https://github.com/rdkit/UGM_2020/blob/master/Notebooks/Landrum_WhatsNew.ipynb
Index ID#: RDKitCB_24
Summary: Enumerate variations of SMILES strings for the same molecule.

    from rdkit import Chem

    # create a mol object
    mol = Chem.MolFromSmiles('CC(N)C1CC1')

    # Generate 100 random SMILES
    smis = []
    for i in range(100):
        smis.append(Chem.MolToSmiles(mol,doRandom=True,canonical=False))

    # remove duplicates
    smis_set = list(set(smis))
    print(smis_set) # output order will be random; doctest skipped

    ['NC(C)C1CC1',
     'C1(C(N)C)CC1',
     'C(N)(C)C1CC1',
     'CC(C1CC1)N',
     'C1C(C(N)C)C1',
     'C1C(C1)C(N)C',
     'C(C1CC1)(C)N',
     'C1(CC1)C(C)N',
     'C1C(C(C)N)C1',
     'C1CC1C(C)N',
     'C(C1CC1)(N)C',
     'C1(C(C)N)CC1',
     'C1C(C1)C(C)N',
     'C(C)(C1CC1)N',
     'C1CC1C(N)C',
     'C1(CC1)C(N)C',
     'C(N)(C1CC1)C',
     'NC(C1CC1)C',
     'CC(N)C1CC1',
     'C(C)(N)C1CC1']

    # If you need the multiple random SMILES strings to be reproducible,
    # the 2020.09 release has an option for this:
    m = Chem.MolFromSmiles('Oc1ncc(OC(CC)C)cc1')
    print(Chem.MolToRandomSmilesVect(m,5))  # output order random; doctest skipped

    ['c1c(cnc(O)c1)OC(CC)C', 'c1c(cnc(c1)O)OC(CC)C', 'c1cc(O)ncc1OC(CC)C', 'O(C(CC)C)c1ccc(nc1)O', 'O(C(C)CC)c1cnc(cc1)O']

    # by default the results are not reproducible:
    print(Chem.MolToRandomSmilesVect(m,5)) # output order random; doctest skipped

    ['c1nc(O)ccc1OC(CC)C', 'n1cc(OC(CC)C)ccc1O', 'c1c(OC(C)CC)ccc(O)n1', 'CCC(Oc1ccc(nc1)O)C', 'O(c1cnc(cc1)O)C(C)CC']

    # But we can provide a random number seed:
    m = Chem.MolFromSmiles('Oc1ncc(OC(CC)C)cc1')
    s1 = Chem.MolToRandomSmilesVect(m,5,randomSeed=0xf00d)
    print(s1)

    ['Oc1ccc(OC(CC)C)cn1', 'CC(CC)Oc1cnc(O)cc1', 'c1(O)ncc(cc1)OC(C)CC', 'c1cc(cnc1O)OC(CC)C', 'c1c(OC(CC)C)cnc(c1)O']

    s2 = Chem.MolToRandomSmilesVect(m,5,randomSeed=0xf00d)
    print(s2 == s1)

    True

### Reorder Atoms

Author: Jeffrey Van Santen and Paolo Tosco
Source: https://sourceforge.net/p/rdkit/mailman/message/37085522/ and https://gist.github.com/ptosco/36574d7f025a932bc1b8db221903a8d2
Index ID#: RDKitCB_28
Summary: Create a canonical order of atoms independent of input.

    from rdkit import Chem
    from rdkit.Chem.Draw import MolsToGridImage

    m = Chem.MolFromSmiles("c1([C@H](C)CC)cccc2ccccc12")
    m1 = Chem.MolFromSmiles("c12ccccc1c(ccc2)[C@H](C)CC")
    print(Chem.MolToSmiles(m) == Chem.MolToSmiles(m1))

    True

    # check if current canonical atom ordering matches
    m_neworder = tuple(zip(*sorted([(j, i) for i, j in enumerate(Chem.CanonicalRankAtoms(m))])))[1]
    m1_neworder = tuple(zip(*sorted([(j, i) for i, j in enumerate(Chem.CanonicalRankAtoms(m1))])))[1]
    print(m_neworder == m1_neworder)

    False

    # add atom numbers in images
    def addAtomIndices(mol):
        for i, a in enumerate(mol.GetAtoms()):
            a.SetAtomMapNum(i)

    addAtomIndices(m)
    addAtomIndices(m1)
    MolsToGridImage((m, m1))

    # renumber atoms with same canonical ordering
    m_renum = Chem.RenumberAtoms(m, m_neworder)
    m1_renum = Chem.RenumberAtoms(m1, m1_neworder)
    addAtomIndices(m_renum)
    addAtomIndices(m1_renum)
    MolsToGridImage((m_renum, m1_renum))

### Conformer Generation with ETKDG

Author: Shuzhe Wang
Source: Direct contribution to Cookbook
Index ID#: RDKitCB_25
Summary: Showcase various tricks for conformer generation with ETKDG

    from rdkit import Chem
    from rdkit.Chem import AllChem

To yield more chemically meaningful conformers, Riniker and Landrum implemented the experimental torsion knowledge distance geometry (ETKDG) method [3] which uses torsion angle preferences from the Cambridge Structural Database (CSD) to correct the conformers after distance geometry has been used to generate them. The configs of various conformer generation options are stored in a EmbedParameter object. To explicitly call the ETKDG EmbedParameter object:

    params = AllChem.ETKDG()

At the moment this is the default conformer generation routine in RDKit. A newer set of torsion angle potentials were published in 2016 [4], to use these instead:

    params = AllChem.ETKDGv2()

In 2020, we devised some improvements to the ETKDG method for sampling small rings and macrocycles [5].

    # this includes addtional small ring torsion potentials
    params = AllChem.srETKDGv3()

    # this includes additional macrocycle ring torsion potentials and macrocycle-specific handles
    params = AllChem.ETKDGv3()

    # to use the two in conjunction, do:
    params = AllChem.ETKDGv3()
    params.useSmallRingTorsions = True

    # a macrocycle attached to a small ring
    mol = Chem.MolFromSmiles("C(OC(CCCCCCC(OCCSC(CCCCCC1)=O)=O)OCCSC1=O)N1CCOCC1")
    mol = Chem.AddHs(mol)
    AllChem.EmbedMultipleConfs(mol, numConfs = 3 , params = params)

One additional tool we used in the paper is changing the bounds matrix of a molecule during distance geometry. The following code modifies the default molecular bounds matrix, with the idea of confining the conformational space of the molecule:

    from rdkit.Chem import rdDistGeom
    import rdkit.DistanceGeometry as DG

    mol = Chem.MolFromSmiles("C1CCC1C")
    mol = Chem.AddHs(mol)
    bm = rdDistGeom.GetMoleculeBoundsMatrix(mol)
    bm[0,3] = 1.21
    bm[3,0] = 1.20
    bm[2,3] = 1.21
    bm[3,2] = 1.20
    bm[4,3] = 1.21
    bm[3,4] = 1.20
    DG.DoTriangleSmoothing(bm)

    params.SetBoundsMat(bm)

Another tool we introduced is setting custom pairwise Coulombic interactions (CPCIs), which mimics additional electrostatic interactions between atom pairs to refine the embedded conformers. The setter takes in a dictionary of integer tuples as keys and reals as values. The following one-liner sets a repulsive (+ve) interaction of strength 0.9 e^2 between the atom indexed 0 and indexed 3, with the idea of keeping these two atoms further apart.

    params.SetCPCI({ (0,3) : 0.9 } )

To use the EmbedParameter for conformer generation:

    params.useRandomCoords = True
    # Note this is only an illustrative example, hydrogens are not added before conformer generation to keep the indices apparant
    AllChem.EmbedMultipleConfs(mol, numConfs = 3 , params = params)

Both of these setters can be used to help sampling all kinds of molecules as the users see fit. Nevertheless, to facilitate using them in conformer generation of macrocycles, we devised the python package github.com/rinikerlab/cpeptools to provide chemcially intuitive bound matrices and CPCIs for macrocycles. Example usage cases are shown in the README.

References

[1]

Hartenfeller, M.; Eberle, M.; Meier,P.; Nieto-Oberhuber, C.; Altmann, K.-H.; Schneider, G.; Jacoby, E.; and Renner, S. A Collection of Robust Organic Synthesis Reactions for In Silico Molecule Design. J. Chem Inf. Model. 2011, 51(12), 3093-3098. https://pubs.acs.org/doi/10.1021/ci200379p

[2]

O’Boyle, N. and Sayle, R. Making a hash of it: the advantage of selectively leaving out structural information. 259th ACS National Meeting Presentation, 2019, San Diego, CA. https://www.nextmovesoftware.com/talks/OBoyle_MolHash_ACS_201908.pdf

[3]

Riniker, S.; Landrum, G. A. Better Informed Distance Geometry: Using What We Know To Improve Conformation Generation. J. Chem. Inf. Model. 2015, 55(12), 2562-74. https://doi.org/10.1021/acs.jcim.5b00654

[4]

Guba, M.; Meyder, A.; Rarrey, M.; Hert, J. Torsion Library Reloaded: A New Version of Expert-Derived SMARTS Rules for Assessing Conformations of Small Molecules. J. Chem. Inf. Model. 2016, 56(1), 1-5. https://pubs.acs.org/doi/10.1021/acs.jcim.5b00522

[5]

Wang, S.; Witek, J.; Landrum, G. A.; Riniker, S. Improving Conformer Generation for Small Rings and Macrocycles Based on Distance Geometry and Experimental Torsional-Angle Preferences. J. Chem. Inf. Model. 2020, 60(4), 2044-2058. https://pubs.acs.org/doi/10.1021/acs.jcim.0c00025

[6]

Schneider, N.; Lowe, D.M.; Sayle, R.A.; Landrum, G. A. Development of a Novel Fingerprint for Chemical Reactions and Its Application to Large-Scale Reaction Classification and Similarity. J. Chem. Inf. Model. 2015, 55(1), 39-53. https://pubs.acs.org/doi/abs/10.1021/ci5006614

## License

This document is copyright (C) 2007-2020 by Greg Landrum and Vincent Scalfani.

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, 543 Howard Street, 5th Floor, San Francisco, California, 94105, USA.

The intent of this license is similar to that of the RDKit itself. In simple words: “Do whatever you want with it, but please give us some credit.”