The full-dimensional TSH dynamics simulations were carried out in exactly the same way as in refs. 26,28 only on different LVC potentials. We propagated the trajectories up to 1 ps using 0.5 fs and 0.005 fs time steps for the nuclear and electronic propagation, respectively. The TSH simulations were performed in in the diagonal basis (which is obtained by full diagonalization of the diabatic LVC potential energy matrix with the off-diagonal elements being the nuclear-coordinate dependent nonadiabatic couplings and the nuclear-coordinate independent SOCs). Initial conditions were sampled from a ground-state Wigner distribution and filtered using a stochastic algorithm based on a 0.1 eV wide energy window centered at the excitation energy to the lowestlying optically-active 1 MLCT state (calculated at the 1 GS equilibrium geometry) and taking into account the singlet oscillator strengths. This 1 GS → 1 MLCT excitation energy was also used to define the number of electronic states included in the LVC models and TSH simulations. For [Fe(bmip) 2 ] 2+ , all states 1 eV above the lowest optically-active 1 MLCT state at the 1 GS geometry were included in the calculations; the application of this criterion led to 14/16 singlets/triplets for RPBE, 12/15 for OPBE, 14/18 for BLYP, 16/18 for B3LYP, 14/21 for PBE0, 16/22 for CAM-B3LYP, and 14/21 for LC-BLYP. For [Fe(terpy) 2 ] 2+ , the utilization of CASPT2 SOCs limited the number of states as further states would require to expand our active space, which would mean infeasibility for CASPT2 calculations. Nevertheless, we could still include all states ca. 0.5 eV above the lowest optically-active 1 MLCT, which should still be sufficient. For [Fe(terpy) 2 ] 2+ , we included the lowest three quintet states for all DFT/TD-DFT methods and 13/18 singlets/triplets for RPBE, OPBE, and BLYP, 15/18 for B3LYP and PBE0, and 15/14 for LC-BLYP. For the triplet states of [Fe(terpy) 2 ] 2+ , calculated by CAM-B3LYP at the ground-state equilibrium geometry, we observed heavy 3 MLCT/ 3 MC mixing, which did not allow assignment of DFT/TD-DFT vs CASPT2 states; therefore, for [Fe(terpy) 2 ] 2+ , we did not carry out TSH simulations on CAM-B3LYP-LVC potentials. We note that in all other cases, the MLCT or MC electronic character could be clearly identified. In the TSH simulations, we employed the decoherence correction scheme of Granucci et al. 57 using a 0.1 a.u. decoherence parameter value. We transformed the adiabatic (spindiabatic) populations obtained from the TSH simulations to the diabatic basis of the given LVC model. All TSH dynamics simulations were carried out using the SHARC2.1 non-adibatic molecular dynamics code. 3 Results