A new variational principle for optimizing thermal density matrices is introduced. As a first application, the variational many body density matrix is written as a determinant of one body density matrices, which are approximated by Gaussians with the mean, width and amplitude as variational parameters. The method is illustrated for the particle in an external field problem, the hydrogen molecule and dense hydrogen where the molecular, the dissociated and the plasma regime are described. Structural and thermodynamic properties (energy, equation of state and shock Hugoniot) are presented.