# The ground state of a spin-crossover molecule calculated by diffusion Monte Carlo

## Abstract

Spin crossover molecules have recently emerged as a family of compounds potentially useful for implementing molecular spintronics devices. The calculations of the electronic properties of such molecules is a formidable theoretical challenge as one has to describe the spin ground state of a transition metal as the legand field changes. The problem is dominated by the interplay between strong electron correlation at the transition metal site and charge delocalization over the ligands, and thus it fits into a class of problems where density functional theory may be inadequate. Furthermore, the crossover activity is extremely sensitive to environmental conditions, which are difficult to fully characterize. Here we discuss the phase transition of a prototypical spin crossover molecule as obtained with diffusion Monte Carlo simulations. We demonstrate that the ground state changes depending on whether the molecule is in the gas or in the solid phase. As our calculation provides a solid benchmark for the theory we then assess the performances of density functional theory. We find that the low spin state is always over-stabilized, not only by the (semi-)local functionals, but even by the most commonly used hybrids (such as B3LYP and PBE0). We then propose that reliable results can be obtained by using hybrid functionals containing about of exact-exchange.

In nature there is a vast class of molecules whose magnetic moment can be altered by an external stimulus. Typical examples of such molecules are the spin-crossover (SC) complexes (1); (2), which, in their most abundant form, contain a Fe ion in octahedral coordination (3) and exhibit a transition from the low spin (LS) (singlet) ground state to a high spin (HS) (quintet) metastable state. Other examples are the cobalt dioxolene molecules (4); (5); (6). These undergo the so-called valence tautomeric interconversion (VTI), namely an interconversion between two redox isomers, which differ in charge distribution and spin configuration. Both the SC transition and the VTI are usually observed for molecules in a single crystal and can be triggered by variations in temperature and pressure or by optical irradiation (7). Furthermore it was also recently suggested that the VTI (8) and the spin ground state of a two-center polar molecule (9) can be controlled by a static electric field .

SC complexes are promising materials candidates for molecular spintronics applications (10); (11). Devices incorporating such molecules are predicted to display drastic changes in the current-voltage curve across the phase transition (12); (13), and several transport experiments have recently achieved encouraging results. Alam et al. (14) were able to distinguish the spin state of a SC molecule placed on graphite by scanning tunnel microscopy, while Prins et al. (15) demonstrated that the temperature dependent conductance of a device incorporating a SC cluster correlates well with the phase transition. In other cases, however, the data are not easy to interpret (16) and the experimental investigations are combined with density functional theory (DFT) simulations. In principle DFT should allow the computation of quantities not easily accessible by experiments and should also provide parameters for effective transport models. However, unfortunately DFT results for SC molecules depend strongly (even qualitatively) on the choice of the exchange-correlation functional used (17); (22); (18); (19); (20); (21) and no standard has yet emerged. This essentially means that DFT is still not a predictive theory for this problem.

Since most of the local and semilocal functionals underestimate the exchange energy, they tend to favor the LS state against the HS one (17); (22). This shortcome often leads to such large errors that even stable HS molecules are described as LS (18); (19). In contrast, the most commonly used hybrid functionals are believed to over-stabilize the HS state (23); (24). Besides, several authors (22); (18); (19); (25) have criticized the common practice, that consists in assessing the performances of the various functionals by direct comparison with the experimental data. In fact, while experiments are usually performed for molecules in the condensed phase, DFT results refer to molecules in the gas phase. Since the properties of SC complexes depend strongly on environmental conditions (counter-ions, interaction between molecules, “strain effects” etc…) (26); (27) their ground state may not be the same in these two phases. The question then becomes: can we produce a robust benchmark for DFT against the problem of predicting the Physics of SC compounds? In order to answer to this question ab-initio methods more accurate than DFT have to be considered.

In the past wave-function based methods were used for this problem (18); (19); (20); (21); (16). However, as the authors themselves pointed out, the results were plagued by systematic errors ascribed to the limited basis set used for Fe and by the fact that the methods themselves neglect dynamic correlation (although this can be partly accounted for through a perturbative treatment). Here we propose an alternative route and perform diffusion Monte Carlo (DMC) (28) calculations for a prototypical Fe spin crossover molecule. As DMC represents one of the most accurate electronic structure method currently available in order to compute ground state energies, our calculations provided a solid benchmark for assessing the performances of DFT.

