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?

  • $\begingroup$ 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. $\endgroup$
    – Andrew
    Commented Sep 20, 2021 at 1:07
  • 1
    $\begingroup$ 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$. $\endgroup$
    – Roger V.
    Commented Sep 20, 2021 at 11:23
  • $\begingroup$ @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. $\endgroup$
    – MadLad
    Commented Sep 20, 2021 at 16:29
  • $\begingroup$ @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)? $\endgroup$
    – Roger V.
    Commented Sep 20, 2021 at 18:38
  • $\begingroup$ 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. $\endgroup$
    – MadLad
    Commented Sep 20, 2021 at 19:03

1 Answer 1


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

  • $\begingroup$ 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. $\endgroup$
    – MadLad
    Commented Sep 20, 2021 at 15:29
  • 1
    $\begingroup$ 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. $\endgroup$
    – Oh My C
    Commented Sep 20, 2021 at 15:56
  • 1
    $\begingroup$ 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. $\endgroup$
    – Oh My C
    Commented Sep 20, 2021 at 15:57
  • $\begingroup$ 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 $\endgroup$
    – MadLad
    Commented Sep 20, 2021 at 16:22
  • $\begingroup$ @MadLad : I added a pseudo code, have a look at at. This is just for understanding. you can write a much better code. $\endgroup$
    – Oh My C
    Commented Sep 20, 2021 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.