Computational approach to the study of thermal spin crossover phenomena

The key parameters associated to the thermally induced spin crossover process have been calculated for a series of Fe(II) complexes with mono-, bi-, and tridentate ligands. Combination of density functional theory calculations for the geometries and for normal vibrational modes, and highly correlated wave function methods for the energies, allows us to accurately compute the entropy variation associated to the spin transition and the zero-point corrected energy difference between the low- and high-spin states. From these values, the transition temperature, T1/2, is estimated for different compounds.


I. INTRODUCTION
In spin crossover (SCO) metal complexes the spin state of the metal ion can be changed from a low-spin (LS) ground state to a high-spin (HS) excited state by a variation in temperature, pressure, or by light irradiation. Since the discovery of the thermally induced spin crossover process by Cambi et al., 1, 2 the phenomenon of spin crossover has been extensively studied. 3 However, the interest in SCO materials increased when it was discovered that conversion between the two spin states can be controlled by light irradiation, thus opening the possibility of using such materials as optically switchable devices. [4][5][6] This phenomenon, called Light-Induced Excited Spin State Trapping (LIESST), initially found in Fe(II) complexes [7][8][9][10][11] and later also observed in systems containing Fe(III), [12][13][14][15] and Ni(II) [16][17][18] has been intensively studied in the last years in order to unravel its mechanism both with experimental techniques 11,[19][20][21][22][23][24][25] and by means of theoretical calculations. [26][27][28][29][30][31][32][33] The most numerous and most studied family of SCO systems involves octahedral Fe(II) complexes in the solid state or in solution. The LS-HS transition in Fe(II) complexes is accompanied by an enlargement of the iron-ligand distances due to the occupation of antibonding e orbitals in the HS state. The variation of the metal-ligand bond lengths influences several properties such as the volume of the unit cell, the magnetic and electric properties, and the vibrational spectra, among others. As a matter of fact, the frequencies of the Fe-ligand stretching modes have been used as indicators of the spin transition since the shorter Fe-ligand distances of the LS state result in higher Fe-ligand stretching frequencies a) Electronic mail: c.sousa@ub.edu compared to the HS state. Internal ligand modes are also influenced by the change of the iron-ligand bond length and can be used as probe of spin conversion as well. 34 Thermal spin crossover is an entropy-driven transition from the populated LS state at low temperatures to the HS state, populated at higher temperatures, and is possible when the zero-point energy difference between the HS and LS states ( H HL ) is small, typically of the order of 0-1000 cm −1 . 35,36 An important parameter to characterize the temperature-driven SCO is the transition temperature (T 1/2 ), which corresponds to the temperature at which the LS and HS states are equally populated. A hysteresis loop can be seen in some systems when a different transition temperature T 1/2 is observed by increasing the temperature and by cooling, when the reverse process takes place. This effect is usually related to strong cooperative interactions 37 and normally the difference between the two transition temperatures is of the order of a few K.
Under the condition of thermodynamic equilibrium, T 1/2 can be approximately estimated by the expression: where the enthalpy difference between the HS and LS states, H HL , can be assumed in good approximation to be temperature independent, whereas S HL (T 1/2 ) corresponds to the variation of entropy between the HS and LS states at the transition temperature.
The increase in entropy ( S HL ) associated to the LS-HS transition changes significantly with temperature and it is recognized as the driving force governing the thermal SCO. Experimental values of S HL for octahedral Fe(II) complexes range from 35 to 80 J K −1 mol −1 . 34 The vibrational contribution is the main proportion of the entropy change due to the downshift of the vibrational frequencies under the spin crossover. 38 The electronic contribution to S HL arises from the change in spin multiplicity from the singlet LS state to the quintet HS state and has a constant value of 13.38 J K −1 mol −1 . Bousseksou et al. 39 reported the vibrational contribution to the entropy change for Fe(phen) 2 (NCS) 2 from Raman spectroscopy, taking into account only the 15 vibrational modes associated to an idealized FeN 6 octahedron. Inclusion of these modes accounts for ∼70%-75% of the vibrational part of S HL and the remaining contribution was ascribed to intermolecular or lattice vibrations. Similar results and conclusions were found by Brehm et al. 40 by combining IR and Raman spectroscopy with Density Functional Theory (DFT) calculations. Recent studies combining various vibrational spectroscopic techniques and DFT calculations on the Fe(phen) 2 (NCS) 2 complex 41, 42 allowed the identification of the vibrational modes that contribute most to the entropy variation under spin crossover. It was concluded that only the lowest 29 normal modes significantly contribute to the entropy difference. These 29 modes include the 15 vibrational modes of an idealized FeN 6 octahedral unit and some ligand librational modes strongly mixed with these.
Throughout the last 10 years, studies based on theoretical calculations have demonstrated to be useful tools to get insight into the SCO behavior, as reviewed by Paulsen et al. 43,44 Many DFT or wave function based studies have been performed to investigate various properties of Fe(II) systems such as equilibrium geometries, vibrational spectra, entropy variations, excited states, or the energy difference between LS and HS states, H HL . It is generally accepted that, irrespective to the functional chosen, DFT is the method of choice to obtain accurate structures and vibrational spectra for both LS and HS states at a reasonable computational cost. However, the relative energy of these electronic states turns out to be critically dependent on the exchange functional chosen. Several studies have been performed to establish the optimal percentage of exact exchange in hybrid functionals, however results depend on the system under study. Nowadays, there is not a definitive DFT functional able to compute H HL with sufficient precision for different series of Fe(II) complexes and, in fact, this issue is currently debated in the literature. [45][46][47][48][49][50][51][52][53][54][55][56][57] Multiconfigurational wave function based methods, on the other hand, have proved to be capable of giving accurate values for the LS-HS energy difference but they are computationally much more expensive, particularly in the calculation of optimized geometries and vibrational frequencies. 26,28,29,51,52,[58][59][60][61][62] The aim of this work is to determine through calculations the key parameters of the thermal SCO process, i.e., the zeropoint corrected energy difference between the LS and HS states, H HL , the entropy change associated to the spin transition, S HL , and an estimation of the transition temperature, T 1/2 , for a set of Fe(II) compounds with ligands of different nature. In order to do that, we have combined DFT calculations on the geometries and vibrational frequencies for the LS and HS states, with multiconfigurational wave function calculations that allow us to compute accurate electronic energy differences. The entropy variation and the zero-point energy correction to the LS-HS energy difference can be extracted from the calculated vibrational frequencies, and, in conjunction with the electronic LS-HS energy difference computed by ab initio wave function methods, the value of the transition temperature can be estimated.
Six different Fe(II) isolated complexes have been analyzed ( Figure 1). First, two iron compounds with monodentate ligands, [Fe(mtz) 6 ] 2 + and [Fe(iso) 6 ] 2 + , where mtz refers to 1-methyl-tetrazole and iso to an isoxazole group. Next, three bidentate Fe(II) complexes have been considered, Fe(phen) 2 (NCS) 2 , [Fe(pic) 3 ] 2 + , and [Fe(bpy) 3 ] 2 + , where phen = 1,10-phenanthroline, pic = 2-picolylamine, and bpy = 2,2 -bipyridine. Finally, a tridentate [Fe(terpy) 2 ] 2 + complex has also been studied, terpy = 2,2 :6 ,2 -terpyridine. In all cases, Fe(II) is surrounded by a N 6 nearly octahedral first coordination sphere. These complexes are found in various SCO materials. The [Fe(mtz) 6 ] 2 + unit is present in the Fe(mtz) 6 (BF 4 ) 2 molecular crystal occupying two nonequivalent crystal positions, site A susceptible to thermal SCO at 78 K and site B, which remains HS down to 10 K. 63 The [Fe(iso) 6 ] 2 + complex is found in two crystals, Fe(iso) 6 (BF 4 ) 2 and Fe(iso) 6 (ClO 4 ) 2 , occupying two non-equivalent sites with a 1:2 ratio. Site A coincides with an inversion center, while site B has a 3-fold symmetry center. The first compound undergoes two reversible SCO transitions at 91 and 192 K, assigned to each one of the two sites, whereas in the second compound both sites undergo a simultaneous spin transition at 213 K. 64 The crystal structure of Fe(phen) 2 (NCS) 2 has been studied showing the existence of a LS-HS transition at around 176 K which is accompanied by an increase of 0.20 and 0.10 Å of the Fe-N(phen) and Fe-N(CS) distances, respectively. 65 The crystal structure of [Fe(pic) 3 ]Cl 2 · EtOH has been characterized 66 and the transition temperature has been measured between 114.0 and 120.7 K. 67 [Fe(bpy) 3 ] 2 + is a LS complex where thermal spin conversion is only possible at very high temperatures. The structure of the LS state of this complex has been resolved both in aqueous solution 19 and in the crystal compound [Fe(bpy) 3 ](PF 6 ) 2 . 68 The tridentate [Fe(terpy) 2 ] 2 + complex, like the [Fe(bpy) 3 ] 2 + , is a LS compound and spin transition is only possible by light irradiation. The structure of this complex in the LS state has been measured in the [Fe(terpy) 2 ](ClO 4 ) · H 2 O crystal 69 showing a tetragonal distortion with two different Fe-N distances, 1.892 Å for the central N of the terpy ligand and 1.988 Å for the distal N atoms.
In the present work, it will be shown that the combination of DFT calculations for geometry and vibrational frequencies and multiconfigurational methods for energy differences between different electronic states, allows to achieve an overall proper description of the thermodynamic and vibrational properties of the thermal SCO process in the systems studied.

