In [1]:

# In Google Colab, install SMSD Pro first; elsewhere this cell does nothing.
import importlib.util, subprocess, sys
if importlib.util.find_spec("smsd") is None:
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "smsd==8.1.0", "rdkit<2027"], check=True)

# In Google Colab, install SMSD Pro first; elsewhere this cell does nothing. import importlib.util, subprocess, sys if importlib.util.find_spec("smsd") is None: subprocess.run([sys.executable, "-m", "pip", "install", "-q", "smsd==8.1.0", "rdkit<2027"], check=True)

Substructure and SMARTS search¶

This notebook uses SMSD 8.1.0 to search a small inline library of common drugs and benzodiazepines with molecule queries and SMARTS patterns. It covers a first search, hydrogen counts, functional-group flags, enumeration, chemistry options, limits, corpus decisions and drawing a hit. Every result keeps its status, so a search that stopped early is never read as an absence. pandas is optional: without it, the table is printed as plain rows.

Setup¶

Select a local Jupyter kernel with SMSD 8.1.0 or later installed (see Installation). The setup cell checks the version; it installs and downloads nothing, and no cell writes a file.

To use Colab, copy the setup cell of feature_search.ipynb to the top of this notebook and run it first. That cell installs an uploaded SMSD 8.1.0 core wheel; for a later wheel, change the version checks in the copy to match. Ignore its notice about NumPy and scikit-learn, which this notebook does not use.

In [2]:

# Uses the SMSD package of this kernel; nothing is installed or downloaded.
import re

import smsd

_release = tuple(int(part) for part in re.findall(r"\d+", smsd.__version__)[:3])
if _release < (8, 1, 0):
    raise RuntimeError("Use SMSD 8.1.0 or later in this notebook's kernel")
print("Using SMSD", smsd.__version__)

try:
    import pandas
except ImportError:
    pandas = None
    print("pandas is not installed; the table is printed as plain rows.")

# Uses the SMSD package of this kernel; nothing is installed or downloaded. import re import smsd _release = tuple(int(part) for part in re.findall(r"\d+", smsd.__version__)[:3]) if _release < (8, 1, 0): raise RuntimeError("Use SMSD 8.1.0 or later in this notebook's kernel") print("Using SMSD", smsd.__version__) try: import pandas except ImportError: pandas = None print("pandas is not installed; the table is printed as plain rows.")

Using SMSD 8.1.0

In [3]:

drugs = {
    "aspirin": "CC(=O)Oc1ccccc1C(=O)O",
    "salicylic acid": "OC(=O)c1ccccc1O",
    "paracetamol": "CC(=O)Nc1ccc(O)cc1",
    "ibuprofen": "CC(C)Cc1ccc(cc1)C(C)C(=O)O",
    "naproxen": "COc1ccc2cc(ccc2c1)C(C)C(=O)O",
    "caffeine": "CN1C=NC2=C1C(=O)N(C(=O)N2C)C",
    "phenol": "Oc1ccccc1",
}
benzodiazepines = {
    "diazepam": "CN1C(=O)CN=C(c2ccccc2)c2cc(Cl)ccc21",
    "nordazepam": "O=C1CN=C(c2ccccc2)c2cc(Cl)ccc2N1",
    "temazepam": "CN1C(=O)C(O)N=C(c2ccccc2)c2cc(Cl)ccc21",
    "oxazepam": "OC1N=C(c2ccccc2)c2cc(Cl)ccc2NC1=O",
    "lorazepam": "OC1N=C(c2ccccc2Cl)c2cc(Cl)ccc2NC1=O",
    "clonazepam": "O=C1CN=C(c2ccccc2Cl)c2cc([N+](=O)[O-])ccc2N1",
    "nitrazepam": "O=C1CN=C(c2ccccc2)c2cc([N+](=O)[O-])ccc2N1",
    "flunitrazepam": "CN1C(=O)CN=C(c2ccccc2F)c2cc([N+](=O)[O-])ccc21",
    "bromazepam": "O=C1CN=C(c2ccccn2)c2cc(Br)ccc2N1",
    "tetrazepam": "CN1C(=O)CN=C(C2=CCCCC2)c2cc(Cl)ccc21",
    "alprazolam": "Cc1nnc2n1-c1ccc(Cl)cc1C(c1ccccc1)=NC2",
    "triazolam": "Cc1nnc2n1-c1ccc(Cl)cc1C(c1ccccc1Cl)=NC2",
    "midazolam": "Cc1ncc2n1-c1ccc(Cl)cc1C(c1ccccc1F)=NC2",
    "clobazam": "CN1C(=O)CC(=O)N(c2ccccc2)c2cc(Cl)ccc21",
}
library = {**drugs, **benzodiazepines}
print(len(drugs), "common drugs and", len(benzodiazepines), "benzodiazepines")

