Benchmarking Density Functional Methods for Calculation of State Energies of First Row Spin-Crossover Molecules.

A systematic study of the performance of several density functional methodologies to study spin-crossover (SCO) on first row transition metal complexes is reported. All functionals have been tested against several mononuclear systems containing first row transition metal complexes and exhibiting spin-crossover. Among the tested functionals, the hybrid meta-GGA functional TPSSh with a triple-ζ basis set including polarization functions on all atoms provides the best results across different metals and oxidation states, and its performance in both predicting the correct ground state and the right energy window for SCO to occur is quite satisfactory. The effect of some additional contributions, such as zero-point energies, relativistic effects, and intramolecular dispersion interactions, has been analyzed. The reported strategy thus expands the use of the TPSSh functional to other metals and oxidation states other than FeII, making it the method of choice to study SCO in first row transition metal complexes. Additionally, the presented results validate the potential use of the TPSSh functional for virtual screening of new molecules with SCO, or its use in the study of the electronic structure of such systems.


Introduction
One of the major outcomes of ligand-field and molecular orbital theories is the fact that a transition metal ion surrounded by a given set of ligands in a certain coordination motif must split their d-based molecular orbitals into different subsets. 1,2 In this context, spin-crossover (SCO) molecules represent a family of transition metal complexes in which the metal has electronic states with different spins but similar electronic energies, therefore allowing the molecule to access two alternative spin-configurations. 3 Since they were first reported, nearly 85 years ago by Cambi and coworkers, 4 spin-crossover molecules have been the focus of intense research in chemistry and physics, mostly due to their intrinsic technological applications. 3,[5][6][7][8][9][10][11][12] In fact, SCO molecules are molecular level switches in which is possible to manipulate the electronic structure of the metal center by means of an external stimuli. It is precisely this behavior what makes SCO systems perfect candidates for molecular memory storage systems, nanoscale and spintronic devices. [13][14][15][16][17][18][19][20] In fact, over the last years SCO molecules have been used in the field of metal-organic frameworks (MOFs) to generate spincrossover frameworks (SCOFs), in which the physical properties of the material can be selectively controlled by different guest molecules, thus making these materials perfect candidates for sensing applications. [21][22][23][24][25][26][27][28][29] More recently it has been shown that SCO molecules can be used as junctions in molecular transport devices, and that the conductivity through the molecule is strongly affected by the spin state of the metal center, turning the SCO molecule into a spin-filter. [30][31][32] Regardless of its obvious technological appeal, the rational design of SCO molecules with tailored properties remains an elusive task. Empirical approaches have been proposed through the years, from early models based on ligand field parameters, 33 to the statistical treatment of the available data for Fe II bidentate systems, 34 which allows to predict the However, DFT is known to struggle when it comes to describe energy differences arising from different spin-state multiplicities with an additional geometry change. 43 Processes that involve major rearrangement of the occupied orbitals (such as spin-crossover) are very sensitive to the description of the exchange correlation energy, [44][45][46] and hybrid functionals usually have a bias towards the high-spin state, while pure functionals tend to overestimate the energy of the low-spin state. [47][48][49] Major attempts to develop a systematic methodology to study spin-state energies in transition metal complexes using DFT methods have been made through the last years. The B3LYP* method, 48,49 derived from the original B3LYP functional but adjusting the amount of Hartree-Fock exchange to 15%, provides with an accurate description of the magnetic properties in the [Fe(NHS4)]L family (NHS4 = 2,2'-bis(2-mercaptophenylthio)diethylamino, L = CO, NO + , PR3, NH3 and N2H4) and has also been successfully used to study the spintransition in the [Fe(phen)2(NCS)2] complex. 48,49 The OPTX 50 exchange functional with the LYP 51 and PBE 52 correlation functionals has also proved to be a useful tool to study spinstate energetics in Fe II systems. Double hybrid functionals also have been used, 53 but their computational cost and strong convergence problems in systems with transition metal atoms makes them less appealing as a potential choice to screen for new SCO systems. Finally, since 2009, several studies have shown that the hybrid meta-GGA functional TPSSh, 54,55 containing 10% of exact exchange Hartree-Fock, provides with quite accurate energy gaps between spin-states for many Fe II transition metal complexes, including dinuclear systems, allowing also for the calculation of the corresponding transition temperature, in very good agreement with the experimental data. [56][57][58][59][60] Very recently, it has also been shown that many DFT methods can improve their performance towards spin-state energetics in iron systems (Fe II and Fe III ) by mixing between 10 and 17% of Hartree-Fock exchange into the functional, but the optimal amount for each individual functional still requires specific parameterizations. 61 From all of the above, major efforts have been made in modeling spin-crossover processes in Fe II metal containing systems (see reference 32 and references within), which is of course the most common situation for spin-crossover to occur, but not the only one. In fact, SCO can also occur for other electron configurations, typically from d 4 to d 7 , and for a given number of d electrons, several metals with different oxidation states may be considered.
Although the hybrid meta-GGA functional TPSSh seems the correct choice when modeling spin-crossover processes on Fe II systems, no systematic study of the performance of this functional towards molecules with other metals and oxidation states with SCO behavior has been done. In this study, the TPSSh functional is used to calculate the spin-state energy gap in a benchmark set of 20  CAM-B3LYP 66 and wB97x 67 ) in order to study its applicability in the modeling of spincrossover processes for first-row transition metal complexes.