II. METHODOLOGY
DFT electronic structure calculations were performed for the six Fe(II) complexes studied. The optimized geometries and the normal vibrational modes were obtained for the singlet LS and quintet HS states. Two different hybrid functionals were employed, B3LYP 70, 71 and PBE0, 72 and two basis sets were applied: split valence plus polarization (SVP) and triple zeta valence plus polarization (TZVP). 73,74 All DFT calculations were performed with the TurboMole 6.3 package. 75 The wave function based calculations were performed applying the CASSCF/CASPT2 76 method, i.e., second-order perturbation theory based on a complete active space selfconsistent field reference wave function, implemented in the MOLCAS 7.4 software. 77,78 Scalar relativistic effects were included using a Douglas-Kroll-Hess Hamiltonian 79, 80 and atomic natural orbital basis sets, specifically designed to include relativistic effects, have been used. 81,82 The contracted Gaussian basis sets applied are: (7s, 6p, 5d, 4f, 3g, 2h) for Fe, (4s, 3p, 1d) for the N atoms bonded to the Fe, (3s, 2p) for the remainder N atoms, for O and C, (4s, 3p) for S, and (2s) for H. The active space used to construct the CASSCF wave functions for the LS and HS states contains ten electrons distributed among 12 orbitals, the five 3d Fe orbitals, two e g -like σ -bonding ligand orbitals, and a second set of diffuse Fe-3d orbitals to account for the large electron correlation effects in the 3d-shell. This active space has been used in previous studies 26,29,51 and aims at a proper description of the different Fe d-states and the ligand to metal charge transfer effects. CASPT2 accounts for the remaining electron correlation by correlating all the electrons except the deep core electrons (1s 2 for N and C and 1s 2 , 2s 2 , 2p 6 for Fe and S). The standard zeroth-order Hamiltonian has been used in the second order perturbative CASPT2 method. 83 This Hamiltonian uses the so-called ionization potential-electron affinity (IPEA) shift of 0.25 a.u. in the definition of the diagonal Fock matrix elements of the active orbitals. Recently, it has been proposed that larger values for the IPEA shift lead to better agreement with experiments 84 and with benchmark coupled cluster calculations. 85 The two approaches are compared in Sec. III.
One-dimensional CASPT2 potential energy surfaces (PES) have been computed for the LS and HS states of the six complexes. Since the largest geometrical variation induced by the spin conversion involves the Fe-N distances, usually the PES are plotted as function of this single variable. Therefore, the reaction coordinate corresponds to the symmetric breathing mode for the Fe(II) systems coordinated to monodentate ligands. However, for complexes with multidentate ligands, the sole variation of the Fe-N distance can result in strains in the structure of the system. Hence, a linear interpolation between the HS and LS B3LYP optimized structures, obtained with the TZVP basis set, was taken as a reaction coordinate for the study of the multidentate complexes. The representation of the CASPT2 potential energy curves along the reaction coordinate permits to obtain an estimate of the CASPT2 Fe-N equilibrium bond distance. However, as the corresponding geometry is not the absolute minimum, calculations of vibrational frequencies are not possible at this level of calculation.

