Assessing the Performance of Approximate Density Functional Theory on 95 Experimentally Characterized Fe(II) Spin Crossover Complexes

: Spin crossover (SCO) complexes, which exhibit changes in spin state in response to external stimuli, have applications in molecular electronics and are challenging materials for computational design. We curate a data set of 95 Fe(II) SCO complexes (SCO-95) from the Cambridge Structural Database that have available low-and high-temperature crystal structures and, in most cases, confirmed experimental spin transition temperatures (T 1/2 ). We study these complexes using density functional theory (DFT) with thirty functionals spanning across multiple rungs of “Jacob’s ladder” to understand the effect of exchange-correlation functional on electronic and Gibbs free energies associated with spin crossover. We specifically assess the effect of varying the Hartree–Fock exchange fraction ( a HF ) in structures and properties within the B3LYP family of functionals. We identify three best-performing functionals, a modified version of B3LYP ( a HF = 0.10), M06-L, and TPSSh, that accurately predict SCO behavior for the majority of the complexes. Contrary to observations from prior studies, double-hybrids with higher a HF values are found to strongly stabilize high-spin states and therefore exhibit poor performance in predicting SCO behavior. Computationally predicted T 1/2 values are consistent among the three functionals but show limited correlation to experimentally reported T 1/2 values. These failures are attributed to the lack of crystal packing effects and counter-anions in the DFT calculations that would be needed to account for phenomena like hysteresis and two-step SCO behavior. The SCO-95 set thus presents opportunities for method development, both in terms of increasing model complexity and method fidelity.

To address this challenge, computational investigations 21,70,71 are essential to provide greater insight into how to tailor ligand designs to fine-tune SCO behavior.A prior study 72 developed accurate benchmarks for spin state ordering of TMCs using CASPT2 which were then used to evaluate optimal DFT functional choice for predicting spin state ordering.Correlated multiconfigurational methods like CASPT2 33,[73][74][75][76][77] and single-reference methods like CCSD(T) 78- 80 can provide reliable energies and geometries 81 of HS and LS states for Fe(II) complexes.It was also shown 82 that extended multistate (XMS) CASPT2 and DLPNO-CCSD(T) provide a good description of Mn(III) SCO complex geometries because they recover dispersion from first principles.However, they are computationally expensive for the large size of most SCO complexes (e.g., 100 atoms or more), particularly for optimizing geometries and computing vibrational frequencies, 46 which limits their applications to small-scale studies.In contrast, density functional theory (DFT) 83 comes with relatively lower computational cost, [84][85][86] but suffers from selfinteraction 87 or delocalization error [88][89][90][91] and static correlation error [91][92][93] .Because improving one of these errors often leads to a worsening of the other, it can be challenging to fine-tune DFT functionals to reproduce TMC properties.For instance, it was shown that predicted spectra and excited state dynamics of photoactive iron complexes obtained by tuning range-separated hybrid parameters to recover exact conditions do not necessarily have improved agreement with experimental data. 77However, results from range-separated hybrids with parameters tuned based on CASPT2 agree better with both CASPT2 and experiments. 77ere is generally not a widespread agreement about which DFT functionals perform best for SCO properties, even if we restrict ourselves to a fixed metal and oxidation state, such as Fe(II).
13 kcal/mol.Furthermore, studies where GGA exchange and meta-GGA exchange were tuned to reproduce spin splitting energies evaluated with CASPT2 showed that the sensitivity of spin state ordering to HF exchange strongly depends on ligand-field strength 105 .Several studies of SCO behavior have been carried out with DFT functionals, and these disagreed on which functionals lead to accurate predictions.Some of these studies found B3LYP 106 or a modified variant (aHF = 0.15) 107,108 to perform best in predicting spin-splitting energies of Fe(II) SCO complexes, while recommending the tuning of the exact exchange parameter to verify the validity of the functional against experiments or correlated wavefunction theory methods 108 .However, other studies working with different test sets comprising small Fe(II/III/IV) and Co(II) complexes, some of which show SCO behavior, found OPBE, 109,110 TPSSh, [111][112][113][114] M06-L, 114 or B2PLYP 115,116 to give the best results.In addition, recent small-scale studies using DFT-based methods on four homoleptic octahedral Fe(II) SCO model complexes with small monodentate ligands (NCH, NH3, H2O, CO) suggested that evaluation of non-empirical functionals on self-interaction-corrected densities 117 and methods such as locally scaled self-interaction correction 118 can mitigate density errors and improve accuracy of spin splitting energies.However, it is unclear if these observations can be extended to larger SCO complexes or SCO complexes with diverse ligands and ligand denticities. 117,118Given that prior studies are mostly small in scale and disagree about what functionals best predict SCO behavior, a significantly larger set of Fe(II) SCO complexes is necessary to determine which functionals make accurate SCO predictions.
Prior studies from our group have addressed the problem of predicting SCO behavior from different angles: using a genetic algorithm (GA) 119 and artificial neural networks (ANNs) trained on different amounts of HF exchange 120 , assessing consensus between a set of 23 density functional approximations (DFAs) 121 , and mining the literature 121 to predict SCO behavior in complexes.We also demonstrated 122 the utility of ANNs trained on hybrid DFT bond lengths, as these ANNs correctly assigned the majority of spin states in a set of 46 Fe(II) SCO complexes data-mined from literature.Following the text-mining protocol from the prior study, 122 we now curate an even larger data set of 95 Fe(II) SCO complexes from the Cambridge Structural Database 123 (CSD) with a diverse set of ligands and confirm the experimental observation of SCO behavior with a combination of literature and structural analysis that we refer to as SCO-95.We evaluate spin splitting energies between HS and LS states and transition temperatures using 23 DFAs across various rungs of "Jacob's ladder" 124 from semi-local generalized-gradient approximations to double hybrids along with seven variants of the B3LYP functional.We evaluate the role of thermodynamic corrections in selecting the best DFA based on spin splitting free energies.Ultimately, T1/2 is the most essential quantity to predict for SCO behavior, and we compare predictions of DFAs to T1/2 values curated from the literature for all our complexes, highlighting the challenges that remain for accurate prediction of SCO behavior by DFT.

