PCCP View Article Online

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

PAPER

Cite this: DOI: 10.1039/c5cp02502h

View Journal

Towards an accurate and computationally-efficient modelling of Fe(II)-based spin crossover materials Sergi Vela,*ab Maria Fumanal,b Jordi Ribas-Arinob and Vincent Roberta The DFT + U methodology is regarded as one of the most-promising strategies to treat the solid state of molecular materials, as it may provide good energetic accuracy at a moderate computational cost. However, a careful parametrization of the U-term is mandatory since the results may be dramatically affected by the selected value. Herein, we benchmarked the Hubbard-like U-term for seven Fe(II)N6-based pseudo-octahedral spin crossover (SCO) compounds, using as a reference an estimation of the electronic enthalpy difference (DHelec) extracted from experimental data (T1/2, DS and DH). The parametrized U-value obtained for each of those seven compounds ranges from 2.37 eV to 2.97 eV, with an average value of U = 2.65 eV. Interestingly, we have found that this average value can be taken as a good starting point since it leads to an unprecedented mean absolute error (MAE) of only 4.3 kJ mol1 in the evaluation of DHelec for the studied compounds. Moreover, by comparing our results on the solid state and the gas phase of the materials, we quantify the influence of the intermolecular interactions on the relative stability of the HS

Received 29th April 2015, Accepted 19th May 2015

and LS states, with an average effect of ca. 5 kJ mol1, whose sign cannot be generalized. Overall, the

DOI: 10.1039/c5cp02502h

crystalline phase of SCO compounds, or the adsorption of individual molecules on organic or metallic

www.rsc.org/pccp

accuracy that is dramatically missing when using bare-DFT functionals.

findings reported in this manuscript pave the way for future studies devoted to understand the surfaces, in which the rational incorporation of the U-term within DFT + U yields the required energetic

1. Introduction Spin crossover (SCO) compounds are molecular materials that can be stable in two electronic states with different magnetic response, and whose relative stability can be interconverted by the application of external stimuli such as light, heat or pressure, among others.1–4 This quality makes them ideal candidates to be exploited in new electronic and spintronic devices5–7 and, consequently, these materials have received a great deal of attention in recent years. The field of SCO compounds has been traditionally dominated by transition-metal compounds based on Fe(II),8–10 Fe(III)11 or Co(II)12 ions coordinated to organic ligands in a pseudo-octahedral fashion. In those materials, the High Spin (HS) and Low Spin (LS) states of the complex are very close in energy due to the equilibrium between different contributions: while enthalpy favors the stability of the LS state, the larger vibrational and electronic entropy of the HS state stabilizes the latter at high temperature. This characteristic can be used to originate a thermally-driven transition between the LS and HS states at a critical temperature a

Laboratoire de Chimie Quantique, Universite´ de Strasbourg, 4 rue Blaise Pascal, F-67000 Strasbourg, France. E-mail: [email protected] b Departament de Quı´mica Fı´sica and IQTCUB, Universitat de Barcelona, Av. Diagonal 645, 08028, Barcelona, Spain

This journal is © the Owner Societies 2015

(T1/2), where the two contributions are equal. Although the SCO phenomenon can be explained at the molecular level, the macroscopic observation is greatly influenced by how the individual transiting units interact in the solid state when the spin transition occurs. The degree of cooperativity of the SCO units defines the shape of the transition, (soft, abrupt, hysteretic)10 and may even originate a region of bistability that is highly appreciated especially if found at room temperature, as it implies a memory effect that is of potential interest in technology.13 The computational modelling of SCO materials is a challenge and, despite some work has been done to rationalize the crystal packing effects,14–16 it has been traditionally limited to the study of the transiting units at the molecular level.17 In our opinion, one of the reasons is the absence of a computational scheme able to simultaneously treat the electron correlation effects (crucial in this kind of systems) and the contributions arising from the crystal packing with enough accuracy and at an affordable cost. The description of electronic correlation can only be properly achieved by using post-HF methods such as CCSD or CASPT2, but the application of these techniques is restricted to systems containing a few atoms; this is clearly at odds with the purpose of describing the solid state of those materials, which requires the use of large cluster models or periodic boundary conditions. Nowadays, the only electronic structure methodology capable of dealing with unit cells with

Phys. Chem. Chem. Phys.

View Article Online

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

Paper