In particular we consider the molecule [FeL](BF) (L=2,6-dypirazol-1-yl-4-hydroxymethylpyridine) (29) (see Fig. 1). We show that its ground state in the gas phase is HS but that a phase transition may exist in the solid state due to a number of crystal-related effects. We then show that the same result can be obtained by DFT hybrid functionals containing approximately of exact exchange, thus confirming early calculations for model molecules (25). This establishes a recipe for the use of DFT for this class of materials, and it opens the opportunity to investigate with confidence the spin crossover transition of molecules in different environments (for instance on surfaces).

DFT calculations are performed with the nwchem code (30). We use several functionals belonging to different classes: 1) the Vosko-Wilk-Nussair local density approximation (LDA) (32); 2) the generalized gradient approximation BP86, which combines the Becke88 exchange functional (33) with the Predew86 correlation one (34); 3) the hybrid functionals B3LYP (35), PBE0 (36); (37) and the Becke-HH (38), which include respectively , and of exact exchange. We also consider a re-parametrization of the B3LYP functional, called B3LYP, which includes only of HF exchange. This was introduced by Reiher and co-workes specifically in order to describe Fe complexes (23); (24). The Ahlrichs triple-zeta polarized basis set (31) is used throughout. DMC calculations are performed by using the casino code (39). The imaginary time evolution of the Schrödinger equation has been performed with the usual short time approximation and time-steps of and a.u. are used. Dirac-Fock pseudopotentials (40); (41) with the “potential localization approximation” (42) have been used. The single-particle orbitals of the trial wave function are obtained through (LDA) DFT calculations performed with the plane-wave (PW) code quantum espresso (43). The same pseudopotentials used for the DMC calculations are employed. The PW cutoff is fixed at Ry and the PW are re-expanded in terms of B-splines (44). The B-spline grid spacing is , where is the length of the largest vector employed in the PW calculations. Periodic boundary conditions are employed for the PW-DFT calculations and supercells as large as Å are considered. In contrast, no periodic boundary conditions are imposed with DMC. Counter ions have been ignored and calculations are presented for the cation (Fig. 1) to which we will generally refer as the molecule.

The crucial quantity for understanding the spin crossover transition is the potential energy surface, schematically displayed in Fig. 2. This is typically plotted for the two different spin configurations as a function of a collective reaction coordinate , which interpolates the molecule geometry along the LS to HS phase transition. In our case DMC and DFT are used to compute the “adiabatic energy gap” (22) defined as

(1) |

where () and [] represent respectively the geometry and the total energy of the LS-singlet (HS-quintet) state. When studying SC molecules at zero temperature is the central quantity, as it indicates whether the molecule ground state is LS () or HS (). We also calculate the “vertical energy gaps”(22)

(2) | |||

(3) |

where [] is the energy of the quintet (singlet) state for the () geometry (see Fig. 2).

Spin state | (Å) | |
---|---|---|

BP86 | LS | |

BP86 | HS | |

B3LYP* | LS | , |

B3LYP* | HS | |

B3LYP | LS | |

B3LYP | HS | |

PBE0 | LS | |

PBE0 | HS | |

Exp. | LS | |

Exp. | HS |

DMC calculations were first carried out by using the molecular geometries optimized with DFT for the molecule in the gas phase. The Fe-ligand bond lengths computed with the various functionals are listed in Tab. 1. As the DMC energy differences between the geometries calculated from BP86, B3LYP*, B3LYP and PBE0 are of the same order of magnitude as the Monte Carlo statistical error, we have not been able to firmly establish which functional produces the best structure. We have then decided to present results only for the structures relaxed with B3LYP, keeping in mind that the same are essentially valid also for BP86 and PBE0.

Method | (eV) | (eV) | (eV) |
---|---|---|---|

BP86 | |||

B3LYP* | |||

B3LYP | |||

PBE0 | |||

Becke-HH | |||

DMC ( a.u) |

