MD simulations rely on a critical albeit simple approximation: The nuclear motion is commonly propagated using Newtonian dynamics, with their driving forces being the negative gradient of the potential energy surface (PES) if no dissipative effects are considered. A PES is represented as a collection of solutions of the electronic Schrödinger equation for different nuclear conformations, resulting from the Born-Oppenheimer approximation. In on-the-fly (or direct) MD, which avoids the cumbersome precomputation of multidimensional surfaces, once the nuclear coordinates are updated, a new prediction of the PES and its gradient is necessary to drive the dynamics.