# Basis sets & AO representation

The [last chapter](molecular-input.md) defined *which* molecule to compute. This one is about the single
biggest lever on the **cost–accuracy trade-off** of a calculation: the **basis set** — the set of functions
in which qc-rs expands the molecular orbitals. Choose it with the `ao=` argument. Get a feel for it here and
you will make good, deliberate choices instead of copying `cc-pvdz` out of habit.

## Theory: why a basis, and why Gaussians

Recall the **LCAO** idea from [Part II](../10-foundations/variational-lcao-basis.md): each molecular orbital
is written as a linear combination of fixed, atom-centred **basis functions** $\{\phi_\mu\}$,

$$
\psi_i(\mathbf r) = \sum_{\mu=1}^{K} C_{\mu i}\,\phi_\mu(\mathbf r).
$$

The set $\{\phi_\mu\}$ *is* the basis, and its size $K$ is the dimension of the Roothaan matrix problem
$\mathbf{FC}=\mathbf{SC}\boldsymbol\varepsilon$. A bigger, more flexible basis lets the SCF describe the true
orbitals more faithfully — at a cost that grows steeply with $K$ (formally $\mathcal O(K^4)$ for the
integrals). Understanding the basis is understanding where your accuracy and your runtime both come from.

What should the $\phi_\mu$ look like? The exact atomic orbitals decay like $e^{-\zeta r}$ — a **Slater-type
orbital** (STO), with a sharp cusp at the nucleus and the right long-range tail. Unfortunately the
two-electron integrals over STOs are prohibitively expensive. The trick that made quantum chemistry
practical is to use **Gaussian** functions $e^{-\alpha r^2}$ instead: their integrals have closed forms (the
Gaussian product theorem turns a product of two Gaussians on different centres into one Gaussian), so they
are *fast*. A single Gaussian is the wrong shape — no cusp, too-fast decay — so we fix that by gluing
several together:

:::{prf:definition} Contracted Gaussian basis function
:label: def-cgto

A **contracted Gaussian-type orbital (CGTO)** is a *fixed* linear combination of **primitive** Gaussians
sharing a centre and angular momentum,

$$
\phi_\mu(\mathbf r) = \sum_{k=1}^{n_\mu} d_{k\mu}\; g_k(\mathbf r),
\qquad
g_k(\mathbf r) \propto x^{a} y^{b} z^{c}\, e^{-\alpha_k r^2},
$$

with the **contraction coefficients** $d_{k\mu}$ and **exponents** $\alpha_k$ fixed by the basis-set
designer. Several primitives combine to mimic the near-cusp Slater shape, while the SCF only ever varies the
*one* coefficient $C_{\mu i}$ per contracted function — so a CGTO gives near-STO accuracy at Gaussian speed.
:::

### Primitive normalization

Before the contraction coefficients $d_{k\mu}$ mean what a basis-set table intends, each primitive $g_k$ must
itself be normalized: $\langle g_k|g_k\rangle=1$. For an angular-momentum-$l$ primitive of exponent $\alpha$,
carrying out that Gaussian integral gives

$$
N(l,\alpha) = \left(\frac{2}{\pi}\right)^{3/4} 2^{\,l}\, \alpha^{(2l+3)/4} \Big/ \sqrt{(2l-1)!!}\,,
$$

where $(2l-1)!! = 1\cdot 3\cdot 5\cdots(2l-1)$ is the double factorial (with $(-1)!!:=1$ for $l=0$, so
$N(0,\alpha)=(2/\pi)^{3/4}\alpha^{3/4}$ for an $s$ primitive). qc-rs applies $N(l,\alpha)$ once, when the
basis is parsed (`qc-mol`), and stores **both** the raw literature coefficients and this primitive-normalized
set side by side — you never call a normalization routine yourself, and the contraction sum above is always
taken with the normalized set.

## Anatomy of a basis set

Basis sets are built up in layers, and the names encode which layers are present:

- **Minimal** (e.g. **STO-3G** — 3 Gaussians contracted per function): exactly one function per occupied
  atomic orbital. Cheap, qualitative, rarely quantitative.
