Implementing a program

The first thing we should do is to define a function that calculates our total Hamiltonian. Luckily, we already know how to do this: the interaction is on the diagonal, and the kinetic energy matrix is as before. Hence, our function just takes as input the current density and returns the matrix.

Then we need a way to calculate a new density given a set of orbitals. This is, again, not really difficult: simply pass the orbitals required, and then sum their squares multiplied by 2 (spin up + spin down electron in each orbital).

We also need an initial guess. As I said previously, using the density of the non-interacting system is a good first try. The rest of it is just iterating on this process over and over again, until the difference between the input and output density is sufficiently small. To check for this, we use the following formula:

\begin{align} \Delta \rho = \int _0^L |\rho _{prev} (x) - \rho _{current}(x)|^2\mathrm{d}x . \end{align}

It's easy to understand that this number will be bigger if the difference between the densities is bigger. Of course, you could invent your own criterion for convergence; for instance, you could check that the energy levels of two different iterations are not too far from each other, or that the maximum difference between the densities is no greater than some number. Nevertheless, it is customary to use the squared integral.

We will once again use the same basis functions as previously, the solutions of the potential well. We calculate the matrix elements by the method we previously learned: that is,

\begin{align} H_{ij} = \int _0^L \phi _i (x)V(x) \phi _j(x)\mathrm{d}x \approx \sum _{x_i}\phi _i(x_i)V(x_i)\phi_j(x_j)\Delta x. \end{align}

This means that on the diagonal we get the eigenvalues of the particle in a box problem plus the contribution from the interaction. That's because

\begin{align} -\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{d}x^2} \phi _i(x) = E_i\phi _i(x) \end{align}

since the particle-in-a-box problem only has the kinetic energy part. Hence,

\begin{align} \int _0^L\phi _j(x)\bigg( -\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{d}x^2}\bigg) \phi _i(x)\mathrm{d}x = E_i\delta _{ij} \end{align}

The non-diagonal elements for the kinetic energy are zero, because the functions \(\phi _i \) are orthogonal, meaning

\begin{align} \int _0^L \phi _i(x) \phi _j(x) \mathrm{d}x = 0. \end{align}

This is again analogous to normal geometry, where the inner product of two orthogonal vectors is also zero.

It is striking that how little we have to actually change the code. The logic still works exactly the same: all you need to do is calculate the matrix differently, and then get the density by a different method. Our eigenvalue solution will now not contain the value of \(\psi\) at every grid point, but it will rather give us the values \(c_i\) in the expansion

\begin{align} \sum _i c_i \phi _i (x)= \psi (x). \end{align}

Hence, we have to calculate the density from our knowledge of the \( c_i \)s. This is not difficult: we just multiply each basis function \(\phi _i \) with the corresponding \(c_i\) from our eigenvector solution. Here's the code [github:dftbasis.py]:


import numpy as np
import matplotlib.pyplot as plt

L=8.0
Np=1000
Norb = 5
lamb = 50
maxiter = 200
Nbasis = 100
x = np.linspace(0,L,Np)
dx = x[1]-x[0]
initrho = np.zeros(Np)

for i in range(1,Norb+1):
    initrho[:]+=2*2.0/L*np.sin(i*x[:]*np.pi/L)**2.0
print(np.sum(initrho)*dx)

def integrate(f):
    return np.sum(f)*dx
def basis_function(i):
    return np.sqrt(2.0/L)*np.sin((i+1)*x[:]*np.pi/L)
def normalized_orbital_density(vec):
    phi = np.zeros(Np)
    for i in range(Nbasis):
        phi[:] += vec[i]*basis_function(i)
    orbdens = phi[:]**2
    
    return 2*orbdens
def hamiltonian(rho):
    H = np.zeros((Nbasis,Nbasis))
    for i in range(Nbasis):
        H[i,i] = ((i+1)**2*np.pi**2)/(2*L)

    for i in range(Nbasis):
        for j in range(Nbasis):
            pot = lamb*0.25*rho[:]
            b1 = basis_function(i)
            b2 = basis_function(j)
            H[i,j] += integrate(b1*pot*b2)

    return H
def get_density(vec):
    newdensity = np.zeros(Np)
    for i in range(Norb):
        orbrho=normalized_orbital_density(vec[:,i])
        newdensity+=orbrho
    return newdensity
def rhodiff(prevrho,rho):
    return integrate((prevrho-rho)**2)
rho = initrho
prevrho = initrho
alpha = 0.1


for i in range(maxiter):
    H = hamiltonian(rho)
    #print(H)
    energies,vec=np.linalg.eigh(H)

    newrho = get_density(vec[:,0:Norb])
    
    if rhodiff(prevrho,newrho) < 1e-5:
        print("woo!")
        break
    rho = (1-alpha)*prevrho+alpha*newrho
    prevrho=rho
    print("iteration ",i,energies[0])
print(dx*np.sum(rho))
plt.plot(initrho)
plt.plot(rho)
plt.show()

Here, I've defined a couple of helper functions. For example, basis_function(i) returns the ith basis function and integrate does numerical integration. The real difference is in calculating the Hamiltonian and the normalized orbital density, as you can see --- but it's nothing we haven't already dealt with before.

Now, this code is actually slower than our first one. Why is that? It's because of all the integrals. We're calculating integrals for \(100\cdot 100 = 10 000\) elements.

What can we do to speed this up? I encourage you to try the following:

Having dealt with one-dimensional toy systems, we'll now turn our attention to real systems: atoms. Though single isolated atoms aren't really found in nature, they are an interesting case to study, since you can find experimental data on the energy levels of atoms. It is also an example of a system where analytical pen-and-paper solution is impossible for everything except the hydrogen atom, so we really do need to do it with some sort of approximation method.

Chapter 2.2Index Chapter 3.1