A Gaussian is the wrong shape

Every molecular electronic-structure calculation published in the last sixty years expands its orbitals in Gaussian functions, and the reason is arithmetic rather than physics. The product of two Gaussians centred on different atoms is a Gaussian centred on a third point between them, so the two-electron integrals that dominate the cost of a calculation have closed forms. Nothing else about them is right.

An atomic orbital has a cusp at the nucleus: the wavefunction comes to a corner there, with a slope fixed by the nuclear charge. A Gaussian has a smooth maximum. An atomic orbital dies as $e^{-\zeta r}$; a Gaussian dies as $e^{-\alpha r^2}$, which is faster by an amount that grows without limit.

So the standard method of quantum chemistry is built on functions that are wrong at both ends of the range. This essay is the smallest honest statement of what that costs: one electron, one nucleus, s functions only, every integral in closed form, and an exact answer to check against.

The calculation #

An s-type Gaussian is $g_i(r) = e^{-\alpha_i r^2}$, and the three integrals a one-electron problem needs are elementary:

$$\langle g_i | g_j\rangle = \left(\frac{\pi}{\alpha_i + \alpha_j}\right)^{3/2}, \qquad \langle g_i | -\tfrac12\nabla^2 | g_j\rangle = \frac{3\alpha_i\alpha_j\pi^{3/2}}{(\alpha_i+\alpha_j)^{5/2}}, \qquad \langle g_i | -\tfrac{Z}{r} | g_j\rangle = \frac{-2\pi Z}{\alpha_i + \alpha_j}.$$

Gaussians on one centre are not orthogonal, so what has to be solved is a generalised eigenvalue problem $HC = SCE$. It is turned into an ordinary one by symmetric orthogonalisation — $S^{-1/2}$, computed by the same routine the localisation transformation uses on four bond orbitals, where the non-orthogonality being removed is between hybrids rather than between Gaussians — and handed to an eigensolver checked against closed forms.

The exponents are not quoted. Published basis sets carry values somebody once optimised, and quoting them would put the one thing this essay is about behind a table, so they are searched for here: a deterministic descent on the logarithms of the exponents, starting from a geometric series, moving each in turn. Logarithms because a basis spans orders of magnitude, and a step that suits the tightest function is far too large for the most diffuse.

The energy converges, and it converges from above #

One Gaussian gives $-0.42441$ hartree against an exact $-0.5$: an error of fifteen per cent. Two give $-0.48581$, three $-0.49698$, six $-0.49995$.

Every one of those is above the exact energy, and has to be: the variational principle says the expectation of the Hamiltonian in any trial function is an upper bound on the ground state, so a basis that cannot represent the exact function cannot get below it. Each added function lowers the energy and none of them overshoots, which is the check that the calculation is doing what it claims — an energy below $-0.5$ would mean an error in an integral rather than a good basis.

Those three numbers — $-0.4244$, $-0.4858$, $-0.4970$ — are the ones a chemist will recognise. They are what a minimal basis of one, two and three Gaussians per orbital gives, and the third is why STO-3G was the workhorse of the 1970s: three functions buy 99.4 per cent of the energy of one, at three times the integral count.

What never improves #

The energy is not the wavefunction, and the two places the shape is wrong are wrong at every basis size.

The cusp. Kato’s condition says the exact wavefunction of an electron at a nucleus of charge $Z$ satisfies $\psi'(0)/\psi(0) = -Z$: the kinetic energy has to produce a singularity that cancels the potential’s, and only a function with a corner can. Every Gaussian is smooth at the origin, so every finite sum of them has a derivative of exactly zero there, whatever its exponents and however many of them there are. The error is not small; it is the whole quantity.

The tail. At eight bohr, one Gaussian has $2 \times 10^{-5}$ of the amplitude the exact function has; six have $0.58$ of it. That is a large improvement and it is not convergence: adding functions moves the point where the fit peels away further out, and past that point the decay is still Gaussian rather than exponential. There is no number of Gaussians for which the ratio settles at one.

Why it works anyway #

The obvious reaction to the two paragraphs above is that the method should not work, and the reason it does is worth stating precisely.

Very little of the energy is at either end. The cusp region is small in volume — the $r^2$ weighting that where the electron is turns on means almost none of the probability is within a tenth of a bohr of the nucleus — and the tail is where the density has already fallen by four orders of magnitude. The energy is dominated by the region in between, and a Gaussian sum fits that region well.

