Numerical calculation of density of states

Question

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?

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