A. Equilibrium geometries and vibrational frequencies
The geometry of the six complexes studied in this work has been fully optimized both for the LS and HS states applying two different functionals, B3LYP and PBE0, and two basis sets, SVP and TZVP, as explained in Sec. II. The optimized values of the Fe-N distances obtained with the different functionals and basis sets are reported in Table S1 of the supplementary material. 86 For a given functional, increasing the basis set has only a small effect on the Fe-N distances, consistently less than 0.01 Å. To analyze the performance of the different computational methods employed in this study, Figure 2 shows the Fe-N average distances for LS and HS states for all six compounds. Reported B3LYP and PBE0 optimized distances correspond to those obtained using the TZVP basis sets. CASPT2 values arise from the minimization of the energy of the HS and LS states along a single coordinate that connects the optimized B3LYP (TZVP basis set) geometries of the LS and HS states. Experimental values are available for all compounds, except for the HS geometries of the two LS complexes, [Fe(bpy) 3 ] 2 + and [Fe(terpy) 2 ] 2 + .
From Figure 2 it can be seen that for all systems B3LYP gives the largest Fe-N distances, CASPT2 tends to somewhat underestimate the Fe-N bond lengths, whereas PBE0 gives the closest values to the experimental measurements. Overall, all methods reproduce the experimental data with sufficient accuracy, and particularly, the variation of the Fe-N distance under spin transition, around 0.2 Å for each one of the systems, is properly reproduced by all three methods.
The optimized values of the N-Fe-N bite angles for the bidentate and tridentate complexes reported in Table I  effects can have an important influence on the rotation angle of the ligands and that DFT geometry optimization in the isolated unit tends to unrealistic angles, around 45 • . 62 Taking into account the environment effects with an approximate model, the experimental angles of the [Fe(mtz) 6 ] 2 + complex are properly reproduced. The results reported in this work correspond to a rotation angle of the methyl-tetrazole ligands of 20 • , which is in the range of the experimental values. For the monodentate [Fe(iso) 6 ] 2 + complex we have studied site A, which lies on an inversion center and undergoes a spin transition at 91 K. For this system, DFT calculations on the HS state satisfactorily reproduce the rotation angles measured experimentally for this site, 64 around 50 • . No experimental data are available for the LS state, DFT calculations give a rotation angle of the isoxazole ligands around 22 • . CASPT2 calculations on this system have been performed at this angle.
The results in Figure 2 and Table I confirm that, in general, DFT methods properly describe the geometric structure of Fe(II) spin crossover complexes containing ligands of different nature.
Vibrational frequencies at the LS and HS equilibrium geometries have been computed analytically in the harmonic approximation for the six complexes under study. The values of the calculated vibrational frequencies do not significantly change with the basis set chosen or the functional applied. A detailed report of the analysis of the whole vibrational spectra of these complexes is beyond the scope of this work; comprehensive studies of the vibrational spectrum of Fe(phen) 2 (NCS) 2 , both by DFT calculations and by spectroscopic techniques, have been provided elsewhere. 40,41 For this particular system, the vibrational frequencies obtained here are in good agreement with those reported. In general terms, and for all the systems studied, the lowest frequencies (normally less than 100 cm −1 ) are related to motions of the ligands. For instance, in the [Fe(mtz) 6 ] 2 + complex, the lowest vibrational modes imply rotation of the methyl-tetrazole ligands, while for the Fe(phen) 2 (NCS) 2 system the modes with frequencies below 100 cm −1 are mainly out-of-plane motions of the phenanthroline groups. For all systems, vibrational modes with frequencies larger than, approximately, 500-600 cm −1 can be associated to ligand modes and almost do not change between LS and HS states. However, the modes with frequencies in the region of 100-600 cm −1 involve ironligand vibrations and significantly change when varying the spin state, with the frequencies of the HS state being lower than those of the LS since the Fe-N bond is weakened in the HS state. Consequently, these vibrational modes contribute the most to the entropy variation and to the zero-point energy difference between the HS and LS states.

