Density Functional Theory of Atoms

KS equations for atoms

For atoms, we will once again solve the Kohn-Sham equations. In this case, they read \begin{align} \bigg( -\frac{1}{2}\frac{\mathrm{d^2}}{\mathrm{d}r^2} +V_h(r)+V_{xc}(r)+ \frac{l(l+1)}{2r^2}\bigg)P _{nl} = \epsilon P _{nl} . \label{ksatom} \end{align}

where the full solution of the equation is of the form \(\psi = R(r) Y_{lm}(\theta , \phi )\). In other words, we solve the equation in spherical coordinates of r (distance from atomic nucleus) and the angles \(\theta\) and \(\phi\). The functions \(Y_{lm}\) are called the spherical harmonics, which we won't worry about; you can find all of them tabulated on e.g. Wikipedia. The problem we really have to solve is finding \(R(r) = \frac{P(r)}{r}\). We will focus on the simples case where l=0, meaning the atom has only two s-orbitals. Particularly, we'll solve beryllium, which has 4 electrons.

Once again, we are confronted with an eigenvalue problem. We could apply the methods of the previous two chapters --- calculating the Hamiltonian in terms of either a grid or basis functions. However, as we noticed, this leads to big matrices. Especially in the case of dividing the real line to intervals ("real space grid") it turns out that for this problem we would need to use many thousands of grid points, and thus the resulting matrix would be extremely troublesome to handle. That particular method is also very unstable for this problem, frequently resulting in nonsensical results if great care is not taken.

We will instead use basis functions. But before getting in to the basis functions we'll use, we need to be able to calculate the potential. The most difficult part of it is the Hartree potential \(V_h\), which is calculated from the electrostatic Poisson equation

\begin{align} \nabla ^2 V_{h}(r) = -4\pi \rho (r) \end{align}

to which we now turn our attention.

The Hartree potential

In spherical systems, where the Hartree potential depends only on the radial coordinate, the Poisson equation is

\begin{align} \frac{\mathrm{d}^2 V_h(r)}{\mathrm{d}r^2}+\frac{2}{r}\frac{\mathrm{d}V_h(r)}{\mathrm{d}r}=-4\pi \rho (r) \end{align}

Conveniently, this too is a differential equation, which we can cast in to a first order form

\begin{align} A(r) &= V'_h(r)\\ A'(r) + \frac{2}{r}A(r) &= -4\pi \rho (r). \end{align}

SciPy contains the necessary tools for solving this problem for us. In order to solve differential equations, though, we need to know the initial conditions, namely the value of the Hartree potential at \( r=0 \). It turns out that for atoms, one can prove that

\begin{align} V_h(0) &= \int _0^\infty 4\pi r\rho (r)\mathrm{d}r \label{hartree1}\\ V'_h(0) &= 0 .\label{hartree2} \end{align}

The argument to see this is somewhat complicated and not particularly illuminating. Basically, one expands the density \( \rho \) in terms of polynomials like \(c_2r^2\) and then solves the Poisson equation analytically near the nucleus. We can calculate the first integral with the quad function from SciPy, replacing the infinity with some sufficiently large number (we choose 30).

In order to solve this equation we need the electron density for some system. In the github repository I've included one in a file called "density". This density is for beryllium, which has 4 electrons.

Since the density is given at particular discrete distances from the nucleus, we have to interpolate to get the points in between. Interpolation simply means "guessing"\ the value in between points. For example, if we the value of a function at two points \(r=1\) and \(r=2\), we might guess the value between these points by drawing a straight line through the points. More sophisticated algorithms exist, of course: what matters for us is that once again SciPy comes to the rescue by providing us with an interpolation routine. Interp1d, when called on some data, returns a function that can be called at any point \(r\) that we desire.

Once we have an interpolated density, we need to solve the differential equation. Here we can use SciPy's solve_ivp routine. You can look at SciPy's documentation yourself: essentially, we give SciPy our differential equation in the form of a two-element array, one element containing A and the other V. The code then reads [github:hartree.py]:


import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.integrate import solve_ivp
from scipy.integrate import quad
dens = np.loadtxt("results")

densi =  interp1d(dens[:,0],dens[:,1])
def integrand(x):
    return densi(x)*x*4*np.pi

start = quad(integrand,dens[0,0],50)[0]
def dV_dr(x,V):
    return V[1],-4*np.pi*densi(x)-2/x*V[1]

sol = solve_ivp(dV_dr, [dens[0,0],40],[start,0]\
                ,method="LSODA")

hartree = np.loadtxt("hartree")
harti = interp1d(hartree[:,0],hartree[:,1])
plt.grid()
plt.plot(sol.t,sol.y[0])
plt.plot(sol.t,harti(sol.t))
plt.show()

Here you see that V[0] contains the function \(V_h\) and V[1] its derivative. We use the "LSODA" method; you can read SciPy's documentation to discover more available methods, and try them out if you like. Also, in the github repo I've included the file "hartree" which contains the correct Hartree potential for berylium; this code should produce exactly the same result. The potential looks like this:

Hartree potential
Hartree potential for an atom.

In summary, to get the Hartree potential given a density, we simply have to solve equations \eqref{hartree1} and \eqref{hartree2} using SciPy's solve_ivp function. If you don't understand what's going on with solve_ivp, I recommend reading some online resources on its use; it really isn't difficult once you get the hang of it.

Exchange-correlation

There are many choices we could make for the exchange-correlation potential; indeed, creating new exchange-correlation functionals is practically its own subfield of physics. Since I'm mostly concerned with showing the basics, we'll just use the simplest possible model, the LDA exchange:

\begin{align} V_{xc}(r) = -\bigg( \frac{3}{\pi} \bigg)^{1/3} \rho (r)^{1/3}. \end{align}

This functional is derived by considering a completely homogenous electron gas; for this simplified system, it is possible to calculate the exchange part of the energy exactly. There are also many different equations for getting the correlation of a homogenous electron gas, but since those are a tad more complicated, we'll just stick to the exchange.

Chapter 2.3Index Chapter 3.2