Kekulé’s benzene: enumerating the Kekulé structures of aromatic hydrocarbons

August Kekulé (1865) proposed that benzene’s six carbons form a ring with alternating single and double bonds, and soon after that the two alternating arrangements interconvert, which explained why benzene has only one ortho-disubstituted isomer. A Kekulé structure is a choice of double bonds that gives every carbon exactly one; in graph language, a perfect matching of the carbon skeleton.

This example uses kekule_structures() to list benzene’s two Kekulé structures and naphthalene’s three, and counts them for larger benzenoid hydrocarbons: a linear acene with \(r\) rings has \(r+1\), while angular phenanthrene has five, one reason it is more stable than its linear isomer anthracene (four).

import matplotlib.pyplot as plt
import numpy as np

from chemistrykit.structure.systems.kekule import count_kekule_structures, kekule_structures


def hexagon(cx, cy, start=90.0):
    """Six vertex positions of a unit hexagon centred at (cx, cy)."""
    a = np.radians(start + 60.0 * np.arange(6))
    return np.column_stack([cx + np.cos(a), cy + np.sin(a)])


# Benzene: atoms 0-5 around the ring.
benzene_xy = hexagon(0.0, 0.0)
benzene_bonds = [(i, (i + 1) % 6) for i in range(6)]

# Naphthalene: two hexagons sharing the bond between atoms 4 and 9.
naphthalene_bonds = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 9), (9, 0), (4, 5), (5, 6), (6, 7), (7, 8), (8, 9)]
left = hexagon(-np.sqrt(3) / 2, 0.0, start=30.0)
right = hexagon(np.sqrt(3) / 2, 0.0, start=30.0)
naphthalene_xy = np.array([left[1], left[2], left[3], left[4], left[5], right[4], right[5], right[0], right[1], right[2]])

benzene_k = kekule_structures(6, benzene_bonds)
naphthalene_k = kekule_structures(10, naphthalene_bonds)
print(f"benzene: {len(benzene_k)} Kekule structures -> {benzene_k}")
print(f"naphthalene: {len(naphthalene_k)} Kekule structures")
assert len(benzene_k) == 2 and len(naphthalene_k) == 3

benzene: 2 Kekule structures -> [[(0, 1), (2, 3), (4, 5)], [(0, 5), (1, 2), (3, 4)]]
naphthalene: 3 Kekule structures

Each benzene bond is double in exactly one of the two structures, so averaging them gives every bond the same order, 1.5. That averaging is how Kekulé’s oscillation hypothesis makes the six bonds equivalent.

double_fraction = [sum((min(b), max(b)) in s for s in benzene_k) / len(benzene_k) for b in benzene_bonds]
print(f"fraction of Kekule structures in which each benzene bond is double: {double_fraction}")

fraction of Kekule structures in which each benzene bond is double: [0.5, 0.5, 0.5, 0.5, 0.5, 0.5]

Kekulé structure counts of larger benzenoid hydrocarbons:

def acene(n_rings):
    """Linear acene skeleton with n_rings fused rings (4n+2 carbons)."""
    top = list(range(2 * n_rings + 1))
    bottom = list(range(2 * n_rings + 1, 4 * n_rings + 2))
    bonds = [(a, a + 1) for a in top[:-1]] + [(b, b + 1) for b in bottom[:-1]]
    bonds += [(top[2 * k], bottom[2 * k]) for k in range(n_rings + 1)]
    return 4 * n_rings + 2, bonds


phenanthrene = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 0), (4, 6), (6, 7), (7, 8), (8, 9), (9, 3), (8, 10), (10, 11), (11, 12), (12, 13), (13, 9)]
counts = {name: count_kekule_structures(*acene(r)) for r, name in enumerate(["benzene", "naphthalene", "anthracene", "tetracene", "pentacene"], start=1)}
counts["phenanthrene"] = count_kekule_structures(14, phenanthrene)
for name, k in counts.items():
    print(f"{name:12s}: K = {k}")
assert counts["anthracene"] == 4 and counts["phenanthrene"] == 5

benzene     : K = 2
naphthalene : K = 3
anthracene  : K = 4
tetracene   : K = 5
pentacene   : K = 6
phenanthrene: K = 5

Draw benzene’s two and naphthalene’s three Kekulé structures:

def draw(ax, xy, bonds, doubles, title):
    double_set = {tuple(sorted(b)) for b in doubles}
    for i, j in bonds:
        is_double = tuple(sorted((i, j))) in double_set
        ax.plot(*zip(xy[i], xy[j]), color="crimson" if is_double else "black", linewidth=4.0 if is_double else 1.5)
    ax.set_aspect("equal")
    ax.axis("off")
    ax.set_title(title, fontsize=10)


fig, axes = plt.subplots(1, 5, figsize=(15, 3.2))
for k, s in enumerate(benzene_k):
    draw(axes[k], benzene_xy, benzene_bonds, s, f"benzene, structure {k + 1}")
for k, s in enumerate(naphthalene_k):
    draw(axes[2 + k], naphthalene_xy, naphthalene_bonds, s, f"naphthalene, structure {k + 1}")
fig.suptitle("Kekulé structures (double bonds drawn thick red)")
fig.tight_layout()
plt.show()