B. Computed enthalpies of the spin crossover process
The energy difference between the HS and LS electronic states, H HL , is one of the most important parameters in the thermal SCO process. However, as commented before, it is usually not properly accounted for with standard DFT methods. Table II reports the computed values of the relative energy between the HS and LS states for the complexes studied.   These energy differences are computed from the minimum of the potential energy surfaces of the two states and therefore, vibrational contributions are neglected. PBE0 and B3LYP values in Table II have been calculated using the TZVP basis set. In contrast to the minor influence in the geometrical parameters, the size of the basis set has a significant effect in the HS-LS energy difference, with variations between SVP and TZVP computed values ranging from 300 to 800 cm −1 (see Table S2 of the supplementary material 86 ). Negative values in Table II mean that the relative ordering of the states is inverted, i.e., that the calculation gives the HS more stable than the LS state. As can be seen in the table, both PBE0 and B3LYP functionals overstabilize the HS with respect to the LS state. PBE0 gives the wrong energetic order in all cases. B3LYP is able to reproduce the right sign of the energy difference for only two of the systems, [Fe(bpy) 3 ] 2 + and [Fe(terpy) 2 ] 2 + , which in fact are LS complexes. Alternatively, energy differences computed at the CASPT2 level show the right order for all the complexes and a reasonable quantitative agreement with experimental estimations, when available. One contribution that is commonly disregarded in the computation of the HS-LS energy difference is the vibrational zero-point energy correction on the spin transition. Since normal vibrational frequencies are larger in the LS than in the HS state, the zero-point energy contribution to the enthalpy variation has a negative value, implying a relative stabilization of the HS state. From the calculation of the frequencies at the equilibrium geometry of the LS and HS states, zero-point energy corrections to the HS-LS energy differences have been computed for all six systems. The zero-point energy contribution amounts to an important part of the energy difference between HS and LS states, H HL . Computed values using B3LYP and a TZVP basis set are: −982 cm −1 for  Table S3 of the supplementary material 86 ).
Owing to the fact that higher frequency vibrational modes correspond to ligand motions, and those are similar in the HS and LS states, the largest contribution to the zeropoint energy correction is concentrated in the lowest vibrational modes. For instance, for the Fe(phen) 2 (NCS) 2 system, the 72% of the zero-point energy correction is recovered by taking into account the 29 lowest normal modes (with frequencies below 400 cm −1 ) and 85% of the zero-point energy correction is accounted for by including the lowest 49 normal modes, with frequencies below 600 cm −1 . This behavior is general for all the complexes.
To obtain a proper estimate of H HL , the zero-point energy corrections have to be added to the energy difference between the HS and LS states reported in Table II. Including this correction leads to smaller spin transition enthalpies, which means that PBE0 and B3LYP values become even more negative. Although the zero-point energy has been computed based on DFT calculations, this correction has been applied to the CASPT2 HS-LS energy difference. This is justified by the observation that DFT is able to properly describe the vibrational spectra of these compounds and moreover, as mentioned before, the values of the vibrational frequencies do not strongly depend on the particular functional or basis set. Therefore, one can rely on the estimation of the zero-point energy correction as obtained by DFT. Particularly, as CASPT2 energy differences have been obtained from an interpolation between optimized B3LYP structures using a TZVP basis set, the corresponding zero-point energy correction has been taken into account. The resulting values of H HL are displayed in Table III. It is worth noting that, after inclusion of the zeropoint energy correction, H HL has the right sign for all complexes and compares satisfactorily with the experimental estimates (Table II). Small values of H HL , less than 1000 cm −1 , are found for the four SCO complexes, [Fe(mtz) 6 ] 2 + , [Fe(iso) 6 ] 2 + , Fe(phen) 2 (NCS) 2 , and [Fe(pic) 3 ] 2 + , while large values are found for the LS complexes, [Fe(bpy) 3 ] 2 + and TABLE III. Computed values of S HL (T 1/2 ) (in J K −1 mol −1 ), H HL (in cm −1 ), T 1/2 (in K), and experimental estimates of H HL and T 1/2 . The choice of using the B3LYP geometries for the construction of the CASPT2 PES may seem somewhat counterintuitive given the results shown in Fig. 2, which show that on average the PBE0 geometries compare slightly better to experiment. One should however keep in mind that the DFT results are also used to estimate the zero-point energy correction. Given the fact that the interpolation between the HS and LS optimized geometries lead to practically the same optimal Fe-N distance at the CASPT2 level for both functionals, that the frequencies are nearly identical and that the energetics of B3LYP are slightly better than the PBE0 ones, we opted to take the B3LYP results as reference. We mention that the results are quantitatively very similar when the PBE0 results are used throughout.
Recently, it has been suggested that appropriate values of H HL can be obtained by applying an IPEA parameter larger than the standard value (0.25 a.u.) in the definition of the zeroth-order Hamiltonian of the CASPT2 method. 84 Values between 0.5 and 0.7 a.u. have been suggested in order to compute HS-LS energy differences in Fe(II)N 6 systems. We have calculated the HS-LS energy difference for some of the complexes studied here using four values of the IPEA shift: the standard value 0.25 and 0.50, 0.60, and 0.70 a.u. Increasing the IPEA value favors the LS state and, hence, larger HS-LS energy differences are found. For the systems studied an increase of ∼1400-1500 cm −1 is found when the IPEA shift is varied from 0.25 to 0.50 a.u., similarly as reported for other complexes. 84 However, for the systems studied here, the values obtained with larger IPEA parameters show less good agreement with experimental data than those obtained with the standard definition of the zeroth-order Hamiltonian. Therefore, this strategy appears to be not generally applicable for all Fe(II)N 6 complexes.

