# 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):

```python
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:

```python
# 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

```python
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}
:label: ex-arom

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} ex-arom
:class: dropdown

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](spectra-dos.md), turns from bonding analysis to orbital-energy
spectra — the density of states and the HOMO–LUMO gap.