We develop a systematic approach to construct energy functionals of the one-particle reduced density matrix (1RDM) for equilibrium systems at finite temperature. The starting point of our formulation is the grand potential Ω[G] regarded as variational functional of the Green’s function G based on diagrammatic many-body perturbation theory and for which we consider either the Klein or Luttinger–Ward form. By restricting the input Green’s function to be one-to-one related to a set on one-particle reduced density matrices (1RDM) this functional becomes a functional of the 1RDM. To establish the one-to-one mapping we use that, at any finite temperature and for a given 1RDM γ in a finite basis, there exists a non-interacting system with a spatially non-local potential v[γ] which reproduces the given 1RDM. The corresponding set of non-interacting Green’s functions defines the variational domain of the functional Ω. In the zero temperature limit we obtain an energy functional E[γ] which by minimisation yields an approximate ground state 1RDM and energy. As an application of the formalism we use the Klein and Luttinger–Ward functionals in the GW-approximation to compute the binding curve of a model hydrogen molecule using an extended Hubbard Hamiltonian. We compare further to the case in which we evaluate the functionals on a Hartree–Fock and a Kohn–Sham Green’s function. We find that the Luttinger–Ward version of the functionals performs the best and is able to reproduce energies close to the GW energy which corresponds to the stationary point.