Working with molecular graphs#
A molecule graph is read-only, but you can walk it freely and compute on it.
The vertex and edge ranges expose everything you need:
graph.vertices / graph.edges; a vertex has id, stringLabel, atomId (the
element — .symbol, and int(...) for the atomic number), charge, degree
and incidentEdges; an edge has source, target, stringLabel and
bondType. This section builds a few classic descriptors from those primitives.
Throughout, “ignore hydrogens” is the test v.atomId.symbol == "H" on a neutral
atom — the H-depleted graph is the usual object of interest.
The adjacency list#
Because a bond stores a fixed (source, target), the neighbor of a vertex across
an edge is simply “the other end”. Collecting those gives the adjacency list.
Running it on cocaine — written with its methyl and phenyl groups spelled out,
since a SMILES bracket atom must be a real element — reports 22 heavy atoms and
their connectivity:
cocaine = smiles("CN1C2CCC1C(C(=O)OC)C(C2)OC(=O)c1ccccc1", "cocaine", add=False)
def is_h(v): return v.atomId.symbol == "H" and int(v.charge) == 0
heavy = [v for v in cocaine.vertices if not is_h(v)]
print(cocaine.numVertices, "atoms,", len(heavy), "heavy")
cocaine.print(_p)
43 atoms, 22 heavy
cocaine
"""9. Reading a molecule graph: the adjacency list.
Molecule graphs are read-only, but you can walk them freely:
* ``graph.vertices`` / ``graph.edges`` -- the atoms and bonds;
* ``v.id``, ``v.stringLabel``, ``v.atomId`` (element), ``v.charge``;
* ``v.incidentEdges`` -- the bonds at an atom;
* ``edge.source`` / ``edge.target`` / ``edge.bondType``.
Because a bond stores a fixed (source, target), the neighbor of ``v`` across an
edge is "the other end". A common option is to ignore hydrogens (the H-depleted
skeleton); here that is ``v.atomId.symbol == "H"`` on a neutral atom.
Run: python 09_adjacency_list.py
"""
import _bootstrap
from molfoundry import smiles, GraphPrinter, enable_post
enable_post()
# Cocaine. (A bracket atom in SMILES must be a real element, so the methyl and
# phenyl groups are written out in full rather than as [Me] / [Ph] shorthands.)
cocaine = smiles("CN1C2CCC1C(C(=O)OC)C(C2)OC(=O)c1ccccc1", "cocaine")
cocaine.print(GraphPrinter())
def is_hydrogen(v):
return v.atomId.symbol == "H" and int(v.charge) == 0
def other_end(edge, v):
return edge.target if edge.source.id == v.id else edge.source
def adjacency_list(m, ignore_hs=True):
adj = {}
for v in m.vertices:
if ignore_hs and is_hydrogen(v):
continue
neighbors = {
other_end(e, v).id
for e in v.incidentEdges
if not (ignore_hs and is_hydrogen(other_end(e, v)))
}
adj[v.id] = neighbors
return adj
print(f"cocaine: {cocaine.numVertices} atoms, {cocaine.numEdges} bonds "
f"(heavy atoms: {sum(1 for v in cocaine.vertices if not is_hydrogen(v))})")
print("\n=== atoms (heavy only) ===")
for v in cocaine.vertices:
if not is_hydrogen(v):
print(f" {v.stringLabel:2s} {v.id:2d}")
print("\n=== adjacency list (H-depleted) ===")
for vid, adj in adjacency_list(cocaine).items():
print(f" {vid:2d} => {sorted(adj)}")
Removing hydrogens#
The H-depleted graph — the molecular skeleton — is a common starting point.
MolFoundry has removeHydrogens() built in, and the
script shows two more ways (by hand from the edge list, and with a graph
grammar) that all yield the same graph. For neopentane the skeleton’s SMILES
writes every carbon in brackets, [C]C([C])([C])[C], because a bracket-less C
would imply the hydrogens back:
neopentane = smiles("CC(C)(C)C", "neopentane", add=False)
skeleton = neopentane.removeHydrogens()
print(skeleton.smiles)
skeleton.print(_p)
[C]C([C])([C])[C]
neopentane
The grammar version is a nice first term rule: ruleDFS("[*]1[H]2>>[*]1") strips
one hydrogen from any atom (the * wildcard), applied once per hydrogen.
"""10. Removing hydrogens (the molecular skeleton).
The H-depleted graph -- every atom minus the hydrogens -- is the classic
"molecular skeleton". Three ways to get it, all yielding the same graph:
1. ``graph.removeHydrogens()`` -- the built-in, one call.
2. By hand with graph operations: keep the non-H atoms and the bonds between
them, and rebuild the graph from GML.
3. With a graph grammar: a rule that deletes one H, applied until none are left.
The example uses neopentane, C(C)(C)(C)C. Its skeleton's SMILES writes every
carbon in brackets, ``[C]C([C])([C])[C]``, because a bracket-less ``C`` would
imply the hydrogens back.
Run: python 10_removing_hydrogens.py
"""
import _bootstrap
from molfoundry import (
smiles, graphGMLString, ruleDFS, DG,
LabelSettings, LabelType, LabelRelation, addSubset, repeat, enable_post,
)
enable_post()
neopentane = smiles("CC(C)(C)C", "neopentane")
def is_hydrogen(v):
return v.atomId.symbol == "H" and int(v.charge) == 0
# --- 1. the built-in --------------------------------------------------------
skeleton = neopentane.removeHydrogens()
print("removeHydrogens(): ", skeleton.smiles)
# --- 2. by hand, via GML ----------------------------------------------------
def remove_hs_by_hand(m):
kept_edges = [e for e in m.edges
if not (is_hydrogen(e.source) or is_hydrogen(e.target))]
kept_vertices = {v.id: v for e in kept_edges for v in (e.source, e.target)}
nodes = "\n".join(f' node [ id {v.id} label "{v.stringLabel}" ]'
for v in kept_vertices.values())
edges = "\n".join(f' edge [ source {e.source.id} target {e.target.id} '
f'label "{e.stringLabel}" ]' for e in kept_edges)
return graphGMLString(f"graph [\n{nodes}\n{edges}\n]", add=False)
print("by hand (graph ops): ", remove_hs_by_hand(neopentane).smiles)
# --- 3. with a graph grammar ------------------------------------------------
# A rule that strips one hydrogen from any atom, written in DFS line notation:
# "[*]1[H]2>>[*]1" keeps the wildcard atom (id 1) and drops the bonded H (id 2).
# The "*" label is a wildcard, so the rule is a term rule.
strip_h = ruleDFS("[*]1[H]2>>[*]1", "strip one H", labelType=LabelType.Term, add=False)
n_h = neopentane.vLabelCount("H")
ls = LabelSettings(LabelType.Term, LabelRelation.Unification)
dg = DG(graphDatabase=[neopentane], labelSettings=ls)
dg.build().execute(addSubset(neopentane) >> repeat[n_h]([strip_h]), verbosity=0)
# The fully-stripped graph is the one with no hydrogens left.
stripped = min((v.graph for v in dg.vertices), key=lambda g: g.vLabelCount("H"))
print("with a graph grammar: ", stripped.smiles)
The Wiener index#
The Wiener index is one of the oldest topological descriptors: the sum of all pairwise shortest-path distances in the H-depleted graph,
a measure of compactness that correlates with the boiling points of alkanes. It falls out of three steps on the molecule graph — adjacency matrix, then a Floyd–Warshall distance matrix, then the sum of the upper triangle. For 2-methylpentane the answer is 32:
mol = smiles("CC(C)CCC", "2-methylpentane", add=False)
# (the script's mol2adjMx / adjMx2distMx / wienerIndex helpers give W = 32)
print(mol.smiles, "-> Wiener index 32")
mol.print(_p)
C(CCC)(C)C -> Wiener index 32
2-methylpentane
"""11. The Wiener index.
The Wiener index W is the sum of all pairwise topological distances in the
H-depleted molecular graph -- an old but useful compactness descriptor (it
correlates with boiling points of alkanes). Computed here in three steps on the
molecule graph:
molecule graph -> adjacency matrix (heavy atoms)
-> distance matrix (Floyd-Warshall)
-> W = sum of the upper triangle
For 2-methylpentane, CC(C)CCC, the answer is 32.
Run: python 11_wiener_index.py
"""
import _bootstrap
from molfoundry import smiles, enable_post
enable_post()
INF = float("inf")
def is_hydrogen(v):
return v.atomId.symbol == "H" and int(v.charge) == 0
def adjacency_matrix(m):
heavy = [v.id for v in m.vertices if not is_hydrogen(v)]
index = {vid: i for i, vid in enumerate(heavy)}
n = len(heavy)
a = [[0] * n for _ in range(n)]
for e in m.edges:
if is_hydrogen(e.source) or is_hydrogen(e.target):
continue
i, j = index[e.source.id], index[e.target.id]
a[i][j] = a[j][i] = 1
return a
def distance_matrix(a):
n = len(a)
d = [[0 if i == j else (1 if a[i][j] else INF) for j in range(n)] for i in range(n)]
for k in range(n): # Floyd-Warshall
for i in range(n):
for j in range(n):
if d[i][k] + d[k][j] < d[i][j]:
d[i][j] = d[i][k] + d[k][j]
return d
def wiener_index(d):
return int(sum(d[i][j] for i in range(len(d)) for j in range(i + 1, len(d))))
mol = smiles("CC(C)CCC", "2-methylpentane")
a = adjacency_matrix(mol)
d = distance_matrix(a)
w = wiener_index(d)
print("distance matrix (heavy atoms):")
for row in d:
print(" ", [int(x) for x in row])
print(f"\nWiener index of {mol.smiles} = {w}")
assert w == 32, w
The automorphism group#
A graph automorphism is a relabeling of the atoms that maps the molecule onto
itself — it is the molecule’s symmetry. MolFoundry computes it natively (no
external group-theory package): aut() returns the
automorphism group in cycle notation, with .gens a generating set, and the
permutations are always label-preserving (a hydrogen never maps to a carbon).
benzene = smiles("c1ccccc1", "benzene", add=False)
print("automorphism group:", benzene.aut())
benzene.print(_p)
automorphism group: [(), (1 2)(10 9)(4 5)(7 8), (0 4 2 3 1 5)(10 7 6 8 9 11), (0 5)(1 4)(10 8)(11 9)(2 3)(6 7)]
benzene
"""14. The automorphism group of a molecule.
A graph automorphism is a relabeling of the atoms that maps the molecule onto
itself -- it captures the molecule's symmetry. MolFoundry computes it natively:
``graph.aut()`` returns the automorphism group, printed in cycle notation, with
``.gens`` giving a generating set. (No external group-theory package needed.)
Run: python 14_automorphism_group.py
"""
import _bootstrap
from molfoundry import smiles, GraphPrinter, DrawMode, enable_post
enable_post()
examples = [
("methane", "C"), # the 4 H's are interchangeable (S4)
("tert-butyl skeleton", "[C]([C])([C])[C]"),
("benzene", "c1ccccc1"), # 6-fold ring symmetry (+ reflections)
]
p = GraphPrinter()
p.drawMode = DrawMode.Circles
for name, smi in examples:
g = smiles(smi, name, add=False)
g.print(p)
aut = g.aut()
gens = aut.gens
print(f"{name:22s} {g.smiles}")
print(f" automorphism group: {aut}")
print(f" {len(gens)} generator(s):")
for gen in gens:
# Each generator maps vertex v -> gen[v]; str() prints it as cycles.
print(f" {gen}")
Simple cycles#
Enumerating a molecule’s elementary (simple) cycles is a good illustration of
walking the graph and then running a plain algorithm on the edge list — no
external graph library needed. Building a heavy-atom adjacency and doing a rooted
depth-first search (each cycle found once, smallest atom first) gives, for the
polycyclic cage C1C2C3C4C1C4=C23, its 14 elementary cycles: two 3-, one 4-,
four 5-, five 6- and two 7-membered.
cage = smiles("C1C2C3C4C1C4=C23", "cage", add=False)
cage.print(_p)
cage
"""15. Elementary cycles of a molecular graph.
An elementary (simple) cycle visits each of its atoms once. Enumerating them is
a nice illustration of walking a MolFoundry graph and then running a plain graph
algorithm on the edge list -- no external graph library required.
The polycyclic cage C1C2C3C4C1C4=C23 has 14 elementary cycles: two 3-, one 4-,
four 5-, five 6- and two 7-membered.
Run: python 15_simple_cycles.py
"""
import _bootstrap
from molfoundry import smiles, GraphPrinter, DGPrinter, enable_post
enable_post()
def is_hydrogen(v):
return v.atomId.symbol == "H" and int(v.charge) == 0
def adjacency(m):
"""Heavy-atom adjacency as {id: set(neighbor ids)}."""
adj = {v.id: set() for v in m.vertices if not is_hydrogen(v)}
for e in m.edges:
if is_hydrogen(e.source) or is_hydrogen(e.target):
continue
adj[e.source.id].add(e.target.id)
adj[e.target.id].add(e.source.id)
return adj
def simple_cycles(adj):
"""All elementary cycles, each canonicalized once (smallest atom first, and
the lexicographically smaller of the two traversal directions)."""
found = set()
for start in sorted(adj):
# DFS from `start`, visiting only atoms > start, so `start` is the
# unique minimum of every cycle it roots (each cycle is found once).
stack = [(start, (start,))]
while stack:
v, path = stack.pop()
for w in adj[v]:
if w == start and len(path) >= 3:
forward = path
backward = (path[0],) + path[1:][::-1]
found.add(min(forward, backward))
elif w > start and w not in path:
stack.append((w, path + (w,)))
return found
mol = smiles("C1C2C3C4C1C4=C23", "cage")
mol.print(GraphPrinter())
adj = adjacency(mol)
cycles = simple_cycles(adj)
by_length = sorted(len(c) for c in cycles)
print(f"molecule {mol.smiles}: {len(adj)} heavy atoms")
print(f"elementary cycles: {len(cycles)}")
print(f"length distribution: {by_length}")
for c in sorted(cycles, key=len):
print(f" {len(c)}-cycle: {list(c)}")