hundreds of atoms, even if at the expense of the proper description of electron correlation, is Density Functional Theory (DFT). However, it also has its drawbacks and so far any DFT functional has been found to provide systematic good HS–LS adiabatic energy gaps, a key magnitude in SCO compounds.18,19 In order to improve the performance of DFT, two main strategies have been used: the use of a customized amount of Hartree– Fock exchange within hybrid functionals,20–22 or the incorporation of a Hubbard-like the U-term.23–25 Both empirical corrections may be tuned to provide correct results but the use of hybrid functionals, despite being a more robust scheme, is still too demanding for periodic large-size calculations. Therefore, and while better strategies are being developed, complementing the bare LDA or GGA functionals with a Hubbard-like U-term in the so-called DFT + U approach seems to be the most-promising strategy to treat the solid state of SCO systems with an energetic accuracy comparable to that of wave-function methods (see Section 4.1 for an introduction to DFT + U). However, a careful parametrization of the U-term is required to perform DFT + U calculations since the results may be dramatically affected by the selected value. In order to circumvent this problem, U is usually chosen to reproduce the experimental property of interest, which penalizes the potential predicting capabilities of computational chemistry, and may be used to obtain any beforehand decided result. In fact, it may even lead to a different value for the same compound depending on the studied property. For instance, a Ueff† = 1.55 eV on the iron ion ´ngya ´n and coworkers23,24 to study was empirically selected by A the crystal structure of Fe(phen)2(NCS)2, whereas Gueddida et al.26 used Ueff† = 2.10 eV to study the adsorption of the same compound in a metallic substrate. To eliminate the arbitrariness of this choice, in a recent paper some of us benchmarked the U-term in order to reproduce the CASPT2 adiabatic energy gap between the HS and LS molecules of the [FeII(1-bpp)]2+ compound at the gas phase.25 We then used PBE + U, with a benchmarked value of U = 2.50 eV, to study the crystalline phase of this compound. As can be seen, such a U value was slightly different from those employed in previous studies of similar compounds (1.55 eV and 2.10 eV). At this point, it became clear that it was necessary to further investigate the appropriate parametrization of the U-term to study those compounds, with particular interest in knowing whether it is necessary to benchmark the DFT + U methodology for every SCO system or if a more-approximate, and transferable, U value yields good agreement with the experiment when applied to a variety of SCO compounds. With this aim in mind, in this manuscript we have focused our attention on SCO materials based on Fe(II), displaying pseudo-octahedral six-fold coordination and having N atoms in the first coordination sphere (FeN6 core), probably the largest family of transition-metal-based SCO materials. DFT + U allows us to work on the crystalline phase of the materials, whereas other methodologies are restricted to work † Ueff = U  J, (where J is a Stoner-like exchange parameter), alternatively utilized in other implementations. For the results obtained in present paper, J has been set to 0 and, thus, Ueff = U.

Phys. Chem. Chem. Phys.

PCCP

at the gas phase. Given that our calculations already incorporate the important crystal packing effects, we can parametrize the DFT + U scheme directly with experimental measurements. This implies one clear advantage: our results can be compared with magnitudes that are accessible experimentally, such as DH and T1/2, usually extracted from DSC and SQUID experiments, respectively. In doing so, we eliminate the problem of choosing a reference methodology, since none of the current ones is fully satisfactory to predict the energy of the low-lying spin states of the studied SCO systems: Coupled-Cluster (CC) calculations are often taken as reference but their application is restricted to model systems,19 not to mention that they are limited to compounds displaying strong monoreferential character.27 CASPT2 has been the workhorse of the chemistry community for many years, and has been successful in the quantification of adiabatic28,29 and vertical30–32 energy gaps, but a recent error detected in the Molcas package33 suggests that a new re-parametrization of the IPEA shift is necessary, which is beyond the scope of this paper. Other perturbative treatments have been proposed to complement the CAS wave-function, such as the NEVPT2 method,34,35 but its performance on SCO systems has never been tested before and, thus, it can be hardly taken as a reference. Therefore, the fine-tuning of the DFT + U methodology has been done here taking experimental data as a reference. Consequently, we had to choose SCO materials with (i) a thermally-driven transition at moderate temperatures, (ii) available HS and LS crystal structures, and (iii) measured thermodynamic data. Finally, after inspection of the literature, six candidates were identified, whose chemical structures are represented in Fig. 1 and with chemical formulae: Fe(phen)2(NCS)2 (1),36 [Fe(abpt)2(NCS)2] (2),9 [Fe(abpt)2(NCSe)2] (3),9 Fe(bapbpy)(NCS)2 (4),37 Fe(HB(pz)3)2 (5)38 and Fe[H2B(pz)2]2(bipy) (6).39 Apart from them, we have decided to benchmark the U value for compound [FeII(1-bpp)][BF4]2 (7)40 based on this new strategy, and to confront it to the value of U = 2.50 eV that some of us reported recently.25 The present manuscript is organized as follows: first, in the Methodology section we comment on the computational strategy that we have employed to parametrize the U value for compounds 1–7, with special attention to the vibrational and electronic contributions to enthalpy and entropy. These thermodynamic variables define the spin transition, and allow us to connect our calculations with the experimental measurements. Then, in the Results and Discussion section, we present the U values that provide a better agreement between DFT + U and the experiment, and compare our calculations at the molecular and solid state levels to evaluate the influence of intermolecular interactions on the spin transition of those materials. The results contained in this section demonstrate that the DFT + U scheme is system dependent, even for materials that share many ingredients in their architecture. In particular, for the studied Fe(II)N6-based pseudo-octahedral SCO compounds, U is found within the range of 2.37–2.97 eV, with an average of U = 2.65 eV. However, we also demonstrate that this average value may be taken as a good starting point when a more accurate benchmark is not possible, or when direct comparison with

This journal is © the Owner Societies 2015

View Article Online

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

PCCP

Paper

Fig. 1 High-spin structures of the seven compounds studied herein, namely: Fe(phen)2(NCS)2 (1), [Fe(abpt)2(NCS)2] (2), [Fe(abpt)2(NCSe)2] (3), Fe(bapbpy)(NCS)2 (4), Fe(HB(pz)3)2 (5), Fe[H2B(pz)2]2(bipy) (6) and [FeII(1-bpp)][BF4]2 (7). Color code for atoms: C (black), N (pale blue), H (rose), S (yellow), Se (green), B (magenta), Fe (brown).