Methodology
As indicated above, the spin-crossover phenomenon appears in systems for which states with different spin have similar electronic energies. In such cases, entropy favors the high-spin state and spin-crossover can appear, with the low-spin state being observed at low temperature, and the high-spin state at high temperature. 3 From this point of view, SCO can be described as a thermodynamic equilibrium between the low-spin (LS) and the high-spin (HS) state.
LS D HS [1] Assuming an ideal system, in which the spin-crossover system is isolated, and due to the fact that the pressure-dependent term to the free enthalpy change is usually small, we can write the Gibbs energy change (ΔG) for the equilibrium expression of [1] as, 68,69 [2] where [3] is the Gibbs free energy associated with spin state i at a given temperature. In Equation 3, the enthalpy term (H i ) includes both electronic ( ) and vibrational ( ) contributions.
For molecular complexes, can be properly estimated by using the harmonic approximation, while the term , describing the electronic energy of spin state i, can be obtained directly from ab initio calculations. The entropy contribution to the free energy (S i ) can also be estimated using the harmonic approximation. Since at equilibrium DG vanishes, the transition temperature (T1/2) can be determined from the changes in enthalpy and entropy using Equation 4. [4] The equilibrium condition defined in Equation 2 can also be expressed in terms of the equilibrium constant, Keq, which in turn can be written in terms of the molar fraction of each spin state, [5] Equation 5 can be recast so the molar fraction of high-spin state depends on the change in the free energy at each temperature as, [6] where R is the gas constant and T the temperature. It is worth mentioning that both DH and DS have some temperature dependence that can be, in principle, neglected without losing accuracy. 70,71 From the above expressions, one must observe that the key term in Equation 3 is the electronic energy change ( ) contribution to the enthalpy change. This term represents the difference on electronic energies between the high-spin and low-spin geometries at its equilibrium points, and the choice of ab initio method to calculate this term will be the critical point in the correct prediction of the spin-crossover phenomena with DFT methods.

Computational details
All Density Functional Theory calculations have been carried out with Gaussian 09 (revision D) using a 10 −8 convergence criterion for the density matrix elements. 72 The fully optimized contracted triple-ζ (TZV) all electron Gaussian basis set developed by Ahlrichs and coworkers was employed for all the elements, and a quadruple-ζ basis set with polarization functions (QZVP) was used on the metal center (Basis Set 1, abbreviated BS1). 73 For the TPSSh and OPBE functionals results, additional calculations with a larger basis set, a tripleζ basis set including polarization functions (TZVP) for all atoms, were also done (Basis Set 2, abbreviated BS2). A systematic study on the effect of increasing the size of the basis set on the metal center was done for selected systems (S7, S8 and S12, corresponding to Mn II , Fe II and Co II systems), which shown that there is very little effect of the basis set size on the metal center with the energy differences (see Table S1 in SI). Scalar relativistic effects were included using the Douglas-Kroll-Hess Hamiltonian (DKH) 74-77 on the optimized geometries at the TPSSh/BS2 level.  Table 1).