2a. Data Set Curation.
We curated a data set of Fe(II) SCO complexes from the CSD 123 following a protocol used in prior work 122 .We performed an initial search of the CSD targeted at octahedral Fe(II) SCO systems via the Conquest graphical interface to the CSD, applied to the v5.41 data set with complexes from November 2019, March 2020, and May 2020 updates.We completed a Conquest query to identify structures containing hexa-coordinate Fe with no restrictions on connecting atom elements or bond orders.We applied filters to the molecular formula to ensure that only one Fe atom was present in a selected complex and the chemical name did not contain Fe(0), Fe(III), or Fe(IV) to select for cases with Fe(II).We excluded polymeric molecules and limited the search to those complexes with 3-dimensional (3D) coordinates determined.This initial search produced 13,028 compounds, each with a single unique refcode, i.e., a six-letter code which in some cases is followed by two more digits identifying additional structure determinations (Figure 1 and Supporting Information Table S1).This new set contained all refcodes reported in the SCO data set from prior work 122 .From this set of 13,028 refcodes, we retained 2,494 refcodes that corresponded to multiple (≥ 2) cases of matching six-letter patterns at the beginning of their refcodes (i.e., 828 sets with unique base refcodes), indicating that they were likely captured in multiple spin states (Figure 1 and Supporting Information Table S1).Using the CSD Python application programming interface (API) along with molSimplify 125 , we checked the geometries of these complexes to ensure that they were octahedral, which resulted in the elimination of 671 refcodes, i.e., 284 sets with unique base refcodes (Supporting Information Tables S1-S2).We then reviewed these sets of refcodes to confirm that at least one component in the crystal had a chemical name containing "iron(ii)" to further ensure Fe(II) complexes were selected, leading to removal of 273 refcodes, i.e., 90 sets with unique base refcodes (Supporting Information Table S1).We further examined all the retained refcodes to ensure that structures with identical base refcodes correspond to complexes with identical connectivity based on the determinant of the atomic-number-weighted connectivity matrix 122 .For each base refcode, complexes corresponding to the most frequently observed molecular graph determinant were selected, and refcodes corresponding to these complexes along with their reported temperatures in the CSD were saved (Figure 1 and Supporting Information Table S1).Consequently, all refcodes corresponding to a given base refcode that exhibit entirely different connectivities are filtered out, resulting in the elimination of 9 refcodes, i.e., 4 sets of unique base refcodes (Supporting Information Table S1).were identified.The data set was narrowed down based on the presence of Fe(II) and octahedral geometries, followed by abstract mining.Lastly, complexes were checked for charge assignment, bond length requirements, and presence of trustworthy H atoms. Representative complexes along with their refcodes are shown after each curation step as insets.Atoms in the insets are colored as: H in white, C in gray, N in blue, O in red, S in yellow, and Fe in brown.When two refcodes are discussed, their structures are overlaid and distinguished by coloring the carbon atoms of one structure in gray and the other in cyan.
For a given set of refcodes corresponding to the same complex from the CSD, the lowesttemperature and highest-temperature structures were recorded as the most likely LS and HS structures.For these complexes, the abstracts were mined 122 using the pybliometrics tool 126 /Scopus package and analyzed for the appearance of spin crossover keywords, lack of negative references to spin crossover, and positive VADER 127 sentiment over sentences containing the spin crossover keywords.From 450 sets of complexes with common base refcodes (1,541 total refcodes) that were candidate SCOs, only 258 sets of complexes with common base refcodes (1,099 total refcodes) showed positive sentiment and spin crossover keywords (Figure 1 and Supporting Information Table S1).To ensure that crystal structures in multiple spin states were indeed recorded in the CSD, we filtered the data to require complexes to have high-temperature (high-T) average equatorial bond lengths at least 0.02 Å greater than low-temperature (low-T) average equatorial bond lengths to remove less likely SCO complexes, removing 58 additional complexes, i.e., corresponding to 29 base refcodes.We further removed 196 complexes corresponding to 98 base refcodes without user-labelled charges or hydrogens that appeared incorrectly added (Figure 1 and Supporting Information Table S1).Lastly, when the same complex, as judged by an identical molecular graph determinant, was found with multiple base refcodes (i.e., 23 sets of base refcodes with 58 total refcodes), the base refcode for which that sum of the R factors for the low-T and high-T structures was the lowest was selected for further analysis, resulting in 96 unique Fe(II) SCO complexes (Figure 1 and Supporting Information Table S1).Structures of all 96 Fe(II) SCO complexes in pairs of low-and high-T crystallization are provided in the Supporting Information zip file.
Geometry optimizations were carried out for all candidate SCO complexes using hybrid (i.e., B3LYP [128][129][130] ) DFT in TeraChem 131 v1.9.All calculations were automated with molSimplify 125 and molSimplify automatic design 132 (mAD).We performed HF exchange resampling on all SCO complexes following the protocol in prior work 102 , in which we varied aHF from 0.00 to 0.30 in increments of 0.05 with the LDA/GGA exchange ratio fixed 99,120 .At each exchange fraction above or below the 0.20 value in B3LYP, the optimized geometry and converged wavefunction from the adjacent exchange fraction (e.g., 0.25 from 0.20 and 0.10 from 0.15) were used as starting points for subsequent optimization.All optimizations were carried out using the LACVP* basis set, which corresponds to the LANL2DZ 133 effective core potential for Fe and heavier elements such as I or Br, and the 6-31G* 134 basis for all remaining elements.
Geometry optimizations were carried out in Cartesian coordinates using L-BFGS algorithm, as implemented in DL-FIND 135 .The default thresholds of 1x10 -6 hartree for self-consistent field (SCF) convergence and 4.5x10 -4 hartree/bohr for the maximum gradient were used.Level shifting 136 of 0.25 Ha on all virtual orbitals was employed to aid with SCF convergence.All singlet states were calculated with a restricted formalism while the non-singlet states were calculated with an unrestricted formalism.For 95 out of 96 SCO complexes, calculations were performed in the singlet, triplet, and quintet spin states using low-T crystal structures in all cases.Despite the fact that we do not use the high-T crystal structure as a starting point for geometry optimizations, we make comparison of the final optimized structure bond lengths to the respective low-T and high-T bond lengths, making the curation of both low-T and high-T structures necessary.For one SCO complex (refcode: IXAQEF) with a radical ligand, calculations were performed in the doublet, quartet, and sextet spin states.However, we discard this SCO complex from analysis because the spin did not localize to the radical ligand equivalently in all spin states and refer to the final set of 95 SCO complexes as SCO-95.
We performed checks of the geometry 137 to determine the fidelity of the optimized geometries at different HF exchange fractions (Supporting Information Table S2).We also checked for significant spin contamination, which we define as cases where the expectation value of the S 2 operator deviates from the anticipated eigenvalue (i.e., S(S + 1)) by greater than or equal to 1 µB 2 .The largest spin deviation observed in our systems was small (i.e., 0.23 µB 2 ).We also checked that the spin on the metal in the SCO complexes was close to the expected value (i.e., the number of unpaired electrons on Fe is within one of the total number of unpaired spins) to ensure that the calculation converged to the correct electronic state.For systems with spin not localized to the metal, calculations were reconverged using wavefunctions from calculations with different HF exchange fractions.After this procedure, all calculations on the 95 SCO complexes were successful, with good geometries and suitable electronic structure.All geometry optimization calculations were performed in both the gas phase and the presence of implicit solvent.Implicit solvent calculations were carried out with the conductor-like polarizable continuum 138 solvation model 139,140 model (C-PCM) as implemented [139][140][141][142] in TeraChem 143 , with a dielectric constant of ε = 10.3 chosen to approximately model the effect of the crystal field.The solute cavity was built for all the elements using the defaults available in TeraChem, i.e., 1.2 x Bondi's van der Waals radii 144 .The gas-phase calculations resulted in some unreasonable geometries that failed the fidelity checks, and hence, we only consider structures obtained from geometry optimizations carried out in C-PCM (Supporting Information Table S3 and Figures S1-S3).Following the geometry optimizations, frequency calculations were carried out in TeraChem to obtain free energy corrections for a temperature of T = 300 K using the corresponding optimized geometries and converged wave functions.We chose T = 300 K instead of the temperature of crystallization to simplify the problem due to the large number of distinct crystallization temperatures and the weak dependence of the crossover temperature on the specific temperature chosen to calculate the correction.All initial and optimized structures are provided in the Supporting Information.
We carried out single-point energy calculations using 23 density functional approximations 121 (DFAs) with a larger, triple z def2-TZVP 145 basis set for all the SCO complexes (Supporting Information Table S4).These calculations were automated following a procedure from prior work 121 where molSimplify was interfaced with a developer version of Psi4 (1.4a2.dev723).In this workflow, we first carried out basis set projection from LACVP* to def2-SVP 145 for functionals with different HF exchange fractions, followed by single-point energy calculations with def2-SVP.The converged def2-SVP wavefunction molecular orbital coefficients from TeraChem single-point energy calculations were used as an initial guess for Psi4 def2-TZVP SCF calculations with varying HF exchange fractions in the functional.Self-consistent singlepoint energies for 22 DFAs were obtained by using the converged Psi4 B3LYP/def2-TZVP wavefunction as an initial guess using the recommended grid size 147,148 with 99 radial points and 590 spherical points.These calculations were run with a maximum of 100 SCF iterations to reach SCF convergence with convergence thresholds of 3.0 × 10 −5 Ha for both energy and density, and all calculations were successful.Single points with Psi4 were obtained in the gas phase due to incompatibility of the previously developed workflow 121 in Psi4 with implicit solvent models.We then computed spin splitting energies between HS and LS states using both electronic energies (DEH-L) as well as Gibbs free energies (DGH-L).Corrections to electronic energies to obtain Gibbs free energies were all obtained from modB3LYP (aHF = 0.10)/LACVP* thermodynamic corrections at T = 300 K. Differences were observed between DEH-L electronic energies obtained with LACVP* and def2-TZVP basis sets, which in turn affected the predictions of SCO complexes (Supporting Information Figure S4).Thus, DEH-L values for these 23 DFAs and modified B3LYP (modB3LYP) with varying HF exchange fractions were obtained with the larger def2-TZVP basis set.Nevertheless, we used the smaller basis set for the thermodynamic corrections in combination with def2-TZVP electronic energies due to the prohibitive computational cost of carrying out frequency calculations with the larger def2-TZVP basis set.Predictions of the spin-crossover temperature, T1/2, are obtained by identifying the temperature at which the HS and LS states have the same Gibbs free energies, which results in the expression, T1/2 = DHH-L/DSH-L.Here, the entropic contributions are assumed to be fairly constant with temperature and are evaluated from our T = 300 K thermodynamic corrections.
3. Results and Discussion.

