I find that there already exists an option for allowing fractional occupation number (FON) to accelerate SCF, so it should be handy to extend it to finite temperature DFT by allowing FON for the converged wavefunction.

Specifically, the FONs with Fermi-Dirac distribution smearing Gaussian smearing are most commonly used.