We present an ab-initio approach for grand canonical ensembles in thermal equilibrium with local or nonlocal external potentials based on the one-reduced density matrix. We show that equilibrium properties of a grand canonical ensemble are determined uniquely by the eq-1RDM and establish a variational principle for the grand potential with respect to its one-reduced density matrix. We further prove the existence of a Kohn-Sham system capable of reproducing the one-reduced density matrix of an interacting system at finite temperature. Utilizing this Kohn-Sham system as an unperturbed system, we deduce a many-body approach to iteratively construct approximations to the correlation contribution of the grand potential.