Sensitivity of Structure and Spin-State Ordering to Hartree-Fock Exchange Fraction.
We curated the SCO-95 data set of 95 octahedral mononuclear Fe(II) SCO complexes from the CSD 123 , which are crystallized at low and high temperatures.A CSV file containing information about each SCO complex is provided in the Supporting Information ZIP file.Our data set consists primarily (41%) of bis tridentate ligands followed by hexadentate (ca.20%).Most complexes are homoleptic, and the vast majority (97%) coordinate Fe(II) with nitrogen, while a minority contain oxygen or sulfur coordinating atoms (Supporting Information Table S5 and Figure S5).
For the complexes in the curated data set, we first examine metal-ligand bond lengths in the experimental low-T and high-T X-ray crystal structures to confirm that they capture a change typical of a spin transition (i.e., LS to HS increase of ca.0.1-0.2Å).For this analysis, we focus on scaled bond lengths evaluated relative to the sum of covalent radii of each ligand connecting atom and iron, as done in prior work 122 , to account for the expected differences of Fe-S bond lengths compared to Fe-N and Fe-O bond lengths due to the larger size of S. A spin transition is evident in the crystal structures, with high-T structures typically exhibiting longer scaled bond lengths by ca.0.07 in comparison to low-T structures, i.e., we observe an increase in bond lengths of ca.0.16 Å (Figure 2 and Supporting Information Table S6).There are some differences within each spin state by ligand denticity, with shorter scaled bond lengths (ca.by 0.02) for SCO complexes with ligands of higher denticities, and longer metal-ligand bond lengths (ca.1.02) are observed for some hexakis monodentate complexes (Figure 2 and Supporting Information Table S7).Outliers are also observed predominantly for bis tridentate cases such as a complex (refcode: QACKAJ 149 ) that has relatively similar average metal-ligand bond lengths (ca.2.13 Å and 2.16 Å) in low-T and high-T structures (Figure 2 and Supporting Information Table S7).The bond length filter employed during data set curation does not screen out complexes with relatively comparable bond lengths since the cut off we use for bond length difference is only > 0.02 Å.Nevertheless, low-T and high-T structures are assigned LS and HS states, respectively through magnetic susceptibility experiments and Mössbauer spectroscopy 149 , and thus the rigidity of the tridentate ligands in this example likely prevents significant bond length changes with spin state.We next analyzed the change in scaled metal-ligand bond lengths with change in spin state observed in the DFT-optimized geometries of these SCOs (Figure 2).The B3LYP scaled metal- ligand bond lengths are shorter in LS states than in HS states by ca.0.09 for all SCO compounds (Figure 2 and Supporting Information Table S6).The B3LYP-optimized HS scaled bond lengths are typically longer than their experimental high-T counterparts, and the LS scaled bond lengths are also longer than the low-T experimental bond lengths for most (87) structures, although there are exceptions (Figure 2 and Supporting Information Table S6).For the 8 SCOs where B3LYPpredicted LS scaled bond lengths are shorter than the low-T experimental values, they correspond to cases where low-T experimental bond lengths were longer than typically expected values (i.e., significantly above 0.95 relative bond lengths 122 , Figure 2 and Supporting Information Table S6).
Closer examination reveals that for one complex (refcode: NIGXUY) 150 the observed low-T structure bond lengths at 100 K correspond to an intermediate-spin (IS) state as indicated by the difference between the low-T and high-T X-ray structure bond lengths, which is much smaller than expected between LS and HS states (a little over 0.02 Å vs a typical value of 0.20 Å), despite the lack of rigidifying elements such as multidentate ligands.The protocol used in data set curation assumes that the lowest-temperature structure from the CSD corresponds to the LS state.However, for this complex, the lowest-temperature structure has been obtained at 100 K and actually corresponds to the IS state, although our workflow protocol led us to assumed it was in the LS state.Nevertheless, magnetic susceptibility experiments confirm that the spin transition captured for this complex is from LS to HS. 150 For the other complexes, the low-T structure LS state assignments were validated by magnetic and spectroscopic data (e.g., refcode QACKAJ 149 , described earlier).Finally, in one case (refcode: VIFNAC 151 ), the structure was obtained at an unusually high temperature (200 K), which likely contributes to its longer experimental low-T bond length than those obtained in the 0 K B3LYP geometry optimization.
Next, we analyzed changes in calculated scaled metal-ligand bond lengths as a function of aHF in the modB3LYP functional 99 (Figure 3).Upon increasing aHF, we observe an increase in predicted metal-ligand bond lengths for LS states and, to a lesser degree, HS states (Figure 3 and Supporting Information Table S8 and Figures S6-S8).Because B3LYP (aHF = 0.20) overestimated experimental low-T bond lengths, the aHF = 0.00 structures are in best agreement with experiment (Figure 3 and Supporting Information Table S8 and Figure S6).In contrast, there is no optimal aHF for improving agreement between DFT of HS states and experimental high-T structures (Supporting Information Table S8 and Figures S7-S8).9][90][91] LS states are more delocalized than HS states and more sensitive to HF exchange, consistent with observations from prior studies 99,102,122 .Thus, although cancellation of errors means that aHF = 0.00 can slightly improve upon standard B3LYP for the prediction of structural properties for LS states, we expect non-zero HF exchange admixture to be essential to describe the electronic structure.
Overall, the bond lengths at all HF exchange fractions are either in very good or at least acceptable agreement with experiment except in select outlier cases.-T) crystal structures and modB3LYP low-spin (LS) optimized geometries.The bin width in these plots is 0.01.The solid gray lines represent zero-axes in all the plots.A representative SCO complex (REFCODE: BAXJUI) is shown in the inset.Hydrogen, carbon, nitrogen, and iron are shown in white, gray, blue, and brown, respectively.A red star near a bar in each pane indicates the average scaled bond length difference between low-T and LS structures of the inset complex.
Following our analysis of bond lengths, we studied adiabatic spin splitting energies, (i.e., the differences in energy of DFT-optimized HS and LS states), as determined from the electronic energies, DEH-L, and with thermodynamically-corrected, Gibbs free energies, DGH-L (Figure 4).work 100,105 (Figure 4 and Supporting Information Figures S9-S11).We assign a prediction from a DFA as correct if it predicts the HS-LS splitting of a complex to be between -5.5 and +5.5 kcal/mol as all of these complexes have experimentally observed SCO behavior.Analysis of DEH-L values obtained at different levels of aHF reveals that modB3LYP with aHF = 0.15 correctly predicts 97% of the systems (i.e., 92 out of 95) to be SCO complexes based on this criterion (Supporting Information Figures S9-S10 and Table S9).However, when the more experimentally relevant DGH-L is considered, the modB3LYP functional with a lower aHF (0.10) correctly predicts the highest number, 93% of the systems (i.e., 88 out of 95), to be SCO compounds, at odds with some earlier work on representative SCO complexes 108,154 that favored a higher aHF = 0.15 value (Figure 4 and Supporting Information Figure S11 and Table S9).The computed DGH-L values for the SCOs exhibit a few consistent outliers at all values of aHF (refcodes: ECODIM 155 , WIHQIQ 156 , YAGYUB 157 , and ZERDOS 158 ) of which some (refcodes: WIHQIQ and YAGYUB) have a slightly stronger preference for HS states while the others (refcodes: ECODIM and ZERDOS) tend to prefer LS states (Figure 4 and Supporting Information Figure S11).However, experimental SCO behavior was reported for crystal structures of all these identified outliers [155][156][157][158] .It is worth noting the diverse chemical composition of these outliers: the homoleptic SCO WIHQIQ has six Fe-S bonds and the ligands of homoleptic YAGYUB contain iodine, while both heteroleptic ECODIM and ZERDOS complexes contain boron (Supporting Information Figure S12).These more varied compositions could cause these predictions to benefit less from error cancellation in HF-exchangetuned functionals in comparison to more conventional Fe(II)/N complexes.The bin width in these plots is 2.5 kcal/mol.The solid gray lines represent zeroaxes in all the plots.A representative SCO complex (REFCODE: ACAHUJ) is shown in the inset and the bin that contains its property value is marked with a red star in each pane.Hydrogen, carbon, nitrogen, oxygen, and iron are shown in white, gray, blue, red, and brown, respectively.
Examination of the free energy contributions reveals that average free energy corrections to DEH-L are comparable for all aHF values, ca.-6.0 kcal/mol (Supporting Information Table S9 and Figure S13).Zero-point energy (ZPE), thermal vibrational energy, electronic (Selec), rotational  S10).For the overall free energy corrections, a few outliers are observed with values more than -3 kcal/mol or less than -10 kcal/mol, strongly stabilizing LS and HS states, respectively (Supporting Information Figure S13).The outliers with free energy corrections < -10 kcal/mol are always hexakis monodentate complexes while those with free energy corrections more than -3 kcal/mol mostly have ligands with higher denticities (Supporting Information Figure S13 and Table S11).Free energy corrections to DEH-L are significant enough to alter which complexes are predicted to display SCO behavior, but they are not strongly dependent on the functional or, except for outlier cases, on the molecule.Thus, two potential strategies for rapid SCO screening would be to use predictions of electronic energies at a higher exchange fraction (i.e., aHF = 0.15 108,154 ) or using a lower exchange fraction (i.e., aHF = 0.10) with a heuristic free energy correction.Nevertheless, given the presence of outliers in the free energy corrections, full characterization of the vibrational contribution is warranted to assess candidate SCOs.