Results
To evaluate the general performance of DFT methods towards accurate calculation of spinstate energetics in spin-crossover systems, we retrieved a benchmark group of 20 mononuclear complexes (see Figure 1 and Table 1) that exhibit spin-crossover behavior with electron configurations d 4 to d 7 . Such systems were selected using as criteria the fact that both, crystallographic and magnetic data (i.e, transition temperature T1/2) was accurately determined. The chemical formula, together with the transition temperature and experimental information for the data set can be found in Table 1, and a schematic drawing of all studied systems is shown in Figure 1.
For all the selected systems, geometry optimizations where performed in both accessible spin-states using the B3LYP, B3LYP*, CAM-B3LYP (CB3LYP), OPBE, TPSSh, M06, M06L, and wB97xD functionals using BS1 (see Computational details section). The double hybrid B2PLYP functional, which already has been used to study SCO systems, 53,60,61,97 was also initially considered, but its computational cost and initial performance for some selected systems made us withdraw it from the test set (see Table S2 on SI).
Among all the tested functionals (see Table 2), TPSSh clearly outperforms the other methods given that is the only one that correctly predicts the ground state (low-spin) for all selected systems, regardless of the metal, number of d electrons and oxidation state. It is also worth mentioning that the OPBE functional seems to perform really well for challenging sandwich type systems of general formula [M(Cp R )2] (see Tables 1 and S2 on SI). However, it fails among the selected d 4 and d 5 molecules, and also with some of the Fe II d 6 systems, thus making it less universal in terms of studying spin-crossover in transition metal complexes.
B3LYP* also works fine, particularly when combined with BS1, but similarly, in this case also the problem seems to be the Mn II d 5 complexes and the Cr II d 4 tested molecules, which again constrains a little bit its use to certain systems. All the other functionals seem to clearly over stabilize the high-spin state, thus leading to a situation in which the high-spin state is more stable than the low-spin state, and therefore the transition can never occur. In Figure 2, the energy differences between both spin states (high-spin -low-spin) have been plotted against the amount of exact exchange Hartree-Fock mixed into the corresponding functional. From the graphic and the data presented in Table 1, we can see that pure functionals such as OPBE and M06L don't do a bad job when it comes to compute spinstate energy differences in SCO systems, but their performance is not satisfactory for such property (see Table S5 on SI). By mixing some amount of HF into the functionals, the performance towards the description of electronic energy differences in spin-crossover systems improves. Following this idea, the proposed B3LYP* functional, 48,49 that lowers the amount of exact exchange Hartree-Fock mixed into the functional to 15% was also tested.
The functional does a reasonable job, particularly for Fe II and Co II systems, and its overall performance is similar to OPBE, but seems to be less general than the TPSSh functional (10% of HF) in terms of the first transition metal row. However, mixing too much HF into the functional can lead to the opposite effect. Above 20% of HF mixing, the stabilization of the high-spin state is such that it becomes the ground state for almost all the studied systems, ΔE"(HS&LS)/kcal.mol &1 " %HF OPBE" TPSSh" B3LYP" M06" M06L" B3LYP*" CB3LYP" that being the case of the B3LYP and M06 functionals, including 20% and 27% of exact exchange Hartree-Fock respectively.
Given the fact that TPSSh seems to be the functional that performs the best among the first row of transition metals, we decided to compute the corresponding transition temperatures for each system. To do that, vibrational frequencies for all optimized systems using the TPSSh functional were calculated. Vibrational frequencies can be accurately calculated, and allows us to apply harmonic approximation to thermally correct the thermochemical quantities. In particular, DG(T) can be approximated using equation [7], [7] where G i can be computed using the harmonic frequencies (nj) for the spin state i and the following expression, [8] In [8], kB is the Boltzman constant, E i the electronic energy of the spin-state i, h the Planck constant and T the temperature. By plotting DG(T) vs. T and fitting the data, a good approximation of the transition temperature (T1/2) can be obtained.
At this point, we used two different basis sets (see Computational details), trying to balance accuracy and computational cost. As can be seen from the data presented in Table 3, the combination of the TPSSh functional with the BS2 triple-z basis set with polarization from Ahlrichs and co-workers provides a good description of the ground state for the whole set of systems. The mean absolute error (MAE) between computed and observed transition temperatures using this methodology (TPSSh/BS2) is 3.70 kcal/mol on the overall data set, which is an excellent result for such type of calculations. Obviously, this error translates into 324 K, which makes the TPSSh functional less useful to quantitatively predict transition temperatures, even though for some systems the accuracy is remarkable (see Table 3). Also, although the TPSSh functional has been previously used to study trends in a given family of SCO compounds, [56][57][58] we couldn't find any correlation between the different families grouped by number of d electrons and the experimental and computed T1/2. Recently, it has been reported that the inclusion of zero-point energy correction, dispersion, and relativistic effects can play a role in the accurate calculation of thermochemical quantities in Fe II spin-crossover systems. 60 Motivated by such results, we have systematically studied such effects on the energy differences between spin-states calculated with the TPSSh functional (using BS2) employing the vibrational information (already employed in Table 3 to compute transition temperatures) to introduce the zero-point energy correction, the D3 scheme 98,99 for dispersion contributions and finally, the DKH method [74][75][76][77] to include relativistic scalar effects.
As can be seen from the Table 4, zero-point energy corrections systematically reduce the energy gap, while relativistic corrections have the opposite effect. Furthermore, electronic energy differences after zero-point energy correction for our study set are reported for the three funcionals that provides with the better results, this is, TPSSh, OPBE and B3LYP* using both basis sets (see Tables S8 and S9). Hence, minimal differences (none for the TPSSh functional) in terms of predicting the correct ground state are observed, but in general, there is a decrease on the energy difference between both spin-states of around 1.7 kcal/mol.
Although DH (electronic energy difference + zero-point correction, see Eq. 3) can eventually be overcome by the entropic term, the TPSSh functional with BS2 gives us values in the range of 2.00 to 8.00 kcal/mol (neglecting the vibrational contribution to DH), an energy window that encloses the usual experimental values for DH, and that should be indicative of having a spin-crossover system.
Dispersion corrections have a much broader behavior (see Table 4), and its larger impact is reflected in two systems (S7 and S15). The inclusion of D3 correction leads to an overestimation of the van der Waals interaction resulting in unrealistic geometries for the more compact low-spin states. In general, for Mn III and Fe II /Fe III systems, inclusion of all such effects together has very little effect in the overall energy difference, leading to similar values for the electronic energy differences. However, the Co II systems seems to be more affected, leading to larger energy gaps when all corrections are accounted for, which eventually will translate in even larger errors when attempting to compute T1/2. In many cases (see Table 4), due to the opposite sign of the correction, the inclusion of the three terms leads to similar results than the non-corrected value.