C. Spin crossover entropies and transition temperature
Total entropy variations associated to the LS-HS transition for the complexes under study have been computed by DFT calculations applying both the PBE0 and B3LYP functionals. The total entropy change can be divided into three contributions: electronic, vibrational, and rotational. As commented in the Introduction, the electronic term, related to the change in spin multiplicity from a singlet to a quintet, is constant for all systems and has a value of 13.38 J K −1 mol −1 . The entropy contribution of the rotational degrees of freedom is expected to be small since no large structural changes occur under spin transition for the systems studied here and, consequently, the moments of inertia of the LS and HS states are similar. The largest contribution to the entropy change is due to vibrations because of the significant variation in the normal modes under spin transition. This term can be easily derived at different temperatures from the vibrational frequencies calculated within the harmonic approximation. The total entropy variation S HL is calculated in the range from 10 to 1000 K at intervals of 1 K, and subsequently the T S HL product is compared to H HL to determine the transition temperature following Eq. (1). The values of H HL , S HL (T 1/2 ), and T 1/2 are reported in Table III for the six complexes.
Like the zero-point energy correction, the largest contribution to the vibrational entropy variation is concentrated in the lowest vibrational modes, as has been pointed out in previous studies. [40][41][42] Moreover, the values of the vibrational contribution to the entropy do not appreciably differ when varying the functional or the basis set used. In Table III, B3LYP values of S HL (T 1/2 ) calculated using the TZVP basis set are collected. The total entropy variation ranges from 29 J K −1 mol −1 for the [Fe(mtz) 6 ] 2 + complex to 93 J K −1 mol −1 for the tridentate [Fe(terpy) 2 ] 2 + system. Computed values for the Fe(phen) 2 (NCS) 2 system and the [Fe(pic) 3 ] 2 + complex, 55 and 51 J K −1 mol −1 , respectively, compare well with the experimental values reported, 49 38 and 59.5 J K −1 mol −1 . 88 The Fe(phen) 2 (NCS) 2 system has been extensively investigated by electronic structure calculations. The value of S HL reported in this work is similar to the values obtained using different functionals and basis sets. 41,89 The calculated transition temperatures, T 1/2 , correlate well with the experimental estimates. For the four SCO materials, the lowest value is found for the [Fe(mtz) 6 ] 2 + complex and the largest for the Fe(phen) 2 (NCS) 2 system, as experimentally found. Calculated values of T 1/2 are less than 50 K below the experimental data. Similar differences have been found in studies of other SCO systems. 55 This discrepancy can be explained by the neglect of crystal effects. In the present work, the systems are treated as isolated units, therefore the influence of counterions, intermolecular interactions, and crystal packing effects are not taken into account. These factors can induce (huge) changes in the transition temperature, as shown by Lemercier et al. 90 and Carbonara et al., 91 among others. For the two LS complexes, [Fe(bpy) 3 ] 2 + and [Fe(terpy) 2 ] 2 + , the calculated transition temperature is found to be very high, in accordance with the fact that in these two compounds thermal SCO is only possible at high temperatures.
Overall, results in Table III show a good agreement of the thermodynamic parameters extracted from ab initio calculations with the experimental observations for different types of Fe(II) complexes. The trend in the transition temperature is well reproduced along the different systems, and this approach, although not quantitative enough to compute accurate values of T 1/2 , allows to reasonably estimate the temperature at which the spin transition takes place.
In order to determine the vibrational contribution to the total entropy, a low-frequency approximation has been proposed 39 based on the expression: where p is the number of normal modes to be considered and ν represents an average wavenumber of the frequencies of the lowest p oscillators for both the LS and HS states. Applying this formula for the SCO systems studied here result in a rather erratic behavior. As has been already pointed out, 40 theν LS /ν H S ratio importantly changes with the number of vibrational modes taken into account. The vibrational contribution to the entropy, S vib , has been calculated considering the modes with frequencies below 600 cm −1 (the ones that change the most when varying the spin state) and applying Eq. (2) for the four SCO complexes. For two of the complexes, Fe(phen) 2 (NCS) 2 and [Fe(pic) 3 ] 2 + , Eq. (2) leads to an astonishing agreement with experimental values. For the first one, inclusion of the lowest 49 normal modes gives a vibrational entropy contribution of 36.7 J K −1 mol −1 , and after adding the electronic contribution, a total entropy of 50.1 J K −1 mol −1 is obtained (experimental value 49 J K −1 mol −1 ). For the [Fe(pic) 3 ] 2 + complex inclusion of the 36 modes with frequencies below 600 cm −1 leads to a total entropy of 59.0 J K −1 mol −1 , in excellent agreement with the experimental data, 59.5 J K −1 mol −1 . However, this accord appears to be accidental since for the rest of the systems the low-frequency approximation gives values that largely differ from the computed S HL . Therefore, results obtained by this approximation can be taken only as a rough estimate of the vibrational contribution to the spin transition entropy.