Spin Splitting Energy Predictions from 23 DFAs.
Beyond the HF exchange fraction, other aspects of the exchange-correlation functional can be expected to be important for identifying SCO behavior.Thus, we evaluated both DEH-L and DGH-L using 23 DFAs distributed across the various rungs of "Jacob's ladder" 124 , following the protocol from prior work 121 (see Computational Details).This set includes three generalized gradient approximation (GGA), four hybrid GGA, four meta-GGA, five meta-GGA hybrid, two range-separated, long-range corrected hybrid, and five double-hybrid DFAs (Supporting Information Table S4).Except for GGAs and meta-GGAs, all other DFAs in the 23 DFA set have aHF varying from 0.10 to 0.71, and MP2 correlation is incorporated over a similarly wide range (between 0.125 and 1.000) in the double-hybrid DFAs (Supporting Information Table S4).
Besides the modB3LYP (aHF = 0.10) functional identified in Sec.S9 and S12).Nevertheless, uncertainties due to the lack of explicit treatment of solvent or crystal environments motivate our broader definition.Using the more stringent requirement to predict canonical SCO behavior, we observe that M06-L (i.e., 78 out of 95 SCOs) performs the best in comparison to TPSSh or modB3LYP (aHF = 0.10) because most of the complexes are predicted to be LS at 0 K with an HS state that is stabilized at room temperature (Figure 5 and Supporting Information Table S12 and Figure S23).

