Basis sets & AO representation

Basis sets & AO representation#

The last chapter 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: 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:

Definition 4 (Contracted Gaussian basis function)

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.

qc-rs stores both the raw and the primitive-normalized contraction coefficients and applies the standard Gaussian normalization for you (qc-mol owns this); you never normalize a basis by hand.

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:

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:

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

# 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.

Bigger basis → lower energy (and higher cost)#

Because Hartree–Fock is variational (Part II), 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.

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), or edit/add a bundled .py file. The NWChem text format and the full editing workflow are in the reference chapter Bundled basis sets 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#

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 2

For water, predict whether 6-31g* or 6-31+g* gives the lower RHF energy, and why. Then check it.
Water/cc-pvdz has \(K=24\) spherical basis functions. How many would it have as Cartesian? (Hint: only the \(d\) shell on oxygen differs.)
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 to Exercise 2

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.
25. The single \(d\) shell on O is \(5\) functions spherically but \(6\) Cartesian, so \(24 - 5 + 6 = 25\) (matches ao_rep="cartesian" above).
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.