  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 $\sim$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 $\sim$2.565.