Discussion
The use of Density Functional methods to compute any property involving more than one electronic state is, by definition, a complicated problem. By combining Equations 2 and 3, we can express the free energy change as, DG = DEelec + DHvib -T·DS [9] in which the energy difference between the two electronic states is explicitly written. Getting show that the combination of the TPSSh functional with a triple-z basis set with polarization functions is, in general, a correct strategy to describe such systems, given that it correctly predicts the ground state for the whole dataset, and the predicted energy gap is in the correct range for spin-crossover to occur for most of cases (see Table 3). The accuracy of the method is quite good in terms of electronic energy differences, but even with a tiny error, this translates, sometimes, in large differences when attempting to compute transition temperatures (see Table 3).
In any case, the vibrational correction of the electronic energies allows for a reasonable modeling of the T1/2, allowing for the study of chemical modifications on the molecule over this physical property. Moreover, this method based on TPSSh functional allows for a rational design of new spin-crossover systems, and given that the functional seems usually to provide with the right range of energies for SCO to occur, the presented strategy can be used in the virtual screening of new targets exhibiting spin crossover. These results are in agreement with previously analysis of similar type, but expand the range of energies for SCO to occur to a higher value. 53 Remarkably, the functional seems to perform properly for complex, the high-spin state shows a strong distortion of the SCN ligands towards the equatorial aromatic rings. In terms of exact exchange mixing, our results seem to indicate that 10% mixing Hartree-Fock is an optimal amount for the hybrid meta-GGA functional to describe SCO systems, although a specific parameterization of this value cannot be ruled out. For B3LYP, it is clear than lowering the amount of exact exchange to 15% (B3LYP*) greatly improves its performance towards SCO systems, although an increase in the quality of the basis set leads to really small energy differences between both spin-states (see Table   S1 in SI).
We want to emphasize at this point that computing transition temperatures is a challenging problem, and although the TPSSh functional can be used to model such quantity and its dependence with the chemical environment around the metal center with great accuracy, such methodology is an approximation that neglects cooperativity effects that can occur in the condensed phase, and it has been shown that inclusion of such effects due to crystal packing can have a significant effect when computing physical properties in SCO systems.
It has been shown, for instance, that different counterions (or solvents) in the unit cell may convert a pure low-spin complex into a spin-crossover molecule, or even a high-spin system. [102][103][104] Much in the same way, if measurements are done in solution, a proper solvation model should be used in order to model such effects in the computed electronic energy differences. However, modeling such effects requires of different computational approaches, which falls out of the scope of this work.

