As a consequence of Bloch’s theorem, the numerical computation of the fermionic ground state density matrices and energies of periodic Schrödinger operators involves integrals over the Brillouin zone. These integrals are difficult to compute numerically in metals due to discontinuities in the integrand. We perform an error analysis of several widely-used quadrature rules and smearing methods for Brillouin zone integration. We precisely identify the assumptions implicit in these methods and rigorously prove error bounds. Numerical results for two-dimensional periodic systems are also provided. Our results shed light on the properties of these numerical schemes, and provide guidance as to the appropriate choice of numerical parameters.