The atom does not bring its own orbital

A molecular orbital picture is usually built from atomic orbitals of a stated size. Hydrogen’s 1s has an exponent of one; a screened orbital in a many-electron atom has an exponent from what an electron actually feels, computed from Slater’s rules for a configuration.

Both of those are properties of an atom. The question here is whether they are the right functions for a molecule, and the answer is available exactly for the one system where every integral has a closed form.

Take the simplest molecule there is: two protons and one electron. Build the trial function as the symmetric combination of two 1s functions with exponent $\zeta$, one on each nucleus, and minimise the energy over $\zeta$ as well as over the separation.

Held at hydrogen’s own $\zeta = 1$, the calculation puts the bond at 2.49 bohr and binds by 0.0648 hartree. The true values are 2.00 and 0.1026.

Let the exponent go free and it puts the bond at 2.0033 bohr, at $\zeta = 1.2380$, and binds by 0.0865.

What is being computed, and why it is exact #

There is no approximation in this beyond the choice of trial function, and it is worth setting out the four integrals so that the claim can be checked.

For a 1s function of exponent $\zeta$ on each of two nuclei a distance $R$ apart, with $w = \zeta R$:

$$S = e^{-w}\left(1 + w + \tfrac{w^2}{3}\right)$$

$$j = \frac{1}{R}\left[1 - (1+w)e^{-2w}\right], \qquad k = \zeta(1+w)e^{-w}$$

being the overlap, the attraction of one nucleus for the other’s electron density, and the exchange attraction. The kinetic energy of one such function is $\zeta^2/2$ and its attraction to its own nucleus is $\zeta$. Everything else is algebra, and the total energy of the symmetric combination is

$$E = \frac{H_{aa} + H_{ab}}{1 + S} + \frac{1}{R}$$

$$H_{aa} = \tfrac{\zeta^2}{2} - \zeta - j, \qquad H_{ab} = -\tfrac{\zeta^2 S}{2} + (\zeta-2)k.$$

Not one electron pair appears anywhere in that. There is one electron, so there is no electron–electron repulsion, no exchange between electrons, no self-consistency and nothing to iterate. It is the largest calculation that can be done exactly without becoming a many-electron calculation.

Checked against a quadrature and against the literature #

A closed-form calculation is exactly the sort of thing that can be transcribed wrongly and still produce plausible curves, so the result is checked twice from outside.

Against a numerical quadrature. The overlap can also be computed by a $90^3$ product Gauss rule over a three-dimensional box, and the closed form above is compared against it at three exponents and three separations. The worst disagreement is $7.6\times10^{-5}$. One route is an exponential times a quadratic; the other is a sum over seven hundred and twenty-nine thousand points, and neither knows about the other.

Against the literature. The numbers above have been the answers to this calculation since Finkelstein and Horowitz in 1928. A fixed exponent gives 2.49 bohr and 0.0648 hartree; a free one gives 2.00, $\zeta = 1.238$ and 0.0865. Anyone who has done the problem recognises them, and they are reproduced here to four decimal places.

Why the orbital contracts #

The two terms that fight are visible in the expressions above.

Kinetic energy goes as $\zeta^2$. Squeezing a function into a smaller volume increases the curvature of the wavefunction, and the kinetic energy is the integral of that curvature. Doubling the exponent quadruples it.

Nuclear attraction goes as $\zeta$. Pulling the density closer to a nucleus increases the attraction in proportion to how close, and the expectation of $1/r$ for a 1s function of exponent $\zeta$ is exactly $\zeta$.

For a free hydrogen atom those two balance at $\zeta = 1$, which is the variational statement of the ground state. In a molecule the electron has two nuclei to be attracted to, so the attraction term is larger for the same $\zeta$ — and the balance moves to a larger exponent.

That is the whole mechanism, and it says that the size of an atomic orbital in a molecule is set by how many nuclei are within reach of it. It is not a property the atom carries in.

The dip nobody expects #

The natural expectation is that the best exponent exceeds one at every separation, relaxing to one as the atoms are pulled apart. It fails at six bohr.

What actually happens is more interesting. The best exponent rises to $1.392$ at 1.4 bohr, falls through one at about 5.1, reaches a minimum of 0.9951 near six, and returns to one by twenty. So at long range the molecule’s best single exponent is slightly smaller than the free atom’s — the orbital expands rather than contracting.

The reason is that at long range the electron gains more by reaching towards the distant nucleus than it pays in kinetic energy. A single spherical function centred on one nucleus cannot lean; the only thing it can do to put density in the direction of the other nucleus is to spread out. So a half per cent expansion is the best a one-parameter function can manage at a distance where what is really wanted is polarisation.

That is a small effect and it is worth recording for two reasons. It is a genuine feature rather than a defect of the search — it survives at every resolution — and it is a case where the trial function’s limitation shows up as a wrong-signed variational answer rather than as a bad energy. The energy at six bohr is essentially right; the exponent it chooses is doing something the exponent was not meant to do.

The virial theorem, and the check it makes possible #