experiment is not available. Indeed, using a value of U = 2.65 eV provides excellent agreement with experiment for compounds 1–7, with a mean absolute error (MAE) of only 4.3 kJ mol1 in the enthalpy difference between the HS and LS states.

2. Methodology A thermally-driven spin transition occurs when the higher enthalpy of the LS state is overcome by the favorable entropy contribution of the HS state at T1/2. At this temperature, DG is equal to zero and the transition temperature T1/2 can be obtained from:   DHtot T1=2   (1) T1=2 ¼ DStot T1=2

DStot(T) = DSelec + DSvib(T)

The total DHtot and DStot values can be obtained experimentally, and are known for the systems studied herein. Following the harmonic-oscillator approximation, the expression for Hvib and Svib can be found in eqn (4) and (5), where it can be seen that their value depends on temperature. Those magnitudes (Hvib and Svib) can be evaluated from first principles calculations provided that the frequencies of the vibrational normal modes are obtained (ni). In turn, Selec (eqn (6)) can be considered, as a good approximation, to be temperature-independent with values 13.38 and 0 J K1 mol1 for the HS state (S = 2) and the LS state (S = 0), respectively. Hvib ¼

As defined in statistical thermodynamics, both enthalpy and entropy can be decomposed in multiple terms accounting for electronic, rotational, vibrational and translational differences between the two spin states involved in the transition. In the present study, we have neglected the contribution of the rotational and translational terms, as they are not expected to contribute significantly, and considered that only the electronic and vibrational terms contribute to enthalpy and entropy (eqn (2) and (3)). DHtot(T) = DHelec + DHvib(T)

This journal is © the Owner Societies 2015

(2)

(3)

Svib ¼

Nvib  X hn i i¼1

 Nvib  X 1 hn i ehvi =kB T hn i þ 2 1  ehvi =kB T i¼1 1

T ehn i =kB T  1

    kB ln 1  ehn i =kB T

Selec = R ln(2S + 1)

(4)

(5)

(6)

The fact that the electronic entropy is constant and equal for most systems allows us to evaluate the vibrational contribution

Phys. Chem. Chem. Phys.

View Article Online

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

Paper

PCCP

to entropy from the experimentally accessible total entropy. In contrast, DHelec‡ and DHvib cannot be directly extracted from the experimental estimate of DHtot. Therefore, we had to develop a strategy in order to estimate DHelec from experiments, using the known values of DHtot and DSvib: first, we have evaluated the ni values of the HS- and LS-optimized structures using standard ab initio methodology. Then, those values have been slightly modified to reproduce the experimental DSvib at T1/2 (see Section 4.2). After, those adjusted frequencies were used to estimate DHvib at the same temperature and, finally, DHelec is quantified according to DHelec = DHtot(T1/2)  DHvib(T1/2). Since this quantity is used all along the paper as our reference value, it will be from now on referred to as DHref elec. In turn, the computational strategy to obtain DHelec is the following. First, we perform variable-cell geometry optimizations, starting from the nuclear coordinates of the unit cells reported experimentally for the HS and LS states as initial coordinates. Then, the energy difference between the minima of the two spin states can be calculated using different U values to obtain the evolution of DHelec as a function of U. Finally, using a linear interpolation, one can extract the U value that reproduces DHref elec (see Section 4.3).

3. Results and discussion 3.1.

Estimation of DHref elec from experimental measurements

Following the working strategy defined in the previous section, we have estimated DHvib and DHref elec for compounds 1–7 (see Table 1). As previously mentioned, DHref elec has been extracted from the experimental value of DHtot and our calculated value of DHvib, according to eqn (2). It is worth noticing the importance of the DHvib term at the transition temperatures. As correctly pointed out by Rudavskyi et al.,28 this contribution (eqn (2)) has often been ignored, and implies an important relative stabilization of the HS state that reduces the LS–HS energy gap (for some systems, DHvib is as large as DHtot at T1/2). As can be seen in Fig. 2, DHtot depends considerably on temperature as a consequence of the vibrational enthalpy term (DHvib for compounds 1 and 2 ranges from 0.8 and 7.9 kJ mol1 at 5 K to 5.6 and 12.8 kJ mol1 at 300 K, respectively), which puts into question the approximation of DHvib and, consequently, DHtot, as temperature-independent quantities. The resulting DHref elec value is our best estimate of the adiabatic energy difference between the HS and LS states of compounds 1–7. It must be stressed that it has been calculated from the minimum energy structures of the HS and LS polymorphs in the crystal and, thus, includes the influence of the intermolecular interactions on both the resulting energies and structures. 3.2.

Benchmark of the U parameter to reproduce DHref elec

At this point, once we have estimated the experimental DHref elec values, those are taken as a reference for the subsequent parametrization of the U parameter within the DFT + U methodology. To this purpose, we have performed variable-cell geometry optimizations using a value of U = 2.50 eV to obtain the ‡ Note that DHelec for the HS–LS transition is often referred to as adiabatic energy difference, given that DH = DE at a constant pressure and volume.

Phys. Chem. Chem. Phys.

Table 1 Experimental and calculated thermodynamic data for compounds 1–7. Temperature is given in K, enthalpy in kJ mol1 and entropy in J K1 mol1

Experimental measurements Compound

T1/2

1 2 3 4 5 6 7

176 180 224 183a B350 160 259

DSvibr(T1/2) DHtot(T1/2) 35.6 19.1 24.6 34.6 36.6 70.5 52.8

