The numerical calculation of the TDPES has been carried out as follows: (1) exact vibronic wavepacket dynamics calculations have been performed with the ElVibRot program 54 on the Tully models; (2) the output of ElVibRot has been analyzed by the Exact Factorization Analysis Code (EFAC) to produce the TDPES. The TDPES can be easily expressed in terms of "standard" (a)diabatic nuclear amplitudes (and their spatial and time derivatives) arising from the Born-Huang representation of the molecular wavefunction. We refer the interested reader to Refs. [ 56-58] for a detailed discussion on those expressions. In any case, ElVibRot outputs the necessary information in the diabatic representation at various time steps throughout the propagation, and EFAC reads this information as input in order to reconstruct the TDPES. Additionally, EFAC imposes the gauge condition, which in the present case is simply taken to be A(x, t) = 0 ∀ t. Following from its definition given in Eq. ( ), and expressing the electronic wavefunction as the ratio of the molecular and nuclear wavefunctions, from Eq. ( ), the TDVP reads