- **Split-valence** (the Pople **6-31G** family): the chemically active *valence* orbitals get **two or more**
  functions (a tight "inner" + a loose "outer"), so a bond can expand or contract. The core stays minimal.
- **Polarization functions** (the `*`/`**` in `6-31G*`, or the built-in `d`,`f`,… of the cc-p sets): higher
  angular-momentum functions ($d$ on heavy atoms, $p$ on H) that let an orbital **distort off-centre** as it
  forms a bond. Almost always needed for real accuracy.
- **Diffuse functions** (the `+`/`++` in `6-31+G`, or the `aug-` prefix): broad, slowly-decaying functions for
  electron density that reaches far from the nucleus — **anions, lone pairs, excited states, and weak
  intermolecular interactions**. Add them when the density is diffuse; skip them otherwise (they can slow
  convergence).

### The families qc-rs bundles

qc-rs ships **229 named basis sets** as editable Python data files. The three families you will use most:

| Family | Examples | Character |
|---|---|---|
| **Pople** | `3-21g`, `6-31g`, `6-31g*`, `6-31+g**`, `6-311g**` | split-valence, compact, widely tabulated |
| **Dunning (correlation-consistent)** | `cc-pvdz`, `cc-pvtz`, `cc-pvqz`, `aug-cc-pvdz` | **systematically convergent** toward the complete-basis limit; the standard for correlated methods |
| **Karlsruhe def2** | `def2-svp`, `def2-tzvp`, `def2-qzvp` | balanced accuracy/cost; **built-in ECP** for heavy elements (Rb onward) |

List them from Python:

```python
from qc import basis
print(len(basis.baslib_filename))              # 229
print("cc-pvtz" in basis.baslib_filename)      # True
print(sorted(n for n in basis.baslib_filename if n.startswith("def2")))
```

## Usage: choosing a basis with `ao=`

The common case is a single name applied to every atom:

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

