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. This variational density matrix is used to calculate properties of the hot, dense hydrogen. First results for a density of rs=2 are presented.