drugs = { "aspirin": "CC(=O)Oc1ccccc1C(=O)O", "salicylic acid": "OC(=O)c1ccccc1O", "paracetamol": "CC(=O)Nc1ccc(O)cc1", "ibuprofen": "CC(C)Cc1ccc(cc1)C(C)C(=O)O", "naproxen": "COc1ccc2cc(ccc2c1)C(C)C(=O)O", "caffeine": "CN1C=NC2=C1C(=O)N(C(=O)N2C)C", "phenol": "Oc1ccccc1", } benzodiazepines = { "diazepam": "CN1C(=O)CN=C(c2ccccc2)c2cc(Cl)ccc21", "nordazepam": "O=C1CN=C(c2ccccc2)c2cc(Cl)ccc2N1", "temazepam": "CN1C(=O)C(O)N=C(c2ccccc2)c2cc(Cl)ccc21", "oxazepam": "OC1N=C(c2ccccc2)c2cc(Cl)ccc2NC1=O", "lorazepam": "OC1N=C(c2ccccc2Cl)c2cc(Cl)ccc2NC1=O", "clonazepam": "O=C1CN=C(c2ccccc2Cl)c2cc([N+](=O)[O-])ccc2N1", "nitrazepam": "O=C1CN=C(c2ccccc2)c2cc([N+](=O)[O-])ccc2N1", "flunitrazepam": "CN1C(=O)CN=C(c2ccccc2F)c2cc([N+](=O)[O-])ccc21", "bromazepam": "O=C1CN=C(c2ccccn2)c2cc(Br)ccc2N1", "tetrazepam": "CN1C(=O)CN=C(C2=CCCCC2)c2cc(Cl)ccc21", "alprazolam": "Cc1nnc2n1-c1ccc(Cl)cc1C(c1ccccc1)=NC2", "triazolam": "Cc1nnc2n1-c1ccc(Cl)cc1C(c1ccccc1Cl)=NC2", "midazolam": "Cc1ncc2n1-c1ccc(Cl)cc1C(c1ccccc1F)=NC2", "clobazam": "CN1C(=O)CC(=O)N(c2ccccc2)c2cc(Cl)ccc21", } library = {**drugs, **benzodiazepines} print(len(drugs), "common drugs and", len(benzodiazepines), "benzodiazepines")

7 common drugs and 14 benzodiazepines

In [4]:

# Phenol as a molecule query, searched in 4-tert-butylphenol
query = smsd.Substructure.find("c1ccccc1O")
result = query.match("CC(C)(C)c1ccc(O)cc1")
print("found:", result.found)
print("status:", result.status.name, "optimal:", result.optimal, "complete:", result.complete)
print("query atom -> target atom:", sorted(result.best.atom_mapping.items()))
if not (result.found and result.status == smsd.SearchStatus.COMPLETED):
    raise RuntimeError("Phenol should be found in 4-tert-butylphenol")

# Phenol as a molecule query, searched in 4-tert-butylphenol query = smsd.Substructure.find("c1ccccc1O") result = query.match("CC(C)(C)c1ccc(O)cc1") print("found:", result.found) print("status:", result.status.name, "optimal:", result.optimal, "complete:", result.complete) print("query atom -> target atom:", sorted(result.best.atom_mapping.items())) if not (result.found and result.status == smsd.SearchStatus.COMPLETED): raise RuntimeError("Phenol should be found in 4-tert-butylphenol")

found: True
status: COMPLETED optimal: True complete: True
query atom -> target atom: [(0, 6), (1, 5), (2, 4), (3, 10), (4, 9), (5, 7), (6, 8)]

Hydrogen counts¶

A molecule query does not compare hydrogen counts, so phenol as a molecule query also matches the ester oxygen of aspirin. A SMARTS H primitive, as in c[OX2H], fixes the count. to_smarts(..., include_h_count=True) fixes the count on every atom, ring positions included, so its pattern matches unsubstituted phenol only.

In [5]:

aspirin = drugs["aspirin"]
print("phenol molecule query in aspirin:", smsd.Substructure.find("c1ccccc1O").match(aspirin).found)
print("c[OX2H] in aspirin:", smsd.Substructure.smarts("c[OX2H]").match(aspirin).found)

