Reactions#

A chemical reaction redistributes the valence electrons that hold a molecule together — some bonds break, others form, charges appear and vanish. Crucially a reaction conserves mass: no atom is created or destroyed, so there is always a one-to-one atom map between the educts and the products. A MolFoundry rule is exactly that map together with the changes layered on top of it, written in the double-pushout style as a span

\[ L \xleftarrow{\;} K \xrightarrow{\;} R \]

left (L) holds the features that are present before but not after — the bonds and labels that vanish;
right (R) holds the features present after but not before — the ones that emerge;
context (K) is everything untouched, shared by both sides.

Atoms are tied across the three blocks by their id. The atoms and bonds that actually change form the reaction center; everything else is context. How much context you spell out sets how specific the rule is: with no context a rule is general and fires on many molecules, while spelling out the entire educt makes a rule that matches exactly one.

Writing a rule in GML#

The recipe for turning a balanced reaction into a rule is mechanical:

Extract the reaction center — the atoms and bonds involved in bond making/breaking.
Number the corresponding atoms on both sides (that is the atom map).
List what vanishes under left, what emerges under right, and what is unchanged under context.

Take the heterolytic ionization of a C–O single bond into a carbocation and an anion. Only the bond and the two neutral labels vanish; two charged labels emerge; nothing is unchanged, so context is empty:

ionization = ruleGMLString('''rule [
    ruleID "ionization reaction"
    left  [ node [ id 1 label "C" ] node [ id 2 label "O" ]
            edge [ source 1 target 2 label "-" ] ]
    context [ ]
    right [ node [ id 1 label "C+" ] node [ id 2 label "O-" ] ]
]''', "ionization", add=False)
ionization.print(rulep)

ionization

The depiction shows the three graphs L, K, R side by side; DrawMode.Circles draws atoms as labeled circles so the shared ids are easy to follow. A rule where only a bond order changes keeps both atoms in context and moves just the edge — for the keto→enol tautomerization the two carbons, the oxygen and the migrating hydrogen all persist:

keto_enol = ruleGMLString('''rule [
    ruleID "keto to enol"
    left  [ edge [ source 1 target 4 label "-" ] edge [ source 1 target 2 label "-" ]
            edge [ source 2 target 3 label "=" ] ]
    context [ node [ id 1 label "C" ] node [ id 2 label "C" ]
              node [ id 3 label "O" ] node [ id 4 label "H" ] ]
    right [ edge [ source 1 target 2 label "=" ] edge [ source 2 target 3 label "-" ]
            edge [ source 3 target 4 label "-" ] ]
]''', "keto to enol", add=False)
keto_enol.print(rulep)

keto to enol

makeInverse() flips a rule — the reverse reaction for free — and ruleGML() loads a rule from a file (with an optional invert=True).

The DFS line notation#

Spelling out every node and edge is precise but verbose. GraphDFS (and RuleDFS for rules) is a compact, SMILES-like alternative: a depth-first walk of the graph in which

a bracketed […] vertex is an arbitrary label (so [C] is a bare carbon with no implicit hydrogens, and [CoA] is perfectly legal), while a bare organic atom C, O, … gets implicit hydrogens exactly as in SMILES;
{…} is an arbitrary edge label, the shorthands are - : = #, an omitted edge is a single bond, and . is a non-edge (start a new, disconnected fragment);
a number written right after a vertex gives it an id; the same number reappearing after an edge is a ring bond back to that vertex, and across a rule’s >> the ids relate the two sides.

The keto→enol rule above is a one-liner, and it parses to the very same rule:

keto_enol_dfs = ruleDFS("[H]4[C]1[C]2=[O]3>>[C]1=[C]2[O]3[H]4", add=False)
print("same rule as the GML above:", keto_enol_dfs.canonKey == keto_enol.canonKey)
graphDFS("CC(=O)C", add=False).print()   # a molecule, hydrogens filled like SMILES

same rule as the GML above: True

g_{0}

Molecules load with graphDFS() and rules with ruleDFS() (pass labelType=LabelType.Term when the labels are term variables — see below).

04_rules_in_gml.py#

