Should all q-points in the first Brillouin zone be sampled in phonon density of states calculations?

Question

For electronic density of states (DoS) calculations with VASP, it is recommended that VASP's internal symmetry routine is turned off by setting ISYM = 0. This way, all k-points in the Brillouin zone are sampled. Should the same be done with q-points while calculating phonon DoS?

Answer 1

Background. The phonon density of states $g$ 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$ is the energy and $\omega_{\mathbf{q}\nu}$ the energy of a phonon of wave vector $\mathbf{q}$ and branch $\nu$. In the first equality the integral is over the full Brillouin zone. The second approximate equality gives a practical expression for the density of states used in numerical calculations, where the integral over the Brillouin zone is replaced by a sum over a discrete set of $N_{\mathbf{q}}$ $\mathbf{q}$-points, and these points should cover the Brillouin zone uniformly. The $\Delta$ function is a narrow function that peaks at $\omega_{\mathbf{q}\nu}$ (e.g. a Gaussian) and replaces the Dirac $\delta$ function in the numerical calculation.

Brillouin zone sampling. The integral or sum cover the entire Brillouin zone, so you need $\mathbf{q}$-points across the entire zone. However, you only need to calculate explicitly the $\mathbf{q}$-points in the irreducible Brillouin zone (whose size depends on the symmetry of the system), because all other points can be constructed by applying the symmetry operations without the need to perform additional calculations. Having said this, the Fourier interpolation that is used to construct the dynamical matrix at arbitrary $\mathbf{q}$-points after the coarse $\mathbf{q}$-point grid has been fully converged is very cheap in computational terms. Therefore, it would not be a problem to brute-force sample the entire Brillouin zone at the Fourier interpolation step.

VASP. I use my own code to calculate phonons starting from the forces calculated from VASP (my code works in a similar fashion to Phonopy but exploits nondiagonal supercell to dramatically reduce the computational cost). Symmetry can be fully expoited in the VASP calculations for the forces, and on top of that only the phonons at the $\mathbf{q}$-points in the irreducible Brillouin zone need to be calculated. I am writing this about my code to show that in principle one can use full symmetry when using VASP as an engine for the force calculations in phonon studies. However, I cannot say what other phonon codes like Phonopy do in practice and whether in that case you can also use full symmetry.