What a contracted Cartesian Gaussian is, and how to normalize one

A molecular orbital is built from atomic orbitals, and in a Gaussian basis each atomic orbital is a fixed sum of Gaussians centered on one atom:

phi(r) = Nc * sum_p cn[p] * x^lx y^ly z^lz * exp(-a_p r^2)

where x, y, z are measured from the atom. The exponents a_p and the raw coefficients come from the basis-set file; they are fixed before the SCF starts and never optimized. All the Gaussians in one sum share the same angular part x^lx y^ly z^lz, and l = lx + ly + lz is the shell's angular momentum. A shell is the set of all components with the same l: one for s, three for p (x, y, z), six Cartesian ones for d (xx, xy, xz, yy, yz, zz), ten for f.

Why Gaussians, when the true atomic orbitals fall off as exp(-r)? Because the product of two Gaussians on different centers is a single Gaussian on a point between them. Every integral the SCF needs then reduces to closed forms and one special function (Boys). A contraction of three to six Gaussians imitates the cusp and the tail well enough.

Normalization. This library works with one convention: every Cartesian component has self-overlap exactly 1. The Gaussian moment

integral x^(2k) exp(-2a x^2) dx = (2k-1)!! / (4a)^k * sqrt(pi / (2a))

fixes the norm of one primitive component as a radial factor times a component factor:

Nr(a, l)          = (2a/pi)^(3/4) (4a)^(l/2)
Nc(lx, ly, lz)    = 1 / sqrt((2lx-1)!! (2ly-1)!! (2lz-1)!!)

Nc cancels the double factorials the moments produce, so the overlap of two primitives of the same shell does not depend on which component they are:

<a|b> = (2 sqrt(ab) / (a+b))^(l + 3/2)

That is why one contracted coefficient per primitive, cn[p], serves every component of a shell. bas_contract multiplies each raw coefficient by Nr(a_p, l) and rescales the whole shell so the contracted self-overlap is

1. The rescale matters because the tabulated coefficients are rounded to 8

figures or so; PySCF does the same renormalization.

PySCF's convention differs for l >= 2. PySCF's Cartesian d and f functions are not unit-normalized: its function equals ours times

F(lx, ly, lz)^2 = 4 pi (2lx-1)!! (2ly-1)!! (2lz-1)!! / (2l+1)!!

and F = 1 for s and p. On the diagonal of PySCF's int1e_ovlp_cart, xx gives 2.513274122871835 and xy gives 0.8377580409572782; the tests below check bas_fc against both.

Component order is PySCF's cart_labels: lx counts down from l, then ly counts down from l - lx, and lz takes the rest.

The packs (atoms, shells, prims) and their column constants are frozen.

Used with import hf.basis.