The Gaussian wavepacket propagation method of Hellsing et al. [Chem. Phys. Lett. 122 (1985) 303] for the computation of equilibrium density matrices, is revisited and modified. The variational principle applied to the 'imaginary time' Schrodinger equation provides the equations of motion for Gaussians in a resolution of (ρ) over cap (T), described by their width matrix, center and scale factor, all treated as dynamical variables. The method is computationally very inexpensive, has favorable scaling with the system size and, with the current implementation, is surprisingly accurate in a wide temperature range, even for cases involving quantum tunneling. Incorporation of symmetry constraints, such as reflection or particle statistics, is discussed as well. (C) 2003 Published by Elsevier B.V.