The DMC adiabatic energy gap is reported in Tab. 2. Our result indicates that the molecule in the gas phase is in its high-spin state, in contrast to the common belief and to the experimental result for the single crystal. Such ground state is indeed quite robust as DMC gives us an adiabatic energy gap of -1.20 eV. Since DMC provides a unequivocal assignment of the molecule ground state, it essentially establishes that no spin crossover transition is expected for [FeL] in the gas phase. Hence, in order to account for the experimentally observed SC transition, one needs to understand how the embedding of the molecule in a crystal is able to reverse the relative order of the HS and the LS states at zero-temperature, i.e. to change the sign of .

(eV) | (eV) | (eV) | |
---|---|---|---|

DMC ( a.u) | |||

DMC ( a.u) |

The argument proceeds then as follows. Firstly, one has to repeat the calculations using the experimental geometries measured for the molecule in the crystal form (29). These are less symmetric and present shorter metal-ligand bond-lengths than those optimized in the vacuum (see Tab. 1). The result is that the DMC-calculated gets smaller, although it maintains the negative sign (compare Tab. 3 with Tab. 2). Secondly, the electrostatic potential felt by the molecule in the crystal and due to both the counter-ions and the other molecules needs to be taken into account. This produces a relative shift of the HS and LS potential energy surfaces. The magnitude of this effect, which tends to stabilize the LS state, has been recently estimated (27) to be of the order of eV. When the effect of the geometry and the electrostatic corrections are both included, DMC allows us to estimate a for the condensed phase of about 0.2 eV. This is now positive, i.e. the ground state is LS, and very close to the typical values of the adiabatic energy gap inferred from experimental data (45).

We finally turn our attention to the assessment of the performaces of the various exchange-correlation functionals. Table 2 displays the vertical and adiabatic energy gaps calculated with DFT. We note that BP86 underestimates the exchange so significantly that the molecule is predicted to be stable in the LS state (). Furthermore the absolute value of [] is much larger than the corresponding one computed with DMC. This means that the standard local density approximation predicts a very stable low spin ground state. B3LYP and PBE0 improve only slightly the accuracy of the calculated gaps and the LS state still remains massively over-stabilized. In contrast, as in the case of small Fe model complexes (25), HH is found to be the functional that performs better, yielding a fair agreement with the DMC gaps.

Importantly our analysis demonstrates that the assessment of the performances of a given DFT functional can be completely erroneous, if one insists on comparing the total energies calculated for the gas phase directly to experiments. If this was done with our DFT data, we would have concluded, as other authors did (23); (24), that B3LYP* was the best functional. Our analysis instead demonstrates that the correct assignment needs to be done against a reliable benchmark, with the result that the best suited functional must carry a fraction of exact exchange near to 50%.

In conclusion, we have shown that for a typical SC molecule in the gas-phase the ground state is high spin and discussed how the this becomes singlet in the condensed phase. We have also assessed the performaces of DFT against this problem, demonstrating that the HH hybrid functional, including of exact exchange, is able to provide a quite accurate estimate of the energetic of the molecule. Our results then finally shed light on the long-standing issue of establishing the ground state of SC complexes.

acknowledgments: We thank N. Baadji for useful discussions. This work is sponsored by EU (Hints project). Computational resources have been provided by TCHPC and ICHEC.

### References

