How to calcualte density of state from computed phonon dispersion?

Question

I have numerically computed phonon dispersion using small displacement method by solving dynamical matrix along high symmetry direction [100], [110] and [111] for f.c.c lattice. Now I want to calculate density of state i.e $D(\omega)$ and plot $D(\omega)$ vs. $\omega$ in a graph. I have managed to do that by generating random $q$ vectors then computed phonon frequencies using my code and constructed a histogram using bin_size = 150 (though not sure what should be the range of $q$ vectors and the normalization of count). There exist a direct method to compute density of state for 1D phonon dispersion but don't know any existing method for 3D lattice. Can anyone please enlighten me in this matter, any reference or detail explanation will be very helpful.

Thanks

Answer 1

I will present the result for the one-dimensional case only. It is readily generalized to your three-dimensional situation. Knowing the phonon dispersion $\omega(q)$, the DOS can be defined as

$$D(\omega) = \frac{1}{2\pi}\int\mathrm{d}q\,\delta(\omega-\omega(q)).$$

Using the formal identity

$$\delta(g(x)) = \sum_{i}\frac{\delta(x-x_i)}{|g'(x_i)|}$$

where $x_i$ are the zeros of $g$ and $g'$ is its derivative, one can write:

$$D(\omega) = \frac{1}{2\pi}\int\mathrm dq \frac{1}{\left|\frac{\mathrm d\omega}{\mathrm dq}\right|}\sum_{\substack{i\\\omega(q_i)=\omega}}\delta(q-q_i) = \frac{1}{2\pi}\sum_{\substack{i\\\omega(q_i)=\omega}}\left|\frac{\mathrm d\omega}{\mathrm dq}(q_i)\right|^{-1}$$

The sum is over all the points where the dispersion reaches the given frequency $\omega$. The flater the dispersion at those points, the more states contribute!