Over the last several decades, molecular dynamics (MD) has become one of the most important and commonly used approaches for studying condensed phase systems. MD calculations generally serve two often complementary purposes. First, an MD simulation can be used to study the dynamics of a system starting from particular initial conditions. Second, MD can be employed as a means of generating a collection of classical microscopic configurations in a particular equilibrium ensemble. The latter of these uses shows that MD is intimately connected with statistical mechanics and can serve as a computational tool for solving statistical mechanical problems. Indeed, even when MD is used to study a system’s dynamics, one never uses just a single trajectory (generated from a single initial condition). Dynamical properties in the linear response regime, computed according to the rules of statistical mechanics from time correlation functions, require an ensemble of trajectories starting from an equilibrium distribution of initial conditions. These points underscore the importance of having efficient and rigorous techniques capable of generating equilibrium distributions. Indeed while the problem of producing classical trajectories from a distribution of initial conditions is relatively straightforward — one simply integrates Hamilton’s equations of motion — the problem of generating the equilibrium distribution for a complex system is an immense challenge for which advanced sampling techniques are often required.