Exact propagation has been carried out with the Chebychev scheme for which the evolution operator is expanded on Chebychev polynomials 60 with a time step 0.1 fs. The Hamiltonian has been normalized in order ensure its spectral range lies in the interval [-1, 1]. The grid size and the number of grid points have been carefully optimized for each Tully model depending on the value of the initial momentum k 0 , in order to avoid reflections at the boundaries and to accurately capture the fine features of the nuclear wavepacket everywhere in space, even in very delocalized situations. The computational parameters for each calculation are listed in Table . ElVibRot provides as output the nuclear wavepackets in the diabatic basis, as well as their spatial and time derivatives. This information is required by EFAC to reconstruct the TDPES (7) in the gauge where the TDVP is zero.