The variational principle

On the previous page we computed the energy of the H-atom trial wavefunction at five values of $\zeta$, and the $\zeta =1$ row came out lowest at exactly the textbook ground-state energy. That wasn't a coincidence. It was the variational principle, the single most important idea in computational quantum chemistry, doing what it always does. This page makes the principle the program's job — instead of sampling a few values of $\zeta$ by hand, the program finds the optimum.

The principle

For any normalized trial wavefunction ${\psi }_{\text{trial}}$, the expectation value of the Hamiltonian is an upper bound on the exact ground-state energy:

The proof is one line. Expand ${\psi }_{\text{trial}}$ in the eigenstates of $\hat{H}$: ${\psi }_{\text{trial}}={\sum }_n{c}_n{\varphi }_n$, with $\hat{H}{\varphi }_n=E_n{\varphi }_n$ and $E_0\le E_1\le E_2\le \dots$. Then

where the last step uses normalization ${\sum }_n|c_n{|}^2=1$. Equality holds only when $c_n=0$ for every $n>0$, i.e., when the trial wavefunction is the exact ground state. Otherwise the energy is strictly above $E_0$.

The implication is the engine of every method we build. Pick a family of trial wavefunctions parametrized by some knobs. Evaluate the energy for any setting of the knobs and you get an upper bound on the ground-state energy. Find the setting that minimizes the energy and you get the tightest upper bound the family can produce. Make the family richer, you get a tighter bound. The method is correct in the sense that more effort never makes the answer worse.

For our 1s Slater family

The trial family on this page is parametrized by a single number: the orbital exponent $\zeta$. The energy as a function of $\zeta$ is the parabola

Differentiate, set the derivative to zero, solve:

The minimum is at $\zeta =1$ and the energy there is $-\frac{1}{2}$ Hartree, which happens to be the exact ground-state energy of hydrogen. That's because the exact ground state of the H atom is a 1s Slater orbital with $\zeta =1$ — the truth lives inside our trial family. When that happens, the variational minimum doesn't just bound the ground-state energy; it equals it.

This is the trivial case. The non-trivial case is when the trial family doesn't contain the exact ground state. Then the minimum is strictly higher than $E_0$, and the only way to improve is to enrich the family. We'll see exactly that on the next page when we switch from Slater orbitals to a finite-Gaussian basis.

The code (new this page)

On the previous page the program just evaluated $E(\zeta )$ at five hand-picked values. This page hands the minimization to scipy.optimize.minimize_scalar and asks for the optimum.

from scipy.optimize import minimize_scalar


def optimise_hydrogen():
    """Find the zeta that minimizes the trial energy.

    Returns (zeta_opt, E_opt). For the H atom in the 1s Slater family the
    answer is exactly zeta = 1, E = -0.5 Hartree, because the exact ground
    state IS a 1s Slater orbital with zeta = 1.
    """
    result = minimize_scalar(hydrogen_energy, bracket=(0.5, 1.5))
    return result.x, result.fun

scipy.optimize.minimize_scalar takes a callable and a bracket — three points known to enclose the minimum. We hand it hydrogen_energy and a bracket that contains $\zeta =1$. It returns a result object with the optimal $\zeta$ in .x and the minimum energy in .fun. For this one-dimensional convex problem the minimization is trivial; we'll need much heavier machinery once the trial space gets bigger.

Run it

Energy as a function of zeta (the variational principle in action):
  zeta = 0.50    E = -0.375000 Hartree
  zeta = 0.80    E = -0.480000 Hartree
  zeta = 1.00    E = -0.500000 Hartree
  zeta = 1.20    E = -0.480000 Hartree
  zeta = 2.00    E = +0.000000 Hartree

Variational minimum: zeta = 1.000000    E = -0.500000 Hartree
Exact ground state:                     E = -0.500000 Hartree

Two things to confirm. (1) The five-row sweep is unchanged from the previous page — every row is still an upper bound on $E_0=-\frac{1}{2}$. (2) The minimization lands on $\zeta =1.000000$ and $E=-0.500000$, matching the analytic answer. The computer found the optimum the same way we did with calculus.

The full file

"""qc.py — Computational QC from scratch.
Step 2: find the variationally optimal energy by minimizing over zeta.
"""

from scipy.optimize import minimize_scalar


def kinetic_1s(zeta):
    """<1s|T|1s> for a normalized 1s Slater orbital with exponent zeta."""
    return 0.5 * zeta ** 2


def nuclear_attraction_1s(zeta):
    """<1s|V|1s> for the same orbital, V = -1/r centred on the nucleus."""
    return -zeta


def hydrogen_energy(zeta):
    """Total energy E = T + V for the trial wavefunction."""
    return kinetic_1s(zeta) + nuclear_attraction_1s(zeta)


def optimise_hydrogen():
    """Find the zeta that minimizes the trial energy."""
    result = minimize_scalar(hydrogen_energy, bracket=(0.5, 1.5))
    return result.x, result.fun


if __name__ == "__main__":
    print("Energy as a function of zeta (the variational principle in action):")
    for zeta in (0.5, 0.8, 1.0, 1.2, 2.0):
        E = hydrogen_energy(zeta)
        print(f"  zeta = {zeta:>4.2f}    E = {E:>+.6f} Hartree")

    print()

    zeta_opt, E_opt = optimise_hydrogen()
    print(f"Variational minimum: zeta = {zeta_opt:.6f}    E = {E_opt:+.6f} Hartree")
    print(f"Exact ground state:                     E = {-0.5:+.6f} Hartree")