Solving the KS equation

To solve the KS equations, we first need to guess the density (or the potential, which amounts to the same thing). We will use the Thomas-Fermi potential as the initial guess. This comes from an arxiv paper, though the Thomas-Fermi theory itself is older than that. We only use this for the initial guess, so there's no use going through it more carefully. You could also guess the density and try another initial guess - in fact, I recommend that you try!

What about the basis functions? Clearly, our cosine basis functions aren't going to be good for atoms, because the electron density should go to zero far away from the nucleus. Instead, we'll use the Slater functions, which Slater specifically invented for describing atomic densities. They have the following form:

\begin{align} R(r) &= N r^{n-1}e^{-\alpha r} \label{slater}\\ N &= (2\alpha)^n \sqrt{\frac{2\alpha}{(2n)!}} \end{align}

We will only use the case n=1, the functions "closest to the nucleus". This will turn out to be enough for beryllium, which we will investigate here. Beryllium only has 2 occupied energy levels and 4 electrons, so we don't need that many basis functions. Also, since our equation actually involves \( P(r)=rR(r)\), we'll multiply the above with r.

We will calculate the Hartree potential by expanding the density from the basis functions and then calculating the potential as discussed previously in the Hartree section.

Now, the kinetic energy operator \( -\frac{1}{2}\frac{d^2}{dr^2} ^2\) can be calculated analytically. The second derivative of the \eqref{slater} is

\begin{align} -\frac{1}{2}\frac{d^2}{dr^2} P(r) = -\frac{1}{2}N(\alpha^2e^{-\alpha r}r -\alpha e^{-\alpha r} -\alpha e^{-\alpha r} ) \end{align}

We see that these terms will be present on not just the diagonal but off diagonal as well. This is because this time around our basis functions aren't orthogonal (meaning that now \( \int _0^L \phi _n (x) \phi _m(x) \mathrm{d}x \neq 0 \) ). For this reason, we also aren't deling with a straightforward eigenvalue equation; for that, we assumed the orthogonality of the basis functions (go back to the chapter on basis functions if you can't remember what I'm talking about). Instead our problem is of the form

\begin{align} \hat{H}\vec{c} = \epsilon S \vec{c} \end{align}

where S is an overlap matrix:

\begin{align} S_{nm} = \int _0^\infty \phi _n (x) \phi _m(x) \mathrm{d}x \end{align}

In the case we previously calculated (for the potential well) this integral was zero except for the diagonal. Hence, our problem had a simpler form.

Fortunately, this won't cause any problems for us, because SciPy's eigh already allows us to specify just such an overlap matrix. All we need to do is calculate this overlap matrix once and supply it to the eigenvalue algorithm every time.

Finally, it should be noted that the normalization of our density is \( \int _0 ^\infty P(r)^2 \mathrm{d}r = 1 \). If we want to plot the average radial density at a given distance, we have to integrate over the angular coordinates, that is

\begin{align} R_n(r) &= P_n(r) \frac{1}{4\pi r^2} \\ \rho (r) &= \int 4\pi r^2 n(r,\theta, \phi) \mathrm{d}\theta \mathrm{d}\phi\\ \implies \rho (r) &= \frac{1}{4\pi r^2} \sum _n 2P_n(r)^2 \end{align}

That's why at the end of the following code, we have the factor \( 4\pi r^2 \) before plotting.

The final part that we must keep in mind is that the nucleus has the potential \( -Z/r \), in our case \( -4/r \).

We're now ready to look at the code, which is here:


import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.linalg import eigh
from scipy.integrate import quad
from scipy.integrate import solve_ivp
import math
N=100
#We solve up to 30 Bohr away from the nucleus.
endpoint = 30
Npoints = 3000
Norb = 2
r = np.linspace(1e-3,endpoint,Npoints)
dx = r[1]-r[0]
def integrate(f):
    return np.sum(f)*dx
def rhodiff(prevrho,rho):
    return integrate((prevrho-rho)**2)
def slater(order,alpha):
    #Slater functions. We use order = 1 throughout, but you could
    #try using also basis functions of higher orders in addition.
    N = (2*alpha)**order*np.sqrt(2*alpha/math.factorial(2*order))
    func = N*r[:]**(order-1)*np.exp(-alpha*r[:])
    return r[:]*func
def slater_derivative(order,alpha):
    #This calculates the -0.5 * second derivative of the basis functions, as described
    #in the text.
    N = (2*alpha)**order*np.sqrt(2*alpha/math.factorial(2*order))
    dfunc = N*(-alpha*np.exp(-alpha*r[:])*r[:] + np.exp(-alpha*r[:]))
    ddfunc = N*(alpha**2*np.exp(-alpha*r[:])*r[:] -alpha*np.exp(-alpha*r[:]) -alpha*np.exp(-alpha*r[:]) )
    return -0.5*ddfunc
def basis_function(i):
    alphastart = 0.3
    alphastep = 0.05
    if i<100:
        return slater(1,alphastart+i*alphastep)