mychk = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom")
```

### A different basis per atom

Pass a **dict** to give elements — or individual labelled atoms — their own basis. This is how you put a
large basis only where it matters (e.g. a big basis on the reacting atom, a small one on spectators):

```python
# a bigger basis on oxygen than on the hydrogens
mychk = qc.chk.new(atom=water, ao={"O": "cc-pvtz", "H": "cc-pvdz"}, unit="angstrom")
```

Lookup goes **atom label → base element → element symbol**, so a labelled atom (`C1`) or a ghost (`H-Bq`,
which inherits `H`) resolves as you would expect — see [Molecular input](molecular-input.md).

### Bigger basis → lower energy (and higher cost)

Because Hartree–Fock is **variational** ([Part II](../10-foundations/hartree-fock.md)), enlarging the basis
can only *lower* (or hold) the energy — it never makes it worse. Here is water (RHF, at the geometry above)
across families, with the number of basis functions $K$:

| `ao=` | layer | $K$ (spherical) | $E$(RHF) / $E_h$ |
|---|---|---|---|
| `sto-3g` | minimal | 7 | −74.962947 |
| `6-31g` | split-valence | 13 | −75.983993 |
| `6-31g*` | + polarization | 18 | −76.009128 |
| `cc-pvdz` | double-ζ + pol | 24 | −76.026794 |
| `cc-pvtz` | triple-ζ + pol | 58 | −76.057161 |

The energy drops steadily and $K$ — hence cost — climbs. It converges toward the **Hartree–Fock
complete-basis-set limit**, *not* to the exact energy: the basis controls how well you solve HF, while the
missing **correlation** is a matter of *method* (DFT or post-HF), not basis. A common workhorse compromise is
a polarized double- or triple-ζ set (`cc-pvdz`/`def2-svp` for quick work, `cc-pvtz`/`def2-tzvp` for
production).

## AO representation: spherical vs Cartesian

There is a second, subtler choice: how the angular part of each shell is represented. Set it with
**`ao_rep=`** (default `"spherical"`; aliases `"sph"`/`"cart"`).

A Cartesian $d$ shell has **six** components $\{x^2, y^2, z^2, xy, xz, yz\}$, but one combination
($x^2+y^2+z^2$) is spherically symmetric — really an $s$ function. The **spherical** (pure) representation
drops it, keeping the **five** true $l=2$ real solid harmonics. So spherical shells have $2l+1$ functions
while Cartesian shells have $(l+1)(l+2)/2$ — equal for $s$ and $p$, but $5$ vs $6$ for $d$, $7$ vs $10$ for
$f$, and so on.

```python
sph  = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom", ao_rep="spherical")  # default
cart = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom", ao_rep="cartesian")
# water/cc-pvdz has one d shell on O: spherical K = 24, Cartesian K = 25
```

**Prefer the default spherical representation** — it is the modern standard, has no redundant functions, and
is what production workflows expect. Reach for `ao_rep="cartesian"` only to match a program or reference that
used Cartesian functions. qc-rs stores the choice canonically in the checkpoint (aliases `"sph"`/`"cart"` are
normalized to `"spherical"`/`"cartesian"`).

:::{tip} Inspecting the basis size
The number of basis functions $K$ is the dimension of the overlap matrix. After materializing the integrals,
`np.asarray(mychk.ints().run().overlap).shape` is `(K, K)`. Watch $K$ — runtime scales steeply with it.
:::

## Custom and heavy-element basis sets

- **Custom / edited basis.** Because the data files are plain Python, you can hand a fully custom basis to
  one element using `qc.basis.nwchem_format(...)` inline in the `ao=` dict (shown in
  [Molecular input](molecular-input.md)), or edit/add a bundled `.py` file. The NWChem text format and the
  full editing workflow are in the reference chapter [Bundled basis sets](../40-reference/basis-list.md) and
  `docs/qc-basis.md`.
- **Heavy elements (ECPs).** For heavy atoms, an **effective core potential** replaces the chemically inert
  core electrons, cutting cost and folding in scalar-relativistic effects. The `def2` sets and the
  correlation-consistent `-pp` sets (e.g. `cc-pvdz-pp`) carry their ECP **inline**, so
  `qc.chk.new(ao="def2-svp")` on a heavy atom automatically runs valence-only — no extra argument. The ECP
  machinery is covered where it is used, in the SCF and heavy-element guidance.

## Worked example & exercise

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

# double-zeta for a quick look, triple-zeta for a better energy
for name in ("cc-pvdz", "cc-pvtz"):
    chk = qc.chk.new(atom=water, ao=name, unit="angstrom").scf(ref="r").run()
    print(f"{name:8}  E(RHF) = {chk.scf.energy:.6f}")
# cc-pvdz   E(RHF) = -76.026794
# cc-pvtz   E(RHF) = -76.057161
```

:::{exercise}
:label: ex-basis

1. For water, predict whether `6-31g*` or `6-31+g*` gives the lower RHF energy, and why. Then check it.
2. Water/`cc-pvdz` has $K=24$ spherical basis functions. How many would it have as **Cartesian**? (Hint:
   only the $d$ shell on oxygen differs.)
3. You are studying the hydroxide anion OH⁻. Which *layer* of functions from the anatomy list is most
   important to add, and which basis name provides it?
:::

:::{solution} ex-basis
:class: dropdown

1. `6-31+g*` is lower (or equal): it adds **diffuse** functions (`+`), enlarging the variational space, and
   HF is variational so more functions cannot raise the energy. Diffuse functions especially help where
   density spreads out.
2. **25**. The single $d$ shell on O is $5$ functions spherically but $6$ Cartesian, so $24 - 5 + 6 = 25$
   (matches `ao_rep="cartesian"` above).
3. **Diffuse** functions — an anion's extra electron occupies a spatially broad orbital that compact bases
   describe poorly. Use an augmented set such as `aug-cc-pvdz` (or a Pople `6-31+g*`).
:::

With a molecule and a basis in hand, the next question is where the SCF **starts** — the initial guess. That
is the [next chapter](initial-guess.md).