IV. SUMMARY AND CONCLUSIONS
The thermodynamic properties of a series of spin crossover Fe(II) complexes with ligands of different nature have been computed by a combination of theoretical methods. Density functional theory methods have demonstrated to be feasible and accurate tools to optimize geometries and calculate vibrational spectra for the LS and HS states involved in the spin conversion. Optimized structures and harmonic vibrational frequencies only slightly vary with the exchangecorrelation functional chosen or the basis set applied. Accordingly, the vibrational contribution to the SCO entropy variation and the zero-point correction energy can be easily deduced from DFT calculations. The lowest normal modes, below approximately 600 cm −1 , contribute the greatest part to the zero-point energy correction and to the vibrational entropy contribution.
A crucial parameter in the spin conversion process is the energy difference between the two spin states. In order to obtain reliable values of the LS-HS electronic energy difference, multiconfigurational wave function methods, like the CASPT2 method employed in this work, are required. However, to quantitatively estimate the energy difference between low-and high-spin states, inclusion of the zero-point correction energy is essential. This zero-point vibrational correction amounts for a significant part of the total LS-HS energy difference. Therefore, the zero-point correction obtained by DFT calculations has been added to the CASPT2 electronic energy difference.
By combining the optimized structure and vibrational frequencies obtained by DFT and the energetics obtained by CASPT2 calculations, the transition temperature of various SCO complexes has been estimated. Calculated transition temperatures for a series of compounds correlate properly with the experimental data.
The mixed approach presented here allows to obtain an overall correct description of the key parameters that characterize the thermal spin crossover process in the systems studied. Although we have focused on the thermal spin crossover in Fe(II) complexes, the approach is equally valid for mononuclear systems with other transition metals provided that the single determinant description is reasonable for both the low-spin and the high-spin state to obtain accurate estimates with DFT of the geometry and the entropy contribution. To apply the strategy to polynuclear complexes one has to face the problem of the size of the active space, which may become prohibitively large. One possible solution may be to rely on a restricted active space (RAS) reference wave function, but this requires additional testing to see whether the calculated RASPT2 energies are as reliable as those obtained within the here-described complete active space approach.

ACKNOWLEDGMENTS
This work is part of the program of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). Financial support has been provided by the Spanish Administration (Projects CTQ2011-23140, FIS2008-02238, and CTQ2012-30751), the Generalitat de Catalunya (Projects 2009SGR462, 2009SGR-1041, and Xarxa d'R+D+I en Química Teòrica i Computacional, XRQTC), and the European Union (COST Action CODECS CM1002). R.W.A.H. acknowledges the Zernike Institute for Advanced Materials of the University of Groningen for financial support ("Dieptestrategie" program).