Gaussian resolutions for equilibrium density matrices

The Gaussian wavepacket propagation method of Hellsing et al. [Chem. Phys. Lett. 122 (1985) 303] for the computation of equilibrium density matrices ρ̂T is revisited and modified. The variational principle applied to the ‘imaginary time’ Schrödinger equation provides the equations of motion for Gaussians in a resolution of ρ̂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.