pattern = smsd.to_smarts("c1ccccc1O", include_h_count=True)
print(pattern)
fixed = smsd.Substructure.smarts(pattern)
for name in ("aspirin", "phenol", "paracetamol"):
    print(f"  {name}: {fixed.match(drugs[name]).found}")

aspirin = drugs["aspirin"] print("phenol molecule query in aspirin:", smsd.Substructure.find("c1ccccc1O").match(aspirin).found) print("c[OX2H] in aspirin:", smsd.Substructure.smarts("c[OX2H]").match(aspirin).found) pattern = smsd.to_smarts("c1ccccc1O", include_h_count=True) print(pattern) fixed = smsd.Substructure.smarts(pattern) for name in ("aspirin", "phenol", "paracetamol"): print(f" {name}: {fixed.match(drugs[name]).found}")

phenol molecule query in aspirin: True
c[OX2H] in aspirin: False
[#6;a;H0]:1(:[#6;a;H1]:[#6;a;H1]:[#6;a;H1]:[#6;a;H1]:[#6;a;H1]:1)-[#8;A;H1]
  aspirin: False
  phenol: True
  paracetamol: False

In [6]:

# Flag four functional groups across the common drugs with one prepared corpus
groups = {
    "carboxylic acid": "C(=O)[OX2H1]",
    "phenol": "c[OX2H]",
    "aryl ester": "cOC(=O)[#6]",
    "amide": "[#6](=O)[#7X3]",
}
corpus = smsd.TargetCorpus.from_smiles(list(drugs.values()))
table = {"name": list(drugs)}
for label, pattern in groups.items():
    results = smsd.Substructure.smarts(pattern).match_many(corpus, num_threads=2)
    if any(result.status != smsd.SearchStatus.COMPLETED for result in results):
        raise RuntimeError(f"The {label} search stopped early, so a False is not proved")
    table[label] = [result.found for result in results]

if pandas is not None:
    print(pandas.DataFrame(table).to_string(index=False))
else:
    print(tuple(table))
    for row in zip(*table.values()):
        print(row)

# Flag four functional groups across the common drugs with one prepared corpus groups = { "carboxylic acid": "C(=O)[OX2H1]", "phenol": "c[OX2H]", "aryl ester": "cOC(=O)[#6]", "amide": "[#6](=O)[#7X3]", } corpus = smsd.TargetCorpus.from_smiles(list(drugs.values())) table = {"name": list(drugs)} for label, pattern in groups.items(): results = smsd.Substructure.smarts(pattern).match_many(corpus, num_threads=2) if any(result.status != smsd.SearchStatus.COMPLETED for result in results): raise RuntimeError(f"The {label} search stopped early, so a False is not proved") table[label] = [result.found for result in results] if pandas is not None: print(pandas.DataFrame(table).to_string(index=False)) else: print(tuple(table)) for row in zip(*table.values()): print(row)

          name  carboxylic acid  phenol  aryl ester  amide
       aspirin             True   False        True  False
salicylic acid             True    True       False  False
   paracetamol            False    True       False   True
     ibuprofen             True   False       False  False
      naproxen             True   False       False  False
      caffeine            False   False       False   True
        phenol            False    True       False  False

In [7]:

# Every embedding of benzene in toluene, with a limit and a time limit
benzene = smsd.Substructure.find("c1ccccc1").with_timeout(1000)
every = benzene.match_all("Cc1ccccc1", limit=100)
occurrences = {frozenset(match.atom_mapping.values()) for match in every.matches}
print("embeddings:", len(every.matches), every.status.name, "complete:", every.complete)
print("distinct ring occurrences:", len(occurrences))

capped = benzene.match_all("Cc1ccccc1", limit=5)
print("limit 5:", len(capped.matches), capped.status.name,
      "optimal:", capped.optimal, "complete:", capped.complete)

# Every embedding of benzene in toluene, with a limit and a time limit benzene = smsd.Substructure.find("c1ccccc1").with_timeout(1000) every = benzene.match_all("Cc1ccccc1", limit=100) occurrences = {frozenset(match.atom_mapping.values()) for match in every.matches} print("embeddings:", len(every.matches), every.status.name, "complete:", every.complete) print("distinct ring occurrences:", len(occurrences)) capped = benzene.match_all("Cc1ccccc1", limit=5) print("limit 5:", len(capped.matches), capped.status.name, "optimal:", capped.optimal, "complete:", capped.complete)

embeddings: 12 COMPLETED complete: True
distinct ring occurrences: 1
limit 5: 5 LIMIT_REACHED optimal: True complete: False

Chemistry options¶

options.chemistry decides which atoms and bonds may match. The default FLEXIBLE aromaticity lets aromatic atoms and bonds of the query match non-aromatic ones, so the benzodiazepine core below also finds tetrazepam, whose phenyl ring is replaced by a cyclohexenyl ring. STRICT aromaticity excludes it. The last two rows show formal-charge and ring-membership matching on small molecules. Options objects take no keyword arguments: set each field as an attribute.

In [8]:

core = "O=C1CN=C(c2ccccc2)c2ccccc2N1"   # 5-phenyl-1,4-benzodiazepin-2-one
strict = smsd.Substructure.Options()
strict.chemistry.aromaticity_mode = smsd.AromaticityMode.STRICT
names, targets = list(benzodiazepines), list(benzodiazepines.values())
hits = {}
for label, search in (("default", smsd.Substructure.find(core)),
                      ("STRICT aromaticity", smsd.Substructure.find(core).with_options(strict))):
    results = search.match_many(targets, num_threads=2)
    if any(result.status != smsd.SearchStatus.COMPLETED for result in results):
        raise RuntimeError("A core search stopped early")
    hits[label] = [name for name, result in zip(names, results) if result.found]
    print(f"{label}: {len(hits[label])} hits")
print("found only with the default:", sorted(set(hits["default"]) - set(hits["STRICT aromaticity"])))


def found(query, target, **fields):
    options = smsd.Substructure.Options()
    for name, value in fields.items():
        setattr(options.chemistry, name, value)
    return smsd.Substructure.find(query).with_options(options).match(target).found


print("C[NH3+] in methylamine, match_formal_charge:",
      found("C[NH3+]", "CN"), "->", found("C[NH3+]", "CN", match_formal_charge=True))
print("propane in cyclohexane, ring_matches_ring_only:",
      found("CCC", "C1CCCCC1"), "->", found("CCC", "C1CCCCC1", ring_matches_ring_only=True))

core = "O=C1CN=C(c2ccccc2)c2ccccc2N1" # 5-phenyl-1,4-benzodiazepin-2-one strict = smsd.Substructure.Options() strict.chemistry.aromaticity_mode = smsd.AromaticityMode.STRICT names, targets = list(benzodiazepines), list(benzodiazepines.values()) hits = {} for label, search in (("default", smsd.Substructure.find(core)), ("STRICT aromaticity", smsd.Substructure.find(core).with_options(strict))): results = search.match_many(targets, num_threads=2) if any(result.status != smsd.SearchStatus.COMPLETED for result in results): raise RuntimeError("A core search stopped early") hits[label] = [name for name, result in zip(names, results) if result.found] print(f"{label}: {len(hits[label])} hits") print("found only with the default:", sorted(set(hits["default"]) - set(hits["STRICT aromaticity"]))) def found(query, target, **fields): options = smsd.Substructure.Options() for name, value in fields.items(): setattr(options.chemistry, name, value) return smsd.Substructure.find(query).with_options(options).match(target).found print("C[NH3+] in methylamine, match_formal_charge:", found("C[NH3+]", "CN"), "->", found("C[NH3+]", "CN", match_formal_charge=True)) print("propane in cyclohexane, ring_matches_ring_only:", found("CCC", "C1CCCCC1"), "->", found("CCC", "C1CCCCC1", ring_matches_ring_only=True))

default: 9 hits
STRICT aromaticity: 8 hits
found only with the default: ['tetrazepam']
C[NH3+] in methylamine, match_formal_charge: True -> False
propane in cyclohexane, ring_matches_ring_only: True -> False

In [9]:

# A stopped search with no matches proves nothing: optimal and complete are False
query = smsd.Substructure.find("c1ccccc1O")
target = "CC(C)(C)c1ccc(O)cc1"


def report(label, result):
    print(f"{label:14}", result.status.name, "found:", result.found,
          "optimal:", result.optimal, "complete:", result.complete)


report("timeout 0", query.with_timeout(0).match(target))

options = smsd.Substructure.Options()
options.node_limit = 3
report("node limit 3", query.with_options(options).match(target))

token = smsd.CancellationToken()
options = smsd.Substructure.Options()
options.cancellation = token
options.timeout_ms = 1000
token.cancel()
report("cancelled", query.with_options(options).match(target))

# A stopped search with no matches proves nothing: optimal and complete are False query = smsd.Substructure.find("c1ccccc1O") target = "CC(C)(C)c1ccc(O)cc1" def report(label, result): print(f"{label:14}", result.status.name, "found:", result.found, "optimal:", result.optimal, "complete:", result.complete) report("timeout 0", query.with_timeout(0).match(target)) options = smsd.Substructure.Options() options.node_limit = 3 report("node limit 3", query.with_options(options).match(target)) token = smsd.CancellationToken() options = smsd.Substructure.Options() options.cancellation = token options.timeout_ms = 1000 token.cancel() report("cancelled", query.with_options(options).match(target))

timeout 0      TIMED_OUT found: False optimal: False complete: False
node limit 3   LIMIT_REACHED found: False optimal: False complete: False
cancelled      CANCELLED found: False optimal: False complete: False

In [10]:

# One decision per slot: targets that cannot contain the core are ruled out
# before matching, with a proof; the rest are searched as usual.
library_corpus = smsd.TargetCorpus.from_smiles(list(library.values()))
decisions = smsd.Substructure.find(core).match_corpus(library_corpus, num_threads=2)
counts = {}
for name, decision in zip(library, decisions):
    if decision.result is None:
        outcome = "ruled out before matching"
    elif decision.threshold_reached:
        outcome = "found"
    elif decision.proven_below_threshold:
        outcome = "searched, proved absent"
    else:
        outcome = "unknown (search stopped)"
    counts[outcome] = counts.get(outcome, 0) + 1
    print(f"{name:15} {outcome}")
print(counts)

# One decision per slot: targets that cannot contain the core are ruled out # before matching, with a proof; the rest are searched as usual. library_corpus = smsd.TargetCorpus.from_smiles(list(library.values())) decisions = smsd.Substructure.find(core).match_corpus(library_corpus, num_threads=2) counts = {} for name, decision in zip(library, decisions): if decision.result is None: outcome = "ruled out before matching" elif decision.threshold_reached: outcome = "found" elif decision.proven_below_threshold: outcome = "searched, proved absent" else: outcome = "unknown (search stopped)" counts[outcome] = counts.get(outcome, 0) + 1 print(f"{name:15} {outcome}") print(counts)

aspirin         ruled out before matching
salicylic acid  ruled out before matching
paracetamol     ruled out before matching
ibuprofen       ruled out before matching
naproxen        ruled out before matching
caffeine        ruled out before matching
phenol          ruled out before matching
diazepam        found
nordazepam      found
temazepam       found
oxazepam        found
lorazepam       found
clonazepam      found
nitrazepam      found
flunitrazepam   found
bromazepam      searched, proved absent
tetrazepam      found
alprazolam      ruled out before matching
triazolam       ruled out before matching
midazolam       ruled out before matching
clobazam        ruled out before matching
{'ruled out before matching': 11, 'found': 9, 'searched, proved absent': 1}

In [11]:

# Draw paracetamol with its aromatic hydroxyl highlighted. depict_mapping takes
# atoms of the drawn molecule as keys, so the search mapping is inverted.
paracetamol = drugs["paracetamol"]
hit = smsd.Substructure.smarts("c[OX2H]").match(paracetamol)
highlight = {target: query for query, target in hit.best.atom_mapping.items()}
print("highlighted target atoms:", sorted(highlight.items()))
svg = smsd.depict_mapping(paracetamol, highlight)

try:
    from IPython.display import SVG, display
except ImportError:
    print("IPython is not available; SVG characters:", len(svg))
else:
    display(SVG(svg))

# Draw paracetamol with its aromatic hydroxyl highlighted. depict_mapping takes # atoms of the drawn molecule as keys, so the search mapping is inverted. paracetamol = drugs["paracetamol"] hit = smsd.Substructure.smarts("c[OX2H]").match(paracetamol) highlight = {target: query for query, target in hit.best.atom_mapping.items()} print("highlighted target atoms:", sorted(highlight.items())) svg = smsd.depict_mapping(paracetamol, highlight) try: from IPython.display import SVG, display except ImportError: print("IPython is not available; SVG characters:", len(svg)) else: display(SVG(svg))

highlighted target atoms: [(7, 0), (8, 1)]

Summary¶

Molecule queries and SMARTS patterns share one prepared-query interface, and every result states whether the search completed. Check the status before reading a missing match as an absence, give enumerations a limit and a time limit, and use SMARTS H primitives when hydrogen counts matter.

For the full reference, see the Substructure search and SMARTS queries chapters of the user guide.