"""4. Writing graph-rewrite rules (GML and DFS).

A reaction rewrites a molecule graph. A rule is a span left <- context -> right:

* ``left``    -- features present in the educts but gone in the products
                 (bonds/labels that vanish);
* ``right``   -- features present in the products but not the educts
                 (bonds/labels that emerge);
* ``context`` -- everything untouched, shared by both sides.

Atoms are related across the three sides by their ``id`` (the atom map). The
"reaction center" is the atoms/bonds that actually change; how much unchanged
context you add controls how specific the rule is.

GML spells all of this out explicitly. There is also a compact **DFS line
notation** (``ruleDFS``, and ``graphDFS`` for molecules): a SMILES-like
depth-first walk where a bracketed ``[..]`` vertex is an arbitrary label, a bare
organic atom gets implicit hydrogens, the trailing id numbers double as ring
bonds and the left<->right atom map, and ``>>`` separates a rule's two sides.

Three rules are loaded from files here (see ``gml/``):
    ionization.gml     C-O single bond  ->  C(+) + O(-)   (empty context)
    keto_enol.gml      keto  ->  enol tautomerization
    aldol_addition.gml enol + carbonyl  ->  aldol adduct

Run:  python 04_rules_in_gml.py
"""
import _bootstrap
from _bootstrap import GML
from molfoundry import ruleGML, ruleDFS, graphDFS, GraphPrinter, DrawMode, enable_post

enable_post()

ionization = ruleGML(str(GML / "ionization.gml"))
keto_enol = ruleGML(str(GML / "keto_enol.gml"))
aldol = ruleGML(str(GML / "aldol_addition.gml"))

# A GraphPrinter tuned for rules: draw atoms as labeled circles so the L/K/R
# atom map (the ids) is easy to follow.
p = GraphPrinter()
p.drawMode = DrawMode.Circles

for r in (ionization, keto_enol, aldol):
    r.print(p)
    print(f"{r.name!r:26s} left={r.numLeftComponents} right={r.numRightComponents}")

# makeInverse() flips a rule -- the reverse reaction, for free (left <-> right).
# Equivalently, ruleGML(path, invert=True) loads a file already inverted.
enol_keto = keto_enol.makeInverse()
enol_keto.print(p)
print(f"\nInverted rule: {enol_keto.name!r}")

# The same keto-enol rule, written compactly in DFS line notation instead of GML
# -- it parses to exactly the same rule (equal canonical keys).
keto_enol_dfs = ruleDFS("[H]4[C]1[C]2=[O]3>>[C]1=[C]2[O]3[H]4", "keto to enol (DFS)")
same = keto_enol_dfs.canonKey == keto_enol.canonKey
print(f"DFS '[H]4[C]1[C]2=[O]3>>[C]1=[C]2[O]3[H]4' == keto_enol.gml: {same}")
# ...and a molecule in GraphDFS (bare atoms get implicit hydrogens, as in SMILES):
print("graphDFS('CC(=O)C') ->", graphDFS("CC(=O)C", add=False).smiles)

print("\n--- ionization.gml as MolFoundry sees it ---")
print(ionization.getGMLString())

Terms: label variables#

Wildcards are everywhere in chemistry — a rule might act on any aldehyde, or any leaving group. Setting labelType "term" turns vertex and edge labels into first-order terms and evaluates them with unification (a Prolog-like logic): a label beginning with _, or a bare *, is a variable that unifies with whatever sits in that position, and a constrainLabelAny block restricts a tuple of variables to an explicit list of allowed value tuples.

The classic example reduces an aldehyde (or thio-aldehyde) with H₂ but leaves ketones alone. The carbonyl partner _X is a variable — one rule covers both C=O and C=S — and the selectivity is a constraint on the carbon’s three neighbors: Ald(_X, _Y, _Z) must unify with Ald(O, H, C) or Ald(S, H, C), i.e. an oxo/thio carbonyl with one hydrogen and one carbon. A ketone would be Ald(O, C, C), which is not in the list, and never matches.

