Thermal Spin Crossover in Fe(II) and Fe(III). Accurate Spin State Energetics at the Solid State

The thermal Spin Crossover (SCO) phenomenon refers to an entropy-driven spin transition in some materials based on d 6 -d 9 transition metal complexes. While its molecular origin is well known, intricate SCO behaviours are increasingly common, in which the spin transition occurs concomitantly to e.g. phase transformations, solvent absorption/desorption, or order-disorder processes. The computational modelling of such cases is challenging, as it requires accurate spin state energies in the solid state. Density Functional Theory (DFT) is the best framework, but most DFT functionals are unable to balance the spin state energies. While few hybrid functionals perform better, they are still too expensive solid-state minima searches in moderate-size systems. The best alternative is to dress cheap local (LDA) or semi-local (GGA) DFT functionals with a Hubbard-type correction (DFT+U). However, the parametrization of U is not straightforward due to the lack of reference values, and because ab initio parametrization methods perform poorly. Moreover, SCO complexes undergo notable structural changes upon transition, so intra- and inter-molecular interactions might play an important role in stabilizing either spin state. As a consequence, the U parameter depends strongly on the dispersion correction scheme that is used. In this paper, we parametrize U for nine reported SCO compounds (five based on Fe II , 1-5 and four based on Fe III when and We Abstract The thermal Spin Crossover (SCO) phenomenon refers to an entropy-driven spin transition in some materials based on d 6 -d 9 transition metal complexes. While its molecular origin is well known, intricate SCO behaviours are increasingly common, in which the spin transition occurs concomitantly to e.g. phase transformations, solvent absorption/desorption, or order-disorder processes. The computational modelling of such cases is challenging, as it requires accurate spin state energies in the solid state. Density Functional Theory (DFT) is the best framework, but most DFT functionals are unable to balance the spin state energies. While few hybrid functionals perform better, they are still too expensive solid-state minima searches in moderate-size systems. The best alternative is to dress cheap local (LDA) or semi-local (GGA) DFT functionals with a Hubbard-type correction (DFT+ U ). However, the parametrization of U is not straightforward due to the lack of reference values, and because ab initio parametrization methods perform poorly. Moreover, SCO complexes undergo notable structural changes upon transition, so intra- and inter-molecular interactions might play an important role in stabilizing either spin state. As a consequence, the U parameter depends strongly on the dispersion correction scheme that is used. In this paper, we parametrize U for nine reported SCO compounds (five based on Fe II , 1 - 5 and four based on Fe III , 6 - 9 ) when using the D3 and D3-BJ dispersion corrections. We analyze the impact of the dispersion correction treatments on the SCO energetics, structure, and the unit cell dimensions. The average U values are different for each type of metal ion (Fe II vs. Fe III ), and dispersion correction scheme (D3 vs. D3-BJ) but they all show excellent transferability, with mean absolute errors (MAE) below chemical accuracy ( i.e. MAE < 4 kJ/mol). This enables a better description of SCO processes and, more generally, of spin state energetics, in materials containing Fe II and Fe III ions.


Introduction
Some organometallic complexes based on d6-d9 transition metal ions undergo a spin-state switch in the presence of an external perturbation. [1][2][3][4][5][6][7][8] These molecules, called Spin Crossover (SCO), are prototypical class of molecular magnetic switches, 9 and are regarded as interesting for a range of technological applications. [10][11][12][13] A key property of SCO materials is the temperature at which the spin transition occurs, called transition temperature ( 1/2 ), as it defines their potential applicability in real-life devices operating at mild temperatures. One of the main energetic contributions controlling 1/2 is the electronic enthalpy difference between the two spin states of the material (∆ ). This term control whether a complex is able to undergo thermal SCO, and at which 1/2 , and is also connected to the thermal relaxation temperature ( ) when the SCO is triggered by light in the so-called LIESST (Light-Induced Excited Spin-State Trapping) effect. [14][15][16][17] In gas-phase conditions, or in the absence of any relevant external influence, ∆ accounts for the ligand-field splitting (∆ ), and intramolecular (e.g. inter-ligand) interactions. In crystals, ∆ is further affected by the crystal packing effects, which might greatly modify its value. 18,19 Consequently, it is rather common to find polymorphs or solvatomorphs of the same SCO complex displaying completely different magnetic transitions. [20][21][22][23] Moreover, other solid state phenomena like phase transitions, solvent absorption/desorption or order disorder processes might occur concomitantly to the spin transition, occasionally leading to multiple SCO pathways, 24 reverse SCO, 25 or a spin state blocking along the entire range of temperatures. 26,27 As a result, the evaluation of ∆ needs to account for the energetic impact of such processes on the SCO. Indeed, most of the important applications of SCO require the description of its spin state energetics in either the solid-state, in surfaces, or in an interface and, hence, it is crucial to evaluate ∆ correctly while accounting for the environment.
There are reliable ways to compute ∆ from electronic structure calculations in gas-, solvent-or crystalline-phase. [28][29][30][31][32][33] The best one to treat the latter phase is DFT+U, which implies the dressing of simple DFT functionals with the Hubbard-like term (U). [34][35][36][37] This scheme has been employed in its GGA (PBE+U) or LDA (LDA+U) flavors to study SCO phenomena since the early 2000's, being pioneered by Letard and coworkers, 38 and Angyan and coworkers. 39 The main advantage of DFT+U is its computational cost, enabling minima searches in any complex environment (e.g. solid state). Initially, the DFT+U scheme had two major disadvantages. The treatment of dispersion interactions was restricted to the underperforming LDA or GGA functionals, and the empirical choice of the U parameter undermined its predictive power. The first issue could be tackled 40 by incorporating cheap dispersioncorrection techniques (such as D2), 41 whereas the latter was solved by our previous benchmark of the PBE+U+D2 scheme. 19 Therein, the adequate value of U able to describe the SCO energetics was determined for a group of seven Fe II -based compounds. The average value of U=2.65 eV showed good transferability, with a mean absolute error (MAE) of 4.4 kJ/mol when applied to all compounds. Such good accuracy prompted us to tackle a number of puzzling phenomena in the context of thermal SCO, including cooperativity. 42 The method has been systematically correct in capturing (i) the relative ordering of the LS and HS states (i.e. correct sign of ∆ ) and (ii) the trends measured experimentally (i.e. correct change of ∆ ). [43][44][45] Moreover, it has provided good quantitative accuracy in comparison to other DFT functionals. However, better dispersion correction schemes have been developed and implemented during this time, and a new benchmark has become necessary. Moreover, the benchmarked U parameter for Fe II is likely not adequate to treat SCO metal complexes based on other metal ions, such as Fe III , and hence an extension was necessary (Fe II -and Fe III -based compounds are the largest families of inorganic SCO systems, 6, 7 together with Co II 8 ).
With this in mind, the motivation of the present paper is, first, to (i) benchmark U under an improved dispersion correction scheme (D3 and D3-BJ) and, second, to (ii) extend the application of the PBE+U scheme to the family of Fe III SCO compounds. To do so, we have benchmarked the U parameter for five Fe II compounds (1)(2)(3)(4)(5), and four Fe III SCO compounds (6)(7)(8)(9). This benchmark is then complemented with gas-phase computations employing a range of DFT functionals. The results mutually reinforce the use of (i) PBE+U with the benchmarked U values when in the solid-state, and of B3LYP* (for Fe II ) and OLYP (for Fe III ) functionals when gas-phase calculations at a higher computational level are feasible. All these schemes achieve chemical accuracy (MAE≤ 4 kJ/mol). PBE+U, due to is low computational cost, is particularly recommended for the application of modern screening methods and accelerated discovery techniques.

Energetics of the SCO transition
In gradual SCO transitions, the temperature at which the SCO occurs ( 1/2 ) corresponds to the ratio between the enthalpy and entropy differences (∆ and ∆ ) involved in the process. Both terms have two major contributions, electronic and vibrational (∆ , ∆ , ∆ and ∆ ). The simplest one, ∆ , accounts for the change in electronic degeneracy between the HS and LS states. The vibrational terms ∆ and ∆ account for the change in the vibrational levels of the molecule, mostly due to the occupation of the antibonding eg orbitals and the concomitant expansion of the coordination sphere upon SCO. In a crystal, these terms also account for the change of lattice phonons. ∆ can be safely modelled using the Harmonic-Oscillator (HO) model, whereas ∆ is better treated combining the HO and Free-Rotor (FR) models. 51 Both models use the vibrational normal modes ( ), whose computational evaluation remains a challenge. The computational cost of accurately evaluating in molecular crystals is large, so calculations usually restrict to isolated molecules. As a result, the evaluation of ∆ and ∆ does not incorporate the effect of lattice vibrations, intermolecular interactions and anharmonicity, which might occasionally have an impact on and, thus, on the vibrational contributions. 52, 53

List of compounds.
Our main criteria to choose the studied compounds has been: (i) available crystal structures, without significant disorder, (ii) complete, single-step spin transition and, when possible, (iii) calorimetry measurements offering an experimental estimate of ∆ and ∆ , which facilitate the assignation of a reference ∆ value (i.e. ) from which we can benchmark our methodology. Moreover, we have favored gradual and, specially, non-hysteretic SCO. The reason is that systems in which cooperativity has a complex energetic fingerprint, 1/2 is no longer defined as the ratio between the enthalpy and entropy differences. 42 The group of Fe II compounds includes (i) the ubiquitous Fe(phen)2(NCS)2 (1) of Sorai and Seki, 54 in the crystal structure reported by Gallois,46 (ii) the X=S (2) and X=Se (3) 26 This set of compounds is almost the same we studied previously in ref. 19 . An exception are compounds [Fe(HB(pz)3]2 and Fe[H2B(pz)2]2(bipy) of our original benchmark, which have not been included in this study. For the former compound, the reason is that its experimental ∆ and ∆ values refer to the first jump (from LS to 1:1 HS:LS) of its two-step transition, whereas our calculations assumed the full LS-to-HS transformation. This was already mentioned in our original paper, but not properly tackled, as it should have never been used as a good reference value. For the latter, the available estimations of ∆ and ∆ were obtained from a fitting of the χT curve with the Slichter-Drickammer model. 56 The fitted value of ∆ (81.9 J/K·mol) is 25% larger than the largest value in our library of Fe II compounds so we did not considered it reliable (see Table 1) and, hence, we disregarded this compound. The list of Fe III compounds has been much more difficult to establish. A comparatively-small amount of works in the literature report DSC measurements, and only few are free from disorder, symmetry breaking, solvent evaporation or hysteresis. Our final list includes three compounds synthesized by Renz and coworkers: [Fe III (L 1 )(NCS)] (6), [Fe III (L Cl )(NCS)] (7) and [Fe III (L Br )(NCSe)] (8), 48,49 in which the thermodynamic quantities were extracted using an Ising-like model developed by Boča and coworkers, which provided very robust values that are in line with the existing literature. 57 The last system is [Fe(qsal-I)2]Otf·MeOH (9) reported by Harding, 50 in Table 2. Calculated thermodynamic data for compounds 1-9, and benchmarked U values under the D3 and D3-BJ schemes. Temperature is given in K, enthalpy in kJ·mol −1 , entropy in J·K −1 ·mol −1 , and the U values in eV. The thermodynamic data is given per molecule. which the thermodynamic parameters were extracted from DSC measurements (see Table 1). The structure of 1-9 is represented in Fig. S1.

D2 vs. D3 vs. D3-BJ dispersion corrections.
The main difference between the D2 and D3 corrections is that, in the former, the C6 coefficients are based exclusively on the atom type, whereas D3 takes into account also their chemical environment through a coordination number. Accordingly, the C6 and C8 coefficients are also different. Finally, each scheme uses a different damping function aimed at correcting the energetic contribution at very close ranges (Fermi function in D2, Chai and Head-Gordon 58 in D3, and Becke and Johnson 59-61 in D3-BJ). Further information about dispersion corrections in general, and the D2 and D3 in particular, can be found in the literature. 62, 63

Estimation of reference values and U benchmark.
Following the working strategy we used in our previous benchmark of the PBE+U+D2 scheme, 19 we have estimated and for compounds 1-9 (see Table 2). The strategy consists in extracting from the experimental (from DSC or fitted) value of and a calculated value of . is computed from slightly-modified vibrational normal modes of the isolated molecule: the eight lowerfrequency of the LS and HS minima are adjusted to reproduce the value that results from subtracting from the experimental . Given that is mostly affected by the high-frequency , this adjustment does not have a significant impact on and, hence, on . Therefore, the resulting value is our best estimate of the adiabatic energy difference between the HS and LS states of compounds 1-9 in the solid state (see Table 2).
Once the reference values have been estimated, the Hubbard-like U parameter under the PBE+U+D3 and PBE+U+D3-BJ schemes is parametrized to reproduce these values. We start by performing variable-cell geometry optimizations to obtain the minima of the HT and LT phases of each compound. Then, at these fixed structures, we perform single-point energy evaluations at different values of U ranging from 1.6 and 3.0 eV. Within this range, the evolution of is mostly linear, so the exact U value reproducing is easily interpolated (see Figure  S2). Notice that, in this benchmark, our philosophy is to employ the same U value for both the LS and HS states of the material. Also, our method is empirical in the sense that the U values are benchmarked to reproduce estimated values. This strategy is typically discouraged in favor of ab-initio benchmarking methods that would result in different U values for each of the LS and HS minima. 37 66 All computations have been performed at the Γ-point of the Brillouin zone. The minimum energy structure has been obtained by performing successive variable-cell geometry relaxations, in which the lattice parameters as well as the atomic positions are optimized simultaneously until the atomic forces are smaller than 1.0·10 −5 atomic units. In these calculations, the number of plane waves has been kept constant at a kinetic energy cutoff of 70 Ry for the wavefunction (ecutwfc) and of 560 Ry for the charge density (ecutrho) throughout the variable-cell relaxations. A final SCF calculation has been done at the final optimized structure with kinetic energy cutoffs of 35 and 280 Ry for ecutwfc and ecutrho, respectively. The spin state of the iron atoms is defined in the initial guess, and maintained along the optimization. We must note here that the LS crystal structure for compounds 2 and 3 has not been reported experimentally. In those cases, and given that the spin transition neither implies a change in the symmetry of the crystal, nor in the orientation of the molecules, the LS minima have been obtained starting from the HS crystals.
PBE+U calculation in the gas phase. To calculate the energy of the isolated molecules, the coordinates of a SCO molecule, whose typical size is ca. 35 Bohr 3 , have been excised from the optimized unit cell and introduced on a cubic cell of 60 Bohr 3 . This fact, together with the application of the Makov-Payne correction within QE, 67 ensures that the molecules are isolated from their virtual counterparts. In these calculations, the geometry has not been optimized in order to preserve the structure found in the variable-cell geometry optimizations. If non-identical SCO units have been found on the optimized crystal, this process has been done for each of those, and the results correspond to an average. In Section 3.3, we have tested several exchange correlation functionals at describing the correct spin-state energy differences, including OLYP, OPBE, B3LYP, B3LYP*, TPSSh, M06L and SCAN. All single point computations have been performed with the quantum chemistry package Q-Chem 5.0. 68 The def2-TZVP 69, 70 basis set was used on all atoms, and a convergence criterion on the energy of 10 -8 was employed for all calculations. Dispersion correction effects were included via Grimme's D3 and D3-BJ schemes. 71

Estimation of and U for 1-9
The protocol described in Section 2.3 to extract the reference values has been applied to compounds 1-9 (see Table 1 . As we will discuss later, the reason is that the D2 correction has a much larger contribution to than D3 and D3-BJ, especially when it comes to the description of the intramolecular interactions within the SCO molecules (see Figure 1). Overall, this results in a net overestabilization of the LS state when using D2 (see Table 3), that needs to be corrected through a larger U parameter.

Role of Intermolecular Interactions
Using the benchmarked U value for each compound and each scheme (D3 and D3-BJ, see Table 2), we have computed for isolated molecules of 1-9 ( ) excised (i.e. without further optimization) from the respective solid-state minima. The difference between and is the energetic influence of the crystal packing to the SCO transition through intermolecular interactions of any kind (dispersion, induction, electrostatic, steric). The crystal-packing effects (CPE) reported herein range from −14.1 to +8.0  Table 3). For most compounds, D2 describes that intermolecular interactions lead to an overall stabilization of the LS state (CPE entry is positive in Table 3), whereas D3 corrections describe the contrary (CPE entry is smaller or negative). One possible reason is that D3 and D3-BJ corrections leads to larger unit cells than the D2 in both spin states (although the unitcell expansion upon SCO is similar, see ∆ in Table  3).Indeed, in some cases the volume of the predicted minima is even larger than the reported X-ray structures, which suggests that D3 and D3-BJ might unrealistically overestimate the unit cell volume. This is consistent with other studies showing that dispersion corrected GGAs like PBE-D3, used in periodic plane-wave DFT computations, systematically overestimate the molecule volume by roughly 2%. The reason might the inclusion of the three-body term, which is repulsive for all solids, 72,73 and of the damping functions used in D3 and D3-BJ. 63 To get better insight, we now trace the D2, D3 and D3-BJ energy contribution to at the same set of PBE+U+D3 solid-state minima (see Figure 1). These comparison enables us to identify to what extend the contribution to shown in Table 3 stems from the construction of the different dispersion correction schemes (see section 2.3). Additionally, we have also computed and displayed its intra-molecular component (dashed bars in Figure 1). The difference between filled and dashed bars corresponds to the inter-molecular component. Similar to what we mentioned earlier, the overall picture suggests that the use of D2 leads to a much larger influence of dispersion interactions to the SCO energetics, favoring the LS state (all red bars show positive values in Figure 1). In turn, D3 and D3-BJ greatly diminish their impact up to the point that, for most systems, it becomes negligible.

D3
D3-BJ Figure 2. Error associated with values obtained with various DFT functionals and the D3 (left) or D3-BJ (right) dispersion correction scheme. The reference is the value obtained with PBE+U using the benchmarked U (different for each compound, see Table 2) at isolated molecules of 1-9 excised from the PBE+U+D3 (left) and PBE+U+D3-BJ (right) solid-state minima. Blue (red) indicates that the LS (HS) is too stable according to the functional. See also Tables S4-S6.
Further comparison between compounds would not meaningful due to the different molecular and unit cell composition. A much broader analysis should be carried out to analyze this issue in detail, 45 which is out of the scope of this work.

Gas phase PBE+U vs. Other DFT Functionals
The PBE+U results in periodic boundary conditions, are now complemented with gas-phase calculations using a range of DFT functionals. These are carried out at the molecular geometries of 1-9 excised from the solid-state minima without further optimization. The main motivation of these calculations is to investigate the performance of DFT functionals when they are confronted to sound reference values ( ). We have chosen the exchange/correlation functionals OLYP, OPBE, TPSSh, B3LYP, B3LYP*, M06L and SCAN, (complemented with the D3 and D3-BJ corrections), which have shown different degrees of success in calculating spin-state energies in SCO systems in the past. 28 The results are summarized on Figure 2, and numerical data can be found in Tables S4 and S5. The functionals TPSSh and SCAN correctly reproduce the ground state, but only as a result of a very large overestabilization of the LS state. Therefore, they do not describe compounds 1-9 as SCO-capable. B3LYP, OPBE and M06L perform slightly better, but the former two often fail at assigning the correct ground state. On the bright side, B3LYP* provides very similar spin-state energies to PBE+U, and OLYP achieves the best mean average error (MAE) for the whole dataset (MAE=6.9 kJ/mol). Further insight can be achieved doing a separate analysis of the results for Fe II and Fe III compounds. If we focus on the Fe II dataset, the MAE associated with B3LYP* is only 0.79 kJ/mol, which is quite remarkable. Indeed, B3LYP* was reparametrized specifically to model spin-state energies in Fe II SCO systems, 74,75 and such tiny error does serve to validate the benchmarked U values for 1-5. For the Fe III systems, the PBE+D3 results are in very good agreement with the ones obtained with the B3LYP and the OPBE functionals (MAE=0.84 and 2.85 kJ/mol, respectively). The latter has already demonstrated its accuracy towards SCO systems, 76,77 and in particular towards d 5 SCO molecules before, 78 which yet again reinforces the benchmarked U values for 6-9. In turn, the good performance of B3LYP is somehow surprising, although the higher ionic character (and hence larger ∆ ) of Fe III , with respect to Fe II , should explain the need for a larger HF-exchange percentage (B3LYP 20% vs. B3LYP* 15%).

Transferability
In section 3.1 we have obtained the U values that reproduce for each system. Eventually, when studying new SCO systems, one would like to skip such benchmark and instead use an approximate value that guarantees the smallest possible error. Logically, one would favor using the average U values proposed in section 3.1 for PBE+U+D3 (2.29 and 2.44 eV) and PBE+U+D3-BJ (2.35 and 2.39 eV). To have an estimation of the error that such approach would entail, single-point computations are performed on very same solid state minima (the geometry is not affected by U) using these U values. The resulting are then compared to and the difference is the error associated with the use of those average values (instead of the best parametrization for each compound). The mean absolute error (MAE) is of 4.5 (Fe II ) and 2.8 kJ/mol (Fe III ) for the PBE+U+D3 minima, and of 4.4 (Fe II ) and 2.3 kJ/mol (Fe III ) with PBE+U+D3-BJ (see Table S7). As expected, the error is much larger in systems whose benchmarked U is further away from the average U (1, 4 and 6). In the case of 1, one possible reason is the reported importance of anharmonic effects in its SCO transition. 52 Finally, we notice that our previous parametrization of PBE+U+D2, with U=2.65 eV, yielded a MAE of 5.1 for compounds 1-5, so we can infer that the new parametrization does improve the transferability of the method.

Conclusions
We have parametrized the U values that must be employed to describe the energetics of the thermal SCO transition of five Fe II (1-5) and four Fe III compounds (6-9), under the PBE+U+D3 and PBE+U+D3-BJ methods. The average U values that result: 2.29 (Fe II ) and 2.44 eV (Fe III ) for the former method, and 2.35 (Fe II ) and 2.39 eV (Fe III ) for the latter, can be used reliably to describe other Fe II and Fe III SCO systems with chemical accuracy. Moreover, and due to its local nature, different U values can be applied to different metal ions in the same computation, which enables the accurate description of mixed Fe II -Fe III systems (i.e. double SCO salts 79 ). Concerning the dispersion corrections, we noticed that the use of D3-based dispersion correction schemes leads to a different landscape of contributions to arising from dispersion interactions than with D2. The latter correction seems to overestimate the importance of dispersion interactions to the SCO transition, with a net stabilization of the LS state that must be compensated by larger U values. Finally, we compared the performance of other DFT functionals at describing the spin state energetics of 1-9 in gas-phase, and the geometry of the solid state minima. Both using D3 and D3-BJ dispersion corrections, we found that Fe II and Fe III complexes are better described using B3LYP* and OLYP, respectively (MAE of 3.5 and 2.9 kJ/mol respectively for D3). This finding complements was is known from the existing literature on the subject, and serves to further validate the benchmark of the PBE+U+D3 and PBE+U+D3-BJ schemes.   Figure S1. HS-minima of Fe II Compounds (1)(2)(3)(4)(5) and Fe III Compounds (6-9)

S2. Linear Response U
In the linear response (LR) approach 1 to obtain U, the response function that is calculated is: = In order to obtain U, the inverse of the self-consistent function is subtracted from the inverse of the bare non-interacting function 0 such as: For the LS and HS states of 1, we obtained 0 and from linear repressions of the relationship between the occupations obtained within a range of potential shifts :