# How to search substructures with explicit hydrogen in RDKit?

For a substructure search I would like to search for structures containing unfused benzyl. The idea was to explicitly add hydrogen. But apparently this does not give the expected results.

So, I assume either my understanding or expectations are wrong or I'm using RDKit not properly. I just learnt (Substructure search with RDKit) that I better use SMARTS expressions with Chem.MolToSmarts() where I get the expected results.

However, if I have much more complex search structures, how can I get from the SMILES to the SMARTS? Do I have always to figure out the SMARTS code myself? Is there maybe a way just drawing a molecule and add explicit hydrogen and search with this somehow? I guess it should be pretty clear (at least to a human) what I want to find (or not find) with structure 2. Are there maybe any tools or code which can assist here?

Code:

from rdkit import Chem

smiles_strings = '''CC1=CC=CC=C1
CC1=C([H])C([H])=C([H])C([H])=C1[H]
C1(C2=CC=CC=C2)=CC=CC=C1
C12=C(CC3=CC=CC=C3C2)C=CC=C1
C12=CC=CC=C1C=C3C(C=CC=C3)=C2
C12=C(C(C(C=C3)=CC=C4)=C4C=C2)C3=CC=C1
'''
smiles_list = smiles_strings.splitlines()

def search_structure(pattern):
found = []
for idx,smiles in enumerate(smiles_list):
m = Chem.MolFromSmiles(smiles)
if m.HasSubstructMatch(pattern):
found.append(idx+1)
print("Structures found: {}".format(found))

smiles_1  = smiles_list[0]
pattern_1 = Chem.MolFromSmiles(smiles_1)
smiles_2a  = smiles_list[1]
pattern_2a = Chem.MolFromSmiles(smiles_2a)
smiles_2b  = Chem.MolToSmiles(pattern_2a)
pattern_2b = Chem.MolFromSmiles(smiles_2b)
smiles_2c = smiles_2b
pattern_2c = Chem.MolFromSmarts(smiles_2c)
smarts_3 = '[cR1]1[cR1][cR1][cR1][cR1][cR1]1-[c,C]'
pattern_3 = Chem.MolFromSmarts(smarts_3)

print("\nSMILES 1 : {}\n# my expectation: [1, 2, 3, 4, 5, 6]".format(smiles_1))
search_structure(pattern_1)
print("\nSMILES 2a: {}\n# my expectation: [1, 2, 3]".format(smiles_2a))
search_structure(pattern_2a)
print("\nSMILES 2b: {} (Mol from SMILES)\n# my expectation: [1, 2, 3]".format(smiles_2b))
search_structure(pattern_2b)
print("\nSMILES 2c: {} (Mol from SMARTS)\n# my expectation: [1, 2, 3]".format(smiles_2c))
search_structure(pattern_2c)
print("\nSMARTS 3 : {}\n# my expectation: [1, 2, 3]".format(smarts_3))
search_structure(pattern_3)


Result:

SMILES 1 : CC1=CC=CC=C1
# my expectation: [1, 2, 3, 4, 5, 6]
Structures found: [1, 2, 3, 4]

SMILES 2a: CC1=C([H])C([H])=C([H])C([H])=C1[H]
# my expectation: [1, 2, 3]
Structures found: [1, 2, 3, 4]

SMILES 2b: Cc1ccccc1 (Mol from SMILES)
# my expectation: [1, 2, 3]
Structures found: [1, 2, 3, 4]

SMILES 2c: Cc1ccccc1 (Mol from SMARTS)
# my expectation: [1, 2, 3]
Structures found: [1, 2, 4]

SMARTS 3 : [cR1]1[cR1][cR1][cR1][cR1][cR1]1-[c,C]
# my expectation: [1, 2, 3]
Structures found: [1, 2, 3]


If you work in a Jupyter Notebook you can visualize your substructure patterns and the search results, so you can see how the patterns work.

