Analytic Hessians are often viewed as essential for the calculation of accurate harmonic frequencies, but the implementation of analytic second derivatives is nontrivial and solution of the requisite coupled‐perturbed equations engenders a sizable memory footprint for large systems, given that these equations are not required for energy and gradient calculations in density functional theory. Here, we benchmark the alternative approach to harmonic frequencies based on finite differences of analytic first derivatives, a procedure that is amenable to large‐scale parallelization. Not only for absolute frequencies but also for isotopic and conformer‐dependent frequency shifts in flexible molecules, we find that the finite‐difference approach exhibits mean errors < 0.1 cm⁻¹ as compared to results based on an analytic Hessian. For very small frequencies corresponding to nonbonded vibrations in noncovalent complexes (for which the harmonic approximation is questionable anyway), the finite‐difference error can be larger, but even in these cases the errors can be reduced below 0.1 cm⁻¹ by judicious choice of the displacement step size and a higher‐order finite‐difference approach. The surprising accuracy and robustness of the finite‐difference results suggests that availability of the analytic Hessian is not so important in today's era of commodity processors that are readily available in large numbers. © 2017 Wiley Periodicals, Inc.