Diffusion Monte Carlo Perspective on the Spin-State Energetics of [Fe

Article pubs.acs.org/JCTC

Diffusion Monte Carlo Perspective on the Spin-State Energetics of [Fe(NCH)6]2+ Maria Fumanal,*,†,‡ Lucas K. Wagner,*,⊥ Stefano Sanvito,*,§ and Andrea Droghetti*,§,∥ †

Institut de Chimie UMR7177 CNRS-Université de Strasbourg, 1 Rue Blaise Pascal BP 296/R8, F-67007 Strasbourg, France Departament de Química Física and IQTCUB, Facultat de Química, Universitat de Barcelona, Av. Diagonal 645, 08028 Barcelona, Spain ⊥ Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, United States § School of Physics, AMBER and CRANN Institute, Trinity College, Dublin 2, Ireland ∥ Nano-Bio Spectroscopy Group and European Theoretical Spectroscopy Facility (ETSF), Universidad del Pais Vasco CFM CSIC-UPV/EHU-MPC and DIPC, Av.Tolosa 72, 20018 San Sebastian, Spain ‡

S Supporting Information *

ABSTRACT: The energy difference between the high spin and the low spin state of the model compound [Fe(NCH)6]2+ is investigated by means of Diffusion Monte Carlo (DMC), where special attention is dedicated to analyzing the effect of the fix node approximation on the accuracy of the results. For this purpose, we compare several Slater−Jastrow and multireference Slater−Jastrow trial wave functions. We found that a Slater−Jastrow trial wave function constructed with the generalized Kohn−Sham orbitals from hybrid DFT represents the optimal choice. This is understood by observing that hybrid functionals account for the subtle balance between exchange and correlation effects and the respective orbitals accurately describe the ligand−metal hybridization as well as the charge reorganization accompanying the spin transition. Finally the DMC results are compared with those obtained by Hartree−Fock, DFT, CASSCF, and CASPT2. While there is no clear reference value for the high spin−low spin energy difference, DMC and high level CCSD(T) calculations agree within around 0.3 eV.

1. INTRODUCTION Spin crossover (SCO) is a phenomenon that occurs in (pseudo)-octahedral 3dn (n = 4−7) transition metal coordination complexes and involves a switching process from a lowspin (LS) to a high-spin (HS) state upon a change of temperature, pressure, or other stimuli such as light irradiation or application of a magnetic field.1,2 A spin state change also takes place in many chemical reactions relevant for biological processes.3 Besides the interesting fundamental aspects, SCO materials have been investigated for several potential applications in sensor, display, and memory devices1 and, more recently, in molecular electronic and spintronics4,5 as the spin state could be read and controlled electrically.6−13 Following the first report of SCO in a Fe(III) dithiocarbamate complex,14 the same phenomenon has also been demonstrated for a multitude of molecules containing Ni(II), Cr(II), Co(III), Mn(II,) and Mn(III) ions.1 However, the most common examples are the Fe(II)-based d6 complexes, where the LS state (spin S = 0) and the HS state (S = 2), have, respectively, the electronic configuration (t2g)6(eg*)0 and (t2g)4(eg*)2. The t2g are the three nonbonding 3dxy, 3dxz, and 3dyz Fe(II) orbitals, while the eg* are the two antibonding orbitals stemming from the hybridization of the 3dx2−y2 and 3dz2 Fe(II) orbitals with the ligand orbitals of the same symmetry. © 2016 American Chemical Society

The ab initio theoretical description of these Fe(II) molecules across the SCO represents a great challenge.15 The relative stability of the HS and LS states at zero temperature depends on the competition among the electronic exchange, the electronic correlation, and the ligand-field energy splitting Δ separating the t2g and the e*g molecular orbitals.16 A large Δ results in a LS state, which is further stabilized by correlation effects. In contrast, when two electrons occupy the eg* antibonding orbitals there is a considerable increase of the metal−ligand distances (about 0.2 Å) and a consequent drastic reduction of the ligand-field energy splitting Δ. This reduction, together with the exchange, promotes the HS state. At finite temperature, entropic effects need also to be considered. Both the vibrons and the spin contribution to the entropy are larger for the HS than for the LS state. Therefore, a molecule with a LS ground state can undergo the thermally driven spincrossover transition when the LS−HS electronic energy difference is balanced by the LS−HS entropy difference and the Gibbs free energy vanishes. Finally the environment largely influences the stability of one spin state over the other. Received: April 1, 2016 Published: August 8, 2016 4233

DOI: 10.1021/acs.jctc.6b00332 J. Chem. Theory Comput. 2016, 12, 4233−4241

Article

Journal of Chemical Theory and Computation

projects the ground state wave function out from a trial wave function (TWF) through a stochastic sampling procedure. In contrast to standard quantum chemistry techniques, DMC scales as ∼N3 for the energy per electron, and almost perfect parallelization on large supercomputers can be achieved.47 As such, DMC has been applied to several molecules containing transition metal ions48−52 as well as solid-state materials, where correlation effects play a crucial role (Mott or charge transfer insulator52,53 and even superconducting cuprates54,55). Because of the Fermion sign problem, DMC calculations for electronic systems rely on the fixed-node approximation (FNA). This means that the nodal surface is constrained to be the same as that of the TWF, which results in an upper bound on the ground state energy. The computed energy equals the ground state energy only if the TWF nodal surface coincides with the ground state wave function nodal surface. Previous applications of DMC to SCO compounds44,45 employed TWFs having the Slater−Jastrow (SJ) functional form, i.e., a single Slater determinant multiplied by a correlation term, the Jastrow factor, which is symmetric with respect to particles exchange. The single-particle orbitals entering into the Slater determinant were obtained from DFT calculations within the local spin density approximation (LSDA), and the spin was fixed by imposing 2S = N↑ − N↓, where N↑ (N↓) is the number of spin-up (-down) electrons. However, such first works did not look in detail at the effects that the FNA has on the accuracy for the spin state energetics, a task that we take on in the present article. Therefore, our goal is to assess whether and how much the zero temperature electronic energies for the HS and LS states computed by DMC depend on the choice of the TWF. This is a key question to find out possible limitations or strengths of DMC for Fe(II) complexes. One of the main advantages of DMC is its variational nature, which implies that any improvements in the TWF leads to a lower electronic total energy so that the relative accuracy of different calculations can be determined. As a test case, we look at the model molecule [Fe(NCH)6]2+ (Figure 1). This has been largely investigated by using different

The theoretical state of art for calculating the zero temperature HS−LS electronic energy difference of gas-phase molecules can be summarized as follows. On one hand, Hartree−Fock (HF), which includes exact exchange but completely neglects the electron correlation energy, stabilizes the HS state over the LS state. The exchange keeps like spin electrons apart so that their electrostatic repulsion energy is reduced with respect to that of unlike spin electrons. On the other hand, density functional theory (DFT) in principle accounts for both electron correlation and exchange exactly. However, in practice, DFT calculations rely on approximate exchange-correlation density functionals and no agreement about the “functional of choice” has been reached yet; see, for example, refs 17−19 for a detailed discussion. Local and semilocal functionals tend to overstabilize the LS state. In contrast, hybrid functionals obtain the HS−LS energy differences that vary almost monotonically as a function of the fraction of exact exchange.20−22 Increasing the challenge, a full assessment of DFT’s performance is hindered by the lack of gas-phase zero temperature data to which large sets of DFT results could be systematically compared. In recent years there have been attempts to account for correlation effects beyond the local functional description by combining DFT with approaches originally developed in solid state physics to describe model systems.23,24 In the DFT+U method25,26 and in the Dynamical Mean Field Theory27,28 the strong Coulomb interaction of the electrons on the metal ion is treated with an additional Hubbard-like term. The main drawback of these methods is that the final results largely depend on the scheme used to select the metal ion Hilbert subspace and on the value of the parameters U and J, which describe the on-site Coulomb and exchange interaction. Among the post-Hartree−Fock (HF) quantum chemistry methods, the complete active space self-consistent field (CASSCF) and its multireference second-order perturbative extension, CASPT2,29 have been commonly considered.30−38 These methods traditionally required a previous knowledge of the system to construct the active space, and the wrong selection of the active space orbitals could lead to erroneous results. However, enormous progress has been made in this regard over the last several years, and a recent study proposed a systematic protocol to assess the more appropriate configurational space.39 Still, in CASPT2 the results may depend on the chosen zeroth-order Hamiltonian because the dynamic electron correlation is accounted for by perturbation theory,40,41 and the energies are not variational. Despite these limitations, Kepenekian et al.42 showed that CASPT2 can provide reliable results, if the active space is constructed in order to reflect the charge reorganization induced by the spin-transition. This was demonstrated in the case of the model molecule [Fe(NCH)6]2+ in the HS state, for which the authors achieved a good quantitative agreement between the results of CASPT2 and the coupled cluster with singles, doubles, and perturbative triples, CCSD(T). Similar conclusions were also reached by Lawson Daku et al.43 for the same system but only after increasing the ionization potential-electron affinity (IPEA) shift defining the zeroth order Hamiltonian to approximately twice the standard value. There is some concern in using CCSD(T) as a benchmark for the LS state, because it has a multiconfigurational character.41 Finally, some of us recently suggested44,45 that fixed node diffusion Monte Carlo (DMC)46 can represent an alternative approach to DFT and CASPT2. DMC is a method that

Figure 1. Model compound [Fe(NCH)6]2+.

computational methods in recent works.41,43,56 Although the ground state has been calculated to be HS so that it is not a spin-crossover compound, its study allows for a deep understanding of the effect that the Fe(II)−ligand hybridization has on the spin-state energetics. The LS state has a more marked multiconfigurational character than the HS state, and the limitations of SJ TWFs for systems with large static correlation (e.g., stretched diatomic molecules) have been sometimes discussed in the literature;57thus, it is an interesting test case for DMC. We consider SJ trial wave functions, where the orbitals are obtained from DFT with various exchange-correlation functionals; this was proposed by one of the authors (L.K.W.)48 and by Korolenč et al.58 as a computationally efficient way to 4234


Article


In this work, we consider either a single Slater determinant, i.e., N = 1 in eq 2, or a linear combination of configuration state functions (CSFs), which are themselves symmetry adapted linear combinations of determinants. A TWF containing a single Slater determinant is labeled 1SJ, while a TWF containing a linear combination of CSFs is labeled nSJ. The considered CSFs were obtained by multireference-configuration interaction (MR-CI) calculations and were selected so that they had larger weights than a certain cutoff. The space of chosen configurations contains the six σ-bonding ligand−Fe orbitals (i.e., eg, t1u, and a1g) and the five Fe-3d-dominated orbitals (i.e., t2g and eg*) with all single and double excitations to the 20 lowest energy virtual orbitals (see Figure 2). We note

optimize the nodal surface. Then we consider multi-SJ wave functions, where the antisymmetric part is expanded in a series of Slater determinants. These two approaches provide enough variational freedom to investigate the effect of the FNA and to gain novel insight on correlation effects in SCO compounds. Morales et al.59 have recently demonstrated the potential of multi-SJ TWFs for improving the DMC methods toward very high levels of accuracy, at least for compounds with light elements. The FNA is the only fundamental approximation in DMC,46 whose effect can be evaluated by extrapolating the results to the zero time-step limit. However, the use of pseudopotentials is necessary for transition metal compounds with more than a few atoms; this may result in another source of systematic errors that cannot be practically distinguished from the FNA. The generation of optimal pseudopotentials for DMC represents by itself an open problem, in spite of important progress in recent years.60−63 In this regard, a more detailed investigation may be warranted in possible future studies. The analysis of the dependence of the DMC energies on the TWF for the model compound [Fe(NCH)6]2+ leads to important conclusions. The neglect of dynamic correlation in CASSCF wave functions affects the one-particle density in such a way that the energy is high in DMC even with a multiple Slater determinant TWF. Hybrid DFT functionals generate the lowest energy nodal surfaces because they obtain the most accurate distribution of charge.

2. COMPUTATIONAL METHODS In this work, we examine the zero temperature electronic energies (ELS, EHS) for the gas phase molecule, and their difference, which is defined as the adiabatic energy gap:18,44,45 ΔEadia = E HS(rHS) − E LS(rLS)

(1) Figure 2. Molecular orbitals energy diagram showing the multiconfigurational space considered for the construction of the DMC trial wave function. This comprises the eg*, t2g, eg, t1u, and a1g molecular orbitals (inside the red-dashed bottom frame) together with all single and double excitations to the 20 lowest energy virtual orbitals, which are schematically represented inside the upper red-dashed frame.

Here rHS(LS) is a collective coordinate, which indicates that EHS(LS) is computed at the HS(LS) equilibrium geometry. A positive (negative) ΔEadia means that LS state has a lower (higher) energy than the HS state. DMC calculations were carried out with the QWalk code.64 The employed TWFs have a multi-Slater−Jastrow form: N

ΨT = e−J ∑ ck det Mk↑ det Mk↓ k=1

that our approach is essentially the same as that used by Morales et al.59,67 and which produced results of very high accuracy for atoms and diatomic molecules as well as a system characterized by large static correlation such as a triradical complex. The Jastrow function and the CSF coefficients were simultaneously energy-optimized with variational Monte Carlo by using the linear method.66 The single-particle orbitals in the Slater part of the TWFs were obtained either by HF and DFT or by CASSCF calculations. These were performed with the GAMESS code.68 No further optimization of the orbitals in the presence of the Jastrow factor was carried out. In the CASSCF calculations an active space consisting of 10 electrons distributed over 7 orbitals was considered, i.e., CASSCF(10,7). These 7 orbitals are the e*g and t2g, which have an Fe-3d character, and the two eg σ-bonding orbitals with large contribution from the lone pairs of the ligands (see Figure 3). Note that we were not able to include additional 3d orbitals to account for the double-shell effect in the GAMESS code. The description of the double-shell effect in CASSCF would ensure a better treatment of radial correlation69,70 and the optimized orbitals may profit from this; the possible

(2)

where J is the Jastrow function. This consists of two parts: the electron−ion and electron−electron correlation terms that are functions of the distance of electron−ion and electron− electron pairs, respectively (see ref 64 for details). ck is the weight of the kth configuration determinant and the matrix Mσk (σ = ↑, ↓) has elements [Mσk ]ij = ϕσ(k)j(ri), where ϕσ(k)j(ri) is a single particle orbital. The Jastrow satisfies the electron−electron cusp condition and is symmetric with respect to particle exchange. It adds short-range correlation, typically capturing around 85% of the correlation energy. The Jastrow factor improves the efficiency of DMC, but not the FNA, since it does not change the nodes. The linear combination of determinants gives the antisymmetric part and determines the nodal surface of the TWF. DMC calculations with a simple single-Slater−Jastrow TWF with HF orbitals are able to capture about 90−95% of the correlation energy, while an accurate description of systems with near degeneracies, for instance, the C2 molecule,65,66 may require a TWF with many determinants. 4235


Article


The CASSCF/CASPT2 calculations were performed with the MOLCAS 7.6 package.80 Relativistic correlation-consistent natural atomic orbitals ANO-RCC basis sets were utilized to properly describe the correlation of the semicore electrons.81 In all cases, the applied contracted basis sets are (8s,7p,6d,4f,3g,2h) for the Fe atom, (4s,3p,2d,1f) for the N atoms, (3s,2p,1d) for the C atoms, (3s,2p,1d) for the O atoms, (5s,4p,2d,1f) for the S atoms, and (2s) for the H atoms. The Cholesky decomposition82 was used to reduce the computational cost associated with the calculation of the two-electron integrals, and scalar relativistic effects were included by using the standard second-order Douglas−Kross−Hess Hamiltonian.83−85 Unlike the CASSCF calculations performed to obtain the orbitals for the DMC TWFs, we properly treat the double-shell effect here. Therefore, we considered an active space that contains 10 electrons distributed over 12 orbitals (see Figure 3): the five t2g and e*g orbitals, the two eg σ-bonding orbitals, and other additional five Fe-3d orbitals. This CASSCF(10,12) active space has been proven to provide a balanced description of all the important nondynamic correlation effects (see, for example, refs 31, 32, 37, and 42). Then, CASPT2 accounts for the remaining electron correlation by correlating all electrons, expect the deep core ones, within the second order in perturbation theory. The standard zeroorder Hamiltonian40 with IPEA shift of 0.25 au was applied in the definition of the diagonal Fock matrix elements of the active orbitals. Here, we do not investigate the effect of rescaling the IPEA shift as proposed by Lawson Daku et al.;43 however, we will note that our results could be easily extrapolated in order to agree with those in ref 43. In order to systematically compare the electronic energies across the different theories and methods (PBE, B3LYP, PBE0, CASPT2, and DMC), we considered the same HS and LS molecular structures in all calculations. These were computed with PBE0 in vacuum without symmetry constraints, and they were confirmed as minima by vibrational frequency calculations. Previous studies have demonstrated the reliability of PBE0 to yield good results for the geometrical parameters of the Fe-based SCO complexes.86 In particular, the optimization of the LS state leads to an Oh equilibrium geometry, where all six Fe−N distances are equal to 1.95 Å. In contrast, the HS minimum energy geometry has D4h symmetry where the Fe−N distances are 2.18 Å in the equatorial plane and 2.19 Å along the axial axis. These results, which are in agreement with those of previous studies,43 reflect the weak Jahn−Teller effect on an octahedral d6 transition metal complex in its HS state.

Figure 3. Molecular orbitals energy diagram displaying the active space used in the CASSCF calculations (within the red dashed frames). In CASSCF(10,7) the e*g , t2g, and eg molecular orbitals are considered. In CASSCF(10,12) five additional 3d states were introduced in order to account for the double-shell effect. The inset shows the e*g , t2g, and eg molecular orbital isosurface.

implications for the DMC results will be discussed more in detail in the following sections. In case of HF and DFT, the restricted (R) and restricted open-shell (RO) formalism,71 R(O)DFT and R(O)HF, were applied for the LS and HS states, respectively. In DFT, we considered the PBE generalized gradient approximation (GGA)72,73 to the exchange-correlation functional and the two hybrid functionals B3LYP74,75 and PBE0,76 which implement 20% and 25% of exact exchange. The orbitals were obtained by solving the Kohn−Sham and the generalized Kohn−Sham equations, respectively, for PBE and B3LYP/PBE0. Finally, we also considered natural orbitals obtained from the MR-CI calculations performed using PBE0 single-particle orbitals. The scalar-relativistic energy-consistent Hartree−Fock pseudopotentials by Burkatzki, Filippi, and Dolg (BFD)62 with the correspondent polarized valence quadruple-ζ basis sets (pVQZ) were utilized for DMC as well as for the HF, DFT, CASSCF(10,7), and MR-CI calculations. All DMC calculations were carried out with t-moves77 in order to guarantee the validity of the variational principle in the presence of pseudopotentials. The imaginary-time evolution of the Schödinger equation was performed with the short-time approximation with time steps Δτ of 0.01 au. In some cases, the dependence of the results on the finite time step has been double-checked by repeating the calculation for Δτ = 0.005 au; negligible numerical bias has been found. Although the main focus of this work is on DMC, we also present a detailed comparison with the results of HF, DFT, and CASPT2 calculations performed according to the standard prescriptions suggested in the literature to deal with SCO compounds. R(O)HF and R(O)DFT calculations were carried out with the Gaussian09 package.78 We performed all-electron calculations with the triple-ζ valence plus polarization (TZVP) basis set by Schäfer et al.79 as well as effective core calculations with the same BFD pseudopotentials used in DMC. We considered the same exchange-correlation functionals employed to obtain the orbitals for the DMC TWFs, namely, PBE, B3LYP, and PBE0.

3. RESULTS 3.1. DMC Results with 1SJ TWFs. We begin by computing the adiabatic energy gap obtained by DMC with 1SJ TWFs, where HF, PBE, B3LYP, PBE0, and CASSCF(10,7) singleparticle orbitals enter into the Slater determinant. The results are summarized in Table 1. The HS state is computed to be more stable than the LS state in all cases (i.e., ΔEadia < 0). The estimates obtained with CASSCF and HF orbitals are consistent with each other within error bars and they are respectively equal to −1.57(4) eV and −1.54(4) eV, and so are those obtained with B3LYP and PBE0 orbitals, i.e., adia ΔEadia DMC/B3LYP = −0.93(5) eV and ΔEDMC/PBE0 = −0.92(5). Finally, the result for PBE orbitals, ΔEadia DMC/PBE = −1.05(5), looks in quite close agreement with those for B3LYP/PBE0 orbitals, with absolute differences of approximately 0.1 eV. The origin of the numerical discrepancies between the various 1SJ4236


Article

Journal of Chemical Theory and Computation Table 1. LS, HS Energies and Adiabatic Energy Gap of [Fe(NCH)6]2+ Computed with 1SJ-DMC by Considering Different Single-Particle Orbitals for 1SJ Trial Wave Functionsa orbitals

ELS (a.u.)

EHS (a.u.)

ΔEadia DMC (eV)

HF CASSCF PBE B3LYP PBE0

−220.666(1) −220.663(1) −220.693(1) −220.700(1) −220.704(1)

−220.722(1) −220.720(1) −220.732(1) −220.734(1) −220.738(1)

−1.54(4) −1.57(4) −1.05(5) −0.93(5) −0.92(5)

The results are for the imaginary time steps Δτ = 0.01 a.u.−1. Statistical errors are shown in parentheses. a

DMC results can be understood by considering the different properties of the orbitals in the TWFs and the physics that they are or are not able to capture. In this regard, it is very useful to inspect separately the total energies of the LS and HS state, ELS and EHS (see again Table 1). These energies are about 0.9 eV higher for HF and CASSCF than for DFT orbitals. Furthermore, the energies for PBE0 and B3LYP orbitals are the lowest among all the TWFs we considered. Since DMC is a variational method, the lowest energies represent the best upper bound to the true ground state energy. The results indicate that correlation is important already at the level of an effective single-particle theory, and it has important effects on the nodal surface. The better performance of hybrid DFT orbitals for transition metal compounds has already been pointed out by other authors50,58,87,88 and has been ascribed to the good description of the d−p hybridization, which is affected by correlation. This can be appreciated at the qualitative level by looking at the charge isosurface for the eg, t2g, and eg* orbitals computed respectively with CASSCF and PBE0 (Figures 4 and 5). The difference between the results is more marked for the LS than for the HS state. While the CASSCF orbitals appear very similar for both the HS and LS states, PBE0 describes the charge reorganization accompanying the spin transition with the eg and t2g orbitals that becomes much more delocalized over the ligands in the LS state. Conversely, the orbitals that are fully localized on the ligands in CASSCF are instead slightly spread over the Fe in PBE0 (not shown). In an attempt to gain a quantitative understanding of these effects we have evaluated the pair correlation function ρ(r) ∝ ⟨∑ i

Figure 4. Electron density isosurface for the CASSCF eg, t2g, and eg* orbitals.

HS and LS states. In the LS state, the results for CASSCF orbitals is changed by the Jastrow more than that for PBE0 orbitals. Next, we analyzed the difference between the pair correlation function for the 1S(J) TWF with PBE0 and CASSCF orbitals, ρPBE0(r) − ρCAS(r). The results are shown in Figure 8; again there is this difference at around 5 Bohr, which is about the distance from the Fe atom to the N−C bond. This suggests that the one-body electron density is smaller there for PBE0 orbitals than for CASSCF orbitals. Indeed, as seen in Figure 9, the Slater determinant of CASSCF orbitals has more electron density in the N−C bond than the Slater determinant of PBE0 orbitals and less density on the Fe atom. The Jastrow factor attempts to correctly redistribute the electronic charge across the molecule, but it is not entirely able to do so in the LS state due to the limited variational freedom. The underestimation of the electron density at the Fe atom in the LS state and the consequent overestimation at the N−C bond is due to lack of short-range electron correlation in CASSCF. This correlation would favor multiple occupation of the Fe atom, since unlike spin electrons are more able to avoid one another. Finally we note that Kepenekian et al.42 suggested a phenomenological picture for the molecule properties where, in the LS state, the reference Fe d6 configuration is mixed with formal d7 and d8 configurations thus reflecting a substantial

∑ δ(r − rij)⟩ j; j ≠ i

(3)

This is the probability to find two electrons at a distance r from one another. As a first step, we calculated ρ(r) for the 1SJTWFs as well as for simpler Slater TWFs (1S TWFs), where the Jastrow factor in eq 2 has been removed. This provides some hints on how much short-range correlation can be recovered just from the orbitals in the TWFs. The computed ρ(r) are displayed in Figure 6. The large bars on the horizontal axis represent the LS and HS Fe−N and Fe−C distances, respectively (LS is closer than HS). The most notable features are a shoulder at the Fe−N distance and a peak at the Fe−C distance, roughly, for all TWFs. In order to highlight the changes, we took the difference between the pair correlation function of each 1SJ TWF and the respective 1S TWF, ρSJ(r) − ρS(r). This is shown in Figure 7 for the case of PBE0 and CASSCF orbitals. We can see that there is a depletion at short distances due to the short-range correlation. However, there is also a depletion of ρSJ(r) − ρS(r) in the ∼5 Bohr range for both 4237


Article


Figure 7. Difference between the DMC pair correlation function for a 1SJ-TWF and the respective 1S-TWF, ρSJ(r) − ρS(r). The results are shown for the case of CASSCF(10,7) and PBE0 orbitals.

Figure 8. Difference between the DMC pair correlation function for a 1SJ-TWF (1S-TWF) with PBE0 and CASSCF orbitals, ρPBE0(r) − ρCAS(r).

Figure 5. Electron density isosurface for the PBE0 eg, t2g, and eg* orbitals.

Figure 9. Difference between the electron density obtained from the Slater determinants with CASSCF orbitals and PBE0 orbitals (yellow = positive values, cyan = negative values).

just shown to be so important. Therefore, it may be that the orbitals obtained from CASSCF(10,12) would also be more suited for DMC than the CASSCF(10,7) orbitals utilized here. 3.2. DMC Results with nSJ TWFs. A systematic way to improve the DMC results consists in employing nSJ TWFs with a progressively larger number of CSFs (see eq 2). We first consider CSFs constructed with CASSCF single-particle orbitals. The results are presented in Figure 10. The energy of the LS state monotonically decreases as a function of the number of CSFs, while the energy of the HS state is constant within the achieved precision. This ultimately results in an overall reduction of the adiabatic gap that shrinks from −1.57(4) eV for a 1SJ TWF to −0.89(4) eV for a nSJ TWF with a CI weight cutoff of 0.005. This cutoff corresponds to considering 89 CSFs and 170 CSFs, respectively, for the HS

Figure 6. DMC pair correlation function for several 1SJ and 1S-TWFs with different orbitals. The lines along the horizontal axis mark the Fe−N and Fe−C distances (∼4 and ∼6 Bohr) for HS and LS (LS closer than HS).

ligand-to-metal charge transfer. In order to correctly account for this by CASPT2, the authors pointed out the necessity of including additional 3d orbitals entering the active space. These orbitals account for the double-shell effect and indeed lead to an improved description of short-range correlation that we have 4238


Article


determinant with DFT orbitals is not just much better than that of a single Slater determinant with CASSCF orbitals, but it is even superior to that of a multideterminant wave function. Having learned that, it becomes interesting to examine whether the accuracy of the DMC results can be further enhanced by using a nSJ TWF with DFT orbitals instead of CASSCF. However, when we do so, we find that the results do not depend on the number of considered CSFs and the HS and LS total energy does not decrease within the precision of our calculations (see Figure 10). Finally, we also tested n-SJ TWF using natural orbitals (NOs), which have been suggested in previous work as an efficient way to obtain better nodal surfaces and very accurate DMC energies.89,90 We extracted the NOs from a MR-CI calculation performed with the PBE0 single-particle orbitals. Then, a second MR-CI calculation over the Slater determinants with this just-computed set of NOs is utilized in order to determine the CSFs for the nSJ-DMC calculations. Once again the computed HS and LS total energies are not lower than those obtained for a 1SJ-TWF with PBE0 orbitals. The PBE0 orbitals therefore are still the best choice energetically. 3.3. Comparison to Other Techniques. The results for the adiabatic energy gap of HF, DFT, CASSCF, and CASPT2 are summarized in Figure 11 together with those of DMC. In

Figure 10. HS, LS total energies and adiabatic energy gap as a function of the weight cutoff in the multi-Slater determinant expansion of the TWF, eq 2. Total energies are shifted by 6000 eV to clarify the graph. The smallest considered determinant weight cutoff is equal to 0.005, while a determinant weight cutoff equal to 0.9 corresponds to one single Slater determinant. The error bars are smaller than the symbols.

and LS cases. These are the largest number of CSFs that we could include at an affordable computational cost. The results clearly demonstrate the different energy contributions responsible for the stabilization of the HS and LS states. On the one hand, the HS is promoted by the exchange and the correlation effects are already fully captured with a 1SJ TWF. On the other hand, in the LS state, shortrange correlation effects play a crucial role in counterbalancing the exchange and in allowing a larger ligand-to-metal hybridization. As seen before, these effects are poorly accounted for by means of a 1SJ TWF with CASSCF orbitals. Their description really benefit from a multireference treatment of the TWF since the expansion coefficients constitute a large set of variationally optimizable degrees of freedom. Although the application of a large number of CSFs has allowed for the systematic improvement of the DMC results with CASSCF orbitals, we must highlight that these results are still worse in terms of total energies than those obtained with a simple 1SJ TWF with DFT orbitals (see Figure 10). In fact, the 1SJ-DMC total energies with either B3LYP or PBE0 orbitals are about 0.3 eV lower than the best nSJ-DMC estimates with CASSCF orbitals. This means the nodal surface of a Slater

Figure 11. Summary of the all results for the adiabatic energy gap, with (BFD) and without (none) pseudopotentials.

this figure, we have plotted all DMC values that are below the “low energy” line in Figure 10, to indicate the uncertainty in values due to the TWF. We have also assessed the errors due to the pseudopotential, which are small for all theories except CASPT2, which we ascribe to a difference in the basis set for the all-electron and pseudopotential calculations. Our CASPT2 results are in good agreement with those by other authors,42 who selected the same active space, with small deviations that can be attributed to the differences in the considered molecular geometries. 4239



■

Starting from the least accurate methods, HF overestimates the relative HS stability by a large amount, while DFT(PBE) underestimates the HS stability by a slightly lesser amount. CASSCF(10,12) without perturbative corrections for the dynamic correlation corrects HF somewhat but is still about 1 eV away from the set of “high accuracy” methods which include DMC, CASPT2, and the hybrid functionals. This highlights the importance of short-range correlation which we also saw in the orbital choice for DMC. Besides, we note that all methods except DFT(PBE) predict a very stable HS ground state and therefore the molecule would not undergo spincrossover. This is consistent with previous studies.41,43,56 The hybrid functionals, DMC, and CASPT2, vary over approximately 0.5 eV for the adiabatic energy gap. This is not yet ideal and, moreover, it is not clear what the reference value should be. Probably the most reliable estimate is the CCSD(T) calculation performed by Lawson Daku et al.43 at a slightly different geometry than we used for our large study. They obtained a HS−LS energy difference of −0.25 eV, neglecting relativistic corrections. On the same geometry as theirs, DMC returns an energy difference of −0.55(6) eV, a discrepancy of 0.30(6) eV.

AUTHOR INFORMATION

Corresponding Authors

*(M.F.) E-mail: [email protected]. *(L.K.W.) E-mail: [email protected]. *(S.S.) E-mail: [email protected]. *(A.D.) E-mail: [email protected]. Funding

M.F. acknowledges the Spanish Government for financial support (Projects MAT2011-25972 and MAT2014-54025-P), the University of Barcelona, and the Fundació Universitaria Agusti ́ Pedro i Pons for funding her stay in Dublin, where this work was started. L.K.W. was supported by NSF DMR1206242. S.S. was supported by the European Research Council (QUEST project). A.D. was supported by the EU project ACMOL (FP7-FET GA618082). Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Gutlich, P., Goodwin, H., Eds. Spin Crossover in Transition Metal Compounds, 1st ed.; Springer-Verlag: Berlin, 2004. (2) Gutlich, P.; Garcia, Y.; Goodwin, H. A. Chem. Soc. Rev. 2000, 29, 419. (3) Harvey, J. N.; Poli, R.; Smith, K. M. Coord. Chem. Rev. 2003, 238, 347. (4) Rocha, A. R.; Garcia-suarez, V. M.; Bailey, S. W.; Lambert, C. J.; Ferrer, J.; Sanvito, S. Nat. Mater. 2005, 4, 335. (5) Sanvito, S. Chem. Soc. Rev. 2011, 40, 3336. (6) Droghetti, A.; Sanvito, S. Phys. Rev. Lett. 2011, 107, 047201. (7) Baadji, N.; Sanvito, S. Phys. Rev. Lett. 2012, 108, 217201. (8) Aravena, D.; Ruiz, E. J. Am. Chem. Soc. 2012, 134, 777. (9) Alam, M.; Stocker, M.; Gieb, K.; Müller, P.; Haryono, M.; Student, K.; Grohmann, A. Angew. Chem., Int. Ed. 2010, 49, 1159. (10) Osorio, E. A.; Moth-Poulsen, K.; van der Zant, H. S. J.; Paaske, J.; Hedegård, P.; Flensberg, K.; Bendix, J.; Bjørnholm, T. Nano Lett. 2010, 10, 105. (11) Meded, V.; Bagrets, A.; Fink, K.; Chandrasekar, R.; Ruben, M.; Evers, F.; Bernand-Mantel, A.; Seldenthuis, J. S.; Beukman, A.; van der Zant, H. S. J. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 245415. (12) Prins, F.; Monrabal-Capilla, M.; Osorio, E. A.; Coronado, E.; van der Zant, H. S. J. Adv. Mater. 2011, 23, 1545. (13) Miyamachi, T.; Gruber, M.; Davesne, V.; Bowen, M.; Boukari, S.; Joly, L.; Scheurer, F.; Rogez, G.; Yamada, T. K.; Beaurepaire, P. O. E.; Wulfhekel, W.; Ohresser, P. Nat. Commun. 2012, 3, 938. (14) Cambi, L.; Szegö, L. Ber. Dtsch. Chem. Ges. B 1931, 64, 2591. (15) Paulsen, H.; Schünemann, V.; Wolny, J. A. Eur. J. Inorg. Chem. 2013, 2013, 628−641. (16) Launay, J.-P.; Verdaguer, M. Electrons in Molecules: From Basic Principles to Molecular Electronics, 1st ed.; Oxford University Press: Oxford, 2014; pp 109−113. (17) Swart, M.; Groenhof, A. R.; Ehlers, A. W.; Lammertsma, K. J. Phys. Chem. A 2004, 108, 5479. (18) Zein, S.; Borshch, S. A.; Fleurat-Lessard, P.; Casida, M. E.; Chermette, H. J. Chem. Phys. 2007, 126, 014105. (19) Lawson Daku, L. M.; Vargas, A.; Hauser, A.; Fouqueau, A.; Casida, M. E. ChemPhysChem 2005, 6, 1393. (20) Reiher, M. Inorg. Chem. 2002, 41, 6928. (21) Reiher, M.; Salomon, O.; Artur Hess, B. Theor. Chem. Acc. 2001, 107, 48. (22) Salomon, O.; Reiher, M.; Hess, B. A. J. Chem. Phys. 2002, 117, 4729−4737. (23) Chen, J.; Millis, A. J.; Marianetti, C. A. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 241111. (24) Vela, S.; Fumanal, M.; Ribas-Arino, J.; Robert, V. Phys. Chem. Chem. Phys. 2015, 17, 16306.

4. CONCLUSIONS We have presented a systematic study of the HS and LS electronic energies by DMC for the prototypical case of the [Fe(NCH)6]2+ model compound. In particular, we have analyzed the effect of using several different TWFs and therefore of the fixed node approximation. Single Slater− Jastrow TWFs with either HF or CASSCF orbitals provide a poor description of the nodal surface likely because these orbitals do not reflect the correct metal−ligand hybridization and charge transfer, which is affected strongly by short-range correlation. The results are then systematically improved by replacing the single Slater determinant with a linear combination of CSFs. In contrast, a single Slater−Jastrow TWF with hybrid DFT orbitals obtains similar or better total energies than a multi-Jastrow TWF with CASSCF or HF orbitals. Overall, the DMC results return an estimate of the adiabatic energy gap of [Fe(NCH)6]2+ for our test geometry which ranges between −0.95 and −0.85 eV. On a slightly different geometry, we compared this method to CCSD(T) results, which differ by around 0.3 eV. Our results highlight the importance of short-range correlation in determining the best nodal structure for the system. Fixed node DMC with PBE0 orbitals seems to be very close to the optimal nodal surface, with minimal improvements even with many added determinants. The DMC results are quite robust to the nodal surface, obtaining the same HS−LS adiabatic energy gap for many reasonable nodal surfaces and an energy difference within a few tenths of an eV of CCSD(T) results. This particular molecule appears to be an interesting test case for future studies to further tighten the constraints on the correct HS−LS adiabatic energy gap.

■

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00332. Molecular geometry for the low-spin (TXT) Molecular geometry for the high-spin (TXT) Energy data (XLSX) 4240


Article

Journal of Chemical Theory and Computation (25) Anisimov, V. I.; Zaanen, J.; Andersen, O. K. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 44, 943. (26) Anisimov, V. I.; Solovyev, I. V.; Korotin, M. A.; Czyżyk, M. T.; Sawatzky, G. A. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 16929. (27) Georges, A.; Kotliar, G.; Krauth, W.; Rozenberg, M. J. Rev. Mod. Phys. 1996, 68, 13. (28) Kotliar, G.; Savrasov, S. Y.; Haule, K.; Oudovenko, V. S.; Parcollet, O.; Marianetti, C. A. Rev. Mod. Phys. 2006, 78, 865. (29) Andersson, K.; Malmqvist, P.; Roos, B. O. J. Chem. Phys. 1992, 96, 1218. (30) Bolvin, H. J. Phys. Chem. A 1998, 102, 7525. (31) Pierloot, K.; Vancoillie, S. J. Chem. Phys. 2006, 125, 124303. (32) Pierloot, K.; Vancoillie, S. J. Chem. Phys. 2008, 128, 034104. (33) Ordejon, B.; de Graaf, C.; Sousa, C. J. Am. Chem. Soc. 2008, 130, 13961. (34) Vancoillie, S.; Zhao, H.; Radoń, M.; Pierloot, K. J. Chem. Theory Comput. 2010, 6, 576. (35) Graaf, C. D.; Sousa, C. Int. J. Quantum Chem. 2011, 111, 3385. (36) Rudavskyi, A.; Havenith, R. W. A.; Broer, R.; de Graaf, C.; Sousa, C. Dalton Trans. 2013, 42, 14702. (37) Rudavskyi, A.; Sousa, C.; de Graaf, C.; Havenith, R. W. A.; Broer, R. J. Chem. Phys. 2014, 140, 184318. (38) Vela, S.; Fumanal, M.; Ribas-Ariño, J.; Robert, V. J. Comput. Chem. 2016, 37, 947−953. (39) Stein, C. J.; Reiher, M. J. Chem. Theory Comput. 2016, 12, 1760−1771. (40) Ghigo, G.; Roos, B. O.; Malmqvist, P.-A. Chem. Phys. Lett. 2004, 396, 142. (41) Kepenekian, M.; Robert, V.; Le Guennic, B. J. Chem. Phys. 2009, 131, 114702. (42) Kepenekian, M.; Robert, V.; Le Guennic, B.; De Graaf, C. J. Comput. Chem. 2009, 30, 2327. (43) Lawson Daku, L. M.; Aquilante, F.; Robinson, T. W.; Hauser, A. J. Chem. Theory Comput. 2012, 8, 4216. (44) Droghetti, A.; Alfè, D.; Sanvito, S. J. Chem. Phys. 2012, 137, 124303. (45) Droghetti, A.; Alfè, D.; Sanvito, S. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 205114. (46) Foulkes, W. M. C.; Mitas, L.; Needs, R. J.; Rajagopal, G. Rev. Mod. Phys. 2001, 73, 33. (47) Gillan, M. J.; Towler, M. D.; Alfe, D. Psi-k Highlight of the Month; 2011; Vol. 103. (48) Wagner, L.; Mitas, L. Chem. Phys. Lett. 2003, 370, 412. (49) Koseki, J.; Maezono, R.; Tachikawa, M.; Towler, M. D.; Needs, R. J. J. Chem. Phys. 2008, 129, 085103. (50) Horváthová, L.; Dubecký, M.; Mitas, L.; Štich, I. Phys. Rev. Lett. 2012, 109, 053001. (51) Granatier, J.; Dubecký, M.; Lazar, P.; Otyepka, M.; Hobza, P. J. Chem. Theory Comput. 2013, 9, 1461. (52) Wagner, L. K. J. Phys.: Condens. Matter 2007, 19, 343201. (53) Kolorenč, J. c. v.; Mitas, L. Phys. Rev. Lett. 2008, 101, 185502. (54) Wagner, L. K.; Abbamonte, P. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 125129. (55) Foyevtsova, K.; Krogel, J. T.; Kim, J.; Kent, P. R. C.; Dagotto, E.; Reboredo, F. A. Phys. Rev. X 2014, 4, 031003. (56) Kepenekian, M.; Le Guennic, B.; Robert, V. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 094428. (57) Toulouse, J.; Umrigar, C. J. J. Chem. Phys. 2008, 128, 174101. (58) Kolorenč, J. c. v.; Hu, S.; Mitas, L. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 115108. (59) Morales, M. A.; McMinis, J.; Clark, B. K.; Kim, J.; Scuseria, G. E. J. Chem. Theory Comput. 2012, 8, 2181. (60) Trail, J. R.; Needs, R. J. J. Chem. Phys. 2005, 122, 174109. (61) Trail, J. R.; Needs, R. J. J. Chem. Phys. 2005, 122, 014112. (62) Burkatzki, M.; Filippi, C.; Dolg, M. J. Chem. Phys. 2007, 126, 234105. (63) Trail, J. R.; Needs, R. J. J. Chem. Phys. 2013, 139, 014101.

(64) Wagner, L. K.; Bajdich, M.; Mitas, L. J. Comput. Phys. 2009, 228, 3390. (65) Toulouse, J.; Umrigar, C. J. J. Chem. Phys. 2007, 126, 084102. (66) Umrigar, C. J.; Toulouse, J.; Filippi, C.; Sorella, S.; Hennig, R. G. Phys. Rev. Lett. 2007, 98, 110201. (67) Koziol, L.; Morales, M. M. J. Chem. Phys. 2014, 140, 224316. (68) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (69) Botch, B. H.; Dunning, T. H.; Harrison, J. F. J. Chem. Phys. 1981, 75, 3466. (70) Andersson, K.; Roos, B. O. Chem. Phys. Lett. 1992, 191, 507. (71) Roothaan, C. C. J. Rev. Mod. Phys. 1960, 32, 179. (72) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (73) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396. (74) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (75) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785. (76) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (77) Casula, M. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 161102. (78) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision E.01; Gaussian, Inc.: Wallingford, CT, 2009. (79) Schäfer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829. (80) Karlström, G.; Lindh, R.; Malmqvist, P.-A.; Roos, B. O.; Ryde, U.; Veryazov, V.; Widmark, P.-O.; Cossi, M.; Schimmelpfennig, B.; Neogrady, P.; Seijo, L. Comput. Mater. Sci. 2003, 28, 222 (Proceedings of the Symposium on Software Development for Process and Materials Design). (81) Roos, B. O.; Lindh, R.; Malmqvist, P.-A.; Veryazov, V.; Widmark, P.-O. J. Phys. Chem. A 2004, 108, 2851. (82) Koch, H.; Sánchez de Merás, A.; Pedersen, T. B. J. Chem. Phys. 2003, 118, 9481. (83) Hess, B. A. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 33, 3742. (84) Reiher, M.; Wolf, A. J. Chem. Phys. 2004, 121, 2037−2047. (85) Reiher, M.; Wolf, A. J. Chem. Phys. 2004, 121, 10945−10956. (86) Papai, M.; Vanko, G.; de Graaf, C.; Rozgonyi, T. J. Chem. Theory Comput. 2013, 9, 509. (87) Wagner, L. K. J. Chem. Phys. 2013, 138, 094106. (88) Dubecky, M.; Derian, R.; Jurecka, P.; Mitas, L.; Hobza, P.; Otyepka, M. Phys. Chem. Chem. Phys. 2014, 16, 20915. (89) Grossman, J. C.; Mitás,̌ L. c. v. Phys. Rev. Lett. 1995, 74, 1323. (90) Kent, P. R. C.; Hood, R. Q.; Towler, M. D.; Needs, R. J.; Rajagopal, G. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 57, 15293.

4241


Diffusion Monte Carlo Perspective on the Spin-State Energetics of [Fe

Recommend Documents