Overlap says the same thing more sharply. Three Gaussians reproduce the exact orbital with an overlap of $0.999376$; six give $0.999992$. That is the fair measure of how good a function is, and it is far better than the fifteen-per-cent energy error of one function suggests the family should manage.

And the failures are in quantities nobody usually wants. A property that samples the wavefunction at the nucleus — a hyperfine coupling, an isomer shift — converges badly in a Gaussian basis for exactly the reason above, and the correction is a subject in its own right. A property that samples the tail — a long-range interaction, an anion’s outermost electron, or the overlap at the separations overlap decides is about — needs diffuse functions bolted on, which is what the “+” in a basis-set name means.

So the trade is not a fudge. It is a choice to be wrong where it is cheap to be wrong, made once, and it bought a factor in speed that made molecular calculation possible at all.

The scaling that nothing put in #

One check here is worth more than the others because nothing in the calculation knows about it.

The one-electron problem is exactly scalable in nuclear charge: the Schrödinger equation for charge $Z$ is the one for charge 1 with lengths divided by $Z$ and energies multiplied by $Z^2$. So the optimal exponents for a one-electron ion of charge two should be exactly four times hydrogen’s, and its energy exactly four times as deep.

Optimise a three-Gaussian basis for each independently, from the same starting geometry, and the ratios come out at $4.000$ for the energy and $4.000$ for every exponent. The optimiser is a blind coordinate descent that has never heard of the scaling; the agreement is a check on the integrals, on the eigensolver and on the search all at once.

It also says what an exponent is. An exponent is a length, in disguise: $\sqrt{1/\alpha}$ is a distance, and optimising it is choosing how far out to put a function. A basis set’s exponents are not properties of an element; they are a set of lengths that suited one atom in one environment — the same finding, in a different currency, as the atom does not bring its own orbital, where a hydrogen in a bond wants an exponent of 1.238 rather than the free atom’s 1 — and the whole business of choosing a basis is choosing whether those lengths suit a different one.

What a dependent basis does #

Add a fourth function and there is always somewhere useful to put it. Add a function that duplicates one already there and there is not, and the arithmetic says so loudly.

Three Gaussians with identical exponents span a one-dimensional space. Their overlap matrix is singular, $S^{-1/2}$ does not exist, and the orthogonalisation refuses rather than returning a number — which is the behaviour a near-dependent basis needs, and the reason real calculations report the smallest eigenvalue of $S$ before they report an energy. Three exponents a thousandth apart, on the other hand, are still a basis: they solve, and they give one Gaussian’s energy, because that is what they are worth.

That is the same distinction hybrids are a basis draws between a change of basis and a change of description: a rotation of a basis changes nothing, a shrinking of one changes what can be represented, and the second is the only kind that matters.

What the exponents came out as #

The three optimised exponents for hydrogen are $0.151$, $0.681$ and $4.500$ — a geometric-ish series with a ratio of about four and a half between neighbours, which is what an even-tempered basis assumes and is here an output rather than an input.

The spacing is not arbitrary. Each function is responsible for a range of radii, the ranges have to cover the orbital without overlapping wastefully, and the widths are proportional to the radii, so the exponents fall on a rough geometric progression. That is why basis sets are built the way they are, and why adding a function to an existing set means adding one at an end rather than in the middle.

The tightest exponent, $4.5$, corresponds to a length of about $0.47$ bohr and is doing the work near the nucleus that the cusp needs and cannot get. The most diffuse, $0.151$, corresponds to $2.6$ bohr and is holding up the tail. Between them they cover the region where the electron is locates the density, and neither reaches the end it is nearest.

What a chemist gets out of this in practice #

Three working consequences follow from the two failures, and they are why the trade is worth knowing rather than merely worth noting.

Energies converge faster than anything else, so basis-set convergence tests on an energy are optimistic. A calculation whose energy has settled to a millihartree may still have a dipole moment, a polarisability or a coupling constant moving in the third figure, because those sample parts of the wavefunction the energy barely weights. That is the same relation between an energy and everything else that a better energy is not a better answer sets out in a model where the exact answer is available.

A property at the nucleus needs a special basis or a correction. Hyperfine couplings, Mössbauer isomer shifts and anything else proportional to the density at a nucleus are the cusp failure made visible, and the standard responses — very tight extra functions, or an explicit correction — are both admissions that the shape is wrong there.

