The setup

Before any of the heavy machinery — Slater determinants, Fock matrices, integrals over Gaussians — we need a foothold. A self-contained Python program that does the simplest meaningful quantum-chemical calculation: evaluate the energy of a single trial wavefunction for the only atom whose Schrödinger equation we can solve in closed form. The hydrogen atom.

This is the file that grows. By the end of the series it will be a working Hartree-Fock implementation for H₂. Today it's twenty lines and doesn't even know what an electron-electron interaction is. That's fine. We'll add things one at a time.

The problem

The time-independent Schrödinger equation for an electron in the field of a nucleus of charge $Z=1$ is $$\hat{H}\psi (\mathbf{r})=E\psi (\mathbf{r}),\hat{H}=-\frac{1}{2}{\nabla }^2-\frac{1}{r}.$$ In atomic units (so $\hslash =m_e=e=1$, lengths in Bohr, energies in Hartree) the exact ground-state solution is the 1s orbital $${\psi }_{1s}(\mathbf{r})=\frac{1}{\sqrt{\pi }}{e}^{-r},E_0=-\frac{1}{2}\text{Hartree}.$$ Two known facts. We're going to use them as a sanity check.

The plan is to guess a wavefunction, evaluate its energy, and see how close we are to the exact answer. The guess we'll use is a 1s Slater-type orbital (STO) with an exponent we can vary: $${\varphi }_{\zeta }(\mathbf{r})=\sqrt{\frac{{\zeta }^3}{\pi }}{e}^{-\zeta r}.$$ For $\zeta =1$ this is the exact ground state. For any other value it's wrong — so its energy will be higher. That fact is the variational principle, and it's what the next page is about. For now we just want a function that evaluates the energy.

The energy expectation value

Energy is $⟨\hat{H}⟩=⟨T⟩+⟨V⟩$. For our normalized 1s STO the kinetic and nuclear-attraction expectation values come out in closed form: $$⟨T⟩=\frac{1}{2}{\zeta }^2,⟨V⟩=-\zeta .$$ These are the only integrals we'll need on this page. (See the existing 1s integrals page if you want the derivation.) Both are exact functions of $\zeta$; no numerical quadrature, no basis-set fitting. Just two formulas.

The total energy is $E(\zeta )=\frac{1}{2}{\zeta }^2-\zeta$, a parabola in $\zeta$. Its minimum is at $\zeta =1$, where $E=-\frac{1}{2}$. So if our program is doing the right thing, plugging in $\zeta =1$ should give back the exact ground-state energy. That's our anchor.

The code (new this page)

Three small functions: kinetic energy, nuclear attraction, total energy. Each takes zeta and returns a number. No data structures yet, no basis sets, no SCF. Just three formulas.

def kinetic_1s(zeta):
    return 0.5 * zeta ** 2


def nuclear_attraction_1s(zeta):
    return -zeta


def hydrogen_energy(zeta):
    return kinetic_1s(zeta) + nuclear_attraction_1s(zeta)

Both kinetic_1s and nuclear_attraction_1s are direct transcriptions of the closed-form expectation values. hydrogen_energy sums them. That's the entire computation for the H atom — at this level of approximation, in this basis. It fits in three lines because hydrogen is special. Helium will be much harder.

Run it

Save the file below as qc.py and run python qc.py. You should see:

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

Three things to notice. (1) The $\zeta =1$ row gives exactly $-0.5$, the textbook hydrogen ground-state energy in Hartree. (2) Every other row is higher (less negative). The trial wavefunction is doing worse than the truth, and the energy is reporting that honestly. (3) $\zeta =2$ gives 0: kinetic and potential exactly cancel. That's not the right answer. A good wavefunction has more negative energy than this one.

The fact that the $\zeta =1$ row is the lowest of the five — and that all five are upper bounds on the true ground-state energy — is the variational principle in action. It's the single most important idea in computational quantum chemistry, and the next page makes it the entire job of the program: don't just evaluate $E(\zeta )$, minimize it.

The full file

Here's qc.py at the end of step 1 — the canonical artifact for this page. Subsequent pages will add to it; nothing on this page will be removed.

"""qc.py — Computational QC from scratch.
Step 1: evaluate the energy of a 1s Slater-orbital trial wavefunction
for the hydrogen atom.
"""


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

    Closed form: T = (1/2) * zeta**2.
    """
    return 0.5 * zeta ** 2


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

    Closed form: V = -zeta.
    """
    return -zeta


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


if __name__ == "__main__":
    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")