M06-L TPSSh
with its property value annotated on each histogram by a red star.Hydrogen, carbon, nitrogen, sulfur, and iron are shown in white, gray, blue, yellow, and brown, respectively.
The remaining functionals with low fractions of exchange, all GGAs (i.e., BLYP 161,162 , BP86 163,164 , and PBE 165 ) and all other meta-GGAs (i.e., TPSS 113 and SCAN 166 ) except the Minnesota meta-GGAs (i.e., M06-L 160 and MN15-L 167 ) stabilize low-spin (LS) states with the average of their DGH-L distributions ranging between 10.1 and 15.7 kcal/mol, predicting SCO behavior for only a few (≤ 7) complexes (Figure 5 and Supporting Information Figures S17-S18 and Table S12).There is nothing obviously different about the few complexes predicted to be SCOs by GGAs or meta-GGAs, as most are representative of the full set (i.e., are Fe(II)/N), except for one Fe(II)/S complex (refcode: WIHQIQ 156 ), and all are experimentally verified SCO complexes with either single-step (five cases) or two-step (two cases) SCO behavior (Supporting Information Table S13 and Figure S24).
We previously noted the good performance of M06-L.This functional on average stabilizes HS states slightly (average DGH-L = -2.6 kcal/mol) but predicts SCO behavior for ca.84% (80 out of 95) of the complexes, on par with modB3LYP with aHF = 0.10 (93% or 88 of 95, Figures 4 and   5 and Supporting Information Table S12).In contrast, MN15-L, a more recently developed Minnesota functional, strongly stabilizes HS states with an average DGH-L of -24.8 kcal/mol, thus not predicting SCO behavior for any of the complexes (Supporting Information Figure S18 and Table S12).The difference in predictions between local meta-GGA Minnesota functionals could be a result of the different parametrization of MN15-L and M06-L (58 vs 34 parameters).The M06-L 160 parametrization was done with data that included bond energies of transition metal (TM) diatomics, and MN15-L 167 used additional datasets, including TM datasets for open-shell TM dimer bond energies and TM atomic excitation energies.Neither parameterization dataset was strongly weighted towards coordination complexes with mid-row transition metals, nor is it obvious why the MN15-L parameterization should be so strongly HS biased.
Next, we examine the DGH-L predicted by the hybrid counterparts of GGA and meta-GGA DFAs more closely, despite their inferior performance with regard to SCO behavior.All Becke three-parameter GGA hybrids, i.e., B3LYP [128][129][130] , B3P86 129,163 , and B3PW91 129,168 , which have aHF of 0.20, strongly favor HS states with average DGH-L between -6.7 to -10.8 kcal/mol (Supporting Information Figure S19 and Table S12).This HS bias is indeed largely attributable to HF exchange fraction, as the higher value (0.25) in PBE0 169 causes it to predict SCO behavior for only one complex (Supporting Information Figure S19 and Table S12).Despite this poor picture for Gibbs free energies, with DGH-L values predicting SCO behavior for very few complexes and B3P86 performing best (30% or 29 of 95), DEH-L values of these DFAs predict SCO behavior for up to 92% of the complexes (Supporting Information Table S12).For the meta-GGA hybrids, we had noted good performance of TPSSh (85% SCO or 81 out of 95), and it slightly stabilizes LS states (average DGH-L = 2.0 kcal/mol) on average (Figure 5 and Supporting Information Table S12).Examination of the remaining meta-GGA hybrids reveals that all of them (i.e., SCAN0 170 , M06 171 , M06-2X 171 , and MN15 172 ) strongly favor HS states and predict few (≤14) SCO complexes, likely due to their higher aHF values (between 0.25 and 0.54, Figure 5 and Supporting Information Figure S20 and Table S12).This reinforces the observation 99,108,154 that the aHF value incorporated in a functional has an overriding effect in predicting the SCO behavior of complexes.
The remaining classes of range-separated hybrids and double hybrids are in principle the highest rungs considered on "Jacob's ladder" in this work.Despite their promise for functional tuning to recover exact conditions, examination of range-separated hybrids reveals that both LRC-wPBEh 173 and wB97X 174 favor HS states and predict only a few systems (≤ 2) to be SCO complexes (Supporting Information Figure S21 and Table S12).All double-hybrids strongly favor HS states (average DGH-L between -39.6 to -59.8 kcal/mol) and thus, never predict any systems to be SCO compounds, in contrast to prior work 175 which showed double-hybrids to predict electronic spin splitting energies accurately (Supporting Information Figure S22 and Table S12).In particular, PBE0-DH 176 was shown 175 to predict accurate DEH-L for 16 small homoleptic transition metal complexes (TMCs) that do not show SCO behavior and consist of Fe(II/III) and Co(II/III) metals with ligands of varying strengths: CO, NCH, NH3, and H2O.The accuracy of these gasphase PBE0-DH/def2-TZVPP 177 calculations in this study 175 was determined through comparison with diffusion Monte Carlo studies as opposed to experimental data.However, in our work with 95 Fe(II) SCO complexes, PBE0-DH strongly favors HS states and fails to predict SCO behavior as shown by both DEH-L and DGH-L (Supporting Information Figure S22 and Table S11).The double-hybrids studied here 159,176,178 , like most double-hybrids, have very high aHF values (between 0.50 and 0.71 across DFAs), which could explain the strong preference towards HS states (Supporting Information Figure S22 and Table S12).This highlights the limitations of higher-rung DFAs in predicting SCO behavior and shows that caution must be exercised in determining the appropriate DFAs for SCO predictions.Now, we focus on the consistency of the predictions among the three best-performing functionals.Comparison of DGH-L predicted by modB3LYP (aHF = 0.10), TPSSh, and M06-L functionals reveals a stronger correlation between DGH-L predicted by modB3LYP (aHF = 0.10) and TPSSh than between modB3LYP (aHF = 0.10) and M06-L (Pearson's r of 0.92 vs 0.76; Figure 5).For most complexes where SCO behavior is predicted by modB3LYP (aHF = 0.10) but not by TPSSh or M06-L, the DGH-L predicted by these functionals differ by ca. 3 kcal/mol or less, with TPSSh stabilizing LS states and M06-L HS states slightly more relative to modB3LYP (Supporting Information Table S14).Closer examination reveals that for one complex with a hexadentate nitrogen-coordinating ligand (refcode: OSABOB 179 ), TPSSh strongly favors an LS ground state compared to modB3LYP (aHF = 0.10), i.e., DGH-L of 11.9 kcal/mol vs 2.1 kcal/mol (Supporting Information Table S14 and Figure S25).This complex has been experimentally 179 determined to exhibit SCO behavior with the aid of variable-temperature NMR spectroscopy, UV/Vis spectroscopy, magnetic susceptibility measurements and differential scanning calorimetry, suggesting this is an example of a limitation of TPSSh or the modelling approach.We also observe that all complexes where DGH-L predicted by M06-L and modB3LYP (aHF = 0.10) differ by a larger amount (refcodes: JAFBUO 180 , JIQDET 181 , QACKAJ 149 , SESNEM 49 , TEPVIW 182 , and WIHQIQ 156 ) have ligands that are rigid with higher denticities (half contain tridentates, Supporting Information Table S14 and Figure S26).In these cases, M06-L strongly stabilizes the HS states relative to modB3LYP and fails to predict SCO behavior.All of these complexes are experimentally shown to exhibit SCO behavior through analysis of X-ray crystal structures and/or magnetic susceptibility measurements.Thus, these molecules represent challenging test cases for further DFA tuning and development.
While there are a few cases not predicted to be SCO complexes by one or more of these three functionals, the DGH-L values are still typically within ca. 2 kcal/mol of the defined range of DGH-L between -5.5 and 5.5 kcal/mol (Figures 4 and 5 and Supporting Information Tables S14-S15, Text S1, and Figures S25-S27).For example, two complexes with refcodes QAXQIR 183 and YAGYUB 157 are not predicted to be SCO complexes by all three functionals (Figures 4 and 5 and Supporting Information Table S14 and Figure S25).However, crystal structure analysis combined with results from high-resolution powder-diffraction data demonstrate gradual and complete SCO behavior of YAGYUB 157 .QAXQIR 183 is also experimentally determined to be an SCO complex based on analysis of X-ray crystal structures, Mössbauer studies, UV/Vis, and NMR spectroscopy data.While SCO behavior is predicted by all three functionals for a majority of the complexes (68 out of 95 systems), the predicted DGH-L values vary quite a bit between these functionals (Supporting Information Figure S27 and Table S15).For example, while a complex (REFCODE: IBUVUW 184 ) is predicted to show SCO behavior by both modB3LYP (aHF = 0.10) and TPSSh, they differ significantly in the DGH-L value (i.e., -0.6 kcal/mol vs 4.1 kcal/mol, Supporting Information Figures S25 and S27 and Table S15).Similarly, OSABOB 179 , which was an outlier for TPSSh, is more consistently predicted as an SCO from both modB3LYP (aHF = 0.10) and M06-L, but they have a large (i.e., 6 kcal/mol) difference in DGH-L values (2.1 kcal/mol vs -4.2 kcal/mol; Supporting Information Figures S25 and S27 and Table S15).Thus, even among top-performing functionals that predict a large number of complexes to have SCO behavior, quantitative SCO free energies may differ, and there remain outliers that are experimentally validated as SCO but that cannot be captured in a one-size-fits-all DFA choice.

Transition Temperature (T1/2) of SCO Complexes.
Beyond qualitative identification of SCO complexes using an energy cutoff, we quantify SCO behavior with T1/2 values predicted by the 23 DFAs and modB3LYP functionals in comparison to experimentally reported T1/2 values.The T1/2 values corresponding to single-step SCO behavior are reported only for 76 of the complexes in SCO-95, while the remaining systems either exhibit two-step SCO behavior or a single-step SCO event without a reported T1/2 in the manuscript (Figure 6 and Supporting Information Figure S28 and Table S16).In order to make a comparison between computationally predicted T1/2 values and experimentally reported T1/2 values for more systems, we also attempted to estimate experimental T1/2 as the average of low-T and high-T of crystallization, but the presence of outliers motivates our focus on only the 76 complexes for which T1/2 is directly reported (Figure 6 and Supporting Information Tables S16  Next, we compare T1/2 values predicted from the 23 DFAs and modB3LYP functionals with varying aHF against those reported experimentally.As could be expected, despite experimentally observed SCO behavior, if DGH-L is not predicted to be within 5 kcal/mol, then T1/2 values will not be reasonable.Thus, only modB3LYP (aHF = 0.10), M06-L, and TPSSh predict reasonable T1/2 values (i.e., between 0 and 1000 K, Figure 7 and Supporting Information Table S18).Nevertheless, the T1/2 values predicted by modB3LYP (aHF = 0.10), TPSSh, and M06-L also show limited correlation with experimentally reported values, with Pearson's r of ca.0.36 for all three functionals (Figure 7 and Supporting Information Figures S32-S33).While all three functionals show similar degree of correlation to experimentally reported values, closer examination reveals that the complexes for which DFA and experimental predictions are consistent (i.e., within < 10 K) vary between the functionals (Figure 7 and Supporting Information Figures We next aimed to identify if the errors in DFA prediction of T1/2 could be attributed to any specific ligand chemistry or SCO behavior.We find that there are complexes for which DFApredicted T1/2 values are considerably larger than their specific experimentally reported values (by 100 K or more) but within the range of overall experimentally reported T1/2 values, i.e., 100 K < T1/2 < 450 K (Supporting Information Table S20).For example, there are three hexakis monodentate complexes (refcodes: GOGSAZ 192 , NIGXUY 150 , and YAGYIP 157 ) for which DFApredicted T1/2 values of all three functionals are larger than experimentally reported values, but examination of the manuscripts for these complexes indicates they all exhibit gradual, single-step SCO behavior 150,157,192 (Supporting Information Table S20 and Figure S35).
There are similar complexes for which DFA-predicted T1/2 values are considerably smaller than experimental T1/2 values, i.e., by 100 K or more (Supporting Information Table S20).For another hexakis monodentate complex (refcode: YAGYUB 157 ), all three functionals predict T1/2 values of 80 K or less while the experimentally reported T1/2 is 205 K (Supporting Information Table S20 and Figure S35).In this case, there are intermolecular N•••H-C and I•••H-C HB interactions reported in the experimental manuscript of this complex, where it is suggested that these interactions are weaker at the low T of crystallization (90 K) but become stronger at room temperature 157 , an interaction we cannot capture in modeling of isolated complexes.We also note intermolecular interactions reported for other complexes with refcodes SACTOG 58 (HBs between imidazoles) and YAZJIT 193 (N-H•••p interactions between the TMC cation and its counter anions) for which DFA-predicted T1/2 values by M06-L and TPSSh are both smaller than experimental T1/2 (Supporting Information Table S20 and Figure S35).Steep SCO behavior observed in these complexes is attributed to the presence of intermolecular interactions 193 .Thus, overall, disagreement between DFA-predicted T1/2 values and experimentally reported T1/2 values could be attributed at least in part due to the lack of explicit solvent 53,194 or crystal packing 63,158,195 as well as the absence of noncovalent interactions 196 between the complex and its counterions 63,193,197,198 in the DFT calculations that are present in the crystal structures.Another source of discrepancy would be hysteresis during SCO 196,[199][200][201][202] or two-step SCO behavior 49,53,155,184,199,203,204 , both of which are not accounted for by the modelling protocol based In summary, we curated a set of 95 Fe(II) SCO complexes (SCO-95) from the CSD that had both available low-and high-temperature crystal structures and were validated to be SCOs using sentiment analysis.We studied these complexes using the modB3LYP functional with varying HF exchange fraction as well as 23 DFAs across multiple rungs of "Jacob's ladder" to identify which functionals accurately predict SCO behavior of these complexes.Because SCOs are characterized by bond elongation from LS to HS states, we also validated that the DFT functionals could reproduce this transition.The B3LYP-predicted bond lengths were typically longer than experimental bond lengths, and pure GGA (i.e., aHF = 0.00) bond lengths had the best agreement with crystal structure geometries.Nevertheless, optimized geometries obtained with other aHF fractions were also in acceptable agreement with experiments.
Consistent with previous studies, 100,105 we found that increasing aHF in modB3LYP led to a stabilization of HS states over LS states whether evaluated with electronic or Gibbs free energies.
While analysis of DEH-L indicated modB3LYP (aHF = 0.15) correctly predicts SCO behavior for 97% of the complexes, for the experimentally relevant DGH-L, modB3LYP with a lower aHF (0.10) better predicted SCO behavior, correctly assigning 93% of the complexes.While the Gibbs free energy contributions shift DEH-L to favor HS states by ca.-6 kcal/mol on an average, individual corrections vary significantly for some complexes, emphasizing the importance of evaluating these contributions for every complex.
Of the additional functionals distributed across the rungs of "Jacob's ladder", only M06-L and TPSSh DFAs predicted SCO behavior for most complexes to the same degree as modB3LYP (aHF = 0.10).While TPSSh has the same exchange fraction as in our tuned functional, M06-L likely predicts SCOs accurately thanks to its more extensive parameterization, including on some openshell transition metal molecules.At odds with prior observations 175 , the high HF exchange fraction in the double hybrids studied here led to strong stabilization of HS states and few predicted SCOs.
Analysis of experimentally characterized SCOs allowed us to extract experimental transition temperatures directly from the manuscript.Comparison of computationally predicted T1/2 to experimentally reported T1/2 revealed that while T1/2 predictions from modB3LYP (aHF = 0.10), M06-L, and TPSSh show limited correlation to experimental T1/2, they are largely consistent among the three DFAs.Thus, while these three DFAs predict SCO behavior for most complexes, they do not necessarily predict the same relative SCO temperatures as experiment.We also identified the reasons for the largest disagreement of DFT predictions of SCO behavior with experiments to be the lack of counter-anions and crystal packing effects in the DFT calculations, as well as the lack of packing effects to predict phenomena like hysteresis and two-step SCO behavior.
Overall, our study of Fe(II) SCO complexes with a variety of DFAs highlights that caution must be exercised while using DFT to study SCO complexes because the predictions depend strongly on the functional, with the amount of HF exchange and parametrization of the functional playing a crucial role.Moreover, accurate prediction of SCO behavior, e.g., by predicting experimental T1/2 values curated in our SCO-95 set, presents an outstanding challenge for the development of both more sophisticated modeling protocols to incorporate packing effects as well as for developing improved wavefunction theory and DFT energies to predict single-step SCO.

Corresponding Author
Supporting Information.List of refcodes retained and eliminated during data set curation; geometry criteria for fidelity evaluation of optimized geometries; list of SCO complexes for which gas-phase geometries failed fidelity evaluation; X-ray crystal structures of gas-phase optimized SCO complexes which failed fidelity evaluation; RMSD of optimized SCO complexes which failed fidelity evaluation; summary of 23 density functional approximations; comparison of LACVP* and def2-TZVP spin splitting energies; ligand types, coordinating atoms, and denticities of SCO complexes; representative SCO complexes with N, S, and O-ligand coordinating atoms; average metal-ligand bond lengths across different ligand denticities; average metal-ligand bond lengths of DFT geometries and crystal structures; average metal-ligand bond length differences; histograms of average metal-ligand bond length differences; histograms of average scaled metalligand bond length differences; histograms of average metal-ligand bond length differences; histograms of DEH-L at aHF = 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30; histograms of DGH-L at aHF = 0.05, 0.15, 0.25; average and standard deviation of modB3LYP DEH-L and modB3LYP DGH-L; Xray structures of SCO complexes that are outliers on DGH-L histograms; histograms of DGH-L corrections, ZPE, thermal vibrational energy, and vibrational entropy at aHF = 0.00, 0.10, 0.20, 0.30; average ZPE, thermal vibrational energy, vibrational and rotational entropy; outliers observed in free energy corrections at different aHF values; histograms of DGH-L of GGA, meta-GGA, GGA hybrid, meta-GGA hybrid, range-separated hybrid, ans double hybrid DFAs; average and standard deviation of DEH-L and DGH-L obtained with 23 DFAs; list of complexes with GGA and meta-GGA predicted SCO behavior; X-ray crystal structures of GGA and meta-GGA predicted SCO complexes; histograms of DEH-L obtained with M06-L and TPSSh; SCO complexes not predicted by modB3LYP (aHF = 0.10), M06-L, TPSSh; X-ray crystal structures of complexes with conflicting SCO predictions; percentile ranks of DGH-L to quantify changes in SCO predictions; plots of percentile ranks of DGH-L and T1/2 values; list of SCO complexes with significantly different DGH-L percentile ranks; plots of experimentally reported and estimated experimental T1/2 values; complexes that exhibit two-step SCO behavior or lack experimental T1/2; mean and standard deviation of experimental low-T, high-T, T1/2; histograms of low-T and high-T of crystallization; SCO complexes with unusual low-T and high-T of crystallization; SCO complexes with identical reported and estimated T1/2; T1/2 predictions by modB3LYP fuunctionals and 23 DFAs; plots of experimental T1/2 vs T1/2 predicted by modB3LYP, M06-L, TPSSh; plots of estimated T1/2 vs T1/2 predicted by modB3LYP, M06-L, TPSSh; list of complexes with comparable experimental and DFA-predicted T1/2; SCO complexes with comparable experimental and DFA-predicted T1/2; list of complexes and X-ray structures of complexes where experimental and DFA-predicted T1/2 differ considerably; plots of experimental T1/2 and average DFA-predicted

Figure 1 .
Figure 1.Description of data set curation and filtering steps involved to obtain unique mononuclear Fe(II) octahedral SCO complexes.From an initial data set of 13,028 refcodes obtained from Conquest search, 2,494 refcodes with multiple instances of the same base refcode 96 SCO pairs after removing duplicates

Figure 2 .
Figure 2. Experimental vs B3LYP average scaled bond lengths for 95 Fe(II) spin crossover (SCO) complexes in SCO-95.SCO complexes corresponding to high-temperature (high-T) and high-spin (HS) are shown as red circles while low-temperature (low-T) and low-spin (LS) complexes are shown as blue circles.1D histograms of each property are shown at top and right with bin widths of 0.01.The solid gray line represents the parity line.A representative SCO complex (REFCODE: VIFNAC) is shown in the inset.Hydrogen, carbon, nitrogen, chlorine, and iron are shown in white, gray, blue, green, and brown, respectively.

Figure 3 .
Figure 3. Histograms of average (avg) scaled metal-ligand bond length errors (bl err) at Hartree-Fock exchange (a HF ) values: (a) a HF = 0.00, (b) a HF = 0.10, (c) a HF = 0.20, and (d) a HF = 0.30.Bond length errors are computed between low-temperature (low-T) crystal structures and modB3LYP low-spin (LS) optimized geometries.The bin width in these plots is 0.01.The solid gray lines represent zero-axes in all the plots.A representative SCO complex (REFCODE: BAXJUI) is shown in the inset.Hydrogen, carbon, nitrogen, and iron are shown in white, gray, blue, and brown, respectively.A red star near a bar in each pane indicates the average scaled bond length difference between low-T and LS structures of the inset complex.
Both the DEH-L and DGH-L distributions indicate an overall change in preferred ground state from LS to HS with increasing fraction of aHF in the modB3LYP functional, consistent with prior

Figure 4 .
Figure 4. Histograms of spin splitting Gibbs free energies for 95 SCO complexes curated from the CSD in SCO-95 (DG H-L , in kcal/mol) at different Hartree-Fock exchange (a HF ) values: (a) a HF = 0.00, (b) a HF = 0.10, (c) a HF = 0.20, and (d) a HF = 0.30, in the modB3LYP functional with def2-TZVP basis set.The bin width in these plots is 2.5 kcal/mol.The solid gray lines represent zeroaxes in all the plots.A representative SCO complex (REFCODE: ACAHUJ) is shown in the inset and the bin that contains its property value is marked with a red star in each pane.Hydrogen, carbon, nitrogen, oxygen, and iron are shown in white, gray, blue, red, and brown, respectively.

(
Srot), and vibrational entropy (Svib) components contribute towards DGH-L, with very limited dependence on aHF, and T*Svib contributes the most towards DGH-L at ca. -4.0 kcal/mol (Supporting Information Figures S14-S16 and Table

Figure 6 .
Figure 6.Histogram of experimentally (expt) reported transition temperature (T1/2) for N = 76 complexes.A representative low-T, LS Fe(II) SCO complex (refcode: BOJLUI) is indicated as a red star in the histogram.The bin width of the histogram is 25 K. Hydrogen, carbon, nitrogen, and iron are shown in white, gray, blue, and brown, respectively.

Figure 7 .
Figure 7. Plots of transition temperature (T1/2) reported experimentally vs T1/2 predicted by (a) modB3LYP (aHF = 0.10) and (b) M06-L, and plots of T1/2) predicted by modB3LYP (aHF = 0.10) vs (c) M06-L and (d) TPSSh.T1/2 is experimentally reported for 76 complexes while computational T1/2 predictions are available for all 95 complexes.The number of data points is specified in each plot.Extreme outliers that are truncated here are shown in SI figures.Representative Fe(II) SCO complexes (refcodes: OSABOB and VIFNAC) are colored as green and blue squares.The parity line is shown in solid gray in all the plots.Hydrogen, carbon, nitrogen, sulfur, chlorine, and iron are shown in white, gray, blue, yellow, green, and brown, respectively.
113, only the M06-L160and TPSSh113DFAs predict the majority of the curated complexes to be SCO materials on the basis of spin slitting Gibbs free energies (Figures4 and 5and Supporting Information FiguresS17-S22).While TPSSh has the same exchange fraction as in our tuned functional (aHF = 0.10), M06-L likely predicts SCOs accurately thanks to its more extensive parameterization, including on some open-shell transition metal molecules.It is worth noting that for canonical SCO behavior, we require an LS state be predicted as the ground state and entropic contributions at raised temperatures to stabilize the HS state (Figures4 and 5and Supporting Information Figures S9 andS23 and Tables T1/2 with and without outliers; plots of estimated experimental T1/2 and average DFA-predicted T1/2 with and without outliers; Pearson's r values of experimental T1/2 vs average DFA-predicted T1/2 plots; study of average T1/2 values; plots of modB3LYP-predicted T1/2 vs T1/2 predicted by M06-L and TPSSh; list of complexes where DFA-predicted T1/2 values are extreme outliers; percentile ranks of T1/2 values to quantify changes in SCO predictions (PDF) X-ray crystal structures of Fe(II) SCO complexes crystallized at low-and high-temperatures of crystallization; implicit solvent optimized geometries of Fe(II) SCO complexes in high-spin quintet and low-spin singlet states with modB3LYP functionals; details of SCO complexes including refcodes, ligand denticities, low-and high-temperatures of crystallization, experimentally reported transition temperatures, estimated experimental transition temperatures, charges, CSD DOI, R factors, ligand connecting atoms, and metal-ligand bond lengths; electronic and Gibbs free energy spin splitting energies, and Gibbs free energy corrections of Fe(II) SCO complexes obtained with modB3LYP functionals and 23 DFAs; transition temperatures of Fe(II) SCO complexes predicted by modB3LYP functionals and 23 DFAs (ZIP)