Conclusions
In this work, a systematic use of different Density Functional Theory methods was tested against a benchmark set of first row transition metal complexes exhibiting spin-crossover.
From our calculations, it is clearly shown that the hybrid meta-GGA functional TPSSh combined with triple-z quality basis set with polarization functions is able to correctly model the energy difference between both spin-states for a wide range of metals and oxidation states, ranging from electron configurations d 4 to d 7 (Cr II , Mn III , Mn II , Fe III , Fe II , and Co II ).
Therefore, we propose this density functional method based on TPSSh functional as the current optimal approach to study spin-crossover systems. Moreover, this functional can be used for virtual screening of new first-row transition metal complexes with potential spincrossover behavior. The inclusion of the zero-point energy correction results in a stabilization of the high-spin state of 1-3 kcal/mol. The computed energy range of electronic energy differences between both spin-states, in the range below 10 kcal/mol, should be indicative of SCO behavior. Therefore, systems with energy gaps in that range will be good candidates to exhibit such behavior. Of course, dispersion, zero-point energy and relativistic corrections can have an important contribution in some particular systems, but in general do not play a major role due to their opposite effects, at least in our dataset. Because of its precision in computing electronic energy differences in transition metal complexes, the electronic energies at DFT level obtained with the TPSSh functional can be corrected using the harmonic frequencies to get the thermal dependency of the thermochemical parameters, and therefore, provide with an estimation of the corresponding transition temperature for such systems. In general, the studied methodology provides with transition temperatures that can be qualitatively compared with the experimental values. Therefore, this theoretical tool can be employed to study the effect of chemical modifications in a family of SCO compounds, but its quantitative performance is not accurate enough.

Supporting Information
The energies for all optimized systems with all tested functionals are provided in the SI. Basis set effect on the energy differences and results for the B3LYP functional are also reported in the SI.
Generalitat de Catalunya for an ICREA Academia award. We thankfully acknowledge the computer resources in the Consorci Serveis Universitaris de Catalunya (CSUC).