  We introduce a method to obtain the specific heat of quantum impurity models
via a direct calculation of the impurity internal energy requiring only the
evaluation of local quantities within a single numerical renormalization group
(NRG) calculation for the total system. For the Anderson impurity model, we
show that the impurity internal energy can be expressed as a sum of purely
local static correlation functions and a term that involves also the impurity
Green function. The temperature dependence of the latter can be neglected in
many cases, thereby allowing the impurity specific heat, $C_{\rm imp}$, to be
calculated accurately from local static correlation functions; specifically via
$C_{\rm imp}=\frac{\partial E_{\rm ionic}}{\partial T} + 1/2\frac{\partial
E_{\rm hyb}}{\partial T}$, where $E_{\rm ionic}$ and $E_{\rm hyb}$ are the
energies of the (embedded) impurity and the hybridization energy, respectively.
The term involving the Green function can also be evaluated in cases where its
temperature dependence is non-negligible, adding an extra term to $C_{\rm
imp}$. For the non-degenerate Anderson impurity model, we show by comparison
with exact Bethe ansatz calculations that the results recover accurately both
the Kondo induced peak in the specific heat at low temperatures as well as the
high temperature peak due to the resonant level. The approach applies to
multiorbital and multichannel Anderson impurity models with arbitrary local
Coulomb interactions. An application to the Ohmic two state system and the
anisotropic Kondo model is also given, with comparisons to Bethe ansatz
calculations. The new approach could also be of interest within other impurity
solvers, e.g., within quantum Monte Carlo techniques.