reduce_aldehyde = ruleGMLString('''rule [
    ruleID "reduce aldehyde, not ketone"
    labelType "term"
    left  [ edge [ source 1 target 2 label "=" ] edge [ source 5 target 6 label "-" ] ]
    context [ node [ id 1 label "C" ]  node [ id 2 label "_X" ]
              node [ id 3 label "_Y" ] node [ id 4 label "_Z" ]
              node [ id 5 label "H" ]  node [ id 6 label "H" ]
              edge [ source 1 target 3 label "-" ] edge [ source 1 target 4 label "-" ] ]
    right [ edge [ source 1 target 2 label "-" ] edge [ source 1 target 6 label "-" ]
            edge [ source 2 target 5 label "-" ] ]
    constrainLabelAny [ label "Ald(_X, _Y, _Z)"
                        labels [ label "Ald(O, H, C)" label "Ald(S, H, C)" ] ]
]''', add=False)
reduce_aldehyde.print(rulep)

reduce aldehyde, not ketone

Term rewriting needs the Term / Unification label setting on the derivation graph — LabelSettings(LabelType.Term, LabelRelation.Unification) — and then the one rule reduces propanal, benzaldehyde and thioacetaldehyde but leaves acetone untouched (see the full script).

A constraint can equally pin bond labels. Benzene has no single Lewis structure — its ring is the average of two Kekulé forms — so a rule acting on a benzene ring must match either. Constraining the six ring bonds Tauto(_X, …, _W) to the two alternations (-,=,-,=,-,=) and (=,-,=,-,=,-) does exactly that: on two naphthalene-dialdehydes the rule reduces both aldehydes of the one where both fused rings match a Kekulé pattern, and only one of the other.

05_terms.py#

"""5. Terms: label variables for polymorphic rules.

Setting ``labelType "term"`` turns vertex/edge labels into first-order terms: a
label beginning with ``_`` is a *variable* that unifies with whatever atom (or
bond) sits in that position. One rule can then cover a whole family of
substrates -- the "wildcard" behavior that is everywhere in chemistry.

The rule below reduces an **aldehyde** (or thio-aldehyde) to the corresponding
primary alcohol/thiol with H2, but leaves **ketones** alone. Two ideas combine:

* the carbonyl partner ``_X`` is a *term variable*, so the single rule reduces
  both C=O and C=S carbonyls (and would extend to C=N, ...);
* aldehyde-vs-ketone selectivity is encoded with a **term constraint**. The
  reacting carbon's three neighbors are the variables ``_X`` (the double-bonded
  hetero-atom) and ``_Y`` / ``_Z`` (its two single-bonded substituents). The
  ``constrainLabelAny`` block requires the tuple ``Ald(_X, _Y, _Z)`` to unify
  with one of ``Ald(O, H, C)`` / ``Ald(S, H, C)`` -- i.e. an oxo/thio carbonyl
  whose other substituents are one hydrogen and one carbon. That is exactly an
  aldehyde; a ketone (``Ald(O, C, C)``) never matches.

This mirrors the classic ``term-aldehyd-not-keton`` example: the ``Ald`` function
constrains ``_X`` to oxygen or sulfur, ``_Y`` to hydrogen, and ``_Z`` to carbon.

A second example below constrains **bond** labels instead of atom labels, to
match either resonance (Kekule) structure of a benzene ring.

(Requires the Term / Unification label setting, passed to the DG.)

Run:  python 05_terms.py
"""
import _bootstrap
from molfoundry import (
    smiles, ruleGMLString, DG, LabelSettings, LabelType, LabelRelation,
    addSubset, repeat, inputGraphs, enable_post,
)

enable_post()

reduce_aldehyde = ruleGMLString(
    """rule [
        ruleID "reduce aldehyde, not ketone"
        labelType "term"
        left [
            edge [ source 1 target 2 label "=" ]
            edge [ source 5 target 6 label "-" ]
        ]
        context [
            node [ id 1 label "C" ]
            node [ id 2 label "_X" ]
            node [ id 3 label "_Y" ]
            node [ id 4 label "_Z" ]
            node [ id 5 label "H" ]
            node [ id 6 label "H" ]
            edge [ source 1 target 3 label "-" ]
            edge [ source 1 target 4 label "-" ]
        ]
        right [
            edge [ source 1 target 2 label "-" ]
            edge [ source 1 target 6 label "-" ]
            edge [ source 2 target 5 label "-" ]
        ]
        constrainLabelAny [
            label "Ald(_X, _Y, _Z)"
            labels [ label "Ald(O, H, C)" label "Ald(S, H, C)" ]
        ]
    ]"""
)

