Radial Schrödinger equation for a hydrogenic atom

This notebook is the first part of a two-notebook sequence. Here the radial Schrödinger equation is solved in a prescribed Coulomb potential. We develop the numerical components needed for an atomic calculation: inward integration, shooting, root finding, normalization, orbital occupation, and construction of the electron density.

In the next notebook, 05_Atom_in_LDA_.ipynb, the potential is no longer prescribed. It depends on the electron density through the Hartree and exchange-correlation potentials, so the same radial solver must be placed inside a self-consistency loop:

$$ \rho(r)\longrightarrow V_{\mathrm{eff}}[\rho](r) \longrightarrow \{\varepsilon_{nl},u_{nl}(r)\} \longrightarrow \rho_{\mathrm{new}}(r). $$

The Schrödinger equation is

$$ \left(-\frac{\hbar^2}{2m}\nabla^2- \frac{Z_{\mathrm{eff}}e^2}{4\pi\varepsilon_0r}\right) \psi(\mathbf r)=E\psi(\mathbf r). $$

We write

$$ \psi(\mathbf r)=Y_{lm}(\hat{\mathbf r})\frac{u(r)}{r} $$

and introduce the dimensionless variables

$$ x=\frac{r}{a_0},\qquad \varepsilon=\frac{E}{\mathrm{Ry}}, $$

where

$$ a_0=\frac{4\pi\varepsilon_0\hbar^2}{me^2}\approx0.529\ \text{Å}, \qquad \mathrm{Ry}=\frac{\hbar^2}{2ma_0^2}\approx13.606\ \text{eV}. $$

The radial equation becomes

$$ u''(x)-\left[\frac{l(l+1)}{x^2}-\frac{2Z_{\mathrm{eff}}}{x}-\varepsilon\right]u(x)=0. $$

For hydrogen, $Z_{\mathrm{eff}}=1$. Later we also use this equation as a simple screened starting model for a many-electron calculation. A single $Z_{\mathrm{eff}}$ is not a complete physical description of screening: the true self-consistent potential depends on radius and approaches the full nuclear potential near the nucleus. Its purpose here is to avoid introducing all radial scales of a heavy atom before the LDA calculation.

