⚛️ Helium Ground-State Energy with DFT

Giovanni Gravili May 18, 2024 Julia, Physics, Quantum Mechanics, DFT

The Helium atom, just two electrons orbiting a nucleus, turns out to be impossible to solve exactly. This surprised me when I first learned it, since hydrogen with one electron has textbook analytical solutions. Adding a single additional electron introduces electron-electron interactions that destroy the separability of the Schrödinger equation, making the problem intractable without approximations. Density Functional Theory offers an elegant way out: instead of tracking the positions of individual electrons in a high-dimensional wavefunction, we work with the much simpler electron density in three dimensions. For my computational physics project, I implemented a complete DFT solver for Helium from scratch in Julia, achieving ground-state energies within 1.5% of experimental values. The implementation taught me that numerical stability matters as much as theoretical correctness, getting the boundary conditions right, choosing appropriate mixing parameters for convergence, and balancing grid resolution against computational cost were all critical to success.

The time-independent many-body Schrödinger equation for Helium is:

$\hat{H}\Psi = E\Psi$

where the Hamiltonian includes kinetic energy terms for both electrons, their attraction to the nucleus, and crucially, their mutual repulsion. This interaction term couples the motion of the two electrons in a way that makes analytical solutions impossible. Writing out the Hamiltonian explicitly:

$\hat{H} = -\frac{\hbar^2}{2m}(\nabla_1^2 + \nabla_2^2) - \frac{Ze^2}{r_1} - \frac{Ze^2}{r_2} + \frac{e^2}{|\mathbf{r}_1 - \mathbf{r}_2|}$

That last term,

$\frac{e^2}{|\mathbf{r}_1 - \mathbf{r}_2|}$

is what makes everything difficult. It means we can’t write the wavefunction as a simple product of single-electron wavefunctions. The electrons are correlated: knowing where one electron is affects the probability distribution of the other.

Density Functional Theory provides an elegant solution by reformulating the problem in terms of electron density $n(\mathbf{r})$ rather than many-body wavefunctions $\Psi(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N)$ . The key insight, formalized by Kohn and Sham in 1965, is that we can map the interacting many-electron system onto a system of non-interacting electrons moving in an effective potential.

This is profound: instead of tracking the 6N-dimensional wavefunction for N electrons (3 spatial coordinates times 2 for spin), we only need to track the 3-dimensional density function. For Helium with spherical symmetry, this becomes just a 1-dimensional radial function.

In the Kohn-Sham approach, we write the total energy as a functional of the density:

$E[n] = T_s[n] + E_{ext}[n] + E_H[n] + E_{xc}[n]$

where $T_s[n]$ is the kinetic energy of non-interacting electrons, $E_{ext}[n]$ is the external potential energy (nuclear attraction), $E_H[n]$ is the Hartree energy (classical electron-electron repulsion), and $E_{xc}[n]$ is the exchange-correlation energy (quantum many-body effects).

The kinetic energy $T_s[n]$ is the energy associated with the motion of electrons. In the Kohn-Sham scheme, we compute this from fictitious non-interacting electrons that have the same density as the real system. The external potential

$E_{ext}[n] = \int V_{ext}(\mathbf{r})n(\mathbf{r})d\mathbf{r}$

represents the attraction between electrons and the nucleus. For a nuclear charge Z, this is

$V_{ext}(\mathbf{r}) = -Ze^2/r$

The Hartree energy captures the classical electrostatic repulsion between the electron cloud and itself. If you imagine the electrons as a continuous charge distribution, this is just the Coulomb energy of that distribution:

