Molecular dynamics simulations numerically solve the classical equations of motion. These simulations, though, assume that the forces on the electrons are derived solely by the ground adiabatic surface. Solving the time-dependent Schrödinger equation numerically incorporates all these effects, but is computationally unfeasible when the system has many degrees of freedom. To tackle this issue, one approach is the mean field or Ehrenfest method, where the molecular dynamics is run on the average potential energy surface given by a linear combination of the adiabatic states. This was applied successfully for some applications, but has some important limitations. When the difference between the adiabatic states is large, then the dynamics must be primarily driven by only one surface, and not an average potential. In addition, this method also violates the principle of microscopic reversibility.