from rdkit import Chem
from rdkit.Chem import Draw

smiles = ['CC1=CC=CC=C1','CC1=C([H])C([H])=C([H])C([H])=C1[H]',
'C1(C2=CC=CC=C2)=CC=CC=C1','C12=C(CC3=CC=CC=C3C2)C=CC=C1',
'C12=CC=CC=C1C=C3C(C=CC=C3)=C2','C12=C(C(C(C=C3)=CC=C4)=C4C=C2)C3=CC=C1']

params = Chem.SmilesParserParams()
params.removeHs=False # draw and work with explicit Hs
mols = [Chem.MolFromSmiles(s, params) for s in smiles]
Draw.MolsToGridImage(mols, molsPerRow=3)


This will give you the compounds and especially the one with the explicit Hs.

Now we get the patterns and especially the one with the explicit Hs.

patt1 = Chem.MolFromSmiles('CC1=CC=CC=C1')
patt2 = Chem.MolFromSmiles('CC1=C([H])C([H])=C([H])C([H])=C1[H]', params) # search with explicit Hs
patt3 = Chem.MolFromSmiles('Cc1ccccc1')
patt4 = Chem.MolFromSmarts('Cc1ccccc1')
patt5 = Chem.MolFromSmarts('[cR1]1[cR1][cR1][cR1][cR1][cR1]1-[c,C]')

patts = [patt1,patt2,patt3,patt4,patt5]
leg = ['CC1=CC=CC=C1','CC1=C([H])C([H])=C([H])C([H])=C1[H]','Cc1ccccc1 smiles',
'Cc1ccccc1 smarts','[cR1]1[cR1][cR1][cR1][cR1][cR1]1-[c,C]']
Draw.MolsToGridImage(patts, molsPerRow=3, legends=leg)


Search and highlight the result.

ms = []
patt = []
allsubs = []
for m in range(len(mols)):
for p in range(len(patts)):
sub = mols[m].GetSubstructMatches(patts[p])
if len(sub) > 0:
ms.append(mols[m])
patt.append(leg[p])
allsubs.append(sub[0]) # substructures could be find multiple time - just take the first

Draw.MolsToGridImage(ms, molsPerRow=3, legends=patt, highlightAtomLists=allsubs)


Compounds 5 and 6 have no nonaromatic bonds, so there is nothing to find.

• thank you very much for your code and explanations. It looks like params.removeHs=False is the bit I was missing. I can only test next week and will let you know. – theozh Mar 7 '20 at 5:21
• If I'm using params.removeHs=False for searching with structure 2, then HasSubstructMatch(pattern) only finds structure 2 (as in your example). Why are structures 1 and 3 not found as well? Maybe because of the explicit non-aromatic extra carbon in search structure 2? What do I have to change that by drawing structure 2 with explicit hydrogen and unspecified (is it or not?) extra carbon will also find structures 1 and 3? – theozh Mar 10 '20 at 6:45
• @theozh It is not the extra carbon, it is just that 1 and 3 do not have explicit Hs. If I add explicit Hs to 3 it is found. I thought C[c;H,h]1[c;H,h][c;H,h][c;H,h][c;H,h][c;H,h]1 (see smartsview.zbh.uni-hamburg.de) will work, but it doesn't. – rapelpy Mar 10 '20 at 18:58
• thank you, yes, today I also found out that if you add Hs to all the structures you're searching in, via mol = Chem.AddHs(mol), you will find what is expected. However, there are other cases which I still do not understand. I need more testing. – theozh Mar 10 '20 at 19:16
• It took me a while, but you can find 1, 2 and 3 with *[c]1[c;H,h][c;H,h][c;H,h][c;H,h][c;H,h]1 and 1,2,3,4,5 with *[c]1[c;H,h,*][c;H,h][c;H,h][c;H,h][c;H,h]1 – rapelpy Mar 12 '20 at 19:48