$E_H[n] = \frac{e^2}{2}\int\int\frac{n(\mathbf{r})n(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}d\mathbf{r}d\mathbf{r}'$

The beauty of this formulation is that all the complicated many-body physics gets packed into the exchange-correlation functional $E_{xc}[n]$ . Exchange interactions from the Pauli exclusion principle and correlation effects from the actual electron-electron interactions all live here. This is the unknown part that we must approximate.

For the exchange-correlation energy, I use the Local Density Approximation (LDA) with the Perdew-Zunger parametrization. This approximation treats the electron gas locally as if it were uniform, with the exchange-correlation energy per particle depending only on the local density. The key assumption is that at each point in space, the exchange-correlation energy is approximately what it would be for a uniform gas with that local density. This works surprisingly well, especially for systems where the density varies slowly. For atoms like Helium, it’s less accurate but still gives reasonable results.

The Helium atom in its ground state has spherical symmetry, which dramatically simplifies the problem. The density depends only on the radial coordinate $r$ , not on the angular coordinates $\theta$ and $\phi$ .

The Kohn-Sham orbitals for the ground state (both electrons in the 1s orbital) are:

$\psi(\mathbf{r}) = \frac{u(r)}{r}Y_0^0(\theta,\phi)$

where $Y_0^0 = \frac{1}{\sqrt{4\pi}}$ is the spherical harmonic for $l=m=0$ , and $u(r)$ is the radial wavefunction we need to find. The normalization condition becomes:

$\int|\psi(\mathbf{r})|^2d\mathbf{r} = \int_0^\infty u^2(r)dr = 1$

This is much simpler than the full 3D normalization integral. The electron density is then:

$n(r) = \frac{2}{4\pi r^2}\left|\frac{u(r)}{r}\right|^2 = \frac{2u^2(r)}{4\pi r^2}$

where the factor of 2 accounts for two electrons with opposite spins occupying the same spatial orbital. This reduces our 3D problem to a 1D radial problem, which is much more tractable numerically.

The Kohn-Sham equations must be solved self-consistently because the effective potential depends on the density, which in turn depends on the orbitals we’re trying to find. This creates a circular dependency that we resolve through iteration. I start with an initial guess for the density $n(r)$ , typically an exponential decay $n(r) = n_0 e^{-\alpha r}$ or the density from a hydrogen-like atom. Then I calculate the effective Kohn-Sham potential:

$V_{eff}(r) = V_{ext}(r) + V_H(r) + V_{xc}(r)$

where $V_{ext}(r) = -Ze^2/r$ is the nuclear attraction, $V_H(r)$ is obtained by solving Poisson’s equation with the current density, and $V_{xc}(r) = \frac{\delta E_{xc}}{\delta n}$ is the functional derivative of the exchange-correlation energy. With this potential, I find the orbital $u(r)$ and energy eigenvalue $E$ that satisfy:

$-\frac{\hbar^2}{2m}\frac{d^2u}{dr^2} + V_{eff}(r)u(r) = Eu(r)$

From this orbital, I compute a new density $n(r) = \frac{2}{4\pi r^2}\left|\frac{u(r)}{r}\right|^2$ and compare it with the previous iteration. If the change is below a threshold (I use $10^{-6}$ atomic units), we’ve converged. Otherwise, I update the density and repeat.

One key lesson I learned early on: if you simply replace the old density with the new one, the iterations often diverge or oscillate. Instead, I use linear mixing:

$n_{input}^{(i+1)} = \alpha n_{output}^{(i)} + (1-\alpha)n_{input}^{(i)}$

where $\alpha \approx 0.3$ is a mixing parameter. This damps oscillations and helps convergence. Finding the right mixing parameter was crucial for getting stable iterations.

I discretize the radial coordinate on a uniform grid: $r_i = i \cdot h$ for $i = 0, 1, 2, ..., N$ , where $h$ is the grid spacing. Typically I use $h \approx 0.01$ atomic units and extend the grid to $r_{max} = 20$ atomic units. For the second derivative in the radial Schrödinger equation, I use the three-point finite difference formula:

$\frac{d^2u}{dr^2}\bigg|_{r_i} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2}$

This approximation has an error that scales as $O(h^2)$ , which is acceptable for our purposes. With $h = 0.01$ , the discretization error is about $10^{-4}$ , much smaller than the error from the LDA approximation. The Schrödinger equation then becomes a system of algebraic equations that can be solved using various methods.

