Trajectory-based calculations are performed with the G-CTMQC code , where the algorithms presented in Section III B are implemented. In all cases, an ensemble of N traj = 1000 classical trajectories have been propagated starting from Wigner-sampled initial conditions. The Gaussian Wigner density in positions and momenta is determined as the Wigner-transform of the quantum nuclear density |χ S0 (ϕ, q, {x j }, t = 0)| 2 from Eq. ( ). Electronic and nuclear evolution equations are integrated with the Runge-Kutta-Gill algorithm and the velocity-Verlet algorithm, respectively, with a time step of 0.0024 fs (0.1 a.t.u.). In Section VI, the results are presented for simulations of duration 300 fs (for the study of the numerical costs of independent-and coupled-trajectory procedures when increasing the system size) and 1 ps (for the study of the effect of the environment on the system dynamics).