And a property far out needs diffuse functions. Anions, weak intermolecular interactions and excited states with diffuse character all fail in a basis whose most diffuse exponent is fitted to a neutral atom’s valence shell. The “aug-” and “+” prefixes are that fix, and the numbers in this essay are why they are needed: a fit optimised for the energy puts its most diffuse function where the density is, not where the tail is.

None of that is a criticism of the method. It is a statement of where the method’s error lives, which is the more useful thing to know.

The cusp is a theorem, so its absence is exactly measurable #

Saying a Gaussian has no cusp understates the situation, because the cusp is not a qualitative feature of the exact solution. It is an equality. For any wavefunction that solves the problem exactly, the spherically averaged density obeys

$$\left.\frac{\partial \bar\rho}{\partial r}\right|_{r=0} = -2Z\,\bar\rho(0)$$

at every nucleus, whatever the atom and whatever else is going on around it. It follows from the potential going as $-Z/r$ and the energy staying finite: the kinetic term has to supply a matching divergence, and only a corner in the wavefunction can.

That turns a shape defect into a number. Every Gaussian is smooth at the origin, so a sum of them has zero slope there, and the ratio a calculation returns for that expression is exactly zero against an exact $-2Z$ — a hundred per cent error, identical at one function and at six, and identical at six hundred. It is the sharpest available statement of what “no number of them fixes it” means: not slowly convergent, not convergent to something slightly wrong, but a quantity on which the sequence never starts.

The consequence is a specific list rather than a general caution, and the list is the properties that live at the nucleus. A hyperfine coupling constant is proportional to the electron density at the nucleus; so is a Mössbauer isomer shift; so is the dominant contribution to nuclear spin–spin coupling. Every one of those is being read off the region the basis is worst in.

The repair is instructive because the energy cannot motivate it. Computing such a property well needs tight functions — Gaussians with very large exponents, packed near the nucleus — and the contractions that published basis sets use have to be undone to let them move. Adding those functions barely changes the total energy at all, so a basis chosen by watching the energy converge will never contain them, and a calculation converged by that criterion will be converged on the wrong thing.

Where the model stops #

One electron. There are no two-electron integrals here, and they are the whole of what makes a many-electron calculation hard. No Gaussian-basis self-consistent field is attempted anywhere here — the smallest many-electron calculation works around it by taking a model with repulsion and no basis at all — on the ground that a wrong one producing plausible numbers is exactly the failure worth avoiding. Nothing in this essay changes that; what it does is say precisely what is missing, which is the two-electron integrals.

s functions only. A real basis has p and d functions with angular parts, and the same integrals acquire more indices and no new principles.

Uncontracted. Published basis sets fix linear combinations of Gaussians and vary only their overall coefficients, which is what “contracted” means and is a further approximation on top of everything above.

No nucleus of finite size, and no relativity. Both matter for a heavy atom and neither is here.

One number from the convergence table deserves separating out, because it is the one that decided the history. Going from one Gaussian to three takes the error from $7.6 \times 10^{-2}$ to $3.0 \times 10^{-3}$ — a factor of twenty-five for a factor of three in the number of functions, and the two-electron integral count goes as the fourth power of the basis size.

So three functions per orbital cost eighty-one times as much as one and buy twenty-five times the accuracy. That is a poor exchange in isolation and an excellent one against the alternative, which was Slater functions whose two-centre integrals had no closed form at all. The whole architecture of computational chemistry follows from that comparison, made once in the 1950s, and every consequence in this essay is a consequence of it.

Why a basis deserves essays of its own #

This is where the question of a basis starts, and it is worth saying what the question is.

The other essays treat orbitals as objects with shapes — nodes, contours, sizes, overlaps. Every one of them is computed from functions written down exactly: hydrogenic radial functions, spherical harmonics, and nothing fitted — which is what lets what an orbital is be careful about the difference between a function and a thing. That is a deliberate austerity and it has an edge, and a basis lies on the other side of it.

A basis is a finite list of functions somebody chose, and the answer depends on the list. The energy converges and the shape does not; the exponents are lengths rather than properties; a nearly dependent basis is worth less than its size suggests; and none of that is visible from inside a calculation that reports only a number. Where the essays on orbitals ask what an orbital is, the essays on bases ask what is being used in place of one, and what the substitution costs.

The questions this one makes possible follow: what a contraction throws away, why a variational energy can be excellent while a property is poor, and what “basis-set superposition” is — an error that arises because two atoms brought together lend each other functions, and which is measured by taking the atoms apart and leaving the functions behind.