A self-consistent atom in the local-density approximation¶
This notebook continues 04_Scipy_Hydrogen_atom_with_solution.ipynb. There, a prescribed Coulomb potential produced radial orbitals and a model density. Here the potential depends on the density itself:
$$ \rho(r)\longrightarrow V_{\mathrm{eff}}[\rho](r) \longrightarrow \{\varepsilon_{nl},u_{nl}(r)\} \longrightarrow \rho_{\mathrm{new}}(r). $$
We repeat this cycle until the input and output densities agree. The calculation is tested first for oxygen, and its total energy is compared with the corresponding non-spin-polarized LDA reference value from NIST.
Kohn–Sham equation and units¶
For a spherical atom with nuclear charge $Z_{\mathrm{nuclear}}$, the Kohn–Sham equation is
$$ \left[-\frac{\hbar^2}{2m}\nabla^2- \frac{Z_{\mathrm{nuclear}}e^2}{4\pi\varepsilon_0r} +V_H(r)+V_{xc}(r)\right]\psi_i=\varepsilon_i\psi_i. $$
We use Bohr radii for length and Rydbergs for energy. In these units,
$$ \left[-\nabla^2-\frac{2Z_{\mathrm{nuclear}}}{r} +V_H(r)+V_{xc}(r)-\varepsilon_i\right]\psi_i=0. $$
Writing $\psi_{lm}(\mathbf r)=Y_{lm}(\hat{\mathbf r})u_l(r)/r$ gives
$$ u_l''(r)=\left[ \frac{l(l+1)}{r^2}-\frac{2Z_{\mathrm{nuclear}}}{r} +V_H(r)+V_{xc}(r)-\varepsilon \right]u_l(r). $$
The nuclear charge remains the physical $Z_{\mathrm{nuclear}}$. Screening is generated by $V_H[\rho]$ and $V_{xc}[\rho]$; it is not represented by one constant screened charge.
import numpy as np
import matplotlib.pyplot as plt
from numba import njit
from numpy.polynomial import Polynomial
from scipy import integrate, optimize
from scipy.interpolate import PchipInterpolator
from excor import ExchangeCorrelation
%matplotlib inline
Radial solver inherited from the hydrogen calculation¶
Define
$$ U_{\mathrm{KS}}(r)=U_H(r)-2Z_{\mathrm{nuclear}}+rV_{xc}(r), \qquad U_H(r)=rV_H(r). $$
The radial equation becomes
$$ u_l''(r)=\left[ \frac{l(l+1)}{r^2}+\frac{U_{\mathrm{KS}}(r)}{r}-\varepsilon \right]u_l(r). $$
The following functions are the Numerov, shooting, root-finding, occupation, and density machinery developed in the preceding notebook, generalized to an arbitrary radial potential.
@njit
def numerov_homogeneous(coefficient, x0, slope, step):
"""Solve x'' = coefficient*x on a uniform mesh."""
x = np.zeros(coefficient.size)
x[0] = x0
x[1] = x0 + step*slope
step2 = step**2
step2_over_12 = step2/12.0
w_previous = x[0]*(1.0 - step2_over_12*coefficient[0])
w_current = x[1]*(1.0 - step2_over_12*coefficient[1])
x_current = x[1]
coefficient_current = coefficient[1]
for i in range(2, coefficient.size):
w_next = 2*w_current - w_previous + step2*coefficient_current*x_current
coefficient_current = coefficient[i]
x_current = w_next/(1.0 - step2_over_12*coefficient_current)
x[i] = x_current
w_previous, w_current = w_current, w_next
return x
@njit
def radial_coefficient(energy, ell, radius, U_ks):
return (ell*(ell + 1)/radius + U_ks)/radius - energy
def compute_radial_orbital(energy, ell, radius, U_ks):
coefficient = radial_coefficient(energy, ell, radius, U_ks)
u = numerov_homogeneous(
coefficient[::-1],
x0=0.0,
slope=-1e-10,
step=radius[0] - radius[1],
)[::-1]
norm = integrate.simpson(u**2, x=radius)
return u/np.sqrt(abs(norm))
def shooting_residual(energy, radius, ell, U_ks):
u = compute_radial_orbital(energy, ell, radius, U_ks)/radius**ell
polynomial = Polynomial.fit(radius[:4], u[:4], deg=3)
return polynomial(0.0)
def find_bound_states(radius, ell, n_states, energy_mesh, U_ks):
states = []
residual_left = shooting_residual(energy_mesh[0], radius, ell, U_ks)
for energy_left, energy_right in zip(energy_mesh[:-1], energy_mesh[1:]):
residual_right = shooting_residual(energy_right, radius, ell, U_ks)
if residual_left*residual_right < 0.0:
energy = optimize.brentq(
shooting_residual,
energy_left,
energy_right,
args=(radius, ell, U_ks),
xtol=1e-12,
)
states.append((ell, energy))
if len(states) >= n_states:
break
residual_left = residual_right
return states
def fill_orbitals(states, radius, n_electrons, U_ks, verbose=False):
rho = np.zeros(radius.size)
eigenvalue_sum = 0.0
electron_count = 0
occupied = []
for ell, energy in states:
capacity = 2*(2*ell + 1)
occupation = min(capacity, n_electrons - electron_count)
if occupation <= 0:
break
u = compute_radial_orbital(energy, ell, radius, U_ks)
rho += occupation*u**2/(4*np.pi*radius**2)
eigenvalue_sum += occupation*energy
electron_count += occupation
occupied.append((ell, energy, occupation, capacity))
if verbose:
print(
f"l={ell}, E={energy/2:11.6f} Ha, "
f"occupation={occupation}/{capacity}"
)
if electron_count >= n_electrons:
break
if electron_count != n_electrons:
raise RuntimeError(
f"Only {electron_count} of {n_electrons} electrons could be placed. "
"Increase n_max or improve the energy mesh."
)
return rho, eigenvalue_sum, occupied
For a bare Coulomb potential, states with the same principal quantum number are degenerate. Small discretization errors can change their numerical order. During the first iteration we therefore use $l$ only as a tiny tie-breaker; after screening lifts the degeneracy, the calculated energy determines the order.
def sort_orbitals(states, tie_break=1e-5):
return sorted(states, key=lambda state: state[1] + tie_break*state[0])
def make_energy_mesh(Z_nuclear, n_points=500, energy_shift=0.1):
# The small positive shift gives fine sampling near the continuum edge.
# Positive roots in a finite radial box must not automatically be
# interpreted as physical atomic bound states.
energy_scale = 1.2*Z_nuclear**2 + energy_shift
return (
-np.geomspace(1e-4, energy_scale, n_points)[::-1]
+ energy_shift
)
A screened hydrogenic starting density¶
As in the preceding notebook, first construct a smooth model density in a reduced Coulomb potential. This is only an initial guess: the LDA calculation below restores the full nuclear potential and determines screening self-consistently.
radius_demo = np.linspace(1e-6, 100.0, 2000)
Z_eff = 1.0
n_electrons_demo = 8
n_max_demo = 4
U_model = -2*Z_eff*np.ones(radius_demo.size)
energy_mesh_demo = make_energy_mesh(Z_eff, n_points=400)
model_states = []
for ell in range(n_max_demo):
model_states += find_bound_states(
radius_demo, ell, n_max_demo - ell, energy_mesh_demo, U_model
)
model_states = sort_orbitals(model_states, tie_break=1e-3)
rho_demo, _, occupied_demo = fill_orbitals(
model_states, radius_demo, n_electrons_demo, U_model, verbose=True
)
normalization = integrate.simpson(
4*np.pi*radius_demo**2*rho_demo, x=radius_demo
)
print("Integrated density =", normalization)
plt.plot(
radius_demo,
4*np.pi*radius_demo**2*rho_demo,
label="initial radial density",
)
plt.xlim(0.0, 25.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$4\pi r^2\rho(r)$")
plt.legend()
plt.show()
l=0, E= -0.500000 Ha, occupation=2/2 l=0, E= -0.125000 Ha, occupation=2/2 l=1, E= -0.125000 Ha, occupation=4/6 Integrated density = 8.0
Step 1: Hartree potential¶
In Rydberg units the electrostatic potential generated by the electron density is
$$ V_H(\mathbf r)=2\int d^3r'\, \frac{\rho(\mathbf r')}{|\mathbf r-\mathbf r'|}. $$
The factor of two appears because one Hartree equals two Rydbergs. The Hartree term contains the classical self-interaction; the exchange-correlation approximation cancels it only approximately.
For the spherically averaged density used here,
$$ \nabla^2V_H(r)=-8\pi\rho(r). $$
Writing $V_H(r)=U_H(r)/r$ reduces Poisson's equation to
$$ U_H''(r)=-8\pi r\rho(r). $$
The boundary conditions are
$$ U_H(0)=0, \qquad U_H(\infty)=2N_e. $$
The second condition depends on the number of electrons, not on the nuclear charge. They are equal only for a neutral atom.
Direct ODE solution¶
A particular solution of Poisson's equation may be shifted by the homogeneous solution $\alpha r$. We choose $\alpha$ so that the large-radius boundary condition is satisfied:
$$ U_H(r)=U_p(r)+\alpha r, \qquad \alpha=\frac{2N_e-U_p(R_{\max})}{R_{\max}}. $$
def hartree_rhs(state, radius, rho_interpolator):
U_h, U_h_prime = state
return np.array([
U_h_prime,
-8*np.pi*radius*rho_interpolator(radius),
])
rho_interpolator = PchipInterpolator(radius_demo, rho_demo)
U_h_ode = integrate.odeint(
hartree_rhs,
[0.0, 0.5],
radius_demo,
args=(rho_interpolator,),
)[:, 0]
alpha = (2*n_electrons_demo - U_h_ode[-1])/radius_demo[-1]
U_h_ode += alpha*radius_demo
plt.plot(radius_demo, U_h_ode, label=r"$U_H(r)$")
plt.axhline(2*n_electrons_demo, color="black", linestyle="--", label=r"$2N_e$")
plt.xlim(0.0, 20.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$U_H(r)$")
plt.legend()
plt.show()
Numerov solution of Poisson's equation¶
Poisson's equation has the inhomogeneous form $x''(t)=q(t)$. With
$$ w_i=x_i-\frac{h^2}{12}q_i, $$
the Numerov recurrence is
$$ w_{i+1}=2w_i-w_{i-1}+h^2q_i, \qquad x_i=w_i+\frac{h^2}{12}q_i. $$
@njit
def numerov_inhomogeneous(source, x0, slope, step):
x = np.zeros(source.size)
x[0] = x0
x[1] = x0 + step*slope
step2 = step**2
step2_over_12 = step2/12.0
w_previous = x[0] - step2_over_12*source[0]
w_current = x[1] - step2_over_12*source[1]
source_current = source[1]
for i in range(2, source.size):
w_next = 2*w_current - w_previous + step2*source_current
source_current = source[i]
x[i] = w_next + step2_over_12*source_current
w_previous, w_current = w_current, w_next
return x
def hartree_u(radius, rho, n_electrons):
source = -8*np.pi*radius*rho
U_h = numerov_inhomogeneous(
source, x0=0.0, slope=0.5, step=radius[1] - radius[0]
)
alpha = (2*n_electrons - U_h[-1])/radius[-1]
return U_h + alpha*radius
U_h_numerov = hartree_u(radius_demo, rho_demo, n_electrons_demo)
plt.plot(radius_demo, U_h_ode, ".", label="ODE solver")
plt.plot(radius_demo, U_h_numerov, label="Numerov")
plt.xlim(0.0, 20.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$U_H(r)$")
plt.legend()
plt.show()
Step 2: Local exchange-correlation potential¶
In LDA,
$$ V_{xc}(r)=V_{xc}(\rho(r)). $$
The supplied excor.py module parametrizes the uniform-electron-gas result using
$$ r_s=\left(\frac{3}{4\pi\rho}\right)^{1/3}. $$
The module returns Hartree units, so its energies and potentials are multiplied by two before being combined with our Rydberg-unit Schrödinger equation.
exchange_correlation = ExchangeCorrelation()
@njit
def density_parameter(rho):
if rho < 1e-100:
return 1e100
return (3/(4*np.pi*rho))**(1/3)
def exchange_correlation_potential(rho):
return np.array([
2*(exchange_correlation.Vx(density_parameter(value))
+ exchange_correlation.Vc(density_parameter(value)))
for value in rho
])
def exchange_correlation_correction(rho):
"""Return 2*(epsilon_xc - V_xc), in Rydbergs."""
return np.array([
2*(exchange_correlation.ExVx(density_parameter(value))
+ exchange_correlation.EcVc(density_parameter(value)))
for value in rho
])
Step 3: Effective Kohn–Sham potential¶
For a neutral atom, $U_H(r)\to2Z_{\mathrm{nuclear}}$, so the Hartree and nuclear $1/r$ terms cancel far from the atom. Because a local LDA exchange-correlation potential follows the exponentially small density, it also decays rapidly. This short-ranged LDA tail is not the exact Kohn–Sham asymptotic behavior; it is one manifestation of the LDA self-interaction error.
Z_nuclear_demo = n_electrons_demo
V_xc_demo = exchange_correlation_potential(rho_demo)
U_ks_demo = U_h_numerov - 2*Z_nuclear_demo + radius_demo*V_xc_demo
fig, axes = plt.subplots(2, 1, figsize=(8, 7), sharex=True)
axes[0].plot(radius_demo, U_ks_demo, label=r"$U_{KS}$")
axes[0].plot(
radius_demo,
U_h_numerov - 2*Z_nuclear_demo,
label=r"$U_H-2Z$",
)
axes[0].set(ylabel="scaled potential", xlim=(0.0, 20.0))
axes[0].legend()
axes[0].grid()
axes[1].plot(radius_demo, V_xc_demo, label=r"$V_{xc}$")
axes[1].set(xlabel=r"$r/a_0$", ylabel="potential (Ry)", xlim=(0.0, 20.0))
axes[1].legend()
axes[1].grid()
plt.show()
Steps 4 and 5: New density and mixing¶
For a given $U_{\mathrm{KS}}$, find the bound states, fill the lowest orbitals, and construct $\rho_{\mathrm{new}}$. Direct replacement can oscillate, so use linear mixing:
$$ \rho^{(k+1)}=(1-\alpha)\rho^{(k)}+\alpha\rho_{\mathrm{new}}^{(k)}. $$
Open subshells are filled fractionally over the magnetic and spin degeneracy. This preserves spherical symmetry and should be understood as a spherical average.
Total energy¶
The sum of Kohn–Sham eigenvalues double-counts the Hartree interaction and contains $V_{xc}$ instead of the exchange-correlation energy per particle. In Hartree units,
$$ E_{\mathrm{tot}} =\sum_{i\in\mathrm{occ}} f_i\varepsilon_i -\frac12\int d^3r\,\rho(r)V_H(r) +\int d^3r\,\rho(r)\left[\epsilon_{xc}(r)-V_{xc}(r)\right]. $$
The code performs the internal calculation in Rydbergs and divides the final result by two. Before self-consistency, this expression mixes orbitals and densities from adjacent iterations, so intermediate energies are convergence diagnostics; the final converged value is the physical LDA result.
def total_energy_rydberg(radius, rho, U_h, eigenvalue_sum):
xc_correction = exchange_correlation_correction(rho)
radial_integrand = (
xc_correction*radius - 0.5*U_h
)*radius*rho*4*np.pi
return eigenvalue_sum + integrate.simpson(radial_integrand, x=radius)
Self-consistent LDA solver¶
The complete iteration is now collected in one function so that oxygen, potassium, and heavier atoms use exactly the same implementation.
def solve_atom_lda(
Z_nuclear,
radius,
n_electrons=None,
n_max=4,
mixing=0.5,
energy_points=500,
max_iterations=100,
energy_tolerance=1e-7,
density_tolerance=1e-6,
initial_effective_charge=1.0,
verbose=True,
):
if n_electrons is None:
n_electrons = Z_nuclear
energy_mesh = make_energy_mesh(Z_nuclear, energy_points)
# Begin from the screened hydrogenic model developed in notebook 04.
# This avoids the severe overcontraction produced by filling the bare
# -2*Z_nuclear/r spectrum. The actual LDA potential constructed below
# still contains the full nuclear singularity.
U_initial = -2*initial_effective_charge*np.ones(radius.size)
initial_energy_mesh = make_energy_mesh(
initial_effective_charge, energy_points
)
initial_states = []
for ell in range(n_max):
initial_states += find_bound_states(
radius,
ell,
n_max - ell,
initial_energy_mesh,
U_initial,
)
initial_states = sort_orbitals(initial_states, tie_break=1e-3)
rho, _, _ = fill_orbitals(
initial_states, radius, n_electrons, U_initial
)
U_h = hartree_u(radius, rho, n_electrons)
V_xc = exchange_correlation_potential(rho)
U_ks = U_h - 2*Z_nuclear + radius*V_xc
previous_energy = None
history = []
for iteration in range(max_iterations):
states = []
for ell in range(n_max):
states += find_bound_states(
radius, ell, n_max - ell, energy_mesh, U_ks
)
states = sort_orbitals(states)
rho_output, eigenvalue_sum, occupied = fill_orbitals(
states, radius, n_electrons, U_ks
)
density_change = integrate.simpson(
4*np.pi*radius**2*np.abs(rho_output - rho), x=radius
)
rho_mixed = mixing*rho_output + (1 - mixing)*rho
U_h = hartree_u(radius, rho_mixed, n_electrons)
V_xc = exchange_correlation_potential(rho_mixed)
U_ks_new = U_h - 2*Z_nuclear + radius*V_xc
energy_rydberg = total_energy_rydberg(
radius, rho_mixed, U_h, eigenvalue_sum
)
energy_hartree = energy_rydberg/2
energy_change = (
np.inf if previous_energy is None
else abs(energy_hartree - previous_energy)
)
charge = integrate.simpson(
4*np.pi*radius**2*rho_mixed, x=radius
)
history.append((energy_hartree, energy_change, density_change))
if verbose:
print(
f"iteration {iteration:2d}: "
f"E = {energy_hartree:14.8f} Ha, "
f"dE = {energy_change:9.2e}, "
f"dRho = {density_change:9.2e}, N = {charge:.8f}"
)
rho = rho_mixed
U_ks = U_ks_new
if (
energy_change < energy_tolerance
and density_change < density_tolerance
):
break
previous_energy = energy_hartree
else:
raise RuntimeError("The self-consistency loop did not converge")
return {
"radius": radius,
"rho": rho,
"U_h": U_h,
"U_ks": U_ks,
"states": states,
"occupied": occupied,
"energy_hartree": energy_hartree,
"history": np.asarray(history),
"iterations": iteration + 1,
}
Oxygen¶
For oxygen, use a dense uniform mesh to resolve the compact $1s$ orbital and the valence tail simultaneously.
radius_oxygen = np.linspace(1e-8, 10.0, 2**14 + 1)
oxygen = solve_atom_lda(
Z_nuclear=8,
radius=radius_oxygen,
n_max=4,
mixing=0.5,
)
print(f"Converged in {oxygen['iterations']} iterations")
print(f"Total energy = {oxygen['energy_hartree']:.9f} Ha")
print("NIST LDA value = -74.473077 Ha")
plt.plot(
oxygen["radius"],
4*np.pi*oxygen["radius"]**2*oxygen["rho"],
)
plt.xlim(0.0, 5.0)
plt.xlabel(r"$r/a_0$")
plt.ylabel(r"$4\pi r^2\rho(r)$")
plt.show()
iteration 0: E = -111.39463819 Ha, dE = inf, dRho = 1.35e+01, N = 8.00000000 iteration 1: E = -86.53029567 Ha, dE = 2.49e+01, dRho = 5.91e+00, N = 8.00000000 iteration 2: E = -77.10259492 Ha, dE = 9.43e+00, dRho = 3.29e+00, N = 8.00000000 iteration 3: E = -75.02274729 Ha, dE = 2.08e+00, dRho = 1.82e+00, N = 8.00000000 iteration 4: E = -74.75344125 Ha, dE = 2.69e-01, dRho = 8.74e-01, N = 8.00000000 iteration 5: E = -74.60864563 Ha, dE = 1.45e-01, dRho = 4.22e-01, N = 8.00000000 iteration 6: E = -74.53934656 Ha, dE = 6.93e-02, dRho = 2.05e-01, N = 8.00000000 iteration 7: E = -74.50553287 Ha, dE = 3.38e-02, dRho = 9.96e-02, N = 8.00000000 iteration 8: E = -74.48901699 Ha, dE = 1.65e-02, dRho = 4.88e-02, N = 8.00000000 iteration 9: E = -74.48091813 Ha, dE = 8.10e-03, dRho = 2.40e-02, N = 8.00000000 iteration 10: E = -74.47693825 Ha, dE = 3.98e-03, dRho = 1.18e-02, N = 8.00000000 iteration 11: E = -74.47497939 Ha, dE = 1.96e-03, dRho = 5.87e-03, N = 8.00000000 iteration 12: E = -74.47401435 Ha, dE = 9.65e-04, dRho = 2.92e-03, N = 8.00000000 iteration 13: E = -74.47353865 Ha, dE = 4.76e-04, dRho = 1.45e-03, N = 8.00000000 iteration 14: E = -74.47330408 Ha, dE = 2.35e-04, dRho = 7.27e-04, N = 8.00000000 iteration 15: E = -74.47318835 Ha, dE = 1.16e-04, dRho = 3.64e-04, N = 8.00000000 iteration 16: E = -74.47313123 Ha, dE = 5.71e-05, dRho = 1.82e-04, N = 8.00000000 iteration 17: E = -74.47310301 Ha, dE = 2.82e-05, dRho = 9.16e-05, N = 8.00000000 iteration 18: E = -74.47308905 Ha, dE = 1.40e-05, dRho = 4.61e-05, N = 8.00000000 iteration 19: E = -74.47308214 Ha, dE = 6.91e-06, dRho = 2.32e-05, N = 8.00000000 iteration 20: E = -74.47307872 Ha, dE = 3.42e-06, dRho = 1.17e-05, N = 8.00000000 iteration 21: E = -74.47307702 Ha, dE = 1.70e-06, dRho = 5.91e-06, N = 8.00000000 iteration 22: E = -74.47307617 Ha, dE = 8.50e-07, dRho = 2.98e-06, N = 8.00000000 iteration 23: E = -74.47307575 Ha, dE = 4.19e-07, dRho = 1.51e-06, N = 8.00000000 iteration 24: E = -74.47307555 Ha, dE = 2.07e-07, dRho = 7.64e-07, N = 8.00000000 iteration 25: E = -74.47307544 Ha, dE = 1.05e-07, dRho = 3.87e-07, N = 8.00000000 iteration 26: E = -74.47307539 Ha, dE = 5.39e-08, dRho = 1.96e-07, N = 8.00000000 Converged in 27 iterations Total energy = -74.473075389 Ha NIST LDA value = -74.473077 Ha
Agreement with the NIST value checks that we have reproduced the same non-spin-polarized LDA calculation. It does not mean that LDA is exact; experiment and more accurate many-body calculations describe a different level of approximation.
Homework¶
- Reproduce the oxygen calculation and compare the total energy and orbital eigenvalues with NIST.
- Solve potassium, $Z=19$. Verify that the self-consistent potential places the occupied $4s$ orbital below the unoccupied $3d$ orbital. This ordering is absent from the degenerate hydrogenic model.
- Explore how far through the periodic table the present uniform-grid algorithm remains reliable. When it fails, determine whether the limitation is the radial spacing, $R_{\max}$, energy mesh, number of angular-momentum channels, or density mixing.