We present an ab initio approach for grand-canonical ensembles in thermal equilibrium (eq) with local or nonlocal external potentials based on the one-reduced density matrix (1RDM). 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 1RDM. We further prove the existence of a Kohn-Sham system capable of reproducing the 1RDM 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.