In our previous TSH/LVC works, we employed the hybrid B3LYP* exchangecorrelation functional, which we selected based on our quantum chemical benchmarks and comparison to experimental reference data. We herein test the performance of several DFT/TD-DFT methods representing different rungs of Jacob's ladder: GGAs -RPBE, OPBE, BLYP, global hybrids -B3LYP and PBE0, and range-separated hybrids -CAM-B3LYP and LC-BLYP. In all DFT/TD-DFT computations, we used the TZVP basis set. Solvation effects were not taken into account as they do not affect significantly the MLCT-MC energetics for the investigated Fe complexes. Two-electron integrals were approximated by using resolution of identity approaches. In all TD-DFT calculations, we utilized the Tammm-Dancoff approximation (TDA). The spin-orbit couplings (SOCs) were calculated at the 1 GS equilibrium geometry and were taken to be the same at any nuclear geometry. In the case of [Fe(bmip) 2 ] 2+ , singlet and triplet states were included in the LVC models and the SOCs were calculated with the given DFT/TD-DFT method (at the 1 GS geometry optimized by the same level of DFT). For [Fe(terpy) 2 ] 2+ , however, the situation is more complicated as in addition the singlet and triplet states, quintets also have to be included. We employ our hybrid approach, according to which the LVC PESs are obtained by DFT/TD-DFT (using a restricted singlet reference for singlets and triplets but an unresricted quintet reference for quintet states) but the SOCs are calculated by multiconfigurational second-order perturbation theory (CASPT2). The CASPT2 SOC matrix for [Fe(terpy) 2 ] 2+ was taken from ref. 26 and was thus identical in all simulations (reindexing was though necessary as state ordering often differed for different TD-DFT methods).