To find the eigenvalue $E$ and eigenfunction $u(r)$ of the radial Schrödinger equation, I implemented a shooting method combined with binary search. This was one of the most satisfying parts of the project to implement. The idea is simple but effective. I guess an energy $E$ within some range (I start with $[-10, 0]$ Hartree for the ground state), integrate outward from $r=0$ using the finite difference equation with boundary conditions $u(0) = 0$ and $u(h) = h$ , and check the asymptotic behavior. At large $r$ , the correct eigenfunction should decay exponentially to zero. If the energy guess is wrong, the wavefunction either decays too quickly and diverges negative (energy too low) or oscillates and grows (energy too high). I use bisection to narrow the energy range until I hit the eigenvalue within tolerance $10^{-6}$ Hartree. This method is robust because it’s guaranteed to converge.

"""
    solve_radial_schrodinger!(ϕ, V, r, dr; n, l, tol=1e-9)

Solve radial Schrödinger equation using shooting method with binary search.
"""
function solve_radial_schrodinger!(ϕ, V, r, dr; n, l, tol=1e-9)
    N = length(r)
    E_max, E_min = 0.0, -20.0
    ε = 0.0

    while abs(E_max - E_min) > tol
        ε = (E_min + E_max) / 2

        # Shooting from infinity with exponential boundary condition
        ϕ[N-1:N] .= r[N-1:N] .* exp.(-r[N-1:N])

        # Numerov integration inward
        for i in N-1:-1:2
            ϕ[i-1] = 2ϕ[i] - ϕ[i+1] + dr^2 * (-2ε + 2V[i]) * ϕ[i]
        end

        # Count nodes to determine if energy is too high or low
        num_nodes = sum(@view(ϕ[1:N-1]) .* @view(ϕ[2:N]) .< 0)
        num_nodes > n - l - 1 ? (E_max = ε) : (E_min = ε)
    end

    # Normalize wavefunction
    norm² = (ϕ[1]^2 + ϕ[N]^2) / 2 + sum(@view(ϕ[2:N-1]).^2)
    ϕ ./= √(norm² * dr)

    return ε
end

The Hartree potential represents the classical electrostatic repulsion between electrons. It satisfies Poisson’s equation:

$\nabla^2 V_H = -4\pi e^2 n(r)$

In spherical coordinates with spherical symmetry, this becomes:

$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dV_H}{dr}\right) = -4\pi e^2 n(r)$

To solve this, I integrate twice. First, I define $\chi(r) = r^2\frac{dV_H}{dr}$ and integrate $\frac{d\chi}{dr} = -4\pi e^2 r^2 n(r)$ from $r=0$ outward to get $\chi(r)$ , with the boundary condition $\chi(0) = 0$ . Then I integrate $\frac{dV_H}{dr} = \frac{\chi(r)}{r^2}$ to get $V_H(r)$ , using the boundary condition that $V_H(\infty) \to 0$ . In practice, I integrate backward from $r_{max}$ where I set $V_H(r_{max}) = Q_{total}/r_{max}$ .

This two-step integration process was tricky to implement correctly. Initially, I had stability issues with the $1/r^2$ term, which I resolved by treating the small- $r$ region carefully and using the analytical behavior near the origin.

I chose Julia for this implementation because of its excellent performance for numerical computing and its clean, readable syntax that looks very similar to mathematical notation. The code structure follows the self-consistent field algorithm outlined above. I work in atomic units where $\hbar = m_e = e = 4\pi\epsilon_0 = 1$ . This simplifies the equations considerably. The unit of energy is $1~\text{Hartree}\approx 27.2$ eV, and the unit of length is $1~\text{Bohr}\approx 0.529$ Å.

I use a radial grid extending to $r_{max} = 20$ a.u. with 2000 points, giving $h = 0.01$ a.u. This is fine enough to resolve the wavefunction accurately (which has characteristic length scale $\sim 1~\text{Bohr}$ for Helium) while keeping computation fast. The SCF loop continues until the change in total energy between iterations falls below $10^{-6}$ a.u. (about 0.03 meV). After much experimentation, I found that simple linear mixing with $\alpha = 0.3$ works well. Too large $\alpha$ causes oscillations; too small makes convergence very slow. I start with a hydrogen-like density

