A novel method for simulating the statistical mechanics of molecular systems in which both nuclear and electronic degrees of freedom are treated quantum mechanically is presented. The scheme combines a path integral description of the nuclear variables with a first-principles adiabatic description of the electronic structure. The electronic problem is solved for the ground state within a density functional approach, with the electronic orbitals expanded in a localized (Gaussian) basis set. The discretized path integral is computed by a METROPOLIS Monte Carlo sampling technique on the normal modes of the isomorphic ring polymer. An effective short-lime action correct to order tau(4) is used. The validity and performance of the method are tested by studying two small lithium clusters, namely Li-4 and Li-5. Structural and electronic properties computed within this fully quantum-mechanical scheme are presented and compared to those obtained within the classical nuclei approximation. Quantum delocalization effects turn out to be significant as shown by the fact that quantum simulation results at 50 K approximately correspond to those of classical simulations carried out at 150 K. The scaling factor depends, however, on the specific physical property, thus evidencing the different character of quantum and thermal correlations. Tunneling turns out to be irrelevant in the temperature range investigated (50-200 K).