Part I: Quantum Chemistry | Chapter 4

The Variational Method

The most powerful approximation technique in quantum chemistry — turning unsolvable eigenvalue problems into optimization

1. Introduction: The Variational Principle

The variational principle stands as the single most important approximation method in quantum chemistry. While the Schrodinger equation can be solved exactly only for hydrogen-like atoms and a handful of model systems, the variational method provides a systematic, rigorous route to approximate solutions for atoms, molecules, and solids. Its power lies in a remarkable guarantee: any trial wavefunction, no matter how crude, yields an energy that is an upper bound to the true ground state energy.

The Central Idea

Given a normalized trial wavefunction $|\phi\rangle$ and a Hamiltonian $\hat{H}$ with ground state energy $E_0$, the variational principle states:

$$\boxed{E[\phi] = \frac{\langle \phi | \hat{H} | \phi \rangle}{\langle \phi | \phi \rangle} \geq E_0}$$

The equality holds if and only if $|\phi\rangle$ is the exact ground state. By minimizing $E[\phi]$ over a family of trial functions parameterized by variational parameters, we systematically approach the true ground state energy from above.

This principle transforms quantum mechanics from an eigenvalue problem into an optimization problem. Instead of solving the Schrodinger equation $\hat{H}|\psi\rangle = E|\psi\rangle$ directly, we search through a space of trial functions to find the one that minimizes the energy functional. The method is variational in the sense of the calculus of variations: we seek stationary points of a functional.

The variational method underpins virtually all of computational quantum chemistry: Hartree-Fock theory, density functional theory, configuration interaction, and coupled cluster methods all rely on variational or quasi-variational principles. Even quantum Monte Carlo and machine-learning wavefunctions are, at their core, variational.

2. Proof of the Variational Theorem

Theorem Statement

Let $\hat{H}$ be a Hamiltonian with a complete set of orthonormal eigenstates $\{|n\rangle\}$ and eigenvalues $E_0 \leq E_1 \leq E_2 \leq \cdots$. For any normalizable trial state $|\phi\rangle$:

$$E[\phi] = \frac{\langle \phi | \hat{H} | \phi \rangle}{\langle \phi | \phi \rangle} \geq E_0$$

Step 1: Expand in the Exact Eigenbasis

Since the eigenstates $\{|n\rangle\}$ form a complete basis, any trial function can be expanded:

$$|\phi\rangle = \sum_{n=0}^{\infty} c_n |n\rangle, \quad c_n = \langle n | \phi \rangle$$

The coefficients $c_n$ are complex numbers. We do not need to know them explicitly — the proof works for arbitrary $c_n$.

Step 2: Compute the Expectation Value

Substituting the expansion into the energy functional and using $\hat{H}|n\rangle = E_n|n\rangle$ and orthonormality $\langle m|n\rangle = \delta_{mn}$:

$$\langle \phi | \hat{H} | \phi \rangle = \sum_{m,n} c_m^* c_n \langle m | \hat{H} | n \rangle = \sum_{m,n} c_m^* c_n E_n \delta_{mn} = \sum_{n=0}^{\infty} |c_n|^2 E_n$$

Similarly, the normalization integral is:

$$\langle \phi | \phi \rangle = \sum_{n=0}^{\infty} |c_n|^2$$

Step 3: Establish the Bound

Since $E_n \geq E_0$ for all $n$, we can replace each $E_n$ with $E_0$ to get a lower bound:

$$\langle \phi | \hat{H} | \phi \rangle = \sum_{n=0}^{\infty} |c_n|^2 E_n \geq \sum_{n=0}^{\infty} |c_n|^2 E_0 = E_0 \langle \phi | \phi \rangle$$

Dividing both sides by $\langle \phi | \phi \rangle > 0$:

$$\boxed{E[\phi] = \frac{\langle \phi | \hat{H} | \phi \rangle}{\langle \phi | \phi \rangle} = \frac{\sum_n |c_n|^2 E_n}{\sum_n |c_n|^2} \geq E_0}$$