Computational estimation DHvib(T1/2)

DHref elec

6.4 7.4 5.9 7.1 3.9 7.5 1.9

15.0 13.2 14.5 16.8 21.4 20.9 19.1

8.6 5.8 8.6 9.7 17.5 13.4 17.2

a

This compound undergoes a two-step transition, and the experimental data included herein refers to the first one upon heating.

optimized HS and LS crystal structures. It must be mentioned here that we have observed that the minimum energy structure does not significantly change when obtained using different U values (at least, within the 2.0–3.0 eV range). Therefore, single point calculations using three different U values (2.0, 2.5 and 3.0 eV) were sufficient to obtain the evolution of DHelec as a function of U (see Fig. 3 and Section 4.2 for further details). Table 2 Computed DHelec values using three different U parameters and 1 DHref elec for compounds 1–7 (in kJ mol ). It is also given the U value (in eV) ref that reproduces DHelec

Compound

DHU=2.0 elec

DHU=2.5 elec

DHU=3.0 elec

DHref elec

Benchmarked U

1 2 3 4 5 6 7 Average

26.7 34.7 38.7 36.5 55.0 38.5 40.5

10.3 18.6 22.8 20.4 37.7 22.7 22.4

3.3 3.0 6.7 4.1 20.4 6.2 4.7

15.0 13.2 14.5 16.8 21.4 20.9 19.1

2.37 2.68 2.76 2.61 2.97 2.55 2.60 2.65

The results obtained have been collected in Table 2, together with the target DHref elec values. As can be seen in Fig. 3, the variation of DHelec with respect to a change in U follows an almost-perfect linear dependence (the same has been recently observed when changing the amount of HF exchange in hybrid functionals20,22). This was expected from the expression of EU given in the Computational details section (see eqn (9)). Indeed, the slope of the lines Is shown in Fig. 2 must be associated with the term lIs i (1  li ) in Is eqn (9), where li is the occupation of the orbitals to which the U correction is applied, in our case the five 3d-shell orbitals of the Fe(II) ion. The different, although similar, slopes for the seven compounds studied herein indicate that the occupation of those orbitals differs from one compound to the other. This could be associated with a different hybridization of 3d orbitals of the Fe ion with the 2p orbitals of the coordinating N atoms, suggesting that there may exist some correlation between the structure of the FeN6 core and the slope of these curves and, thus, the proper U of each compound. This issue will be left for future inspections. In any case, this linear dependence allowed us to easily estimate which U value reproduces the corresponding DHref elec for compounds 1–7. Interestingly, all the benchmarked U values using this strategy are found between 2.37 eV for 1 and 2.97 eV for 5,

This journal is © the Owner Societies 2015

View Article Online

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

PCCP

Paper

Fig. 2 Temperature dependence of Htot for the HS (blue) and LS states (red) of compounds 1 and 2, together with the associated DHtot (black), LS calculated as: DHtot = HHS tot  Htot.

compounds with a notable degree of accuracy (MAE = 4.3 kJ mol1). The MAE reported here for the DFT + U methodology using a U value of 2.65 eV is much better than the one reported for the best bare-DFT functional: TPSSh, a hybrid functional with a 10% of HF exchange (4.3 vs. 11 kJ mol1).41 Therefore, one can conclude that this average value can be taken as a good starting point in cases for which a more accurate benchmark is not possible, or for which experimental data are not available. 3.3.

Fig. 3 Dependence of the DHelec term with U, for compounds 1–7. The empty squares represent DHref elec. Note that the orange and red lines, corresponding to compounds 3 and 6, respectively, are overlapped.

Table 3 Computed DHelec values using the average benchmarked U value of 2.65 eV (DHU=2.65 ) together with DHref elec elec and the absolute error

Compound

DHU=2.65 elec

DHref elec

Absolute error

1 2 3 4 5 6 7 Mean absolute error (MAE)

6.8 14.0 17.9 15.5 32.5 17.6 17.2

15.0 13.2 14.5 16.8 21.4 20.9 19.1

8.3 0.8 3.4 1.4 11.1 3.3 1.9 4.3

with an average value of 2.65 eV. It must be mentioned that the proper U value for compound 7 found using this strategy (2.60 eV) is in very good agreement with our previous estimate of 2.50 eV, benchmarked using CASPT2 calculations in the gas phase.25 In contrast, the proper U value for compound 1 is considerably different than the one empirically selected in previous computational work: 1.55 eV23,24 and 2.10 eV.26 At this point, one would like to assess the transferability of this average value U = 2.65 eV in the estimation of DHelec for compounds 1–7. From Fig. 3, the value of DHU=2.65 is readily available and is reported in Table 3, together elec with the MAE associated with those values. It can be seen that our benchmarked value is able to reproduce DHelec for all

This journal is © the Owner Societies 2015

Influence of the intermolecular interactions on DHelec

