  A Gaussian resolution method for the computation of equilibrium density
matrices rho(T) for a general multidimensional quantum problem is presented.
The variational principle applied to the ``imaginary time'' Schroedinger
equation provides the equations of motion for Gaussians in a resolution of
rho(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 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 also discussed.
