“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Variational Method with Rational Trial Function

Quantum Chemistry

The textbook variational example uses $p s i (r) = e^{- a lp ha r}$ on hydrogen and gets lucky — that shape IS the exact ground state up to a scale, so optimizing $a lp ha$ hits $E_{0} = - 0.5$ hartree on the nose. That's misleading about what the variational principle actually promises. The bound $E [p s i] g e q E_{0}$ holds for any normalizable trial; what changes is how loose it gets.

Where rational forms succeed in practice is as part of the Padé-Jastrow factor in quantum Monte Carlo. The trial wavefunction is $P s i = e^{J} D$ (Slater determinant times Jastrow), and the two-body Jastrow exponent is $u (r_{ij}) = o m e g a_{0}, r_{ij}, /, (1 + o m e g a, r_{ij})$ . The numerator coefficient $o m e g a_{0}$ is fixed, not variational — it's set by the electron-electron Kato cusp condition (0.5 antiparallel spins, 0.25 parallel); only $o m e g a$ is optimized. So the rational shape lives inside the exponent. $e x p (e x t r a t i o na l)$ gets correct exponential decay at large $r$ AND a built-in cusp — both pathologies cancel through the wrapping.

This page strips that shape out of the exponent and uses it as $p s i$ directly: $p s i (r) = 1/ (1 + a r + b r^{2})$ . Now there's nothing to cancel. The decay is algebraic $1/ r^{2}$ instead of exponential, and the function is smooth at the origin instead of cusped. The variational principle still gives an upper bound; the question is how much binding the wrong shape costs. See variational principle for the exponential-trial baseline.

Trial Wavefunction

We use a rational form for the trial wavefunction:

ψ (r) = \frac{1}{1 + a r + b r ^{2}}

where $a$ and $b$ are variational parameters to be optimized.

Energy Functional

For the hydrogen atom, the radial Hamiltonian in atomic units is:

\hat{H} = - \frac{1}{2} (\frac{d ^{2}}{d r ^{2}} + \frac{2}{r} \frac{d}{d r}) - \frac{1}{r}

The variational energy is computed as:

E [a, b] = \frac{⟨ ψ ∣ H ^ ∣ ψ ⟩}{⟨ ψ ∣ ψ ⟩} = \frac{\int _{0}^{\infty} ψ H ^ ψ 4 π r ^{2} d r}{\int _{0}^{\infty} ψ ^{2} 4 π r ^{2} d r}

This requires computing the Laplacian $\nabla^{2} ψ$ , which for an s-wave is:

\nabla^{2} ψ = \frac{d ^{2} ψ}{d r ^{2}} + \frac{2}{r} \frac{d ψ}{d r}

Implementation

The following code implements the variational optimization:

import numpy as np
from scipy.integrate import simps
from scipy.optimize import minimize

# ------------------------------------------------------------
# Trial wavefunction: rational form
# ψ(r) = 1 / (1 + a r + b r^2)
# ------------------------------------------------------------
def psi(r, a, b):
    return 1.0 / (1.0 + a*r + b*r*r)

# ------------------------------------------------------------
# Second derivative needed for kinetic energy
# d²ψ/dr² computed analytically
# ------------------------------------------------------------
def psi_second_derivative(r, a, b):
    denom = 1 + a*r + b*r*r
    d1 = -(a + 2*b*r) / denom**2
    d2 = (2*(a + 2*b*r)**2) / denom**3 - (2*b) / denom**2
    return d2

# ------------------------------------------------------------
# Energy functional E[a,b]
# Hydrogen radial Hamiltonian in atomic units:
# H = -½ (d²/dr² + 2/r d/dr) - 1/r
#
# We work with full 3D wavefunction ψ(r), and the energy integral is:
# E = ∫ ψ H ψ 4π r^2 dr / ∫ ψ^2 4π r^2 dr
#
# ------------------------------------------------------------
def energy(params, r):
    a, b = params
    psi_r = psi(r, a, b)

    # derivatives
    psi_rr = psi_second_derivative(r, a, b)

    # kinetic term T = -1/2 ∫ ψ ∇²ψ dV
    # ∇²ψ (for s-wave) = ψ'' + 2ψ'/r
    # Let's do dψ/dr analytically:
    dpsi_dr = -(a + 2*b*r) / (1 + a*r + b*r*r)**2
    laplacian = psi_rr + 2*dpsi_dr/r

    T_density = -0.5 * psi_r * laplacian

    # potential term V = -1/r
    V_density = -psi_r**2 / r

    # total integrand * 4π r^2
    integrand_T = T_density * 4*np.pi * r*r
    integrand_V = V_density * 4*np.pi * r*r
    integrand_N = psi_r**2 * 4*np.pi * r*r

    T = simps(integrand_T, r)
    V = simps(integrand_V, r)
    N = simps(integrand_N, r)

    return (T + V) / N   # normalized energy

# ------------------------------------------------------------
# Main execution: minimize E[a,b]
# ------------------------------------------------------------
if __name__ == "__main__":
    # grid for r
    r = np.linspace(1e-6, 20.0, 20000)

    # initial guess for parameters (close to expected)
    x0 = np.array([1.0, 0.1])

    print("Optimizing variational parameters...")

    res = minimize(lambda p: energy(p, r), x0,
                   method='Nelder-Mead',
                   options={'maxiter': 200, 'xatol':1e-6, 'fatol':1e-6})

    a_opt, b_opt = res.x
    E_opt = energy(res.x, r)

    print("\nOptimal parameters:")
    print("a =", a_opt)
    print("b =", b_opt)
    print("\nVariational energy:", E_opt)
    print("Exact ground-state energy:", -0.5)

Results

The optimization finds the optimal parameters $a$ and $b$ that minimize the energy. The exact ground-state energy for hydrogen is $E_{0} = - 0.5$ atomic units, and the variational principle guarantees that our result will be an upper bound.