Key insight: each term $|c_n|^2(E_n - E_0) \geq 0$ is manifestly non-negative. The energy functional is a weighted average of the exact eigenvalues, weighted by $|c_n|^2$. Since the lowest eigenvalue is $E_0$, the weighted average cannot fall below it. Equality holds when $c_n = 0$ for all $n \neq 0$, i.e., when $|\phi\rangle$ is exactly the ground state.

3. Helium Atom Ground State

The helium atom is the simplest system that cannot be solved exactly, making it the perfect testing ground for the variational method. With two electrons orbiting a nucleus of charge $Z = 2$, the Hamiltonian (in atomic units, $\hbar = m_e = e = 1$, $4\pi\epsilon_0 = 1$) is:

$$\hat{H} = -\frac{1}{2}\nabla_1^2 - \frac{1}{2}\nabla_2^2 - \frac{Z}{r_1} - \frac{Z}{r_2} + \frac{1}{r_{12}}$$

The electron-electron repulsion term $1/r_{12}$ makes exact solution impossible. We apply the variational method with a physically motivated trial wavefunction.

Trial Wavefunction: Screened Hydrogenic

We use a product of hydrogen-like 1s orbitals with an effective nuclear charge $Z_{\text{eff}}$ as the variational parameter:

$$\psi(r_1, r_2) = \frac{Z_{\text{eff}}^3}{\pi a_0^3} \, e^{-Z_{\text{eff}} r_1 / a_0} \, e^{-Z_{\text{eff}} r_2 / a_0}$$

The physical idea: each electron partially screens the nuclear charge from the other. Instead of seeing the full $Z = 2$, each electron sees an effective charge $Z_{\text{eff}} < 2$.

Computing the Energy Components

Since the trial function is a product of hydrogenic 1s functions with $Z \to Z_{\text{eff}}$, we can use known results for hydrogen-like atoms. Working in atomic units ($a_0 = 1$):

Kinetic energy (two electrons):

$$\langle T \rangle = 2 \times \frac{Z_{\text{eff}}^2}{2} = Z_{\text{eff}}^2$$

Electron-nucleus attraction (two electrons, actual charge $Z$):

$$\langle V_{en} \rangle = -2Z \cdot Z_{\text{eff}}$$

Note: the wavefunction is parameterized by $Z_{\text{eff}}$ but the potential involves the true nuclear charge $Z$. Using $\langle 1/r \rangle = Z_{\text{eff}}/a_0$ for a hydrogenic 1s orbital gives $\langle -Z/r_i \rangle = -Z \cdot Z_{\text{eff}}$ per electron.

Electron-electron repulsion (the crucial integral):

$$\langle V_{ee} \rangle = \left\langle \frac{1}{r_{12}} \right\rangle = \frac{5}{8} Z_{\text{eff}}$$

This six-dimensional integral can be evaluated by expanding $1/r_{12}$ in Legendre polynomials. The result $5Z_{\text{eff}}/8$ (in atomic units) is one of the classic results of quantum chemistry.

Total Energy and Minimization

Combining all three terms, the variational energy (in Hartrees) is:

$$E(Z_{\text{eff}}) = Z_{\text{eff}}^2 - 2Z \cdot Z_{\text{eff}} + \frac{5}{8} Z_{\text{eff}}$$

Minimizing with respect to $Z_{\text{eff}}$:

$$\frac{dE}{dZ_{\text{eff}}} = 2Z_{\text{eff}} - 2Z + \frac{5}{8} = 0$$

$$\boxed{Z_{\text{eff}} = Z - \frac{5}{16} = 2 - \frac{5}{16} = \frac{27}{16} = 1.6875}$$

Substituting back:

$$E_{\text{min}} = -\left(Z - \frac{5}{16}\right)^2 = -\left(\frac{27}{16}\right)^2 = -2.8477 \text{ Hartree} = -77.5 \text{ eV}$$

The exact ground state energy of helium is $-79.0$ eV. Our single-parameter variational calculation achieves 98.1% accuracy. The effective charge $Z_{\text{eff}} = 27/16 \approx 1.69$ confirms that each electron screens about 5/16 of a unit charge from the other.

4. The Linear Variational Method

The most powerful and widely used form of the variational method expands the trial wavefunction as a linear combination of known basis functions. This transforms the variational problem into a matrix eigenvalue problem that can be solved by standard linear algebra.

