Faussurier [Phys. Rev. E 65, 016403 (2001)] proposed to use a variational principle relying on Jensen-Feynman (or Gibbs-Bogoliubov) inequality in order to optimize the accounting for two-particle interactions in the calculation of canonical partition functions. It consists of a decomposition into a reference electron system and a first-order correction. The procedure appears to be very efficient in order to evaluate the free energy and the orbital populations. In this work, we present numerical applications of the method and propose to extend it using a reference energy which includes the interaction between two electrons inside a given orbital. This is possible, thanks to our efficient recursion relation for the calculation of partition functions. We also show that a linear reference energy, however, is usually sufficient to achieve a good precision and that the most promising way to improve the approach of Faussurier is to apply Jensen's inequality to a more convenient convex function.