Smearing techniques are widely used in first-principles calculations of metallic and magnetic materials, where they improve the accuracy of Brillouin zone sampling and lessen the impact of level-crossing instabilities. Smearing introduces a fictitious electronic temperature that smooths the discontinuities of the integrands; consequently, a corresponding fictitious entropic term arises and needs to be considered in the total free energy functional. Advanced smearing techniques – such as Methfessel-Paxton and cold smearing – have been introduced to guarantee that the system’s total free energy remains independent of the smearing temperature at least up to the second order. In doing so, they give rise to non-monotonic occupation functions (and, for Methfessel-Paxton, non-positive definite), which can result in the chemical potential not being uniquely defined. We explore this shortcoming in detail and introduce a numerical protocol utilizing Newton’s minimization method that is able to identify the desired Fermi energy. We validate the method by calculating the Fermi energy of ~20,000 materials and comparing it with the results of standard bisection approaches. In passing, we also highlight how traditional approaches, based on Fermi-Dirac or Gaussian smearing, are actually equivalent for all practical purposes, provided the smearing width is appropriately renormalized by a factor ~2.565. This data set contains the AiiDA databases, scripts, and other data necessary for the reproduction of our results.