Setting Up the Expansion

Choose a set of $N$ linearly independent basis functions $\{\phi_1, \phi_2, \ldots, \phi_N\}$ and expand:

$$|\psi\rangle = \sum_{i=1}^{N} c_i |\phi_i\rangle$$

The coefficients $\{c_i\}$ are the variational parameters. Unlike the nonlinear variational method where we optimize parameters inside the wavefunction, here the optimization is over the linear combination coefficients.

The Energy Functional

The energy functional becomes a ratio of quadratic forms in the coefficients:

$$E[\mathbf{c}] = \frac{\langle \psi | \hat{H} | \psi \rangle}{\langle \psi | \psi \rangle} = \frac{\sum_{i,j} c_i^* c_j H_{ij}}{\sum_{i,j} c_i^* c_j S_{ij}}$$

where the Hamiltonian and overlap matrix elements are:

$$H_{ij} = \langle \phi_i | \hat{H} | \phi_j \rangle, \qquad S_{ij} = \langle \phi_i | \phi_j \rangle$$

Deriving the Generalized Eigenvalue Problem

At the stationary point, $\partial E / \partial c_k^* = 0$ for all $k$. Starting from $E \cdot \mathbf{c}^\dagger \mathbf{S} \mathbf{c} = \mathbf{c}^\dagger \mathbf{H} \mathbf{c}$ and differentiating:

$$\frac{\partial}{\partial c_k^*}\left[\sum_{i,j} c_i^* H_{ij} c_j - E \sum_{i,j} c_i^* S_{ij} c_j\right] = 0$$

$$\sum_{j} H_{kj} c_j - E \sum_{j} S_{kj} c_j = 0, \quad k = 1, 2, \ldots, N$$

This is the generalized eigenvalue equation in matrix form:

$$\boxed{\mathbf{H}\mathbf{c} = E\,\mathbf{S}\mathbf{c}}$$

If the basis is orthonormal ($S_{ij} = \delta_{ij}$), this reduces to the standard eigenvalue problem $\mathbf{H}\mathbf{c} = E\mathbf{c}$. In general, the basis functions are not orthogonal, and the overlap matrix $\mathbf{S}$ is non-trivial.

The Secular Determinant

Non-trivial solutions exist only when the secular determinant vanishes:

$$\boxed{\det(\mathbf{H} - E\,\mathbf{S}) = 0}$$

For a 2-function basis $\{\phi_1, \phi_2\}$, this gives the explicit $2 \times 2$ secular equation:

$$\begin{vmatrix} H_{11} - ES_{11} & H_{12} - ES_{12} \\ H_{21} - ES_{21} & H_{22} - ES_{22} \end{vmatrix} = 0$$

This yields a polynomial of degree $N$ in $E$, giving $N$ approximate eigenvalues. By the Hylleraas-Undheim-MacDonald theorem, the $k$-th eigenvalue provides an upper bound to the exact $k$-th energy level. Increasing the basis size $N$ can only lower (or maintain) each eigenvalue estimate.

5. The Hylleraas Variational Wavefunction

While the simple screened-hydrogenic trial function gives 98% accuracy for helium, achieving chemical accuracy (error $< 1$ kcal/mol $\approx 0.04$ eV) requires a fundamentally different approach. In 1929, Egil Hylleraas made a breakthrough by including the electron-electron distance $r_{12}$ explicitly in the trial wavefunction.

Hylleraas Coordinates

Instead of the independent-particle coordinates $(r_1, r_2)$, Hylleraas introduced:

$$s = r_1 + r_2, \quad t = r_1 - r_2, \quad u = r_{12} = |\mathbf{r}_1 - \mathbf{r}_2|$$

The coordinate $u = r_{12}$ directly encodes the electron-electron distance, allowing the wavefunction to describe electron correlation explicitly. This is the key physical insight: the two-electron cusp condition requires $\psi$ to have a specific linear dependence on $r_{12}$ as $r_{12} \to 0$.

The Trial Function

The Hylleraas trial function takes the form:

$$\psi(s, t, u) = e^{-ks/2} \sum_{l,m,n} c_{lmn} \, s^l \, t^{2m} \, u^n$$

