How to find Density of states numerically?

[PRB147]

If one has already known the dispersion dataset,
for example, for 2D crystal, we know the 1000*2000 dataset for E=E(k_x,k_y),
How to find the density of states numerically?

 

[PhysTech]

Equation (8.63) from Ashcroft and Mermin will be helpful. I have listed it below for convenience:

$g_n(\mathcal{E})=\int_{S_n(\mathcal{E})} \frac{dS}{4\pi^3}\frac{1}{|\nabla\mathcal{E}_n(k_x,k_y,k_z)|}$

I know you want to know how to compute density of states numerically. I was referring to this formula merely to give an insight behind the procedure involved in doing this computation. So what is happening in the above formula is this: in k-space you have found a surface with equal energy $\mathcal{E}$ and you are integrating the quantity $|\nabla\mathcal{E}_n(\textbf{k})|^{-1}$ over the surface. In other words, this integral is "counting" the number of states that lie on this surface of equal energy. This is the density of states by definition. Now, in your data set, you pretty much have to: (1) pick an energy $\mathcal{E}$ (2) determine how many data points, which are indexed by $(k_x,k_y)$, have an energy $\mathcal{E}(k_x,k_y)$ such that $\mathcal{E} \le \mathcal{E}(k_x,k_y) \le \mathcal{E} + \Delta \mathcal{E}$. Then your density of states will be i.e. $g(\mathcal{E}) = N(\mathcal{E})/\Delta \mathcal{E}$, where $N(\mathcal{E})$ is the number of data points are counted.

 

[PRB147]

Thank you very much, PhysTech. Thank you!

 

[daveyrocket]

The easy way to do this numerically is to use the equation g(\epsilon) = \sum_{k} \delta(\epsilon - \epsilon_k)
and use an approximation for the delta function. Basically you setup a grid for your energies E, initialize your density of states function g(E) on those points to zero, and then loop over your momentum (and band index) and add the delta function approximation to your g(E).

 

[PRB147]

The easy way to do this numerically is to use the equation g(\epsilon) = \sum_{k} \delta(\epsilon - \epsilon_k)
and use an approximation for the delta function. Basically you setup a grid for your energies E, initialize your density of states function g(E) on those points to zero, and then loop over your momentum (and band index) and add the delta function approximation to your g(E).

Numerically, it seems that the dos is calculated from the integration over BZ through the so called tetrahedron methods. But I have no experience on that.

 

[PRB147]

Thank Daveyrocket and PhysTech very much, discussion is always helpful

 

[OhYoungLions]

I've used daveyrocket's method with a Lorentzian of small width in the past. Not sure if this is the most efficient approximation for the delta function, but it's real easy to work with.