All spectra presented in this work were calculated using the cumulant approach as implemented in our opensource spectroscopy Python package. The vertical excitation energies computed along the MD trajectories were used to construct classical autocorrelation functions, that, following Eqns. 11 and 8, yield the linear response function and thus the absorption spectrum for a given electronic transition. Nonadiabatic effects were assumed to be negligible (see Sec. II A), such that the full absorption spectrum was computed as the sum of the S 1 and the S 2 spectrum respectively.