Only even powers of $t$ appear because the spatial wavefunction for the singlet ground state must be symmetric under electron exchange ($t \to -t$). The parameter $k$ is an additional nonlinear variational parameter (related to $Z_{\text{eff}}$).

Systematic Convergence

The convergence with the number of terms is remarkable:

Terms	Form	Energy (eV)	Error (eV)
1	$e^{-ks/2}$	-77.5	1.5
3	$e^{-ks/2}(1 + c_1 u + c_2 t^2)$	-78.7	0.3
6	Through $s^2, t^2, u^2, su, \ldots$	-78.95	0.05
10	Higher powers	-78.999	0.001
Modern (thousands)	Extended Hylleraas + log	-79.0051	$< 10^{-6}$

The inclusion of just the linear $r_{12}$ term reduces the error by a factor of 5. Modern Hylleraas calculations by Drake and others use thousands of terms and achieve agreement with experiment to better than $10^{-6}$ eV — one of the most precise theory-experiment comparisons in all of physics.

Why $r_{12}$ is So Important

The electron-electron cusp condition (due to Kato, 1957) requires that the exact wavefunction satisfy:

$$\left.\frac{\partial \psi}{\partial r_{12}}\right|_{r_{12}=0} = \frac{1}{2}\psi(r_{12}=0)$$

Any trial function built from products of single-electron orbitals (like $e^{-\alpha r_1} e^{-\beta r_2}$) cannot satisfy this condition — it has zero derivative with respect to $r_{12}$ everywhere. The Hylleraas functions, by including $r_{12}$ explicitly, can capture this cusp and the associated correlation energy.

6. Modern Applications

Quantum Monte Carlo

Variational Monte Carlo (VMC) evaluates $E[\phi]$ by stochastic sampling rather than analytic integration. The trial function is parameterized by hundreds or thousands of parameters (Jastrow factors, backflow corrections, neural network wavefunctions). Diffusion Monte Carlo (DMC) then projects out the exact ground state by imaginary time propagation, using the VMC wavefunction as a guide.

Basis Set Optimization

The Gaussian basis sets (STO-3G, 6-31G*, cc-pVDZ, etc.) used in all quantum chemistry codes are optimized variationally. The exponents and contraction coefficients of Gaussian primitives are determined by minimizing atomic energies. Dunning's correlation-consistent basis sets are specifically designed for systematic convergence of the correlation energy.

Excited States

The variational principle extends to excited states through orthogonality constraints. If we restrict the trial function to be orthogonal to the ground state, $\langle \psi_0 | \phi \rangle = 0$, then $E[\phi] \geq E_1$. This is the basis for state-averaged CASSCF, equation-of-motion coupled cluster (EOM-CC), and the Hylleraas-Undheim-MacDonald theorem for the linear variational method.

Neural Network Wavefunctions

Recent breakthroughs use deep neural networks as variational wavefunctions. FermiNet (Pfau et al., 2020) and PauliNet (Hermann et al., 2020) parameterize antisymmetric many-electron wavefunctions with neural networks, achieving chemical accuracy for small molecules. The variational principle provides the loss function for training.

7. Historical Context

1870s

Lord Rayleigh established the variational principle for classical vibrating systems, showing that the frequency of any trial mode shape is always greater than or equal to the fundamental frequency. This classical result directly inspired the quantum mechanical version.

1909

Walther Ritz formalized the variational method for differential equations, introducing the linear expansion approach (the "Rayleigh-Ritz method"). He showed that expanding in a finite basis and solving the resulting secular equation gives systematic upper bounds to eigenvalues.

1928

Kellner applied the screened-charge variational method to helium, obtaining $Z_{\text{eff}} = 27/16$ and $E = -77.5$ eV — the first quantitative success of approximate quantum mechanics for a many-electron system.

1929

Egil Hylleraas introduced the explicitly correlated wavefunction with coordinates $(s, t, u)$. Using a 6-term expansion computed by hand, he achieved an accuracy of 0.01 eV — a stunning triumph that validated quantum mechanics for many-body systems.

1930

