Figure : Flowchart of the simulation protocol used in this work. BSE@GW computations (1) are used as input to generate an LVC model for the entire diabatic PES (2). Wave packets are propagated on this model in TDH approximation (3) to compute the correlation matrix of coordinate expectation values (4). The clustering of the resulting adjacency matrix ( ) generates an ML-tree (6), which serves as the input for executing ML-MCTDH dynamics simulations to obtain diabatic state populations (7). the preceding layer. The resulting structure of the wave packets can be visualized using so-called ML-trees. Equations of motion for SPFs and coefficient vectors can be obtained by applying the Dirac-Frenkel time-dependent variational principle. Convergence can be controlled by inspecting the so-called natural orbital populations related to the SPFs. Below we will define thresholds for the largest population of the least occupied natural orbital, p NO .