The exact factorization offers a way to decompose the dynamical problem represented by the molecular time-dependent Schrödinger equation into coupled nuclear and electronic equations for the marginal and for the conditional amplitudes, respectively. The electrons effectively act on the nuclei as time-dependent scalar and vector potentials that completely determine the evolution of the nuclear wavefunction via a nuclear Hamiltonian, thus in the form of a (nuclear) time-dependent Schrödinger equation. These time-dependent potentials incorporate nonadiabatic effects beyond the Born-Oppenheimer approximation , and thus involve the ground and excited electronic states which are coupled through the nuclear motion. While only few attempts have been made to solve exactly the coupled equations for the marginal and for the conditional amplitudes, the presence of the time-dependent potentials in the nuclear Hamiltonian allows one to easily introduce a picture of the nuclear dynamics in terms of trajectories. Ciccotti and coworkers rigorously reformulated the nuclear equation using the method of characteristics, that can be interpreted as quantum trajectories in a Bohmian sense , but the main computational advantage has been achieved with a more classical-like treatment of the nuclear dynamics . Algorithms such as the coupled-trajectory mixed quantum-classical (CTMQC) scheme , along with its energyconserving variation CTMQC-E , and coupled-trajectory Tully surface hopping (CTTSH) allowed us to simulate nonadiabatic processes of various molecular models and molecular systems using an on-the-fly approach .