In this work, we present a super-polynomial improvement in the precision scaling of quantum simulations for coupled classical-quantum systems. Such systems are found, for example, in molecular dynamics simulations within the Born-Oppenheimer approximation. By employing a framework based on the Koopman-von Neumann formulation of classical mechanics, we express the Liouville equation of motion as unitary dynamics and utilize phase kickback from a dynamical quantum simulation to calculate the quantum forces acting on classical particles. This approach allows us to simulate the dynamics of these classical particles without the overheads associated with measuring gradients and solving the equations of motion on a classical computer, resulting in a super-polynomial advantage at the price of increased space complexity. We demonstrate that these simulations can be performed in both microcanonical and canonical ensembles, enabling the estimation of thermodynamic properties from the prepared probability density.