# Four substrates spanning {aldehyde, ketone} x {oxo, thio}. The constraint
# reduces every carbonyl carbon that carries a hydrogen (an aldehyde), whether
# its double-bond partner is O or S, and leaves the fully substituted ones
# (ketones / thio-ketones) untouched.
smiles("[H][H]", "H2")
smiles("O=CC=S", "oxal-dial")        # aldehyde + thio-aldehyde -> both reduce
smiles("CC(=O)C(=S)", "keto-thial")  # ketone + thio-aldehyde  -> only the thial
smiles("O=CC(=S)C", "al-thioketone") # aldehyde + thio-ketone  -> only the aldehyde
smiles("CC(=O)C(=S)C", "di-ketone")  # ketone + thio-ketone    -> no reaction

# Term rewriting requires the Unification label relation.
ls = LabelSettings(LabelType.Term, LabelRelation.Unification)
dg = DG(graphDatabase=inputGraphs, labelSettings=ls)
dg.build().execute(addSubset(inputGraphs) >> repeat[2]([reduce_aldehyde]), verbosity=0)
dg.print()

print(f"DG: {dg.numVertices} molecules, {dg.numEdges} reactions")
print("Reductions (note: di-ketone, with no aldehyde carbon, produces no reaction):")
for e in dg.edges:
    educts = " + ".join(sorted(v.graph.name for v in e.sources))
    products = " + ".join(sorted(v.graph.smiles for v in e.targets))
    print(f"    {educts:26s} ->  {products}")


# -----------------------------------------------------------------------------
# A second example: a term constraint on BOND labels (mesomerism / resonance).
# -----------------------------------------------------------------------------
# Benzene has no single Lewis structure -- its ring is the average of two Kekule
# forms whose single/double bonds alternate the opposite way. This rule reduces
# an aldehyde attached to a benzene ring, but only when the ring's six bonds
# match one of those two alternations, stated by constraining the tuple of bond
# variables Tauto(_X, ..., _W).
match_kekule = ruleGMLString("""
rule [
    ruleID "reduce aldehyde on a benzene ring (either Kekule form)"
    labelType "term"
    left [
        edge [ source 1 target 9 label "=" ]
        edge [ source 10 target 11 label "-" ]
    ]
    context [
        node [ id 1 label "C" ] node [ id 2 label "C" ] node [ id 3 label "C" ]
        node [ id 4 label "C" ] node [ id 5 label "C" ] node [ id 6 label "C" ]
        node [ id 7 label "C" ] node [ id 9 label "O" ]
        node [ id 10 label "H" ] node [ id 11 label "H" ]
        edge [ source 1 target 2 label "-" ]
        edge [ source 2 target 3 label "_X" ] edge [ source 3 target 4 label "_U" ]
        edge [ source 4 target 5 label "_Y" ] edge [ source 5 target 6 label "_V" ]
        edge [ source 6 target 7 label "_Z" ] edge [ source 7 target 2 label "_W" ]
    ]
    right [
        edge [ source 1 target 10 label "-" ]
        edge [ source 1 target 9 label "-" ]
        edge [ source 9 target 11 label "-" ]
    ]
    constrainLabelAny [
        label "Tauto(_X,_U,_Y,_V,_Z,_W)"
        labels [ label "Tauto(-, =, -, =, -, =)" label "Tauto(=, -, =, -, =, -)" ]
    ]
]
""", add=False)

