The TSH approach belongs to the mixed quantum/classical class of methods and is one of the computationally most expedient way to include non-adiabatic effects. Nonetheless, TSH simulations require the evaluation of the nuclear forces from first-principles simulations, i.e., by solving the time-independent electronic Schrödinger equation, and those forces have to be evaluated for each nuclear geometry along the trajectory. Furthermore, due to its stochastic form (see section II A), the TSH approach requires to run a statistical ensemble of trajectories, which means that, in practice, several hundred thousands of geometries have to be considered. Solving the electronic Schrödinger equation is very computationally demanding and often have a very steep scaling with the number of electrons considered. These considerations limit considerably the applications of the TSH method. In particular, the three main limitations arise from: (i) the size of the systems that can be treated, (ii) the number of trajectories required to recover the proper statistical properties, and (iii) the duration of the physical processes that can be studied, i.e. the total length of the simulations.