An ab initio variational grand-canonical electronic structure mean-field method, based on the Gibbs–Peierls–Bogoliubov minimum principle for the Gibbs free energy, is applied to the di-lithium (Li+Li) system at temperatures around T \approx 10,000 K and electronic chemical potential of μ \approx -0.1Eh. The method is an extension of the Hartree–Fock approach to finite temperatures. We first study the Li2 molecule at a frozen inter-nuclear distance of R = 3 \AA as a function of temperature. The mean-field electronic structure changes smoothly as temperature increases, up to 104 K, where a sharp spontaneous spin-polarization emerges as the variational mean-field solution. Further increase in the temperature extinguishes this polarization. We analyze the mean-field behavior using a correlated single-site Hubbard model and show it arises from an attempt of the mean-field to mimic the polarization of the spin–spin correlation function of the exact solution. Next, we keep constant the temperature at 104 K and examine the electronic structure as a function of inter-nuclear distance R. At R = 3.7 \AA, a crossing between two free energy states occurs: One state is “spin-unpolarized” (becomes lower in energy when R \ge 3.7 \AA), while the other is “spin polarized”. This crossing causes near-discontinuous jumps in calculated properties of the system and is associated with using the noninteracting electron character of our mean-field approach. Such problems will likely plague FT-DFT calculations as well. We use second-order perturbation theory (PT2) to study effects of electron correlation on the potential of mean force between the two colliding Li atoms. We find that PT2 correlation free energy at 104 K is larger than at 0 K and tends to restore the spin-polarized state as the lowest free energy solution.