Based on a generalization of Hohenberg-Kohns theorem, we propose a ground state theory for bosonic quantum systems. Since it involves the one-particle reduced density matrix $gamma$ as a natural variable but still recovers quantum correlations in an exact way it is particularly well-suited for the accurate description of Bose-Einstein condensates. As a proof of principle we study the building block of optical lattices. The solution of the underlying $v$-representability problem is found and its peculiar form identifies the constrained search formalism as the ideal starting point for constructing accurate functional approximations: The exact functionals for this $N$-boson Hubbard dimer and general Bogoliubov-approximated systems are determined. The respective gradient forces are found to diverge in the regime of Bose-Einstein condensation, $ abla_{gamma} mathcal{F} propto 1/sqrt{1-N_{mathrm{BEC}}/N}$, providing a natural explanation for the absence of complete BEC in nature.