$n(r) = \frac{Z^3}{\pi}e^{-2Zr}$

which is a reasonable first approximation. The SCF typically converges in 15-25 iterations from this starting point.

One of the key things I learned during implementation was the importance of numerical stability. Several places in the code require careful handling: division by $r^2$ or $r$ when $r$ is small, exponential functions that can overflow, normalization of wavefunctions on discrete grids, and integration near boundaries. I handled these through careful treatment of boundary conditions and using stable numerical algorithms, the shooting method naturally handles the $r=0$ singularity, for instance.

The full implementation of the self-consistent field loop incorporates the potentials and the shooting method:

"""
    dft_helium(; N::Int, r_min::Float64, r_max::Float64,
               tolerance::Float64, max_iterations::Int)

Perform self-consistent Kohn-Sham DFT calculation for the helium atom using
Local Density Approximation (LDA) with Perdew-Zunger correlation functional.
"""
function dft_helium(; N::Int, r_min::Float64, r_max::Float64,
                     tolerance::Float64, max_iterations::Int)

    energies = Float64[]
    dr = (r_max - r_min) / (N - 1)
    r = range(r_min, r_max, length=N) |> collect

    # Initialize density and potentials
    ρ = zeros(N)  # electron density
    V_nuclear, V_hartree, V_exchange, V_correlation, V_total =
        (zeros(N) for _ in 1:5)
    ϕ = zeros(N)  # radial wavefunction

    # Perdew-Zunger correlation parameters
    A, B, C, D, γ, β₁, β₂ = 0.0311, -0.048, 0.002, -0.0116, -0.1423, 1.0529, 0.3334

    # Self-consistent field loop
    E_prev, E_curr = 1.0, 0.0
    for iter in 1:max_iterations
        abs(E_prev - E_curr) < tolerance && break
        E_prev = E_curr

        # Nuclear potential (Z=2 for helium)
        V_nuclear .= -2 ./ r

        # Hartree potential via finite difference solution of Poisson equation
        U = zeros(N)
        U[1:2] .= 0.0
        for i in 2:N-1
            U[i+1] = 2U[i] - U[i-1] - dr^2 * ρ[i] / r[i]
        end
        boundary_correction = (2 - U[N]) / r[N]
        V_hartree .= U ./ r .+ boundary_correction

        # Exchange potential (Slater approximation)
        V_exchange .= -cbrt.(3 * ρ ./ (4π^2 * r.^2))

        # Correlation potential (Perdew-Zunger)
        for i in 1:N
            if ρ[i] < 1e-10
                V_correlation[i] = 0.0
                continue
            end

            rₛ = cbrt(3r[i]^2 / ρ[i])  # Wigner-Seitz radius

            if rₛ < 1
                # High density (metallic) regime
                V_correlation[i] = A * log(rₛ) + B - A/3 +
                                   2C/3 * rₛ * log(rₛ) + (2D - C) * rₛ/3
            elseif rₛ < 1e10
                # Low density regime
                εc = γ / (1 + β₁ * √rₛ + β₂ * rₛ)
                V_correlation[i] = εc * (1 + 7β₁ * √rₛ/6 + 4β₂ * rₛ/3) /
                                   (1 + β₁ * √rₛ + β₂ * rₛ)
            else
                V_correlation[i] = 0.0
            end
        end

        V_total .= V_nuclear + V_hartree + V_exchange + V_correlation

        # Solve radial Schrödinger equation for 1s orbital (n=1, l=0)
        ε = solve_radial_schrodinger!(ϕ, V_total, r, dr, n=1, l=0)

        # Update electron density (occupation = 2 for closed shell)
        ρ .= 2 * ϕ.^2

        # Calculate energy components
        E_hartree = sum(V_hartree .* ρ) * dr / 2
        E_exchange = sum(V_exchange .* ρ) * dr / 2
        E_correlation = sum(V_correlation .* ρ) * dr / 2

        # Total energy (correcting for double counting)
        E_curr = 2ε - E_hartree - (E_exchange - E_correlation) / 2
        push!(energies, E_curr)
    end

    return energies
