A detailed comparison has been made of the performance of molecular dynamics and hybrid Monte Carlo simulation algorithms for calculating thermodynamic properties of 2D Lennard‐Jonesium. The hybrid Monte Carlo simulation required an order of magnitude fewer steps than the molecular dynamics simulation to calculate reproducible values of the specific heat. The ergodicity of the two algorithms was compared via the use of intermediate scattering functions. For classical systems the intermediate scattering functions should be real; however, a simple analysis demonstrates that this function will have a significant imaginary component when ergodicity breaks down. For q vectors near the zone boundary, the scattering functions are real for both algorithms. However, for q vectors near the zone center (i.e., harmonic, weakly coupled modes), the scattering function calculated via molecular dynamics had a significantly larger imaginary component than that calculated using hybrid Monte Carlo. Therefore, the hybrid Monte Carlo algorithm is more ergodic and samples phase space more efficiently than molecular dynamics for simulations of 2D Lennard‐Jonesium.