There is a third check on a calculation of this kind, and it is the only one here that the calculation cannot arrange to pass by getting the answer right.

The virial theorem for a molecule at its equilibrium separation says that the total kinetic energy is minus the total energy and the potential energy is twice it: $\langle T\rangle = -E$ and $\langle V\rangle = 2E$. Away from equilibrium the statement is not that, it is

$$2\langle T\rangle + \langle V\rangle = -R\,\frac{dE}{dR}$$

and the familiar form is this one at the separation where the derivative vanishes. What makes it a theorem about this calculation rather than a quotation is where it comes from. Scale every length in the trial function by $\zeta$: the kinetic energy becomes $\zeta^2$ times a function of $\zeta R$ and the potential becomes $\zeta$ times another, because kinetic energy has the dimensions of an inverse length squared and a Coulomb potential of an inverse length. Differentiating the total with respect to $\zeta$ at fixed $R$ and setting it to zero gives the identity above. A trial function whose scale is freely optimised satisfies it; one whose scale is held fixed has no reason to.

And nothing on the way to the energy computes either side of it. The four integrals are assembled into $H_{aa}$ and $H_{ab}$, those into one number, and the division of that number between kinetic and potential energy is never formed — so $2\langle T\rangle + \langle V\rangle$ is a quantity the arithmetic could get wrong with the energy unmoved. That is what a check has to be, and it is why this one is worth the algebra.

Which turns out to be four lines rather than the page it looks like. For a 1s function of exponent $\zeta$, $-\tfrac12\nabla^2$ acting on it returns $\left(-\zeta^2/2 + \zeta/r\right)$ times it, so

$$\langle a|T|a\rangle = \tfrac{\zeta^2}{2}, \qquad \langle a|T|b\rangle = -\tfrac{\zeta^2 S}{2} + \zeta k$$

and the potential parts are whatever is left of $H_{aa}$ and $H_{ab}$. The exchange term $\zeta k$ in the second is the one that is easy to lose, and losing it does not look like an error: the residual then comes out at $-0.31$ hartree at the optimum, which reads as a broken theorem rather than as a missing term.

Run, it holds. At the equilibrium separation of 2.0033 bohr and the optimal exponent of 1.2380 the residual is $-8\times10^{-10}$ hartree, and the kinetic energy is $0.5865065014$ against a total energy of $-0.5865065022$ — nine figures, which is the precision of the search rather than of the theorem. Away from that separation the residual is not zero and is not supposed to be: at 1.6, 2.6 and 3.4 bohr it comes out at $+0.1074$, $-0.0711$ and $-0.0959$, and a numerical derivative of the optimised energy curve gives $-R\,dE/dR$ as the same three numbers to nine figures. One of those quantities is read off the wavefunction and the other off the slope of a curve, and they have no reason to agree if either is wrong.

Held at $\zeta = 1$ the theorem fails, and by a chemical amount. At that curve’s own minimum of 2.4928 bohr the residual is $-0.182$ hartree, and the kinetic energy is $0.3827$ where minus the total energy is $0.5648$ — the ratio is 0.68 where the theorem requires 1. That is a completely independent statement of what the extra parameter buys: not merely a lower energy, but a wavefunction with the right balance between kinetic and potential energy, which the fixed one does not have at any separation. A calculation can be handed a better number and still be describing the wrong object, and this is the measurement that tells the two apart.

What one parameter is worth #

The comparison between the fixed and free calculations is worth stating in fractions, because it says something about where effort goes in a real calculation.

Bond length. Fixed: 2.4928 bohr, 24.6 per cent too long. Free: 2.0033, wrong by 0.16 per cent. One parameter took essentially all of the error.

Binding energy. Fixed: 0.0648 hartree, 63 per cent of the exact value. Free: 0.0865, 84 per cent. One parameter took a third of the remaining error.

The asymmetry between those two is exactly the second-order behaviour of a better energy is not a better answer, read from the other side. The bond length is the position of a minimum, so it is a first-derivative property and it responds strongly to an improvement in the wavefunction. The binding energy is the depth of the minimum, a second-order quantity, and it responds less.

So the calculation gets the geometry nearly right while getting the energy substantially wrong, and that is the expected pattern rather than a surprise. It is also the pattern real quantum chemistry shows: geometries converge with basis-set size long before energies do.

What this does to other orbital pictures #

Every molecular orbital drawn from atomic functions at an atomic exponent is, as a consequence, drawn slightly too large.

The size of the effect is knowable for the one case computed: a quarter, at the equilibrium separation of the simplest molecule. For heavier atoms with more shielding the fractional effect is smaller, because a valence orbital is already large and diffuse and the extra nucleus is a smaller perturbation on it — but it is in the same direction everywhere.

Two things follow.

No stated contour changes. Every orbital figure states the fraction of its own density it encloses, and that statement is a property of the function drawn. Scaling the exponent scales the contour radius and leaves the enclosed fraction exactly where it was, so say what it encloses is unaffected and none of the levels in the orbitals index moves.

