  Based on a generalization of Hohenberg-Kohn's 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, $\nabla_{\gamma} \mathcal{F} \propto
1/\sqrt{1-N_{\mathrm{BEC}}/N}$, providing a natural explanation for the absence
of complete BEC in nature.