Carl Eckart proved rigorous upper and lower bounds and extended the variational method to excited states. He also demonstrated the connection between variational principles and perturbation theory.

1933

Hylleraas & Undheim, MacDonald independently proved that each eigenvalue of the secular equation provides an upper bound to the corresponding exact eigenvalue — the "HUM theorem" that justifies computing excited states from the linear variational method.

8. Computational Laboratory: Python

The Python simulation below performs two variational calculations. First, it optimizes the effective nuclear charge for the helium atom, plotting energy versus $Z_{\text{eff}}$ and finding the analytical minimum at $27/16$. Second, it implements the full linear variational method for a particle in a box using polynomial basis functions, solving the generalized eigenvalue problem $\mathbf{H}\mathbf{c} = E\mathbf{S}\mathbf{c}$ and comparing the resulting wavefunctions to exact solutions.

Variational Method: Helium Atom & Linear Variational Particle-in-a-Box

Python

Optimizes Z_eff for helium ground state energy and implements the linear variational method with polynomial basis for particle-in-a-box. Solves generalized eigenvalue problem Hc = ESc. Uses numpy only.

script.py151 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# Part 1: Variational Optimization for Helium Atom
# Trial wavefunction: psi = exp(-Z_eff * r1 / a0) * exp(-Z_eff * r2 / a0)
# Energy functional: E(Z_eff) = (Z_eff^2 - 2*Z*Z_eff + 5*Z_eff/8) * E_h
# where Z = 2 for helium, E_h = 27.2114 eV (Hartree)
# ============================================================

E_h = 27.2114  # Hartree in eV
Z = 2  # nuclear charge for helium

# Energy as function of Z_eff (in Hartrees)
def E_var(Z_eff):
    # Kinetic: Z_eff^2 (in Hartrees, for two electrons)
    # Nuclear attraction: -2 * Z * Z_eff
    # Electron repulsion: 5 * Z_eff / 8
    return Z_eff**2 - 2 * Z * Z_eff + 5.0 * Z_eff / 8.0

Z_eff_range = np.linspace(1.0, 2.5, 500)
E_range = E_var(Z_eff_range) * E_h  # convert to eV

# Analytical minimum: dE/dZ_eff = 2*Z_eff - 2*Z + 5/8 = 0
Z_opt = Z - 5.0 / 16.0  # = 27/16 = 1.6875
E_opt = E_var(Z_opt) * E_h
E_exact = -79.0  # experimental value in eV

print("=" * 55)
print("  VARIATIONAL METHOD: HELIUM ATOM GROUND STATE")
print("=" * 55)
print(f"  Optimal Z_eff = {Z_opt:.4f} (= 27/16)")
print(f"  Variational energy = {E_opt:.2f} eV")
print(f"  Exact energy       = {E_exact:.1f} eV")
print(f"  Error              = {abs(E_opt - E_exact):.2f} eV ({abs((E_opt - E_exact)/E_exact)*100:.2f}%)")
print(f"  Screening: Z_eff = {Z_opt:.4f} < Z = {Z}")
print(f"  Each electron screens {Z - Z_opt:.4f} of nuclear charge")
print("=" * 55)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
fig.patch.set_facecolor('#0a0a0a')

# Plot 1: Energy vs Z_eff
ax1 = axes[0]
ax1.set_facecolor('#111111')
ax1.plot(Z_eff_range, E_range, color='#34d399', linewidth=2.5, label='E(Z_eff)')
ax1.axhline(y=E_exact, color='#f87171', linestyle='--', linewidth=1.5, label=f'Exact = {E_exact:.1f} eV')
ax1.plot(Z_opt, E_opt, 'o', color='#fbbf24', markersize=10, zorder=5,
         label=f'Min: Z_eff={Z_opt:.4f}, E={E_opt:.1f} eV')
ax1.set_xlabel('Z_eff', color='white', fontsize=12)
ax1.set_ylabel('Energy (eV)', color='white', fontsize=12)
ax1.set_title('Helium Variational Energy', color='white', fontsize=14)
ax1.legend(fontsize=9, facecolor='#1a1a1a', edgecolor='#333', labelcolor='white')
ax1.tick_params(colors='white')
ax1.grid(True, alpha=0.2)
for spine in ax1.spines.values():
    spine.set_color('#333')

