For condensed phase systems, a reasonable approximation of quantum time correlation functions is obtained by running classical trajectories sampled initially on the multidimensional equilibrium Wigner distribution. In this context, we present a new efficient method to generate such a distribution using Gaussian phase packets propagation in imaginary time starting from a high-temperature classical Maxwell-Boltzman distribution. The method is illustrated for the Wigner equilibrium distribution and the vibrational correlation function of a proton in a model H-bond. (c) 2006 Elsevier B.V. All rights reserved.