How to search substructures with explicit hydrogen in RDKit?

Question

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]

Answer 1

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.