In the following, we will quantify the effect of the intermolecular interactions on the DHelec term of the studied compounds. To this purpose, DFT + U provides an excellent framework given its capacity to describe the crystalline phase of the materials. The goal of this section is not to identify what kind of interactions affects this value, but to quantify them as a whole. A deep study on the effect of the crystal packing in SCO compounds can be alternatively found in other references,16,23 ref. 25 being focused on compound 7. Regarding the accuracy of those values, it must be mentioned that the contribution of the intermolecular interactions may depend on the chosen U value, as it modifies the electron density in the SCO molecules. However, this effect cannot be much relevant within the U = 2–3 eV range and, in any case, the results shown in this section have been obtained using the benchmarked U values for each compound. Instead, we believe that the values reported essentially depend on the quality of the Grimme-D2 dispersion correction coupled with PBE. Although some problems have been reported in the application of this empirical correction when dealing with metallic surfaces,42,43 it has recently produced excellent results for structures and cohesion energies of molecular crystals.44,45 In Table 4 we show a comparison between the DHelec values calculated in the gas phase (i.e. isolated units) and in the crystalline phase. The difference between both values stems from the intermolecular interactions between SCO molecules (or SCO-counterion interactions) in the crystal, absent when the molecules are isolated. The influence of those interactions on the relative stability of the HS and LS states ranges from 8.1 kJ mol1 in 1 to 6.5 kJ mol1 in 4, the average influence being ca. 5 kJ mol1. Indeed, for compounds 1–3 this contribution accounts for ca. 50% of this

Phys. Chem. Chem. Phys.

View Article Online

Paper

PCCP

Table 4 Computed DHelec values for compounds 1–7 at the crystalline and gas phase levels, together with the difference between those values, which must be ascribed to the influence of the intermolecular interactions. All values are given in kJ mol1

package,51 in which the correction to the energy can be finally written as:

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

EU Compound

DHelec crystal

DHelec gas phase

Diff.

1 2 3 4 5 6 7

15.0 13.2 14.5 16.8 21.4 20.9 19.1

6.9 5.9 7.6 23.3 22.9 25.9 15.9

8.1 7.3 6.9 6.5 1.5 5.0 3.2

value at the crystalline phase. Therefore, at this point one may wonder whether the development of more accurate computational frameworks, limited to describe isolated molecules, is relevant in the study of properties that are severely influenced by crystal packing effects. It is important to note here that intermolecular interactions have an impact whose sign cannot be generalized. This puts into question the strategy employed in previous computational work consisting of adding an empirical contribution of ca. 50 kJ mol1 to the value obtained at the gas phase, in order to account for the intermolecular interactions.46,47 The adoption of this large value is based on previous work of Kepenekian et al., who calculated the contribution of the Madelung field in model systems.15

4. Computational details 4.1.

Introduction to the DFT + U methodology

In LDA and GGA DFT-functionals there is an incomplete cancellation of the electronic self-interaction.48 That is, the interaction of an electron with the exchange and correlation potential generated by its own charge, which results in an unrealistic delocalization of orbitals. This fact is manifested by the partial occupation of nearly-degenerated orbitals, such as those of the ‘d’ shell of Fe. In order to cure this defect, the bare DFT Hamiltonian is complemented by a Hubbard-like U-term (EU, eqn (7)) that introduces a strong intra-atomic interaction that penalizes partial occupations. EDFT+U[r(r)] = EDFT[r(r)] + EU[{nIs mm}]

(7)

nIs mm

are generalized atomic where r(r) is the electronic density, orbital occupations with spin s associated with the I atom, and EDFT[r(r)] is the standard LDA or GGA functional. Since EDFT[r(r)] already contains an approximate correlation contribution, a term intended to model such a contribution, EDC[{nIs}], must be subtracted to avoid double counting, DC (eqn (8)).   EDFTþU ½rðrÞ ¼ EDFT ½rðrÞ þ EHUB nIs mm0  Is   EDC n (8) sum of the occupations corresponding to all where nIs  is the    P Is orbitals represents the ‘‘correct’’ nmm0 and EHUB nIs mm0 m

on-site correlation energy. Along this manuscript, this methodology has been used following the formulation of Liechtenstein et al.49 and Dudarev et al.50 as implemented in the Quantum Espresso

Phys. Chem. Chem. Phys.

 Is  U X X Is   nmm0 ¼ li 1  lIs i 2 I;s i

(9)

Eqn (9) already includes the double counting correction, and imposes a penalty (mediated by U) for fractional occupations, thus favoring the stability of either fully occupied (l = 1) or empty (l = 0) orbitals. One must note that, under the definition of Liechtenstein and Dudarev, the term U corresponds to the difference U  J (where J is a Stoner-like exchange parameter), alternatively utilized in other implementations. When referring to other computational papers in which this definition has been employed, we have defined an effective U value (Ueff), which corresponds to U  J. For the results obtained in present paper, J has been set to 0 and, thus, Ueff = U. Finally, it must be mentioned again that we have only applied the U correction to the ‘d’ orbitals of Fe. 4.2 Adjusting the vibrational normal modes to reproduce DSvib(T1/2). Given that the computation of analytical frequencies at the crystalline phase would be prohibitive (even using the finitedifference technique would involve thousands of single-point evaluations), we computed the minimum energy structures of HS and LS at the molecular level with Gaussian09,52 for which we obtained the frequencies of the vibrational normal modes (ni). To this purpose, we have systematically used four different computational schemes: PBE/def2TZVP, PBE/6-31+g(d), B3LYP/ TZVP and BP86/TZVP, successfully employed in previous computational work to compute ni values, and the resulting frequencies have been used to evaluate the vibrational entropy at the transition temperature using eqn (5) (note that for those calculations, U = 0). All methods provided extremely unsatisfactory results, often overestimating DSvib(T1/2) by 100%. The reason for such a disagreement probably relies on the ni values of the soft modes, since those are the most affected by the approximations involved (isolated molecule, harmonic approximation), and those small errors may induce large changes in the computed DSvib values. In fact, the value of those frequencies (few cm1) lies at the edge of the accuracy of any computational methodology, and changes significantly with the functional and the basis set used. For all those reasons, one believes that there is no special reason to stick to the computed ni values for the soft vibrational modes (i.e. o100 cm1). Therefore, and given the inability to obtain reliable frequencies, we decided to iteratively modify the ni values obtained with PBE/def2TZVP for the first eight vibrational modes, until the whole set of adjusted frequencies provided the experimental DSvib(T1/2) (see Table 5). As mentioned before, those softer modes are the ones that suffer from the largest computational error and, at the same time, they concentrate the largest impact on DSvib.53,54 This means that adjusting them seems to be a good strategy and, by doing so, one can rapidly achieve the desired DSvib(T1/2) value while keeping the other key quantity,DHvib(T1/2), almost unchanged (see below).