Writing $y=(u,u')$ gives the first-order system

$$ \frac{dy}{dx}=(u',u''). $$

The bound-state boundary conditions are

$$ u(0)=0,\qquad u(\infty)=0. $$

Because these conditions are imposed at opposite ends of the interval, this is a two-point boundary-value problem. We solve it using a shooting method.

Shooting algorithm

Suppose the boundary conditions are $u(a)=u(b)=0$.

1. Choose a trial energy $\varepsilon$ and an arbitrary nonzero normalization for the derivative at one end.
2. Integrate the differential equation to the other end and evaluate the boundary mismatch.
3. Vary $\varepsilon$ and use a root-finding method to make that mismatch vanish.
4. Continue scanning the energy to find additional bound states.

The overall amplitude is arbitrary because the Schrödinger equation is linear; the wavefunction is normalized after the correct energy has been found.

Numerical remarks

• Outward integration of the radial Schrödinger equation is often unstable because a tiny admixture of the exponentially growing solution eventually dominates. Inward integration suppresses this problem.
• Evaluating the solution on a logarithmic radial mesh resolves the region near the nucleus and the long-distance tail efficiently. odeint remains adaptive internally; the logarithmic array specifies where the solution is returned and subsequently integrated.
• Starting at a finite $R_{\max}$ with $u(R_{\max})=0$ and a small nonzero derivative approximates the boundary condition at infinity. Results must be checked for convergence with respect to $R_{\max}$ and the radial grid.

Plan

1. Use scipy.integrate.odeint to demonstrate inward and outward integration. Distances are measured in Bohr radii and energies in Rydbergs.
2. Define a shooting residual at the origin and use scipy.optimize.brentq to locate bound-state energies.
3. Introduce Numerov's method, which exploits the absence of a first-derivative term in the radial equation.
4. Normalize the radial orbitals and construct a spherically averaged density by filling each $l$ subshell with up to $2(2l+1)$ electrons.
5. Study convergence with respect to the radial mesh and the extrapolation to $r=0$.

We use different names for different physical roles:

• Z_eff: effective charge in the preliminary hydrogenic potential;
• Z_nuclear: physical nuclear charge used in the LDA calculation;
• n_electrons: number of electrons used to fill the orbitals.

Recall:

\begin{eqnarray} && y = (u(r), u'(r) )\\ && dy/dr = (u'(r), u''(r)) \end{eqnarray}

\begin{eqnarray} u''(r)= \left(\frac{l(l+1)}{r^2}-\frac{2Z}{r}-\varepsilon\right)u(r) \end{eqnarray}

import numpy as np
import matplotlib.pyplot as plt

from IPython.display import Image
from numba import njit
from numpy.polynomial import Polynomial
from scipy import integrate, optimize

%matplotlib inline


@njit
def schroedinger_derivative(state, radius, ell, energy, Z_eff=1.0):
    """Return (u', u'') for the radial Schrödinger equation."""
    u, u_prime = state
    coefficient = ell*(ell + 1)/radius**2 - 2*Z_eff/radius - energy
    return np.array([u_prime, coefficient*u])

First use a linear output mesh and integrate outward. For hydrogen, $Z_{\mathrm{eff}}=1$ and the ground-state energy is $E_{1s}=-1$ Ry. The exact reduced radial wavefunction is proportional to $r e^{-r}$.

R = np.linspace(1e-10, 20.0, 500)
ell = 0
energy_1s = -1.0
Z_eff = 1.0

outward_solution = integrate.odeint(schroedinger_derivative,[0.0, 1.0],R,args=(ell, energy_1s, Z_eff))

plt.plot(R, outward_solution[:, 0], label="outward integration")
plt.plot(R, R*np.exp(-R), label="exact shape")
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$, arbitrary normalization")
plt.legend()
plt.show()

Recall the geometrical ideas behind Euler's method and Runge–Kutta methods. The issue below is not that odeint is inaccurate; it is that the unwanted solution of this boundary-value problem grows exponentially during outward integration.

Image(url="https://upload.wikimedia.org/wikipedia/commons/thumb/7/7e/Runge-Kutta_slopes.svg/1920px-Runge-Kutta_slopes.svg.png")

The outward integration eventually becomes dominated by the unwanted growing solution. We therefore integrate inward from a sufficiently large radius.

R = np.linspace(1e-10, 20.0, 500)
ell = 0
energy_1s = -1.0
Z_eff = 1.0

inward_solution = integrate.odeint(schroedinger_derivative, [0.0, -1e-7], R[::-1], args=(ell, energy_1s, Z_eff))
u = inward_solution[:, 0][::-1]
u /= np.sqrt(integrate.simpson(u**2, x=R))

plt.plot(R, u, label="inward integration")
plt.plot(R, 2*R*np.exp(-R), ".", label="exact normalized solution")
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$")
plt.legend()
plt.show()

plt.plot(R, u, "-", label="numerical")
plt.plot(R, 2*R*np.exp(-R), "o", label="exact")
plt.xlim(0.0, 0.05)
plt.ylim(0.0, 0.10)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$")
plt.legend()
plt.show()

Inward integration is stable for this example. Excited orbitals extend to much larger radii, so a logarithmic output mesh can represent both the near-origin region and the long-distance tail efficiently. We first package the calculation into a function.

def solve_schroedinger_ode(energy, ell, R, Z_eff=1.0):
    solution = integrate.odeint(schroedinger_derivative,[0.0, -1e-7],R[::-1],args=(ell, energy, Z_eff))
    u = solution[:, 0][::-1]
    u /= np.sqrt(integrate.simpson(u**2, x=R))
    return u

ell = 1
n = 3
energy = -Z_eff**2/n**2

R_linear = np.linspace(1e-6, 100.0, 1000)
u_linear = solve_schroedinger_ode(energy, ell, R_linear, Z_eff)

plt.plot(R_linear, u_linear, "o-", label="linear output mesh")
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$")
plt.legend()
plt.show()

R_log = np.logspace(-6, 2.0, 100)
u_log = solve_schroedinger_ode(energy, ell, R_log, Z_eff)

plt.plot(R_linear, u_linear/R_linear, "s-", label="linear output mesh")
plt.plot(R_log, u_log/R_log, "o-", label="logarithmic output mesh")
plt.xlim(0.0, 2.0)
plt.ylim(-0.1, 0.4)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)/r$")
plt.legend()
plt.show()

Turning the boundary condition into a root

For a trial energy, inward integration gives a solution that generally fails to have the required regular behavior at the origin. Since

$$u(r)\propto r^{l+1},$$

the quantity $u(r)/r^l$ should approach zero linearly as $r\rightarrow0$. We extrapolate it to the origin and use that value as the shooting residual. A bound-state energy is a zero of this residual.

def shoot_ode(energy, R, ell, Z_eff=1.0):
    u = solve_schroedinger_ode(energy, ell, R, Z_eff)
    regularized_u = u/R**ell
    value_0, value_1 = regularized_u[:2]
    return value_0 + (value_1 - value_0)*(0.0 - R[0])/(R[1] - R[0])

R = np.logspace(-6, 2.2, 500)

shoot_ode(-1.0, R, ell=0, Z_eff=1.0), shoot_ode(-1/3**2, R, ell=1, Z_eff=1.0)

(np.float64(-9.516227259926638e-09), np.float64(4076.9722451440994))

Bracketing all bound states

We evaluate the shooting residual on an energy mesh. Every sign change supplies a bracket for optimize.brentq, which then determines the zero accurately. The scan must be fine enough that no roots are skipped.

def find_bound_states_ode(R, ell, n_states, energy_mesh, Z_eff=1.0):
    """Find at most `n_states` radial bound states for one ell."""
    bound_states = []
    residual_left = shoot_ode(energy_mesh[0], R, ell, Z_eff)

    for energy_left, energy_right in zip(energy_mesh[:-1], energy_mesh[1:]):
        residual_right = shoot_ode(energy_right, R, ell, Z_eff)

        if residual_left*residual_right < 0.0:
            energy = optimize.brentq(shoot_ode,energy_left,energy_right,args=(R, ell, Z_eff),xtol=1e-14)
            bound_states.append((ell, energy))
            principal_n = ell + len(bound_states)
            exact_energy = -Z_eff**2/principal_n**2
            print(f"E = {energy:14.9f}, exact = {exact_energy:14.9f} l = {ell}")
            if len(bound_states) >= n_states:
                break
        residual_left = residual_right

    return bound_states

energy_mesh = -1.2/np.arange(1.0, 20.0, 0.2)**2

ell = 0
shooting_residual = np.array([shoot_ode(energy, R, ell, Z_eff) for energy in energy_mesh])

plt.plot(energy_mesh, shooting_residual, "o-")
plt.axhline(0.0, color="black", linewidth=1)
plt.xlabel("trial energy (Ry)")
plt.ylabel("shooting residual at the origin")
plt.show()

R = np.logspace(-6, 2.2, 500)
n_max = 7

bound_states = []
for ell in range(n_max-1):
    bound_states += find_bound_states_ode(R, ell, n_max - ell, energy_mesh, Z_eff)

E =   -1.000000018, exact =   -1.000000000 l = 0
E =   -0.249999997, exact =   -0.250000000 l = 0
E =   -0.111111111, exact =   -0.111111111 l = 0
E =   -0.062500000, exact =   -0.062500000 l = 0
E =   -0.040000002, exact =   -0.040000000 l = 0
E =   -0.027777773, exact =   -0.027777778 l = 0
E =   -0.020409865, exact =   -0.020408163 l = 0
E =   -0.249999997, exact =   -0.250000000 l = 1
E =   -0.111111112, exact =   -0.111111111 l = 1
E =   -0.062500001, exact =   -0.062500000 l = 1
E =   -0.040000003, exact =   -0.040000000 l = 1
E =   -0.027777889, exact =   -0.027777778 l = 1
E =   -0.020421586, exact =   -0.020408163 l = 1
E =   -0.111111111, exact =   -0.111111111 l = 2
E =   -0.062500001, exact =   -0.062500000 l = 2
E =   -0.040000004, exact =   -0.040000000 l = 2
E =   -0.027778516, exact =   -0.027777778 l = 2
E =   -0.020414271, exact =   -0.020408163 l = 2
E =   -0.062500000, exact =   -0.062500000 l = 3
E =   -0.040000001, exact =   -0.040000000 l = 3
E =   -0.027778126, exact =   -0.027777778 l = 3
E =   -0.020405512, exact =   -0.020408163 l = 3
E =   -0.040000000, exact =   -0.040000000 l = 4
E =   -0.027777799, exact =   -0.027777778 l = 4
E =   -0.020410583, exact =   -0.020408163 l = 4
E =   -0.027777787, exact =   -0.027777778 l = 5
E =   -0.020406833, exact =   -0.020408163 l = 5

bound_states

[(0, -1.0000000175237345),
 (0, -0.24999999672054443),
 (0, -0.11111111124087078),
 (0, -0.06250000042506298),
 (0, -0.040000002241628384),
 (0, -0.027777772745669085),
 (0, -0.020409864997276442),
 (1, -0.24999999669248937),
 (1, -0.11111111161933616),
 (1, -0.06250000145389889),
 (1, -0.04000000294999185),
 (1, -0.02777788871543995),
 (1, -0.02042158629447343),
 (2, -0.11111111126358558),
 (2, -0.06250000053822075),
 (2, -0.04000000419609433),
 (2, -0.02777851578620133),
 (2, -0.020414270711033227),
 (3, -0.06249999980977969),
 (3, -0.040000001233071264),
 (3, -0.02777812596328541),
 (3, -0.020405511668884154),
 (4, -0.03999999995488807),
 (4, -0.027777798500991942),
 (4, -0.02041058292625117),
 (5, -0.02777778742556402),
 (5, -0.020406832575415756)]

In a pure Coulomb potential, the energy depends only on the principal quantum number $n$, so states with different $l$ can be exactly degenerate. Energy sorting alone therefore does not define a unique order within each degenerate hydrogenic shell.

sorted(bound_states, key=lambda state: state[1])

[(0, -1.0000000175237345),
 (0, -0.24999999672054443),
 (1, -0.24999999669248937),
 (1, -0.11111111161933616),
 (2, -0.11111111126358558),
 (0, -0.11111111124087078),
 (1, -0.06250000145389889),
 (2, -0.06250000053822075),
 (0, -0.06250000042506298),
 (3, -0.06249999980977969),
 (2, -0.04000000419609433),
 (1, -0.04000000294999185),
 (0, -0.040000002241628384),
 (3, -0.040000001233071264),
 (4, -0.03999999995488807),
 (2, -0.02777851578620133),
 (3, -0.02777812596328541),
 (1, -0.02777788871543995),
 (4, -0.027777798500991942),
 (5, -0.02777778742556402),
 (0, -0.027777772745669085),
 (1, -0.02042158629447343),
 (2, -0.020414270711033227),
 (4, -0.02041058292625117),
 (0, -0.020409864997276442),
 (5, -0.020406832575415756),
 (3, -0.020405511668884154)]

For constructing a deterministic spherically averaged model density, we use $l$ only as a tie-breaker between numerically degenerate states. This is a bookkeeping convention, not a derivation of the Aufbau ordering of real many-electron atoms. In the LDA notebook, the self-consistent potential lifts the hydrogenic $l$ degeneracy and the orbitals should be filled according to their calculated energies.

def bound_state_sort_key(state):
    ell, energy = state
    principal_n = round(Z_eff/np.sqrt(-energy))
    return principal_n, ell


bound_states = sorted(bound_states, key=bound_state_sort_key)

bound_states

[(0, -1.0000000175237345),
 (0, -0.24999999672054443),
 (1, -0.24999999669248937),
 (0, -0.11111111124087078),
 (1, -0.11111111161933616),
 (2, -0.11111111126358558),
 (0, -0.06250000042506298),
 (1, -0.06250000145389889),
 (2, -0.06250000053822075),
 (3, -0.06249999980977969),
 (0, -0.040000002241628384),
 (1, -0.04000000294999185),
 (2, -0.04000000419609433),
 (3, -0.040000001233071264),
 (4, -0.03999999995488807),
 (0, -0.027777772745669085),
 (1, -0.02777788871543995),
 (2, -0.02777851578620133),
 (3, -0.02777812596328541),
 (4, -0.027777798500991942),
 (5, -0.02777778742556402),
 (0, -0.020409864997276442),
 (1, -0.02042158629447343),
 (2, -0.020414270711033227),
 (3, -0.020405511668884154),
 (4, -0.02041058292625117),
 (5, -0.020406832575415756)]

We can now fill these effective hydrogenic orbitals to construct a model density. The model is useful for testing normalization and occupation counting, but it is not yet a self-consistent many-electron atom.

n_electrons = 28
electron_count = 0
rho = np.zeros(R.size)

for ell, energy in bound_states:
    u = solve_schroedinger_ode(energy, ell, R, Z_eff)
    capacity = 2*(2*ell + 1)
    occupation = min(capacity, n_electrons - electron_count)

    rho += occupation*u**2/(4*np.pi*R**2)
    electron_count += occupation

    print(f"adding (l={ell}, E={energy:14.9f}) with occupation {occupation}/{capacity}; N={electron_count}")

    if electron_count >= n_electrons:
        break

adding (l=0, E=  -1.000000018) with occupation 2/2; N=2
adding (l=0, E=  -0.249999997) with occupation 2/2; N=4
adding (l=1, E=  -0.249999997) with occupation 6/6; N=10
adding (l=0, E=  -0.111111111) with occupation 2/2; N=12
adding (l=1, E=  -0.111111112) with occupation 6/6; N=18
adding (l=2, E=  -0.111111111) with occupation 10/10; N=28

Spherically averaged model density

The quantity $4\pi r^2\rho(r)$ is the radial electron distribution. Its integral must equal n_electrons.

plt.plot(R, 4*np.pi*R**2*rho, label="radial electron distribution")
plt.xlim(0.0, 60.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$4\pi r^2\rho(r)$")
plt.legend()
plt.show()

integrate.simpson(4*np.pi*R**2*rho, x=R)

np.float64(28.000000000000004)

ell = 3
energy = -1/4**2
R_test = np.logspace(-2, 2.0, 1000)
u_test = solve_schroedinger_ode(energy, ell, R_test, Z_eff)

plt.plot(R_test, u_test)
plt.xlim(0.0, 100.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$")
plt.show()

For large $l$, accumulated numerical error prevents accurate inward integration all the way to extremely small radii. Starting farther from the origin can hide the instability, but does not solve the underlying problem. We therefore introduce a method tailored to second-order equations without a first-derivative term.

Numerov algorithm

General-purpose ODE solvers do not exploit the special structure of the radial Schrödinger equation. Numerov's method is designed for equations of the form $x''=f(t)x+q(t)$ and achieves high accuracy using a three-point recurrence.

Numerov's method applies to a second-order linear differential equation of the form

$$ x''(t)=f(t)x(t)+q(t). $$

Expand the solution one step forward and one step backward:

$$ \begin{aligned} x(t+h)&=x(t)+hx'(t)+\frac{h^2}{2}x''(t)+\frac{h^3}{3!}x^{(3)}(t) +\frac{h^4}{4!}x^{(4)}(t)+\cdots,\\ x(t-h)&=x(t)-hx'(t)+\frac{h^2}{2}x''(t)-\frac{h^3}{3!}x^{(3)}(t) +\frac{h^4}{4!}x^{(4)}(t)+\cdots. \end{aligned} $$

Adding these expressions cancels every odd derivative:

$$ x_{i+1}-2x_i+x_{i-1} =h^2x_i''+\frac{h^4}{12}x_i^{(4)}+O(h^6). $$

If the fourth-derivative term were discarded, this would reduce to the ordinary centered second-difference method, with local truncation error $O(h^4)$ and global error $O(h^2)$. Numerov retains an approximation to the fourth-derivative contribution and thereby improves the local truncation error to $O(h^6)$; accumulated over a fixed interval, the usual global error is $O(h^4)$.

Let

$$ g_i=x_i''=f_i x_i+q_i. $$

The fourth derivative of $x$ is the second derivative of $g$. Using its centered approximation,

$$ x_i^{(4)}=g_i'' =\frac{g_{i+1}-2g_i+g_{i-1}}{h^2}+O(h^2), $$

gives the Numerov relation

$$ x_{i+1}-2x_i+x_{i-1} =\frac{h^2}{12}\left(g_{i+1}+10g_i+g_{i-1}\right)+O(h^6). $$

It is convenient to introduce

$$ w_i=x_i-\frac{h^2}{12}g_i =x_i\left(1-\frac{h^2}{12}f_i\right)-\frac{h^2}{12}q_i. $$

The recurrence then becomes

$$ w_{i+1}=2w_i-w_{i-1}+h^2(f_i x_i+q_i), $$

and the physical solution is recovered at each point from

$$ x_i=\frac{w_i+\frac{h^2}{12}q_i} {1-\frac{h^2}{12}f_i}. $$

For the radial Schrödinger equation used below, $q_i=0$.

@njit
def numerov(f, x0, dx, step):
    """Solve x'' = f(t)x on a uniform mesh using Numerov's method."""
    x = np.zeros(f.size)
    x[0] = x0
    x[1] = x0 + step*dx

    step2 = step**2
    step2_over_12 = step2/12.0
    w_previous = x[0]*(1.0 - step2_over_12*f[0])
    w_current = x[1]*(1.0 - step2_over_12*f[1])
    x_current = x[1]
    f_current = f[1]

    for i in range(2, f.size):
        w_next = 2*w_current - w_previous + step2*f_current*x_current
        f_current = f[i]
        x_current = w_next/(1.0 - step2_over_12*f_current)
        x[i] = x_current
        w_previous, w_current = w_current, w_next

    return x

For Numerov algorithm we can evaluate derivative part $f(r)$ for all points at once:

Because Schroedinger Eq is: \begin{eqnarray} u''(r)= \left(\frac{l(l+1)}{r^2}-\frac{2Z}{r}-\varepsilon\right)u(r) \end{eqnarray} the $f$ function is \begin{eqnarray} f(r)= \left(\frac{l(l+1)}{r^2}-\frac{2Z}{r}-\varepsilon\right) \end{eqnarray}

def radial_coefficient(energy, ell, R, Z_eff=1.0):
    return ell*(ell + 1)/R**2 - 2*Z_eff/R - energy

This implementation of Numerov's method requires a uniformly spaced mesh. Numerov can be reformulated after a change of variables, but it cannot be applied unchanged to a logarithmic mesh.

The SciPy tools used in this notebook are integrate.odeint for the initial demonstration, integrate.simpson for normalization and density integrals, and optimize.brentq for bracketed root finding.

R_linear = np.linspace(1e-6, 100.0, 1000)
ell, energy = 0, -1.0
coefficient = radial_coefficient(energy, ell, R_linear[::-1], Z_eff)
u = numerov(coefficient,x0=0.0,dx=1e-7,step=R_linear[1] - R_linear[0])[::-1]
u /= np.sqrt(integrate.simpson(u**2, x=R_linear))

plt.plot(R_linear, u, label="Numerov")
plt.plot(R_linear, 2*R_linear*np.exp(-R_linear), label="exact")
plt.xlim(0.0, 10.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$")
plt.legend()
plt.show()

For this equation and uniform mesh, Numerov is both faster and more accurate than the general-purpose integration used above. Accuracy must still be checked by changing the step size, $R_{\max}$, and the treatment of the origin.

R_linear = np.linspace(1e-6, 100.0, 1000)
ell, energy = 3, -1/4**2
coefficient = radial_coefficient(energy, ell, R_linear[::-1], Z_eff)
u = numerov(coefficient, 0.0, 1e-7, R_linear[1] - R_linear[0])[::-1]
u /= np.sqrt(integrate.simpson(u**2, x=R_linear))

plt.plot(R_linear, u)
plt.xlim(0.0, 100.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$")
plt.grid()
plt.show()

R_linear = np.linspace(1e-6, 100.0, 1000)
ell, energy = 1, -1/3**2
coefficient = radial_coefficient(energy, ell, R_linear[::-1], Z_eff)
u = numerov(coefficient, 0.0, -1e-7, R_linear[0] - R_linear[1])[::-1]
u /= np.sqrt(integrate.simpson(u**2, x=R_linear))

plt.plot(R_linear, u, label=r"$u(r)$")
plt.plot(R_linear, u/R_linear, label=r"$u(r)/r$")
plt.xlabel(r"$r/a_0$")
plt.legend()
plt.show()

Put it all together

def radial_coefficient(energy, ell, R, Z_eff=1.0):
    return ell*(ell + 1)/R**2 - 2*Z_eff/R - energy


def compute_schroedinger(energy, R, ell, Z_eff=1.0):
    coefficient = radial_coefficient(energy, ell, R[::-1], Z_eff)
    u = numerov(
        coefficient,
        x0=0.0,
        dx=-1e-7,
        step=R[0] - R[1],
    )[::-1]
    norm = integrate.simpson(u**2, x=R)
    return u/np.sqrt(abs(norm))


def shoot(energy, R, ell, Z_eff=1.0):
    u = compute_schroedinger(energy, R, ell, Z_eff)/R**ell
    value_0, value_1 = u[:2]
    return value_0 + (value_1 - value_0)*(0.0 - R[0])/(R[1] - R[0])


def find_bound_states(R, ell, n_states, energy_mesh, Z_eff=1.0):
    bound_states = []
    residual_left = shoot(energy_mesh[0], R, ell, Z_eff)

    for energy_left, energy_right in zip(energy_mesh[:-1], energy_mesh[1:]):
        residual_right = shoot(energy_right, R, ell, Z_eff)

        if residual_left*residual_right < 0.0:
            energy = optimize.brentq(
                shoot,
                energy_left,
                energy_right,
                args=(R, ell, Z_eff),
                xtol=1e-14,
            )
            bound_states.append((ell, energy))
            principal_n = ell + len(bound_states)
            exact_energy = -Z_eff**2/principal_n**2
            print(
                f"E = {energy:14.9f}, exact = {exact_energy:14.9f}, "
                f"l = {ell}"
            )

            if len(bound_states) >= n_states:
                break

        residual_left = residual_right

    return bound_states


def bound_state_sort_key(state):
    ell, energy = state
    return energy + ell/1000.0

def charge_density(bound_states, R, n_electrons, Z_eff=1.0):
    rho = np.zeros(R.size)
    electron_count = 0

    for ell, energy in bound_states:
        u = compute_schroedinger(energy, R, ell, Z_eff)
        capacity = 2*(2*ell + 1)
        occupation = min(capacity, n_electrons - electron_count)

        rho += occupation*u**2/(4*np.pi*R**2)
        electron_count += occupation

        print(
            f"adding (l={ell}, E={energy:14.9f}) with "
            f"occupation {occupation}/{capacity}"
        )

        if electron_count >= n_electrons:
            break

    return rho

energy_mesh = -1.2/np.arange(1.0, 20.0, 0.2)**2
R = np.linspace(1e-6, 100.0, 2000)
Z_eff = 1.0
n_electrons = 28
n_max = 6

bound_states = []
for ell in range(n_max - 1):
    bound_states += find_bound_states(
        R, ell, n_max - ell, energy_mesh, Z_eff
    )
bound_states = sorted(bound_states, key=bound_state_sort_key)

rho = charge_density(bound_states, R, n_electrons, Z_eff)
integrated_density = integrate.simpson(4*np.pi*R**2*rho, x=R)
print("Integrated density =", integrated_density)

plt.plot(R, 4*np.pi*R**2*rho, label="radial electron distribution")
plt.xlim(0.0, 80.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$4\pi r^2\rho(r)$")
plt.legend()
plt.show()

E =   -0.999921923, exact =   -1.000000000, l = 0
E =   -0.249990167, exact =   -0.250000000, l = 0
E =   -0.111108194, exact =   -0.111111111, l = 0
E =   -0.062498769, exact =   -0.062500000, l = 0
E =   -0.039999312, exact =   -0.040000000, l = 0
E =   -0.027736583, exact =   -0.027777778, l = 0
E =   -0.250000016, exact =   -0.250000000, l = 1
E =   -0.111111117, exact =   -0.111111111, l = 1
E =   -0.062500003, exact =   -0.062500000, l = 1
E =   -0.039999959, exact =   -0.040000000, l = 1
E =   -0.027743636, exact =   -0.027777778, l = 1
E =   -0.111111111, exact =   -0.111111111, l = 2
E =   -0.062500000, exact =   -0.062500000, l = 2
E =   -0.039999977, exact =   -0.040000000, l = 2
E =   -0.027754397, exact =   -0.027777778, l = 2
E =   -0.062500000, exact =   -0.062500000, l = 3
E =   -0.039999992, exact =   -0.040000000, l = 3
E =   -0.027765321, exact =   -0.027777778, l = 3
E =   -0.039999999, exact =   -0.040000000, l = 4
E =   -0.027773150, exact =   -0.027777778, l = 4
adding (l=0, E=  -0.999921923) with occupation 2/2
adding (l=0, E=  -0.249990167) with occupation 2/2
adding (l=1, E=  -0.250000016) with occupation 6/6
adding (l=0, E=  -0.111108194) with occupation 2/2
adding (l=1, E=  -0.111111117) with occupation 6/6
adding (l=2, E=  -0.111111111) with occupation 10/10
Integrated density = 28.0

The Numerov energy for the $1s$ state still has a noticeable discretization error. Increasing the number of grid points helps, but first we should identify which part of the calculation limits the accuracy.

R = np.linspace(1e-6, 100.0, 1000)
energy_1s_numerical = optimize.brentq(shoot, -1.1, -0.99, args=(R, 0, Z_eff), xtol=1e-14)
print("Energy error:", energy_1s_numerical + 1.0)

Energy error: 0.0005816429635072673

The residual changes sign close to the exact energy, so root finding is not the limiting step. The dominant error comes from representing the near-origin boundary condition on a finite grid.

shoot(-1.0, R, ell=0, Z_eff=Z_eff), shoot(energy_1s_numerical, R, ell=0, Z_eff=Z_eff)

(np.float64(1.8982382736093747e-07), np.float64(-1.2303153558868712e-19))

Inspect the numerical solution close to the origin.

u = compute_schroedinger(-1.0, R, ell=0, Z_eff=Z_eff)
plt.plot(R, u, "o-")
plt.xlim(0.0, 0.5)
plt.ylim(0.0, 1.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$u(r)$")
plt.show()

The first few grid points are relatively far apart on the scale over which the $1s$ orbital bends. A linear extrapolation of $u(r)/r^l$ to the origin therefore introduces an avoidable error. We replace it with a cubic extrapolation through the first four points.

Polynomial.fit(x, y, deg=3) constructs a cubic polynomial from the selected near-origin points. Evaluating that polynomial at zero gives the improved shooting residual.

def shoot_cubic(energy, R, ell, Z_eff=1.0):
    regularized_u = compute_schroedinger(energy, R, ell, Z_eff)/R**ell
    polynomial = Polynomial.fit(R[:4], regularized_u[:4], deg=3)
    return polynomial(0.0)

shoot_cubic(-1.0, R, ell=0, Z_eff=Z_eff), shoot_cubic(energy_1s_numerical, R, ell=0, Z_eff=Z_eff)

(np.float64(1.0576086117985994e-09), np.float64(-1.85605589231308e-07))

R = np.linspace(1e-6, 100.0, 2000)

optimize.brentq(shoot_cubic, -1.1, -0.9, args=(R, 0, Z_eff), xtol=1e-14) + 1.0

5.718114559627452e-08

For this grid, cubic extrapolation reduces the energy error from about $8\times10^{-5}$ to about $6\times10^{-8}$. This shows that the treatment of the origin, rather than the root finder, was limiting the result.

A convergence test should vary the grid spacing systematically. In this example, reducing the spacing by a factor of ten reduces the cubic-extrapolation error by approximately a factor of $10^3$, consistent with an observed error close to $O(h^3)$ over this range. This is an empirical convergence rate for the complete calculation, not the formal local order of Numerov's recurrence.

R = np.linspace(1e-6, 100.0, 2000)

linear_error = optimize.brentq(shoot, -1.1, -0.9, args=(R, 0, Z_eff), xtol=1e-14) + 1.0
cubic_error = optimize.brentq(shoot_cubic, -1.1, -0.9, args=(R, 0, Z_eff), xtol=1e-14) + 1.0

print("Error with linear extrapolation:", linear_error)
print("Error with cubic extrapolation: ", cubic_error)

Error with linear extrapolation: 7.807747000920795e-05
Error with cubic extrapolation:  5.718114559627452e-08

Now repeat the calculation using a ten-times denser mesh.

R = np.linspace(1e-8, 100.0, 20000)

linear_error = optimize.brentq(shoot, -1.1, -0.9, args=(R, 0, Z_eff), xtol=1e-14) + 1.0
cubic_error = optimize.brentq(shoot_cubic, -1.1, -0.9, args=(R, 0, Z_eff), xtol=1e-14) + 1.0

print("Error with linear extrapolation:", linear_error)
print("Error with cubic extrapolation: ", cubic_error)

Error with linear extrapolation: 8.297740849982915e-08
Error with cubic extrapolation:  -1.0893730362226961e-11

Connection to the LDA notebook

The reusable output of this notebook is not a realistic heavy atom; it is the numerical machinery:

• solve a radial eigenvalue equation for a prescribed potential;
• locate eigenvalues with shooting and bracketed root finding;
• normalize radial orbitals;
• fill subshells and construct $\rho(r)$;
• verify normalization and convergence.

In 05_Atom_in_LDA_.ipynb, the physical nuclear charge is denoted by Z_nuclear. The full singular nuclear potential is retained near the origin, while screening is generated by the self-consistent Hartree and exchange-correlation potentials rather than by one constant Z_eff.

Homework

Using the effective hydrogenic model, construct spherically averaged densities containing 28, 92, and 94 electrons. Verify that each density integrates to the requested electron number.