Aromaticity indices

Aromaticity — the special stability and delocalization of conjugated rings like benzene — is one of chemistry’s most useful concepts and, like atomic charge, one with no single definition. It shows up in geometry (equalized bond lengths), in electron delocalization, and in magnetic response, so different indices probe different facets. qc-rs collects them under qc.prop.arom; this chapter shows how to compute and compare them, using benzene as the reference aromatic.

Theory: three windows on aromaticity

An aromatic ring differs from a non-aromatic one in measurable ways, and each family of index looks through a different window:

• Geometric — an aromatic ring has equalized bond lengths (no localized single/double alternation). HOMA (harmonic-oscillator model of aromaticity) scores this with a single formula over the \(n\) ordered ring bonds, $\( \text{HOMA} = 1 - \frac{1}{n}\sum_{i=1}^n \alpha\,(R_{\text{opt}}-R_i)^2 , \)\( where \)R_i\( is bond \)i\('s actual length, \)R_{\text{opt}}\( is the tabulated bond length of a *perfectly aromatic* bond of that atom-pair type (Krygowski's reference values, one table entry per bond type — C–C, C–N, …), and \)\alpha\( is a normalization constant chosen so a fully-localized reference ring gives HOMA \)=0\(. A perfectly aromatic ring (\)R_i=R_{\text{opt}}\( everywhere) scores exactly **HOMA \)=1$**; any deviation — bonds too short, too long, or unequal — subtracts a squared penalty, pushing HOMA toward 0 or negative as bonds localize. BIRD is a related bond-order-based index using the same squared-deviation idea applied to Gordy bond orders instead of lengths.

• Electronic / delocalization — aromaticity means electrons are shared around the ring. FLU (aromatic fluctuation index, ≈ 0 when aromatic), PDI (para-delocalization index), MCI (multicenter index), and I_ring measure this delocalization from the density matrix.

• Information-theoretic / other — Shannon-entropy-based and related descriptors (shannon, info_theoretic, lolipop, …).

Because they measure different things, you report several and look for agreement rather than trusting a single number.

Usage

The indices take the ring — the list of atom indices making up the ring, in order. For benzene (carbons at indices 0, 2, 4, 6, 8, 10 in the geometry below):

import qc, math
# build a D6h benzene
R, Rh = 1.39, 1.39 + 1.08
atoms = []
for i in range(6):
    a = math.radians(60 * i)
    atoms.append(f"C {R*math.cos(a):.4f} {R*math.sin(a):.4f} 0")
    atoms.append(f"H {Rh*math.cos(a):.4f} {Rh*math.sin(a):.4f} 0")
benzene = ";".join(atoms)

m = qc.chk.new(atom=benzene, ao="6-31g", unit="angstrom").scf(ref="r").run()   # E = -230.624263
ring = [0, 2, 4, 6, 8, 10]   # the six carbons

qc.prop.arom.homa(m, ring=ring)      # 0.999    geometric — essentially perfect
qc.prop.arom.bird(m, ring=ring)      # 99.99    bond-order geometric index (~100 = aromatic)
qc.prop.arom.indices(m, ring=ring)   # {'pdi': 0.093, 'flu': 0.001, 'iring': 0.034, 'mci': 0.051, ...}

Benzene scores as textbook-aromatic on every index: HOMA 0.999 (≈ 1), BIRD 99.99 (≈ 100), and a near-zero FLU 0.0006 — all three windows agree the ring is fully aromatic. The indices leaf bundles the delocalization measures (PDI, FLU, I_ring, MCI, and the AV1245/AVmin descriptors for larger rings) in one call.

Tip

Getting the ring right The ring= list must be the atoms of the ring, given in connectivity order around the cycle (not arbitrary). Getting the order wrong scrambles the multicenter indices (MCI, I_ring), which depend on going around the ring. The geometric HOMA is more forgiving.

Comparing rings

Aromaticity indices are most useful relative — is this ring more aromatic than that one? Compute the same index on each ring and compare:

# schematic: compare a benzene ring to a non-aromatic ring
homa_benzene = qc.prop.arom.homa(m, ring=ring)      # ~0.999 (aromatic)
# homa_cyclohexane would be far below 1 (localized single bonds)

A HOMA near 1 and a FLU near 0 flag an aromatic ring; a HOMA well below 1 (or negative) and a larger FLU flag a localized, non-aromatic (or anti-aromatic) one.

Worked example

import qc, math
R, Rh = 1.39, 2.47
atoms = []
for i in range(6):
    a = math.radians(60 * i)
    atoms += [f"C {R*math.cos(a):.4f} {R*math.sin(a):.4f} 0",
              f"H {Rh*math.cos(a):.4f} {Rh*math.sin(a):.4f} 0"]
m = qc.chk.new(atom=";".join(atoms), ao="6-31g", unit="angstrom").scf(ref="r").run()
ring = [0, 2, 4, 6, 8, 10]

print("HOMA :", round(qc.prop.arom.homa(m, ring=ring), 3))          # 0.999
print("BIRD :", round(qc.prop.arom.bird(m, ring=ring), 2))          # 99.99
print("FLU  :", round(qc.prop.arom.indices(m, ring=ring)["flu"], 4)) # 0.0006

Exercise 14

1. A ring gives HOMA = 0.30 and FLU = 0.05. Aromatic or not? Do the two indices agree, and what does the low HOMA physically mean?

2. Why is it good practice to report both a geometric index (HOMA) and a delocalization index (FLU/MCI) rather than just one?

3. Your arom.indices MCI value looks wrong for a ring you know is aromatic. What is the most likely input mistake?

Solution to Exercise 14

1. Not aromatic (localized). HOMA = 0.30 is far below the aromatic ≈ 1, and FLU = 0.05 is well above the aromatic ≈ 0 — the two agree. The low HOMA means the ring has alternating/unequal bond lengths (localized single and double bonds) rather than the equalized bonds of an aromatic ring.

2. They probe different facets (geometry vs electron delocalization) that usually but not always agree; requiring both to point the same way guards against a false positive from any single measure.

3. The ring= atom order — the multicenter index MCI depends on traversing the ring in connectivity order; a scrambled list gives a meaningless MCI even though the geometric HOMA still looks fine.

The last property chapter, spectra & DOS, turns from bonding analysis to orbital-energy spectra — the density of states and the HOMO–LUMO gap.