This journal is © the Owner Societies 2015

View Article Online

PCCP

Paper

Table 5 Computed and adjusted ni (in cm1) for the HS and LS minima of compound 1. The computed values have been obtained using PBE/def2TZVP

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

HS

LS

i

ni (computed)

ni (adjusted)

ni (computed)

ni (adjusted)

1 2 3 4 5 6 7 8

12.24 14.57 15.96 17.37 21.90 28.24 29.65 62.88

13.49 16.06 17.59 19.14 24.13 31.12 32.67 69.29

14.57 16.13 17.97 18.99 29.98 37.91 41.50 77.65

13.08 14.48 16.14 17.05 26.92 34.04 37.27 69.73

The adjustment of the original ni values has been done by multiplying those of the HS minima by a factor of (1 + x), and by (1  x) in the case of the LS minima, where x may be positive or negative depending on whether the target DSvib(T1/2) value is larger or smaller than the one obtained with the original frequencies. For the studied compounds 1–7, x has been found to be between 0.1 and 0.4, which typically implies a modification of 1–10 cm1 (see Table 5). As discussed in the main text, we have subsequently used the adjusted vibrational frequencies to compute DHvib(T1/2) and used this value to extract the so-called DHref elec, our reference value for the posterior parametrization of the U-term within DFT + U. We have confirmed that the use of the original set of frequencies or the modified ones does not imply a significant change in the resulting DHvib(T1/2) and, thus, in the benchmarked U values. For instance, in the case of compound 1 it changes less than 0.15%. This difference slightly increases for other compounds for which T1/2 is located at higher temperatures. 4.3.

DFT + U calculations

All the DFT + U calculations have been performed using the Quantum Espresso package (QE),51 using the PBE functional within the spin unrestricted formalism, the D2 correction of Grimme, Vanderbilt pseudopotentials, and a G-point sampling of the Brillouin zone. See below for further details. 4.3.1 DFT + U calculations in the solid state. The minimum energy structure of the HS and LS crystals for all compounds 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  105 atomic units. In these calculations, the number of plane waves has been kept constant at a kinetic energy cutoff of 70 Ry throughout the variable-cell relaxations. A constant number of plane waves imply no Pulay stress but a decreasing precision of the calculation as the volume of the super-cell increases. The large cutoff employed in these calculations ensures that the artifacts arising from this change of precision are negligible. The spin state of the iron atoms is obtained by defining an appropriate initial guess (LS or HS) that is 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 crystal structures have been obtained computationally from

This journal is © the Owner Societies 2015

the HS crystals by defining the spin state of all the SCO units as a singlet, and allowing for the relaxation of the unit cell. 4.3.2 DFT + 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. 25 Bohr3, have been excised from the optimized unit cell and introduced on a cubic cell of 60 Bohr3. This fact, together with the application of the Makov– Payne correction within QE, 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 performed on the corresponding crystal and, thus, only single point evaluations have been done. If non-identical SCO units have been found on the optimized crystal, this process has been done for each of those, and the results shown in the main text correspond to an average, with very small deviations observed between the individual values.

5. Conclusions and outlook In the present paper we have fine-tuned the U-parameter to correctly describe the HS–LS energy difference of seven SCO compounds based on an Fe(II) ion coordinated to six N atoms in a pseudo-octahedral fashion. We have derived a way to extract an estimation of the DHref elec value from accessible experimental data, which has been then taken as a reference in the subsequent parametrization of the Hubbard term. According to our results, the proper U value ranges from 2.37 eV to 2.97 eV, which evinces that this parameter is system-dependent, even for systems that share a common architecture, in this case an FeN6 core. However, we have also demonstrated the excellent transferability of this average U value of 2.65 eV, as it is capable of describing the energetics of the whole set of studied compounds with excellent accuracy. Indeed, the MAE of 4.3 kJ mol1 is the smallest reported up to now in the evaluation of DHelec, a key quantity in defining SCO compounds. Therefore, we believe that this value of U = 2.65 eV can be regarded as a good starting point when performing DFT + U calculations on this class of materials. Furthermore, we have evaluated the importance of the intermolecular interactions in the relative stability of the HS and LS states of the studied SCO compounds. Herein, we demonstrate that those interactions qualify for a contribution that ranges from 8.1 kJ mol1 in 1 to 6.5 kJ mol1 in 4, and influences the value of DHelec by an average of ca. 5 kJ mol1. It is important to note here that intermolecular interactions have an impact whose sign cannot be generalized, and that for compounds 1–3 contributes up to 50% to the DHelec value at the crystalline phase. Therefore, at this point one may wonder whether the development of more accurate computational frameworks, limited to describe isolated molecules, is relevant in the study of properties that are severely influenced by crystal packing effects. The findings reported in the present manuscript pave the way for future studies devoted to understand the crystalline