# ============================================================
# Part 2: Linear Variational Method - Particle in a Box
# Basis: phi_i(x) = x^i * (L - x) for i = 1, 2, ..., N
# These satisfy boundary conditions phi(0) = phi(L) = 0
# ============================================================

L = 1.0  # box length
N_basis = 5  # number of basis functions

# We work in units where hbar^2/(2m) = 1, so exact energies = n^2 * pi^2 / L^2
exact_energies = [(n * np.pi / L)**2 for n in range(1, N_basis + 1)]

# Numerical integration grid
x = np.linspace(0, L, 2000)
dx = x[1] - x[0]

# Basis functions: phi_i(x) = x^i * (L - x)
def basis(i, x):
    return x**i * (L - x)

def d_basis(i, x):
    # d/dx [x^i * (L - x)] = i * x^(i-1) * (L - x) - x^i
    return i * x**(i - 1) * (L - x) - x**i

# Build H and S matrices
H = np.zeros((N_basis, N_basis))
S = np.zeros((N_basis, N_basis))

for i in range(N_basis):
    for j in range(N_basis):
        phi_i = basis(i + 1, x)
        phi_j = basis(j + 1, x)
        dphi_i = d_basis(i + 1, x)
        dphi_j = d_basis(j + 1, x)
        # H_ij = integral of (d phi_i/dx)(d phi_j/dx) dx  [with hbar^2/2m = 1]
        H[i, j] = np.trapezoid(dphi_i * dphi_j, x)
        # S_ij = integral of phi_i * phi_j dx
        S[i, j] = np.trapezoid(phi_i * phi_j, x)

# Solve generalized eigenvalue problem Hc = ESc
# Transform: S^{-1/2} H S^{-1/2} y = E y, c = S^{-1/2} y
evals_S, evecs_S = np.linalg.eigh(S)
S_inv_half = evecs_S @ np.diag(1.0 / np.sqrt(evals_S)) @ evecs_S.T
H_prime = S_inv_half @ H @ S_inv_half
energies, coeffs_y = np.linalg.eigh(H_prime)
coeffs = S_inv_half @ coeffs_y

print()
print("=" * 55)
print("  LINEAR VARIATIONAL METHOD: PARTICLE IN A BOX")
print("  Basis: phi_i(x) = x^i * (L - x), i = 1..5")
print("=" * 55)
print(f"  {'Level':>6s}  {'Variational':>14s}  {'Exact':>12s}  {'Error':>10s}")
print("-" * 55)
for n in range(min(N_basis, 4)):
    err = abs(energies[n] - exact_energies[n]) / exact_energies[n] * 100
    print(f"  {n+1:>5d}   {energies[n]:>14.6f}  {exact_energies[n]:>12.6f}  {err:>9.4f}%")
print("=" * 55)

