An approximate method for integrating the Liouville equation to obtain the time-dependent classical phase space density distribution at constant energy or temperature is presented. The density distribution of each degree of freedom is represented by a single Gaussian phase packet (GPP) whose center and width obey variationally optimized equations of motion. The constant energy dynamics is applied to the calculation of equilibrium thermodynamic averages for a Lennard-Jones cluster and fluid to demonstrate the feasibility and utility of this approximate method for the simulation of many-body condensed phase systems. The rate of kinetic energy equipartitioning is examined for GPP dynamics using a generalization of the ergodic measure and found to be significantly faster than for standard molecular dynamics simulation. A global optimization algorithm is developed based on simulated annealing of the phase space density distribution. This method is applied to the global energy minimization of Lennard-Jones clusters and found to be superior to simulated annealing methods employing classical point particles.