- O. Kahn, Molecular magnetism (VCH, New York, 1993).
- P. Gütlich, H.A. Goodwin, in Spin Crossover in Transition Metal Compounds (eds. P. Gütlich, H.A. Goodwin) (Springer-Verlag, Berlin-Heidelberg, 2004).
- SC molecules incorporating Fe, Co and Mn have also been reported. However these are much less common then those containing Fe ions.
- O. Sato, J. Tao, Y.Z. Zhang, Angew. Chem. Int. Ed. 46, 2152 (2007).
- D.N. Hendrickson and C.G. Pierpont, Top. Curr. Chem. 234, 63 (2004).
- A. Beni, C. Carbonera, A. Dei, J.-F. Létard, R. Righini, C. Sangregorio, L. Sorace, J. Braz. Chem. Soc. 17, 1522 (2006).
- A. Hauser, J. Chem. Phys. 94, 2741 (1991).
- A. Droghetti, S. Sanvito, Phys. Rev. Lett. 107, 047201 (2011)
- N. Baadji, M. Piacenza, T. Tugsuz, F. Della Sala, G. Maruccio and S. Sanvito, Nature Materials 8, 813 (2009).
- L. Bogani, W. Wernsdorfer, Nature Materials 7,179 (2008).
- S. Sanvito, Chem. Soc. Rev. 40, 3336 (2011).
- N. Baadji and S. Sanvito, Phys. Rev. Lett., in press, also arXiv:1201.2028 (2012).
- D. Aravena, E. Ruiz, J. Am. Chem. Soc. 134, 777 (2011).
- M.S. Alam, M. Stocker, K. Gieb, P. Müller, M. Haryono, K. Student and A. Grohmann, Angew. Chem. Int. Ed. 49, 1159 (2010)
- F. Prins, M. Monrabal-Capilla, E.A. Osorio, E. Coronado and H.S.J. van der Zant, Adv. Mater. 23,1545 (2011).
- V. Meded, A. Bagrets, K. Fink, R. Chandrasekar, M. Ruben, F. Evers, A. Bernand-Mantel, J.S. Seldenthuis, A. Beukman and H.S.J. van der Zant, Phys. Rev. B, 83, 245415 (2011).
- M. Swart, A. Groenhof, A.W. Ehlers and K. Lammertsma, J. Phys. Chem. A 108, 5479 (2004).
- A. Fouqueau, S. Mer, M. Casida, L.M.L. Daku, A. Hauser, T. Mineva, F. Neese, J. Chem. Phys. 120, 9473 (2004).
- A. Fouqueau, M. Casida, L.M.L. Daku, A. Hauser and F. Neese, J. Chem. Phys. 122, 044110 (2005).
- K. Pierloot and S. Vancoilie, J. Chem. Phys. 125, 124303 (2006).
- K. Pierloot and S. Vancoilie, J. Chem. Phys. 128, 034104 (2006).
- S. Zein, S.A. Borshch, P. Fleurat-Lessard, M.E. Casida and H. Chermette, J. Chem. Phys. 126, 014105 (2007).
- M. Reiher, Inorg. Chem. 41, 6928 (2002).
- M. Reiher, O. Salomon and B.A. Hess, Theor. Chem. Acc. 107, 48 (2001).
- A. Droghetti, D. Alfè and S. Sanvito, in preparation.
- M. Kepenekian, B. Le Guennic and V. Robert, J. Am. Chem. Soc. 131, 11498 (2009).
- M. Kepenekian, B. Le Guennic and V. Robert, Phys. Rev. B 79, 094428 (2009).
- W.M.C. Foulkes, L. Mitas, R.J. Needs and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
- V.A. Money, J. Elhaïk, M.A. Halcrow and J.A.K. Howard, Dalton Trans. 10, 1516 (2004).
- M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J. van Dam, D. Wang, J. Nieplocha, E. Apra, T.L. Windus and W.A. de Jong, Comput. Phys. Commun. 181, 1477 (2010).
- A. Schäfer, C. Huber and R. Ahlrichs, J. Chem. Phys. 100, 5829 (1994).
- S.J. Vosko, L. Wilk, M. Nusair and Can. J. Phys. 58, 1200s (1980).
- A.D.Becke, Phys. Rev. A 38, 3098 (1988).
- J.P. Perdew, Phys. Rev B 33, 8822 (1986).
- P.J. Stephens, F.J. Devlin, C.F. Chabalowski and M.J. Frisch, J. Phys. Chem. 98, 11623 (1994).
- M. Ernzerhof and G.E. Scuseria, J. Chem. Phys. 110, 5029 (1999).
- C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
- A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
- R.J. Needs, M.D. Towler, N.D. Drummond and P. López Ríos., J. Phys.: Condens. Matter 22, 023201 (2010).
- J.R. Trail and R.J. Needs, J. Chem. Phys. 122, 174109 (2005).
- J.R. Trail and R.J. Needs, J. Chem. Phys. 122, 014112 (2005).
- L. Mitas, E.L. Shirley and D.M. Ceperley, J. Chem. Phys. 95, 3467 (1991).
- P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de 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 21, 395502 (2009).
- D. Alfè and M.J. Gillian, Phys. Rev. B 70, 161101(R) (2004).
- G. Ganzenmüller, N. Berkaïnem A. Fouqueau, M.E. Casida and M. Reiher, J. Chem. Phys. 122, 234321 (2005).