Gaussian Trial Functions for Hydrogen

Introduction

When the variation principle is applied to hydrogen with the trial function e^(−ζr), the optimum exponent turns out to reproduce the exact 1s orbital, so the calculation gives the exact energy. That is reassuring but slightly unrealistic: in real problems we almost never guess the exact form. A far more instructive test is to use a function that is deliberately the wrong shape — a Gaussian, e^(−αr²). The answer shows how forgiving the variation principle is, and it explains why nearly every modern quantum chemistry program is built on exactly these "wrong" functions.

Core explanation

The trial function. Take ψ(r) = N e^(−αr²), where α > 0 is the variational parameter and N is the normalisation constant. Working in atomic units (energies in hartrees, Eₕ; lengths in bohr, a₀), the hydrogen Hamiltonian is Ĥ = −½∇² − 1/r.

The two energy terms. Evaluating the integrals with standard Gaussian results gives:

- mean kinetic energy ⟨T⟩ = 3α/2 - mean potential energy ⟨V⟩ = −2√(2α/π)

so the variational energy is

E(α) = 3α/2 − 2√(2α/π)

The competition. A large α makes the Gaussian compact. The electron then sits close to the nucleus, which lowers the potential energy, but it is squeezed into a small region, which raises the kinetic energy (a curvature effect you met with the particle in a box). A small α spreads the function out: low kinetic energy, but weak attraction. The kinetic term grows linearly in α while the potential term grows only as √α, so there is a single best compromise.

Minimising. Setting dE/dα = 0:

3/2 − √(2/π) α^(−1/2) = 0, so α^(1/2) = (2/3)√(2/π) and α = 8/(9π) ≈ 0.283 a₀⁻².

Substituting back, ⟨T⟩ = 4/(3π) and ⟨V⟩ = −8/(3π), so

E min = −4/(3π) Eₕ ≈ −0.424 Eₕ ≈ −11.5 eV.

The exact ground-state energy is −0.500 Eₕ (−13.6 eV). The Gaussian estimate lies above it, as the variation theorem guarantees, and is about 15% too high.

Where the error comes from. The exact 1s function has a cusp at the nucleus — a sharp point with a finite, non-zero slope — whereas every Gaussian has zero slope at r = 0 and is rounded there. At large r, the exact function decays as e^(−r) but a Gaussian decays as e^(−αr²), far too quickly. A single Gaussian therefore puts too little density both very close to the nucleus and far out in the tail.

Why Gaussians are used anyway. The overwhelming advantage is the Gaussian product theorem: the product of two Gaussians centred on different atoms is itself a single Gaussian centred at a point on the line between them. This turns the difficult three- and four-centre electron-repulsion integrals of a molecule into expressions that can be evaluated analytically and very fast. The shape problem is cured by using several Gaussians with different exponents combined into one contracted function.

Formulae

ψ = N e^(−αr²), with N = (2α/π)^(3/4)

E(α) = 3α/2 − 2√(2α/π) (atomic units)

α opt = 8/(9π) ≈ 0.283 a₀⁻²; E min = −4/(3π) Eₕ ≈ −0.424 Eₕ

1 Eₕ ≈ 27.21 eV

Step-by-step reasoning

1. Choose a normalised trial function containing an adjustable exponent α. 2. Evaluate ⟨T⟩ and ⟨V⟩ as functions of α. 3. Add them to obtain E(α). 4. Differentiate, set dE/dα = 0 and solve for α. 5. Substitute the optimum α back into E(α). 6. Compare with the exact energy: the estimate must lie above it.

Visual explanation

Plot the exact 1s function e^(−r) and the optimised Gaussian on the same axes, both normalised. The exponential has a sharp peak at r = 0 and a long tail. The Gaussian is flat-topped at the nucleus, bulges slightly higher at intermediate distances around 1 to 2 a₀, and then falls away steeply beyond about 3 a₀. The areas under the squared curves are equal; only the distribution differs.

Real-world analogy

A Gaussian is like building a pointed church spire out of rounded Lego bricks. One brick is a poor spire, but stacking several bricks of decreasing size gets surprisingly close to the true shape, and the bricks snap together far more easily than carved stone.

Real-world example

The widely used STO-3G basis represents each Slater-type orbital by three Gaussians with fixed coefficients. For hydrogen alone, a single Gaussian gives −0.424 Eₕ, while a good three-Gaussian fit gives an energy within about 1% of the exact −0.500 Eₕ, and larger expansions approach the exact value very closely.

Why?

Why does a badly shaped function still give an energy within 15%? The variation principle makes the energy error second order in the wavefunction error: a function that is wrong by a small amount everywhere gives an energy wrong by a much smaller amount. The optimisation of α also automatically balances kinetic and potential energy.

Common misconception

"A Gaussian is used because it resembles an atomic orbital." It does not: it has no cusp and the wrong tail. Gaussians are chosen for computational convenience — their integrals are easy — and accuracy is recovered by combining many of them.

Worked example

Question: Show that the optimised Gaussian obeys the virial theorem for a Coulomb potential, ⟨V⟩ = −2⟨T⟩.

Reasoning: At α = 8/(9π), ⟨T⟩ = 3α/2 = 4/(3π) ≈ 0.424 Eₕ. ⟨V⟩ = −2√(2α/π) = −2√(16/(9π²)) = −8/(3π) ≈ −0.849 Eₕ.

Answer: −0.849 = −2 × 0.424, so ⟨V⟩ = −2⟨T⟩ and E = −⟨T⟩ = −0.424 Eₕ. Optimising the scale parameter always enforces the virial theorem.

Quick check

1. Why must the optimised Gaussian energy for hydrogen lie above −0.500 Eₕ? Answer: Because the variation theorem states that any trial function gives an energy at or above the true ground-state energy.

Exam focus

Be ready to differentiate E(α), find α opt and state E min with units. Always compare with the exact value and comment on the sign of the error. Explain both the cusp problem and the fast tail decay, then give the Gaussian product theorem as the reason Gaussians are still used.

Advanced insight

Integrals over Gaussians reduce to error functions and Boys functions, which programs compute to machine precision in nanoseconds. A molecule with a few hundred basis functions requires billions of such integrals, so this speed decides what is computable. Some specialised codes use Slater functions directly, but they need numerical quadrature and are much slower for large molecules.

Summary

A Gaussian trial function e^(−αr²) gives E(α) = 3α/2 − 2√(2α/π) for hydrogen. Minimising gives α = 8/(9π) and E = −4/(3π) ≈ −0.424 Eₕ, about 15% above the exact −0.500 Eₕ. The error arises from the missing nuclear cusp and a tail that decays too quickly. Gaussians remain the standard building blocks because the product of two Gaussians is another Gaussian, which makes molecular integrals fast.

Practice questions

1. Give the two physical reasons a single Gaussian is a poor description of the hydrogen 1s orbital. Answer: It has zero slope at the nucleus instead of a cusp, and it decays as e^(−αr²), much faster than the true e^(−r) tail. 2. Convert the optimised Gaussian energy, −0.424 Eₕ, into electronvolts and find the percentage error. Answer: −0.424 × 27.21 ≈ −11.5 eV; the error is (0.500 − 0.424)/0.500 ≈ 15%. 3. If α were chosen larger than the optimum, which energy term would increase in magnitude most quickly? Answer: The kinetic energy, which rises linearly with α, whereas the potential energy magnitude grows only as √α. 4. State the Gaussian product theorem and explain why it matters for molecules. Answer: The product of two Gaussians on different centres is a single Gaussian on an intermediate centre, so multi-centre integrals collapse to simpler ones that can be computed analytically and quickly.