def basis_derivative(i):
    alphastart = 0.3
    alphastep = 0.05
    if i<100:
        return slater_derivative(1,alphastart+i*alphastep)

def normalized_orbital_density(vec):
    phi = np.zeros(Npoints)
    for i in range(N):
        phi[:] += vec[i]*basis_function(i)
    orbdens = phi[:]**2
    normorbdens = orbdens/(integrate(orbdens))
    return 2*normorbdens
def get_density(vec):
    rho = np.zeros(Npoints)
    for i in range(Norb):
        rho +=  normalized_orbital_density(vec[:,i])
    return rho/(4*np.pi*r[:]**2)
def get_hartree(rho):
    densi = interp1d(r,rho)
    def integrand(x):
        return densi(x)*x*4*np.pi
    
    start = quad(integrand,r[0],endpoint)[0]
    def dV_dr(x,V):
        return V[1],-4*np.pi*densi(x)-2/x*V[1]
    sol = solve_ivp(dV_dr, [r[0],endpoint],[start,0]\
                ,t_eval = r,method="LSODA")
    return sol.y[0]

def get_xc(rho):
    #The commented part is the correlation potential. You can uncomment it
    #and see what happens.
    '''y0 = -0.10498
    b = 3.72744
    c = 12.9352
    A = 0.0621814

    rs = (3/(4*np.pi*rho[:]))**(1.0/3)
    Q = np.sqrt(4*c - b**2)
    y = np.sqrt(rs)
    def get_Y(y,b,c):
        return y**2 + b*y + c
    ec = A/2 * (np.log(y**2/get_Y(y, b, c)) + 2*b/Q * np.arctan(Q/(2*y+b))  \
   - b*y0/get_Y(y0, b, c) * ( \
            np.log((y-y0)**2 / get_Y(y, b, c)) \
            + 2*(b+2*y0) / Q * np.arctan(Q/(2*y+b)) \
          ) )
    Vc =  ec - A/6 * (c*(y-y0)-b*y0*y)/((y-y0)*get_Y(y, b, c))'''
    return -3/(4*np.pi) * (3*np.pi**2*rho[:])**(1.0/3)# + Vc
def hamiltonian(rho):
    H=np.zeros((N,N))
    Vh = get_hartree(rho)
    Vxc = get_xc(rho)
    Va = -4/r[:]
    vks = Vh+Vxc+Va
    for i in range(N):
        for j in range(i,N):
            H[i,j]+=integrate(basis_function(i)*vks[:]*basis_function(j))
            H[i,j]+= integrate(basis_function(i)*basis_derivative(j))
            H[j,i] = H[i,j]
    return H
def get_overlap():
    S = np.zeros((N,N))
    for i in range(N):
        for j in range(i,N):
            if i==j:
                S[i,i] = 1
            else:
                S[i,j] = integrate(basis_function(i)*basis_function(j))
                S[j,i] = S[i,j]
    return S
def get_initial_potential():
    Z = 4
    x = r[:] * (128*Z/(9*np.pi**2)) ** (1.0/3)
    alpha = 0.7280642371
    beta = -0.5430794693
    gamma = 0.3612163121
    Z_eff = Z * (1 + alpha*np.sqrt(x) + beta*x*np.exp(-gamma*np.sqrt(x)))**2 * \
        np.exp(-2*alpha*np.sqrt(x))
    for i in range(len(Z_eff)):
        if Z_eff[i] < 1:
            Z_eff[i] = 1
    return  -Z_eff / r[:]
#Initialization. We first calculate the overlap and get the initial TF potential.
overlap = get_overlap()

H = np.zeros((N,N))
vtf = get_initial_potential()
for i in range(N):
    for j in range(i,N):
        H[i,j]+=integrate(basis_function(i)*vtf*basis_function(j))
        H[i,j]+= integrate(basis_function(i)*basis_derivative(j))
        H[j,i] = H[i,j]

energies,vec=eigh(H,b=overlap)
initrho = get_density(vec)
rho = initrho
prevrho = initrho
maxiter = 200
alpha = 0.2
#The standard KS iteration
for i in range(maxiter):
    H = hamiltonian(rho)
    energies,vec=eigh(H,b=overlap)
    newrho = get_density(vec[:,0:Norb])

    rhod = rhodiff(prevrho,newrho)
    if rhodiff(prevrho,newrho) < 1e-5:
        print("woo!")
        break

    rho = (1-alpha)*prevrho+alpha*newrho
    prevrho=rho
    print("iteration ",i,rhod,energies[0:2])


plt.plot(r,4*np.pi*r**2*rho)
plt.show()

The lowest two energy eigenvalues we get from this process are around -3.55 and -0.13. The reference values from NIST are around -3.85 and -0.20, so our result is not too bad for 130 lines of code! See if you can improve the accuracy by starting the grid closer to 0 and using also the correlation potential.

That's it for now. I'm presently considering some simple-to-calculate 3D models which I might add to this little booklet later.

Again, don't hesitate to contact me on Twitter if you have questions.

at the contact form on the front page.

Chapter 2.3Index