Phys. Chem. Chem. Phys.

View Article Online

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

Paper

phase of SCO compounds, where hundreds of atoms are present, as well as the adsorption of individual molecules on organic or metallic surfaces. Both scenarios require the use of large unitcells and periodic boundary conditions, which limits the available computational tools to DFT-based approaches. In this sense, the rational incorporation of the U-term within DFT + U may provide the energetic accuracy that is dramatically missing when using bare-DFT functionals.

Abbreviations phen 1,10-Phenanthroline abpt 4-Amino-3,5-bis(pyridin-2-yl)-1,2,4-triazole babpby N-(6-(6-(Pyridin-2-ylamino)pyridin-2-yl)pyridin-2-yl)pyridin-2-amine pz Pyrazine bipy 2,2 0 -Bipyridine 1-bpp 2,6-Di(pyrazol-1-yl)pyridine

PCCP

12 13

14 15 16

17 18 19

Acknowledgements

20

S.V., M.F. and J.R.-A. acknowledge the Spanish Government for ´n financial support (Project MAT2011-25972) and a ‘‘Ramo y Cajal’’ fellowship to J.R-A. We also acknowledge Universitat de Barcelona for a PhD grant to M.F., the BSC and CSUC for the allocation of massive computer time and the Catalan DURSI (Grant 2014SGR1422). S.V. is also thankful to the LabEx-Chemistry of Complex Systems for a post-doctoral grant. (ANR-10-LABX-0026CSC).

21 22 23 24 25

Notes and references 1 M. A. Halcrow, Spin-Crossover Materials: Properties and Applications, Wiley, 2013. ´tard, P. Guionneau and ¨tlich, H. A. Goodwin, J.-F. Le 2 P. Gu L. Goux-Capes, Spin Crossover in Transition Metal Compounds Vol. I-III, Springer, Berlin/Heidelberg, 2004, vol. 235, pp. 1–19. 3 A. Bousseksou, G. Molnar, L. Salmon and W. Nicolazzi, Chem. Soc. Rev., 2011, 40, 3313–3335. ¨tlich, A. Hauser and H. Spiering, Angew. Chem., Int. Ed. 4 P. Gu Engl., 1994, 33, 2024–2054. 5 O. Kahn and C. J. Martinez, Science, 1998, 279, 44–48. 6 T. Mahfoud, G. Molnar, S. Cobo, L. Salmon, C. Thibault, C. Vieu, P. Demont and A. Bousseksou, Appl. Phys. Lett., 2011, 99, 053307. 7 T. Miyamachi, M. Gruber, V. Davesne, M. Bowen, S. Boukari, L. Joly, F. Scheurer, G. Rogez, T. K. Yamada, P. Ohresser, E. Beaurepaire and W. Wulfhekel, Nat. Commun., 2012, 3, 938. 8 P. Gutlich, Y. Garcia and H. A. Goodwin, Chem. Soc. Rev., 2000, 29, 419–427. 9 J. A. Real, A. B. Gaspar and M. C. Munoz, Dalton Trans., 2005, 2062–2079, DOI: 10.1039/b501491c. 10 G. A. Craig, O. Roubeau and G. Aromı´, Coord. Chem. Rev., 2014, 269, 13–31. 11 P. Koningsbruggen, Y. Maeda and H. Oshio, in Spin Cross¨tlich and over in Transition Metal Compounds I, ed. P. Gu

Phys. Chem. Chem. Phys.

26 27 28 29 30 31 32 33

34 35 36 37

