Periodic Density Functional Calculations in Order to Assess the Cooperativity of the Spin Transition in Fe ( phen ) 2 ( NCS ) 2

Periodic density functional calculations combined with the Hubbard model (DFT+U) have been performed for the archetype spin crossover complex Fe(phen)2(NCS)2 with phen = 1,2-phenanthroline. The relative energies of the 16 different configurations of two possible spin states for each of the four molecules in the unit cell have been calculated in order to determine from first principles the phenomenological interaction parameter Γ of the Slichter-Drickamer model. These kind of calculations may help to predict important spin crossover characteristics like the abruptness or hysteresis of the transition.


Introduction
A change from an electronic high spin (HS) to a low spin (LS) ground state, commonly called spin crossover (SCO), can be observed in certain transition metal complexes when the temperature is decreased or external pressure is applied.The investigation of the spin crossover phenomenon on its own is a fascinating field of basic research (a renowned review has been edited by Gütlich and Goodwin [1][2][3]).The discovery that some SCO compounds can be reversibly switched between metastable HS and LS states by irradiation with light (light induced excited spin state trapping, called LIESST effect [4]) made SCO complexes promising materials for memory devices with extremely high capacities.
Spin crossover is an intrinsically molecular phenomenon-even in solutions it is possible to switch the spin state by irradiation with light [5].However, it is well known that inter-molecular interactions can sensitively influence the character of the spin crossover.For example, the choice of the counterion of a SCO complex could be used to tune the transition temperature T 1/2 , i.e., the temperature where half of the complexes are in either spin state.Spin transition curves-as the temperature dependent HS-fraction γ HS (T)-exhibit certain features like abruptness, steps or hysteresis only in case of solid samples.
Slichter and Drickamer [6] presented a theoretical model (the so called Slichter-Drickamer model or SD model for short) with only a few empirical parameters, from which many features of the spin transition can be derived.This model belongs to the class of mean field theories since the intermolecular interactions are described by a term that contains only the average magnetization of the crystal (given by the HS-fraction γ HS ) and a phenomenological interaction parameter Γ.The latter is positive if neighboring molecules tend to have the same spin-and negative in the opposite case.In a crystal with diluted SCO complexes, and, therefore, no correlation between the spin states of different molecules, Γ should vanish.Within the SD model the free enthalpy per molecule can be written as where the mixing entropy S mix is defined by and ∆E HL and ∆S HL are the differences of the total electronic energy and entropy, respectively, between the HS and LS states.The contribution of volume expansion, p∆V, can be neglected at ambient pressure (usually it amounts to less than 1 J mol −1 [7]).The free enthalpy of the pure LS phase has been arbitrarily set to zero: G(0) = 0. Minimization of G(γ HS ) with respect to γ HS allows to determine the HS-fraction as a function of temperature.
The essential parameters of the SD model-∆E HL , ∆S HL and Γ-may be obtained by fitting experimental data.In this way the SD model may help to interpret and analyze the experiment.If, on the other hand, it would be possible to calculate these parameters from first principles, the SD model might be used to aid the design of new materials, e.g., optical storage devices.Since the pioneering study of Bolvin [8] who investigated the spin state splitting of an FeN 6 octahedron with Hartree Fock (HF) calculations a large number of studies have been published which report about electronic structure calculations on SCO complexes (see [7,[9][10][11] for an overview).Wavefunction based methods have the advantage that in principle any order of accuracy can be reached.However, since the spin crossover is the result of a delicate balance of different contributions to the electronic energy, reasonably accurate results for ∆E HL can only be reached with sophisticated post-HF methods, which can currently only be applied to systems with a few atoms [12].For this reason practically all electronic structure calculations on realistic SCO systems have been performed with methods utilizing density functional theory (DFT), either with pure density functionals or with hybrid functionals which include some amount of HF exchange like B3LYP* [13].In the first years these calculations were restricted to SCO complexes in the gas phase.Since the entropy difference ∆S HL is almost purely of molecular origin, this parameter can in some cases be estimated with acceptable accuracy by calculations for the SCO molecule in the gas phase (see for example [14][15][16][17]).
The calculation of ∆E HL turned out to be more difficult.Tuning the amount of HF exchange in hybrid functionals as well as detailed investigations of the performance of different density functionals with regard to the accuracy of ∆E HL [18][19][20][21][22][23][24][25][26] led to results that are at least in the right order of magnitude.This is especially true if the effect of small ligand modifications should be predicted [27][28][29].All methodological progress, however, cannot change the fact that due to inter-molecular interactions the electronic energy difference in the crystal might significantly deviate from the corresponding value for a complex in the gas phase.The phenomenological parameter Γ is of course by definition out of the reach of a calculation for a molecule in vacuo.
In order to assess also cooperative effects two approaches have been chosen: The first one is given by molecular calculations for fragments of polymeric SCO materials [30][31][32][33] as well as calculations for polynuclear SCO complexes [34][35][36][37], whereas the second one consists of electronic structure calculations with periodic boundary conditions.The latter approach was feasible for SCO materials only in the last decade, a first study was published in 2006 by Lemercier et al. [38], who used the local density approximation (LDA) and atom centered basis functions.While part of the periodic DFT studies are concerned with networks of SCO complexes (see for instance [19,[39][40][41]) others are focussed on crystals of monomeric SCO complexes and many of them [12,17,[42][43][44][45] investigated the complex Fe(phen) 2 (NCS) 2 with phen = 1,2-phenanthroline (Figure 1), a complex which-like the majority of the SCO complexes known so far-contains a central iron(II) ion, octahedrally coordinated by six nitrogen atoms.This complex was one of the first SCO complexes discovered at all [46,47] and belongs to the SCO complexes most thoroughly studied, experimentally as well as theoretically, and is for this reason often denominated as an archetype SCO complex.It undergoes an abrupt and complete transition from the S = 2 HS state to the S = 0 LS state when the temperature is decreased below the transition temperature T 1/2 ≈ 176 K.The attention this complex has attracted is also due to the fact that it was the first SCO complex for which the crystal structure was available for the HS and for the LS phase [48].A more recent study by Legrand et al. [49] presents the crystal structure of the LS and of the (metastable) HS phase at 15 K. Fe(phen) 2 (NCS) 2 crystallizes in the space group Pbcn, with four complexes in the unit cell.At 15 K the volume of the unit cell increases by about 61.3 Å 3 upon transition from the LS to the HS phase which corresponds to an increase of 2.8%, a typical value for the volume expansion upon spin crossover [1].Fe(phen) 2 (NCS) 2 is also one of the few SCO complexes where the entropy difference ∆S HL , which is built up out of two components, an electronic one, ∆S el , and a vibrational one, ∆S vib , has been most accurately determined, experimentally as well as theoretically [14][15][16]50,51].The electronic component is due to the higher multiplicity of the HS state and can be easily calculated exactly (electronic excitations can be neglected in good approximation), while the vibrational component stems from the fact that loosened metal ligand bonds in the HS state lead to lower vibrational frequencies and-at elevated temperatures-to an increase of entropy.
While, at least in the case of Fe(phen) 2 (NCS) 2 , ∆E HL and ∆S HL are well accessible by molecular and periodic electronic structure calculations, Γ, the parameter describing cooperativity, is more difficult to access.In a recent study from Sinitskiy et al. [52] Γ has been derived from atom-atom potentials.In this study a different approach was been chosen.Periodic LDA+U calculations were performed for different configurations of the spin states of the four molecules of the unit cell of Fe(phen) 2 (NCS) 2 and the resulting energies were used to estimate Γ as described below.At zero temperature the entropic terms in this equation vanish and it can be rewritten as where the total electronic energy per molecule in the pure LS phase (γ HS = 0) has been arbitrarily set to zero.The index c denotes the configuration of the crystal with respect to the spin states (HS or LS) of the individual SCO molecules.In general, for a crystal with N SCO molecules there are 2 N possible configurations and the HS-fraction γ HS could be written as a function of the configuration c.While there is only one configuration that corresponds to the pure LS phase (γ HS = 0) and the same is valid for the pure HS phase, there are ( N N/2 ) different configurations for the critical HS-fraction γ HS = 1/2.In the latter case not all configurations c are energetically favorable to the same degree and the electronic energy E(γ HS , c) is a function of c.Since the HS-fraction γ HS (c) is a function of c too, we can omit γ HS as an argument and write the electronic energy shortly as E(c).Since the left side of Equation ( 3) is varying with the configuration c even at constant γ HS , the same must be true for the right side of this equation resulting in a configuration dependent parameter for configurations c with 0 < γ HS < 1 (if γ HS = 0 or γ HS = 1 the term γ HS (γ HS − 1)Γ in Equations ( 1) and (3) vanishes and Γ becomes irrelevant).The macroscopic interaction parameter Γ that is used in the SD model as defined in Equation ( 1) can now be written as the weighted sum where the weights w(c) = f (c)/Z are given by the Boltzmann factor normalized by the partition function Z = ∑ c f (c) (the numerator of the fraction in the preceding equation reflects the fact that at the critical temperature both spin states are in equilibrium and the Boltzmann weight of the pure LS phase should be the same as the one of the pure HS phase).
For sufficiently steep transitions it is well justified to replace the temperature T in the exponent of the Boltzmann factor by the critical temperature T 1/2 as it is done here, since we have γ HS ≈ 0 for T T 1/2 and γ HS ≈ 1 for T T 1/2 (in case of a pronounced gradual transition the temperature T has to be written as a function of the HS-fraction leading to a system of equations that has to be solved self consistently).In this way the cooperativity of a SCO complex can be predicted by first principles calculations without need of empirical potentials like in the approach from Sinitskiy et al. [52].The disadvantage however is, that the approach presented here is computationally much more demanding and for this reason restricted to the unit cell or at maximum to small super cells, which might reduce the accuracy.

Results and Discussion
From variable cell geometry optimization with U = 2.5, 4.0, 5.6 and 6.0 eV for the pure LS and the pure HS phase an almost perfectly linear dependence between ∆E HL and U was obtained.For U = 5.6 eV the total energy difference amounts to ∆E HL ≈ 9.77 kJ mol −1 fitting well to the experimentally determined 9 kJ mol −1 from Sorai and Seki [53].A similar linear dependence between ∆E HL and U was found after including the Grimme-D2 dispersion correction [54].The obtained Hubbard U of 5.6 eV is at the upper limit of the range of values (2.5 to 5.0 eV) that were reported [12,17,[42][43][44][45] for this parameter in case of calculations for Fe(phen) 2 (NCS) 2 .However, all studies on this complex that employed a small Hubbard parameter of about 2.5 eV [12,17,42] were performed with gradient corrected density functionals which should better account for the localization of the strongly correlated iron 3d electrons and thus have less need for the Hubbard correction.
The volume expansion upon complete spin transition was calculated to 65 Å 3 and 76 Å 3 with and without Grimme-D2 correction, respectively, which is in reasonable agreement with the experimental result of 79 Å 3 [49].However, with dispersion correction the volume of the HS unit cell results to 1832 Å 3 , whereas without this correction 1996 Å 3 were obtained which is closer to the experimental value of 2248 Å 3 .This finding might be traced back to the notorious underestimation of bond lengths by LDA methods [55].Since dispersion corrections usually lead to smaller calculated volumes, neglecting dispersion in LDA calculations might lead to a partial cancellation of errors with respect to the volume.The calculations discussed in the following were performed without dispersion correction.
In order to assess the cooperativity of spin transitions using first principles methods the 16 possible spin configurations of the four molecules in the unit cell of Fe(phen) 2 (NCS) 2 have been studied (see also Figure 2).For each configuration c (labeled with a four letter code like LLHL for a configuration with three LS molecules and one HS molecule in a given order) a complete variable-cell geometry optimization has been performed.The resulting total energies per molecule (Table 1) clearly do not depend linearly on the HS-fraction γ HS (in other words the number of HS molecules divided by the number of molecules in the cell).If the energy needed to switch a molecule from the LS to the HS state would be independent from the spin state of the neighboring molecules, the total energy per molecule of configurations with γ HS = 0.5 (like LLHH or LHLH where the four letters stand for the spin states of the four molecules of the unit cell, L indicating a molecule in the low spin state and H indicating a molecule in the high spin state) should be ∆E HL /2 ≈ 4.9 kJ mol −1 .Instead we found significantly larger values of in average 7.0 kJ mol −1 , indicating that different spin states for neighboring molecules are energetically not favorable.In the SD model this should be reflected by a positive parameter Γ.
For the four configurations that correspond to γ HS = 0.75 (the two rightmost columns in Table 1) as well as for the configurations LHLH and HLHL (called group 1 in the following ) values for Γ(c) are obtained with the help of Equation ( 4) that range from 4.9 to 5.7 kJ mol −1 as compared to values between 8.5 and 12.1 kJ mol −1 for the remaining eight configurations (group 2).The configurations of group 1 have in common that one yz-plane of the unit cell is filled with HS molecules and that, therefore, the size of the unit cell is close to the size in the pure HS phase.In this way there is sufficient space for the HS molecules and the electronic energy per HS molecule is in average 11.6 kJ mol −1 , roughly 1.8 kJ mol −1 more than ∆E HL .In each configuration of group 2 on the other hand there is at least one LS molecule in each of the two yz-planes of the unit cell and the volume of the cell is less than in the pure HS phase.Therefore, there is insufficient space for the HS molecules and their electronic energy is in average 16.6 kJ mol −1 , which is about 6.8 kJ mol −1 more than ∆E HL .As a result the correlation between the spin states of neighboring molecules is less in group 1 than in group 2 (no correlation would correspond to an electronic energy per HS molecule of exactly ∆E HL ) and consequently the interaction parameter Γ(c) for the configurations of group 1 is smaller than for the configurations of group 2.  As an unweighted average value we would get Γ ≈ 8.4 kJ mol −1 , which is more than twice as large as the experimental value of 3.0 kJ mol −1 from Sinitskiy et al. [52] but in the same order of magnitude.A weighted average according to Equation ( 5) gives a value of 7.3 kJ mol −1 that fits somewhat better to the experimental value but is still more than twice as large.One quite obvious reason for this deviation is the small cell with only four SCO molecules that has been used for the present study.Studies that employ the Ising model together with Monte Carlo simulations on large supercells indicate that large clusters of equal spins are energetically favorable to configurations with neighboring molecules with different spin states [56].Finally, the Γ parameter predicted by Sinitskiy et al. with a completely different theoretical approach (between 5 and 8 kJ mol −1 ) is larger than the experimental value as well-thus leaving room for the speculation that the SD model does not perfectly describe the spin transition of Fe(phen) 2 (NCS) 2 .A possible reason for this may be the influence of inter-and intra-molecular mode coupling on the cooperative behavior of the spin transition observed by Mebs et al. [57].

Computational Details
All calculations were performed with the software suite Quantum ESPRESSO (Version 5.2.0, QE Foundation, Trieste, Italy, 2015) [58].Valence electrons were represented explicitly by plane waves with cutoffs of 50 Ry for the waves and 400 Ry for the charge density.Self consistent wavefunctions and forces were converged to an accuracy of 10 −7 and 10 −5 atomic units, respectively.Less strict cutoffs for the forces influenced adversely the accuracy of the calculated interaction parameter Γ.The Brillouin zone sampling was restricted to the Γ point.The interactions between the valence electrons were described within the local density approximation with the help of the exchange and correlation functional from Perdew and Zunger [59].Ultrasoft pseudopotentials [60] were employed [61] to treat the core electrons implicitly by replacing the all-electron Hamiltonian by a pseudo-Hamiltonian that describes the interactions between the ions and the valence electrons.In order to describe better the localization of the strongly correlated iron 3d electrons the DFT+U model with the simplified scheme of Cococcioni et al. [62] was used with a Hubbard U of 5.6 eV if not explicitly stated otherwise.For some calculations in addition the Grimme-D2 correction [54] was applied to include the influence of the otherwise poorly described London forces.
Since the electronic structure calculations were performed for the state at zero temperature the X-ray structures from Legrand et al. [49] was used as starting points, where the structure of the LS-and of the metastable HS-phase have been determined at T = 15 K.

Conclusions
From first principles calculations for different spin configurations of the unit cell of a SCO complex a reasonable estimate for the phenomenological interaction parameter Γ was obtained that allows to describe the cooperativity of the spin transition within the SD model.Future calculations will have to show whether a better agreement between the theoretical and the experimental value for Γ can be reached, for instance by choosing a gradient corrected density functional together with a dispersion correction or by performing calculations for supercells.In any case such kind of calculations might help to predict important characteristics like the abruptness and the hysteresis of the spin transition.Beyond the validity of the SD model the type of calculations presented here may also be useful to deliver input parameters to more sophisticated models that describe the spin transition in SCO materials.

Figure 1 .
Figure 1.Calculated structure of the low spin (LS) phase of Fe(phen) 2 (NCS) 2 with crystallographic axes and unit cell.One molecule is represented by balls and sticks (color scheme: Fe (red), S (yellow), N (blue) and C (gray); hydrogens are omitted for clarity), the other molecules are represented by wireframes in different colors.For the two molecules with x coordinates equal to zero for the iron centers also the symmetrically equivalent molecules with x = 1 are shown, resulting in Z + 2 = 6 depicted molecules in total.

Figure 2 .
Figure 2. 16 different spin configurations of a unit cell: (a) one configuration with γ HS = 0; (b)-(e) four configurations with γ HS = 0.25; (f)-(k) six configurations with γ HS = 0.5; (l)-(o) four configurations with γ HS = 0.75 and (p) one configuration with γ HS = 1.Low spin (LS) and high spin (HS) molecules are colored in blue and red, respectively, and only the central FeN 6 octahedron is shown.As in Figure 2 Z + 2 molecules are shown in total for each configuration.The four letter codes (like LLHH or LHLH) stand for the spin states of the four molecules of the unit cell, L indicating a molecule in the low spin state and H indicating a molecule in the high spin state.

Table 1 .
Total electronic energies per molecule, E(c), in kJ mol −1 for different spin configurations c given by four letter codes (like LLHH or LHLH) as defined in Figure2.The energy of the LLLL configuration (γ HS = 0) is set to zero.