Obtaining phonon density of states from Density Functional Theory

Question

From Statistical Mechanics obtain an expression for the heat capacity of a solid is given by¹:

\begin{equation} \tag{1} {C_V}\left( T \right) = k\int\limits_0^\infty {\frac{{{{\left( {uv} \right)}^2}{e^{uv}}}}{{{{\left( {{e^{uv}} - 1} \right)}^2}}}g\left( v \right)dv} \end{equation}

with $u=h/kT$ ($h$: Plank constant, $k$ Boltzmann constant, $T$: temperature) and $g(v)$ is the phonon density of states distribution.

My question is: Is it possible to calculate $g(v)$ using Density Functional Theory?

1. C. A. Tavares, et al. Solving ill-posed problems faster using fractional-order Hopfield neural network. J. Comp. Appl. Math. 381 112984 (2021) (DOI: 10.1016/j.cam.2020.112984)

Answer 1

Short answer. Yes, it is possible to calculate the phonon density of states using density functional theory. You can calculate the phonon frequencies on an arbitrarily large $\mathbf{q}$-point grid to construct the density of states, and most DFT codes will have the functionality to do this.

Longer answer. The density of states is given by:

$$ \tag{1} g(\omega)=\sum_{\nu}\int\frac{d\mathbf{q}}{(2\pi)^3}\delta(\omega-\omega_{\mathbf{q}\nu})\approx\frac{1}{N_{\mathbf{q}}}\sum_{\nu}\sum_{\mathbf{q}}\Delta(\omega-\omega_{\mathbf{q}\nu}), $$

where $\omega_{\mathbf{q}\nu}$ is the phonon frequency at wave vector $\mathbf{q}$ and branch $\nu$. The first equality is the analytical definition, and the second approximate equality is a practical expression in which the integral over the Brillouin zone is replaced by a sum over a discrete grid of $\mathbf{q}$-points, and the delta function is replaced by a function $\Delta$ of some width (e.g. a Gaussian).

Therefore, all you need for a calculation of the density of states is to have the phonon frequencies $\omega_{\mathbf{q}\nu}$ at many different $\mathbf{q}$-points for the sum to provide a converged result. The practical question then becomes calculating the phonon frequencies using DFT. To do that, you need to diagonalize the dynamical matrix $D(\mathbf{q})$ at that $\mathbf{q}$-point, whose eigenvalues are the squares of the frequencies. In turn, the dynamical matrix is the Fourier transform of the matrix of force constants, which roughly speaking measures the force that an atom feels when another atom moves. The entries in the matrix of force constants decay away as the distance between atoms increases, so in a simple real-space picture, if you have a large enough supercell, then the matrix of force constants is converged. Once you have a converged matrix of force constants, you can calculate the Fourier transform to build the dynamical matrix at any $\mathbf{q}$-point you want. I went over some detail about how to do this in this answer.

This procedure is relatively well-established and it is relatively easy to obtain well-converged results. Major DFT codes support the calculation of phonons, and as a simple post-processing step they also support the calculation of phonon densities of states.