H. A. Goodwin, Springer, Berlin Heidelberg, 2004, ch. 10, vol. 233, pp. 259–324. S. Hayami, Y. Komatsu, T. Shimizu, H. Kamihata and Y. H. Lee, Coord. Chem. Rev., 2011, 255, 1981–1990. ´tard, P. Guionneau and L. Goux-Capes, Spin Crossover J.-F. Le in Transition Metal Compounds III, Springer, Berlin Heidelberg, 2004, ch. 10, vol. 235, pp. 221–249. B. Le Guennic, S. Borshch and V. Robert, Inorg. Chem., 2007, 46, 11106–11111. M. Kepenekian, B. Le Guennic and V. Robert, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 094428. ´fuel, S. Shova, J. A. Wolny, F. Dahan, G. Lemercier, N. Bre M. Verelst, H. Paulsen, A. X. Trautwein and J.-P. Tuchagues, Chem. – Eur. J., 2006, 12, 7421–7432. ¨nemann and J. A. Wolny, Eur. J. Inorg. H. Paulsen, V. Schu Chem., 2013, 628–641. L. M. Lawson Daku, A. Vargas, A. Hauser, A. Fouqueau and M. E. Casida, ChemPhysChem, 2005, 6, 1393–1410. L. M. Lawson Daku, F. Aquilante, T. W. Robinson and A. Hauser, J. Chem. Theory Comput., 2012, 8, 4216–4231. M. Reiher, O. Salomon and B. Artur Hess, Theor. Chem. Acc., 2001, 107, 48–55. M. Reiher, Inorg. Chem., 2002, 41, 6928–6935. A. Slimani, X. Yu, A. Muraoka, K. Boukheddaden and K. Yamashita, J. Phys. Chem. A, 2014, 118, 9005–9012. T. Bucko, J. Hafner, S. Lebegue and J. G. Angyan, Phys. Chem. Chem. Phys., 2012, 14, 5389–5396. ´ngya `gue, S. Pillet and J. G. A ´n, Phys. Rev. B: Condens. S. Lebe Matter Mater. Phys., 2008, 78, 024433. S. Vela, J. J. Novoa and J. Ribas-Arino, Phys. Chem. Chem. Phys., 2014, 16, 27012–27024. S. Gueddida and M. Alouani, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 144413. M. Kepenekian, V. Robert, B. Le Guennic and C. De Graaf, J. Comput. Chem., 2009, 30, 2327–2333. A. Rudavskyi, C. Sousa, C. de Graaf, R. W. A. Havenith and R. Broer, J. Chem. Phys., 2014, 140, 184318. A. Rudavskyi, R. W. A. Havenith, R. Broer, C. de Graaf and C. Sousa, Dalton Trans., 2013, 42, 14702–14709. S. Saureu and C. de Graaf, Chem. Phys., 2014, 428, 59–66. C. de Graaf and C. Sousa, Chem. – Eur. J., 2010, 16, 4550–4556. ´n, C. de Graaf and C. Sousa, J. Am. Chem. Soc., B. n. Ordejo 2008, 130, 13961–13968. An error has been found in the contraction of the ANO-RCC basis set for C atoms. Although it is now corrected, this error has been present in Molcas versions from 6.4 to 8.0 that span from 2006 to 2014. C. Angeli, R. Cimiraglia and J.-P. Malrieu, J. Chem. Phys., 2002, 117, 9138–9153. C. Angeli, R. Cimiraglia and J.-P. Malrieu, Chem. Phys. Lett., 2001, 350, 297–305. B. Gallois, J. A. Real, C. Hauw and J. Zarembowitch, Inorg. Chem., 1990, 29, 1152–1158. S. Bonnet, M. A. Siegler, J. S. Costa, G. Molnar, A. Bousseksou, A. L. Spek, P. Gamez and J. Reedijk, Chem. Commun., 2008, 5619–5621, DOI: 10.1039/B811746B.

This journal is © the Owner Societies 2015

View Article Online

Published on 21 May 2015. Downloaded by UNIVERSITY OF NEW ORLEANS on 05/06/2015 08:25:32.

PCCP

38 L. Salmon, G. Molnar, S. Cobo, P. Oulie, M. Etienne, T. Mahfoud, P. Demont, A. Eguchi, H. Watanabe, K. Tanaka and A. Bousseksou, New J. Chem., 2009, 33, 1283–1289. ˜ oz, J. Faus and X. Solans, Inorg. Chem., 39 J. A. Real, M. C. Mun 1997, 36, 3008–3013. 40 J. M. Holland, J. A. McAllister, Z. Lu, C. A. Kilner, M. Thornton-Pett and M. A. Halcrow, Chem. Commun., 2001, 577–578, DOI: 10.1039/b100995h. 41 K. P. Jensen and J. Cirera, J. Phys. Chem. A, 2009, 113, 10033–10039. 42 G. Mercurio, E. R. McNellis, I. Martin, S. Hagen, F. Leyssner, S. Soubatch, J. Meyer, M. Wolf, P. Tegeder, F. S. Tautz and K. Reuter, Phys. Rev. Lett., 2010, 104, 036102. ´ ski and N. Lo ´pez, 43 N. Almora-Barrios, G. Carchini, P. Błon J. Chem. Theory Comput., 2014, 10, 5002–5009. 44 T. Bucko, J. Hafner, S. Lebegue and J. G. Angyan, J. Phys. Chem. A, 2010, 114, 11814–11824. 45 S. Vela, F. Mota, M. Deumal, R. Suizu, Y. Shuku, A. Mizuno, K. Awaga, M. Shiga, J. J. Novoa and J. Ribas-Arino, Nat. Commun., 2014, 5, 4411. ` and S. Sanvito, Phys. Rev. B: Condens. 46 A. Droghetti, D. Alfe Matter Mater. Phys., 2013, 87, 205114. `guerie and 47 N. Suaud, M.-L. Bonnet, C. Boilleau, P. Labe ´ry, J. Am. Chem. Soc., 2009, 131, 715–722. N. Guihe 48 J. P. Perdew and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079. 49 A. I. Liechtenstein, V. I. Anisimov and J. Zaanen, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 52, R5467–R5470. 50 S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505.

This journal is © the Owner Societies 2015

Paper

51 P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502. 52 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. ¨ . Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and Daniels, O D. J. FoxI. Gaussian, Wallingford CT, 2009. 53 A. Bousseksou, J. J. McGarvey, F. Varret, J. A. Real, J.-P. Tuchagues, A. C. Dennis and M. L. Boillot, Chem. Phys. Lett., 2000, 318, 409–416. 54 K. L. Ronayne, H. Paulsen, A. Hofer, A. C. Dennis, J. A. Wolny, A. I. Chumakov, V. Schunemann, H. Winkler, H. Spiering, A. Bousseksou, P. Gutlich, A. X. Trautwein and J. J. McGarvey, Phys. Chem. Chem. Phys., 2006, 8, 4685–4693.

Phys. Chem. Chem. Phys.