The phonon density of states tells us how many phonon modes there are at every frequency. This can be calculated by choosing a uniform grid of $\vec{k}$ states in the first Brillouin zone and calculating the frequencies of the phonon modes for each $\vec{k}$ state. A histogram is then constructed out of these frequencies to plot the density of states.

For a simple cubic lattice, the Brillouin zone is a cube with $-\frac{\pi}{a} < k_x,\,k_y,\,k_z < \frac{\pi}{a}$. This cube can be divided into six pyramids, where the base of each pyramid is a face of the cube and the peak of the pyramid is the center of the cube. By symmetry, it is clear that the distribution of phonon modes in each of the pyramids must be the same; therefore, we only need to calculate the frequencies of the $\vec{k}$ states in one of the pyramids. The pyramids themselves have 4-fold symmetry, so we only need to calculate the frequencies in a quarter of a pyramid. The quarters of the pyramids have a remaining reflection symmetry about a plane that passes through a corner of the cube, the middle of the face of the cube, and the center of the cube. This is a polyhedron with four triangular faces and corners at $\Gamma$, $X$, $M$, and $R$. The frequencies only need to be calculated in the range $0 < k_x < \frac{\pi}{a}$, $0 < k_y < k_x$, $0 < k_z < k_y$. The code from the previous page was used to generate the density of states seen below.

$\frac{C_2}{C_1} = $  

$D(\omega)$
	$\large \sqrt{\frac{M}{C_1}}\omega$

There is still some numerical noise visible in the density of states because the $\vec{k}$ state grid is not fine enough. To get a smooth density of states, it would be necessary to calculate more points than is possible in the time it takes to load a web page. Notice that the number of numerical errors is listed. The code that calculates the eigenvalues of the matrix sometimes fails. This is common for simple routines that calculate eigenvalues. As long as there are not too many errors compared to the number of $\vec{k}$ states calculated, it should not have a large influence on the result.

Matlab code for performing this calculation more precisely can be found here.