h2 = smiles("[H][H]", "H2", add=False)
# Two naphthalene-dialdehydes: in TAU-1 both fused six-rings match a Kekule
# pattern (so both aldehydes reduce); in TAU-2 only one ring does.
print("\nMesomerism -- one rule matching either Kekule structure:")
for name, smi in [("TAU-1", "O=CC(=C1)C2=C(C=C1)C(C=O)=CC=C2C"),
                  ("TAU-2", "O=CC1=C2C(=CC=C1)C(C=O)=CC=C2C")]:
    mol = smiles(smi, name, add=False)
    dg2 = DG(graphDatabase=[h2, mol], labelSettings=ls)
    dg2.build().execute(addSubset([h2, mol]) >> repeat[3]([match_kekule]), verbosity=0)
    print(f"    {name}: {dg2.numEdges} reduction(s)  "
          f"({'both rings match' if dg2.numEdges > 1 else 'only one ring matches'})")

Rule composition#

If a reaction sequence A —r₁→ B —r₂→ C keeps recurring, you can fuse r₁ and r₂ into a single rule r₁∘r₂ that does A → C in one step. This is the DPO concurrency construction: glue the product side of the first rule onto the educt side of the second along an overlap, and read off the composite — with the atom map carried through. MolFoundry exposes it through RCEvaluator and four operators; rcCommon enumerates every common overlap and returns one composite per overlap.

Two staples of terpene cation chemistry — an allyl-cation shift and a proton loss — written in DFS and composed:

allyl_shift = ruleDFS("[C+]1[C]2=[C]3>>[C]1=[C]2[C+]3", "allyl-shift", add=False)
h_loss      = ruleDFS("[C+]1[C]2[H]3>>[C]1=[C]2.[H+]3", "H-loss", add=False)
composites  = RCEvaluator().eval(allyl_shift * rcCommon * h_loss)
print(f"{len(composites)} overlaps ->",
      len({r.canonKey for r in composites}), "distinct composed rules")
composites[0].print(rulep)

4 overlaps -> 4 distinct composed rules

common

Each result is an ordinary Rule you can print, apply, or compose again.

12_rule_composition.py#

"""12. Rule composition.

If a reaction sequence A --r1--> B --r2--> C recurs, you can fuse r1 and r2 into
a single rule r1 o r2 that does A --> C in one step (composition preserves the
atom map across the sequence). MolFoundry composes rules with an *evaluator* and
composition *operators*:

    RCEvaluator().eval(r1 * rcCommon * r2)

``rcCommon`` enumerates every way the product side of r1 can overlap the educt
side of r2 -- i.e. every subgraph they can share -- and returns one composed
rule per overlap. (Other operators: ``rcSuper`` / ``rcSub`` require one side to
fully contain the other; ``rcParallel`` places the rules side by side.)

Here r1 is an allyl-cation shift and r2 loses a proton to make a double bond --
two staples of terpene cation chemistry -- both written in the compact DFS line
notation (the ids carry the atom map across the ``>>``).

Run:  python 12_rule_composition.py
"""
import _bootstrap
from molfoundry import ruleDFS, GraphPrinter, DrawMode, enable_post
from molfoundry.rulecomposition import RCEvaluator, rcCommon

enable_post()

# C(+)-C=C  ->  C=C-C(+)   and   C(+)-C-H  ->  C=C  +  H(+)
allyl_shift = ruleDFS("[C+]1[C]2=[C]3>>[C]1=[C]2[C+]3", "allyl-shift")
h_loss = ruleDFS("[C+]1[C]2[H]3>>[C]1=[C]2.[H+]3", "H-loss")

p = GraphPrinter()
p.drawMode = DrawMode.Circles
allyl_shift.print(p)
h_loss.print(p)

# Compose. Deduplicate the results by canonical key: distinct overlaps can lead
# to the same composed rule.
composites = RCEvaluator().eval(allyl_shift * rcCommon * h_loss)
unique = {}
for r in composites:
    unique.setdefault(r.canonKey, r)

print(f"allyl-shift (o) H-loss: {len(composites)} overlap(s), "
      f"{len(unique)} distinct composed rule(s)")
for i, r in enumerate(unique.values()):
    r.print(p)
    print(f"\n--- composed rule #{i} "
          f"(left components: {r.numLeftComponents}, right: {r.numRightComponents}) ---")
    print(r.getGMLString())