# Numerical calculation of density of states

I am trying to figure out the numerical interpretation of density of states for a fermionic system under a periodic potential.

The equation for the density of states reads

$$DOS(E) = \sum_{k \in BZ, n} \delta(E-E_n(k)),$$

where $$E_n(k)$$ are the eigenvalues of the particular Hamiltonian matrix I am solving. I would like to use the Cauchy/Lorentzian approximation of the Delta function such that the first equation now becomes

$$DOS(E) = \frac{1}{\pi} \lim_{\epsilon\to 0} \sum_{k \in BZ, n} \frac{\epsilon}{(E-E_n(k))^2+\epsilon^2}.$$

From here on, I am confused on how to interpret the second equation numerically. I have the respective eigenvalues of the Hamiltonian, but I don't know how to obtain the DOS using $$E$$. How do I include $$E$$ in my code? Discretizing E to me means that I grab a certain energy window around a certain value $$E$$, but I don't know how to structure it, or if needs to be an array, a grid...or something else. If E should be a grid, should it be a grid between the minimum and maximum values of the energy eigenvalues?

EDIT: Hey all. After pondering about Murali's answer I have come up with a pseudo code that is rather bad but I would like to know if I am going in the right direction.

I basically coded a function for the Lorentzian broadened delta function like so:


def delta_l(x):
return (1/np.pi)*(epsilon/(epsilon**2 + x**2))

def dos(Egrid,Eigen):
DOS = np.zeros((AllK,1))

for j in range(Allk):
DOS[j] = (1/AllK)*sum([delta(Egrid[j]-Eigen[i]) for i in range(np.shape(Eigen)[0])])

return DOS



Here epsilon is given a value of 0.1 just to test. The eigenvectors of the Hamiltonian were obtained by inputting points of the FBZ:


AllK = len(np.arange(0, 1, 0.01)) * len(np.arange(0, 1, 0.01))

E  = np.zeros((AllK,4*n), float)

count = 0
for m in np.arange(0, 1, 0.01):
for f in np.arange(0, 1, 0.01):
kx = (m-f) * np.sqrt(3)/2
ky = (m+f) * 3.0/2 - 1
E[count] = Hamiltonian(kx*Kmag, ky*Kmag)
count = count + 1

import pandas as pd
EinBZ = E.flatten()


So then I get all the eigenvalues of the FBZ in this array. Am I going in the right direction?

• Naively I would have thought you'd want to define a function DOS which takes as input $E$ and returns the density of states at that value of $E$. Then you can do things like plot the density of states vs energy by passing an array of $E$ values to this function. Commented Sep 20, 2021 at 1:07
• You need to perform numerical summation/integration over $k$ and $n$ for each value of $E$. In case where $k$ is actually a continuous variable, you might be better off first carrying integration analytically in terms of the roots $E_n(k^*)=E$ and then simply summing the contributions - thsi would avoid any uncertainty due to the value of $\epsilon$. Commented Sep 20, 2021 at 11:23
• @RogerVadim k in this case is not continuous. What I am doing is choosing a grid in momentum space such that I have around 10,000 K points that cover a piece of the Brillouin zone. Commented Sep 20, 2021 at 16:29
• @MadLad I am not sure I understand you: you discretize it in order to do computer calculations, but initially it is continuous (because the crystal is infinite)? Commented Sep 20, 2021 at 18:38
• Right, sorry about that @RogerVadim. We start from an irreducible cell in momentum space, and then we would naturally do an integral in k but since I am doing the numerics I just grab some arrays in k that cover the area of this cell. Commented Sep 20, 2021 at 19:03

I think you did not correctly understand what density of states (DOS) mean. DOS is a probability density function (PDF). As Andrew pointed out, it takes energy as input and returns the number of states for a given energy.

You cannot discretize $$E$$ as they are not eigenvalues of any observable. It is the input parameter, and discretizing it simply does not make any sense. The values of $$E_n(k)$$ are discretized as they are eigenvalues of the electronic hamiltonian.

If you consider the 1st equation in your question,

$$DOS(E) = \sum_{k \in BZ, n} \delta(E-E_n(k)),$$ For Energies $$E\neq E_n(k), DOS(E) = 0$$. However, this is not what you observe in experiments. Typically we see some finite DOS for energies not equal to eigenstates due to Heisenberg uncertainty. To account for this, we add a small finite electronic broadening parameter ($$\epsilon$$) as shown in the second equation of your question.

To compute the DOS, you take E as a parameter which can take any value and fix $$\epsilon$$, then compute the double summation over Brillouin zone (BZ) and bands (n) which are the eigenvalues of your hamiltonian. Summing over BZ is simply summation over k points in the Brillouin zone and divide the obtained sum by the total number of k points. Choose a reasonable k point grid and make sure it is converged. Have a look at the following link (http://www.iiserpune.ac.in/~smr2626/hands_on/week1/july1/bzsums_mastani.pdf) if you don't have any idea about BZ summation

Pseudo code:

def delta_l(x):
return delta function(x)

def E(k):
return Eigen values for Each k

def dos(E): (let us compute for some E value. This is very inefficient way. just writing for your understanding)
sum = 0
for i in kpoints:
for j in total_number_of_bands:
sum = sum + delta_l(E-E(i)[j]) #where E(i)[j] is $$j^{th}$$ eigen value of $$i^{th}$$ kpoint
return sum/N # N is total number of kpoints


• I am aware that I need to sum over K points in the Brillouin zone, and also over the bands. "You take E as a parameter which can take any value". This is what I am mostly confused with. So should E be an array in increments going from the minimum energy to the maximum energy in the system? The Lorentzian in that case would give me values close to 1 for each energy that is close to a band energy, and so DOS will be higher if bands at specific Kx,Ky points have energy close to that E value. Commented Sep 20, 2021 at 15:29
• The Lorentzian is not 1 at E = E(k). It blows up (1/epsilon). Yes, E can take any value, not just between Min and Max. I will give an analogue. Consider a 1D gaussian distribution ( ~e^-{kx^2}). no matter what your x is, you get some value for the distribution (what you are finding is probability density at x). It is distribution, which is a function of x. cont. Commented Sep 20, 2021 at 15:56
• As you go farther to the mean, your values falls off very rapidly. Similarly, you want number of states for a given energy. If there are no states with that energy, your DOS is 0 (For Eg. if you choose E in the band gap, your DOS = 0, as there are no states in band gap). Of course, your dos will be higher if they match with the band energies. Commented Sep 20, 2021 at 15:57
• I think you gave me a better idea on how to think of it and compute it. I will try some stuff, thank you! I will comment anything if I have doubts. @Murali Commented Sep 20, 2021 at 16:22
• @MadLad : I added a pseudo code, have a look at at. This is just for understanding. you can write a much better code. Commented Sep 20, 2021 at 19:07