# Plot wavefunctions
ax2 = axes[1]
ax2.set_facecolor('#111111')
colors = ['#34d399', '#22d3ee', '#a78bfa', '#fb923c']
for n in range(min(3, N_basis)):
    psi = np.zeros_like(x)
    for i in range(N_basis):
        psi += coeffs[i, n] * basis(i + 1, x)
    # Normalize
    norm = np.sqrt(np.trapezoid(psi**2, x))
    psi /= norm
    if psi[len(x)//2] < 0:
        psi = -psi
    # Exact
    psi_exact = np.sqrt(2.0 / L) * np.sin((n + 1) * np.pi * x / L)
    ax2.plot(x, psi, color=colors[n], linewidth=2.5, label=f'Var. n={n+1}')
    ax2.plot(x, psi_exact, color=colors[n], linewidth=1.5, linestyle='--', alpha=0.6)

ax2.set_xlabel('x / L', color='white', fontsize=12)
ax2.set_ylabel('Wavefunction', color='white', fontsize=12)
ax2.set_title('Linear Variational Wavefunctions', color='white', fontsize=14)
ax2.legend(fontsize=9, facecolor='#1a1a1a', edgecolor='#333', labelcolor='white')
ax2.tick_params(colors='white')
ax2.grid(True, alpha=0.2)
for spine in ax2.spines.values():
    spine.set_color('#333')

plt.tight_layout()
plt.savefig('variational_method.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print()
print("Plot saved: variational_method.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Computational Laboratory: Fortran

The Fortran code below computes variational energies for hydrogen-like atoms ($Z = 1, 2, 3$) using three different trial wavefunctions: the exact exponential form $e^{-\alpha r}$, a Gaussian $e^{-\alpha r^2}$, and a rational function $(1 + \alpha r)^{-2}$. For each, it sweeps the variational parameter $\alpha$ and finds the energy minimum, demonstrating how the choice of functional form affects the variational bound.

Variational Energies for Hydrogen-Like Atoms

Fortran

Computes variational energies for Z=1,2,3 hydrogen-like atoms using exponential, Gaussian, and rational trial wavefunctions. Numerical radial integration with parameter optimization.

variational_hydrogen.f90135 lines

program variational_hydrogen
  implicit none
  ! Variational energies for hydrogen-like atoms
  ! Compare different trial wavefunctions against exact results
  !
  ! Trial 1: psi = exp(-alpha * r)  (exact form)
  ! Trial 2: psi = exp(-alpha * r^2) (Gaussian)
  ! Trial 3: psi = (1 + alpha*r)^(-2) (Lorentzian-like)
  !
  ! Integration by Simpson rule on radial grid

integer, parameter :: N = 10000
  real(8), parameter :: pi = 3.14159265358979323846d0
  real(8), parameter :: r_max = 40.0d0
  real(8) :: dr, r, alpha, E_best, alpha_best
  real(8) :: T_int, V_int, norm_int, psi_val, dpsi_val, E_trial
  real(8) :: exact_E
  integer :: i, j, Z, trial
  real(8) :: alpha_arr(200), E_arr(200)
  integer :: n_alpha
  character(len=30) :: trial_name

dr = r_max / dble(N)

print *, "======================================================"
  print *, "  VARIATIONAL METHOD: HYDROGEN-LIKE ATOMS"
  print *, "  Trial wavefunctions with radial optimization"
  print *, "======================================================"
  print *, ""

do Z = 1, 3
    exact_E = -0.5d0 * dble(Z*Z)
    print '(A,I2,A,F12.6,A)', "  Z = ", Z, "   Exact E = ", exact_E, " Hartree"
    print *, "  ------------------------------------------------"

do trial = 1, 3
      if (trial == 1) trial_name = "Exponential exp(-a*r)"
      if (trial == 2) trial_name = "Gaussian exp(-a*r^2)"
      if (trial == 3) trial_name = "Rational (1+a*r)^-2"

E_best = 1.0d10
      alpha_best = 0.0d0
      n_alpha = 200

do j = 1, n_alpha
        alpha = 0.05d0 + dble(j - 1) * 0.05d0
        alpha_arr(j) = alpha

T_int = 0.0d0
        V_int = 0.0d0
        norm_int = 0.0d0

do i = 1, N
          r = dble(i) * dr
          if (r < 1.0d-12) cycle

if (trial == 1) then
            psi_val = exp(-alpha * r)
            ! -d2/dr2 - (2/r)d/dr acting on R(r)
            ! For psi=e^{-ar}: T = alpha^2/2 per electron (3D radial)
            dpsi_val = alpha * alpha * psi_val - 2.0d0 * alpha * psi_val / r
          else if (trial == 2) then
            psi_val = exp(-alpha * r * r)
            dpsi_val = (4.0d0*alpha*alpha*r*r - 6.0d0*alpha) * psi_val
          else
            psi_val = 1.0d0 / (1.0d0 + alpha * r)**2
            dpsi_val = 6.0d0 * alpha * alpha / (1.0d0 + alpha*r)**4 &
                     - 4.0d0 * alpha / (r * (1.0d0 + alpha*r)**3)
          end if

! Simpson weight
          T_int = T_int + 0.5d0 * psi_val * dpsi_val * r * r * dr
          V_int = V_int - dble(Z) * psi_val * psi_val * r * dr
          norm_int = norm_int + psi_val * psi_val * r * r * dr
        end do

T_int = T_int * 4.0d0 * pi
        V_int = V_int * 4.0d0 * pi
        norm_int = norm_int * 4.0d0 * pi

if (norm_int > 1.0d-20) then
          E_trial = (T_int + V_int) / norm_int
        else
          E_trial = 1.0d10
        end if
        E_arr(j) = E_trial

if (E_trial < E_best) then
          E_best = E_trial
          alpha_best = alpha
        end if
      end do

print '(A,A30,A,F8.4,A,F12.6,A,F8.4,A)', &
        "    ", trial_name, "  a=", alpha_best, &
        "  E=", E_best, "  err=", &
        abs(E_best - exact_E) / abs(exact_E) * 100.0d0, "%"
    end do

print *, ""
  end do

! Print table for Z=1 exponential scan (for plotting)
  print *, "======================================================"
  print *, "  Z=1 Exponential: Energy vs alpha"
  print *, "======================================================"
  print *, "  alpha        Energy(Hartree)"
  print *, "  ---------------------------"
  alpha = 0.0d0
  do j = 1, 40
    alpha = 0.2d0 + dble(j - 1) * 0.05d0
    T_int = 0.0d0
    V_int = 0.0d0
    norm_int = 0.0d0
    do i = 1, N
      r = dble(i) * dr
      if (r < 1.0d-12) cycle
      psi_val = exp(-alpha * r)
      dpsi_val = alpha * alpha * psi_val - 2.0d0 * alpha * psi_val / r
      T_int = T_int + 0.5d0 * psi_val * dpsi_val * r * r * dr
      V_int = V_int - psi_val * psi_val * r * dr
      norm_int = norm_int + psi_val * psi_val * r * r * dr
    end do
    T_int = T_int * 4.0d0 * pi
    V_int = V_int * 4.0d0 * pi
    norm_int = norm_int * 4.0d0 * pi
    if (norm_int > 1.0d-20) then
      E_trial = (T_int + V_int) / norm_int
    else
      E_trial = 0.0d0
    end if
    print '(F8.4, F16.8)', alpha, E_trial
  end do

end program variational_hydrogen

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Key Equations Summary

Variational Principle

$$E[\phi] = \frac{\langle \phi|\hat{H}|\phi\rangle}{\langle \phi|\phi\rangle} \geq E_0$$

Helium Optimal Charge

$$Z_{\text{eff}} = Z - \frac{5}{16} = \frac{27}{16}$$

Secular Determinant

$$\det(\mathbf{H} - E\,\mathbf{S}) = 0$$

Electron Repulsion Integral

$$\left\langle\frac{1}{r_{12}}\right\rangle = \frac{5}{8}Z_{\text{eff}}$$

Kato Cusp Condition

$$\left.\frac{\partial\psi}{\partial r_{12}}\right|_{r_{12}=0} = \frac{1}{2}\psi(0)$$

Generalized Eigenvalue

$$\mathbf{H}\mathbf{c} = E\,\mathbf{S}\mathbf{c}$$

The Physical Picture

1. The variational theorem guarantees that any trial wavefunction gives an upper bound to $E_0$, because the energy functional is a weighted average of exact eigenvalues with non-negative weights $|c_n|^2$.

2. For helium, the screened-charge trial function $e^{-Z_{\text{eff}}(r_1 + r_2)/a_0}$ with $Z_{\text{eff}} = 27/16$ captures electron screening and yields 98% of the exact energy.

3. The electron-electron repulsion integral $\langle 1/r_{12}\rangle = 5Z_{\text{eff}}/8$ is one of the most important integrals in quantum chemistry.

4. The linear variational method reduces quantum mechanics to linear algebra: the secular determinant $\det(\mathbf{H} - E\mathbf{S}) = 0$ gives upper bounds to all eigenvalues simultaneously.

5. Hylleraas's inclusion of $r_{12}$ in the wavefunction was transformative — satisfying the electron-electron cusp condition and converging rapidly to chemical accuracy.

6. From Hartree-Fock to neural network wavefunctions, virtually all of quantum chemistry rests on the variational principle as its foundation.

← Chapter 3: Hartree-Fock Method Chapter 5: Perturbation Theory →