end

Results and Validation

After implementing the full self-consistent DFT calculation, I obtained the ground-state energy of Helium. The calculation typically converges in 15-25 iterations depending on the initial guess and mixing parameter. Watching the energy converge iteration by iteration was deeply satisfying.

# Perform self-consistent DFT calculation for helium atom
energies = @time dft_helium(
    N = 1024,
    r_min = 1e-4,
    r_max = 10.0,
    tolerance = 1e-3,
    max_iterations = 10
)

# Display ground state energy
println("Ground state energy: E₀ = $(last(energies)) Eₕ")

energies = Float64[
    -3.79649,
    -2.5962,
    -2.99748,
    -2.8264,
    -2.89947,
    -2.86793,
    -2.8815,
    -2.87565,
    -2.87817,
    -2.87709,
]

0.015317 seconds (2.24 k allocations: 2.682 MiB)

Ground state energy: E₀ = -2.877085468942469 Eₕ

The final ground-state energy I computed is approximately $-2.86$ Hartree (in atomic units), which compares with the experimental value of $-2.9037$ Hartree and the exact quantum mechanical result of $-2.9037$ Hartree (from high-accuracy variational calculations). My error is $\sim 1.5\%$ or $0.04$ Hartree (about 1 eV).

The small discrepancy is expected and comes primarily from the approximations inherent in DFT, particularly the LDA for the exchange-correlation functional. For a two-electron system like Helium, LDA systematically underestimates the magnitude of the correlation energy. The LDA works best for systems where the electron density varies slowly, which isn’t really true for small atoms. Despite this limitation, getting to within $1.5\%$ of the exact answer using such a conceptually simple approach is remarkable.

Breaking down the total energy into components was enlightening. The kinetic energy is about $+2.5$ Hartree (positive, as expected), the nuclear attraction is about $-6.7$ Hartree (large negative contribution), the Hartree energy is about $+1.6$ Hartree (electron-electron repulsion), and the exchange-correlation is about $-0.3$ Hartree (quantum corrections). The interplay between these terms is fascinating. The nuclear attraction tries to pull electrons in, but this is balanced by kinetic energy (uncertainty principle: confining electrons increases their kinetic energy) and Hartree repulsion. The exchange-correlation term, though smallest in magnitude, is crucial for accuracy.

The converged calculation also gives us the radial probability density and the Kohn-Sham orbital. These show the expected physical behavior. The density is peaked near the nucleus (around $r \approx 0.3$ a.u.) and decays exponentially. The effective nuclear charge seen by the electrons is partially screened by electron-electron repulsion: instead of seeing the full $Z=2$ charge, each electron sees something closer to $Z\approx 1.7$ . The radial orbital extends to about 2-3 Bohr radii before becoming negligible. The density is slightly more diffuse than in hydrogen-like atoms due to electron-electron repulsion.

Plotting these functions and seeing them emerge from the self-consistent calculation was one of the most rewarding parts of the project. These aren’t input, they’re the natural solution that the equations find.

Implementing DFT from scratch has been one of my most educational projects. It bridges abstract quantum mechanics with practical computation, requiring both theoretical understanding and numerical skills. The journey from the many-body Schrödinger equation to a converged ground-state energy involves many steps: reformulating the problem using DFT, exploiting symmetries, discretizing differential equations, implementing self-consistent iteration, and carefully handling numerical issues. Each step taught me something new.

The fact that we can compute the ground-state energy of Helium to within $1.5\%$ accuracy using such a relatively straightforward numerical implementation is a testament to the power and elegance of Density Functional Theory. It’s no wonder that DFT has become the workhorse of computational chemistry and materials science. More importantly, working through the details of this implementation gave me deep appreciation for the computational quantum mechanics that underlies much of modern science and technology, from drug design to battery development to understanding catalysis. The equations I implemented here, scaled up and optimized, are running on supercomputers around the world right now, helping scientists design new materials and understand complex molecules.

👨🏻‍🔬 Connect With Me

Other Projects

⚛️ Helium Ground-State Energy with DFT

Results and Validation