Overlaps computed at atomic exponents are too large, not too small — which is the opposite of what one would guess and is the sharper half of the result. At 2.0 bohr, $S$ falls from 0.5865 to 0.4641 when $\zeta$ goes from 1 to 1.238, because two smaller orbitals reach into each other less than two large ones do.

So the molecule chooses a configuration with less overlap and a lower energy. That is not a paradox and it is exactly the point overlap is not interaction makes from the other direction: overlap is one factor in a bond and not the quantity being minimised. What the contraction buys is nuclear attraction — each nucleus holds the density near it more tightly — and it pays for that in both kinetic energy and overlap. The balance comes out in favour, by 0.033 hartree.

The comparison at the two equilibrium geometries is the tidiest form of it. At $\zeta = 1$ the molecule sits at 2.4928 bohr with $S = 0.4600$; at $\zeta = 1.238$ it sits at 2.0033 bohr with $S = 0.4631$. Almost the same overlap, a bond a quarter shorter, and a third more binding. Anything that read the overlap as the measure of bond strength would find the two structures indistinguishable.

The exponent as a report on the environment #

Because the exponent is chosen by the molecule, it can be read as a measurement of what the electron’s surroundings are like — and the reading is worth having beside the site’s other measure of the same thing.

What an electron actually feels computes an effective nuclear charge from Slater’s rules: a screened charge that says what the nucleus looks like from where the electron is, after the other electrons have been allowed for. That is a number about a free atom in a stated configuration, and it comes from a fitted set of rules.

The exponent here is the same kind of number computed a completely different way: not from rules, but by minimising an energy, and not for an atom but for a molecule. For hydrogen the free-atom value is one and the molecular value at equilibrium is 1.238, so the electron in $\mathrm{H}_2^+$ behaves as though it were bound to a nucleus of charge 1.238.

That is a physically sensible number. The electron sees two protons, so it might naively be expected to feel something up to two — and it feels 1.238 because it is spread between them and only ever close to one at a time. Halfway between “one proton” and “two protons” would be 1.5; the answer is below that, which says the sharing is not efficient enough to feel like a single doubly charged centre.

The pair of numbers also says something about which route is more general. Slater’s rules apply to any atom and give a number in a second; the variational route applies to any system and needs a calculation. Where both are available they answer slightly different questions — “what does this electron feel in this atom” against “what size of function best describes this electron here” — and the second is the one that changes when a bond forms.

Where this stops #

The boundary is the usual one for exact models, and it is closer here than usual.

One electron is the limit. Adding a second electron to this molecule adds a repulsion integral, and the honest calculation then needs the two-electron integrals over these functions, which are not computed here. The reason for stopping has not changed: a wrong two-electron calculation produces numbers that look right, and the whole method here is arranged against exactly that.

What is not deferred is what a single parameter shows. That an atom’s orbital is not the molecule’s is a statement this calculation makes exactly, in the one system where nothing is approximated, and the size of the effect it gives — a quarter at a bond — is the number every larger calculation has to reproduce.

And the same variational freedom exists in every basis set. A modern calculation does not optimise an exponent; it supplies several functions of different exponents and lets the coefficients do the work. That is the same freedom bought differently, and it is why a minimal basis with atomic exponents gives poor geometries and a double-zeta basis gives good ones. The parameter computed here is what the second zeta is for.

What a modern basis does instead of optimising one number #

The exponent optimised here is a single parameter let loose, and it is worth saying how the same freedom is supplied in practice — because the modern arrangement gets it without optimising anything at the time of the calculation.

Letting a molecule choose the size of its atomic functions means re-optimising a non-linear parameter for every geometry, which is expensive and awkward. The standard alternative is to supply two functions of different size for each valence orbital — one compact and one diffuse — and let the molecule mix them with ordinary linear coefficients, which the eigenvalue problem determines for nothing.

A mixture of a tight function and a loose one behaves very nearly like a single function of intermediate size, and the mixture is adjustable continuously and for free. So a basis with two functions per valence orbital reproduces the effect of an optimised exponent at every geometry, without any exponent ever being varied.

That is exactly what a split-valence basis is for, and it is why the smallest basis anybody recommends for a real calculation has two valence functions rather than one. A minimal basis, with one function per orbital and a fixed exponent, is the case measured here — and its errors are the ones measured here: a bond a quarter too long and a binding a third too weak.

So the 1.238 matters to basis sets as much as to orbitals. The contraction a molecule wants is real and it is not small, and a basis that cannot supply it is a basis that gets the geometry wrong before any question about electron correlation has been asked.

What the free exponent adds #

The usual account of orbitals covers what an orbital is, where the electron is in it, how the shells sit, and what screening does to its size. All of it treats the atomic orbital as the given object.

This essay says the given object is not given. A molecule chooses the size of the functions it is described in, by a variational condition that has nothing to do with the atom, and for hydrogen the choice is a contraction of a quarter that corrects almost the whole error in the bond length.

The many-electron version is the open question. What can be said now is that the one-electron version is exact, that its numbers agree with a calculation from 1928 to four decimals, and that its overlap agrees with a three-dimensional quadrature to eight parts in a hundred thousand.