Detailed Pair Natural Orbital-Based Coupled Cluster Studies of Spin Crossover Energetics.

In this work, a detailed study of spin-state splittings in three spin crossover model compounds with DLPNO-CCSD(T) is presented. The performance in comparison to canonical CCSD(T) is assessed in detail. It was found that spin-state splittings with chemical accuracy, compared to the canonical results, are achieved when the full iterative triples (T1) scheme and TightPNO settings are applied and relativistic effects are taken into account. Having established the level of accuracy that can be reached relative to the canonical results, we have undertaken a detailed basis set study in the second part of the study. The slow convergence of the results of correlated calculations with respect to basis set extension is particularly acute for spin-state splittings for reasons discussed in detail in this Article. In fact, for some of the studied systems, 5Z basis sets are necessary in order to come close to the basis set limit that is estimated here by basis set extrapolation. Finally, the results of the present work are compared to available literature. In general, acceptable agreement with previous CCSD(T) results is found, although notable deviations stemming from differences in methodology and basis sets are noted. It is noted that the published CASPT2 numbers are far away from the extrapolated CCSD(T) numbers. In addition, dynamic quantum Monte Carlo results differ by several tens of kcal/mol from the CCSD(T) numbers. A comparison to DFT results produced with a range of popular density functionals shows the expected scattering of results and showcases the difficulty of applying DFT to spin-state energies.

Introduction

The relative energies of different spin states in transition metal complexes, also called spin-state splittings, are of central importance for fields like spin crossover, magnetism, and bioinorganic chemistry.[1−3] Here, multistate reactivity and catalysis pose formidable challenges for established computational methods such as density functional theory (DFT) which has been evaluated in numerous studies.[4−6] On the other hand, wave function theory (WFT) methods, though systematically improvable and thus, in principle, as accurate as one desires, often suffer from their complexity, nonblackbox character, and, most of all, the often overwhelming computational expense.[7−9] It is thus not at all trivial to decide which method one should apply when faced with the task of computing accurate spin-state energies of a given compound.[10,11]

We would therefore like to briefly reiterate the underlying physics of the problem at hand. The correlation energy, defined as the energy difference between the full configuration interaction (FCI) limit and the Hartree–Fock (HF) level, is made up of two parts: static and dynamic electron correlation. Static electron correlation arises from energetically close-lying or degenerate electronic configurations such that more than one Slater determinant is needed to provide a qualitatively correct description of the wave function. This can, for example, be treated by multiconfigurational methods like complete active space self-consistent field (CASSCF) and its variants.

Dynamic electron correlation, on the other hand, stems from the instantaneous interaction of two or more electrons and is not properly described by the mean-field nature of the HF method which does not account for the Coulomb cusp. This is remedied by including many excited Slater determinants in the wave function thus providing the electrons with more freedom to evade each other. This can be realized in a number of ways, e.g., through the methods of configuration interaction (CI), many-body perturbation theory (usually truncated at second order, MP2), or coupled cluster theory (CC). Especially the latter, in the variant including single, double, and perturbative triple excitations (CCSD(T)), is often called the gold standard of quantum chemistry.[12,13] Dynamical correlation in systems with substantial multiconfigurational characters are usually treated with a second-order perturbation method on top of a CASSCF wave function such as CASPT2[14] and NEVPT2.[15,16]

When investigating spin-state energetics, one is concerned with two or more states that differ significantly in their electronic structure. This is due to the interelectronic repulsion depending on the orbital occupancy and the relative spin alignment of the respective electrons. In spin crossover phenomena, the orbital occupancy changes for the involved spin states, e.g., from t2g6eg0 (LS) to t2g4eg2 (HS) for Fe(II). Given different numbers of unpaired electrons and thus different amounts of interelectronic repulsion in spin-crossover phenomena, the differential dynamic correlation energy of the constituting states can become very large. In particular, the number of antiparallel spin pairs is higher in the low spin states than in the high spin states. In the d-shell of Fe(II), there are six antiparallel spin pairs in the low spin (1A1g) state, three of which are intraorbital while there are only five antiparallel spin pairs with one intraorbital pair in the high spin 5T2g state. Since the electron repulsion is generally larger for antiparallel intraorbital electron pairs than for antiparallel interorbital pairs and much larger than for parallel interorbital pairs, this is also the order in which the dynamic correlation contributions increase. Hence, generally, there will be an increased amount of dynamic correlation in the low spin states. Errors in this differential dynamic correlation energy will not cancel out upon taking the energy difference between the high and low spin states. Consequently, one needs to recover nearly 100% of the physical correlation energy in order to calculate spin-state splittings accurately, not just the basis set correlation energy of a given finite basis. In other words, when attempting to compute spin-state splittings, one cannot hope for systematic error cancellation but rather has to attempt to recover the full correlation energy of both states. If a method is chosen that recovers significantly less than 100% of the correlation energy, a systematical bias will be introduced into the resulting spin-state splitting.

This clearly constitutes a formidable methodological challenge given the notoriously slow convergence of the dynamical electron correlation with respect to the orbital basis set size.[17,18] Very large basis sets of quintuple zeta quality or even beyond together with basis set extrapolation or explicit correlation techniques (F12-methods) have to be employed. Unfortunately, the latter have not yet found their way into the chemistry of open-shell transition metal ions. On the other hand, extended basis set extrapolation is not feasible when employing methods such as CCSD(T) given its steep N7 scaling. Obviously, explicitly correlated methods are a viable alternative here but have yet to be developed to be fully applicable to transition metal systems.[19,20]

Herein, we investigate the application of a local coupled cluster approximation to the problem of spin-state splittings, namely, domain-based pair natural orbital coupled cluster with singles, doubles, and perturbative triples (DLPNO-CCSD(T)).[21−24] This method has previously been shown to provide results of near-CCSD(T) quality at a fraction of the cost and linear scaling with respect to system size.[25] Moreover, it has been supplemented by the recent implementation of iterative triples (denoted (T1)) further enhancing its capabilities to mirror the accuracy of canonical CCSD(T).[24] We present a detailed analysis of the capability of DLPNO-CCSD(T) to recover the canonical correlation energy, both with the semicanonical (T0) and the full iterative triples (T1). It is shown that the previous difficulties of the DLPNO method with spin-state energetics can be traced back to the (T0) approximation which is analyzed and compared to the full iterative triples in detail. We also compare our results with the available literature and DFT results and discuss the sources of arising discrepancies from the results presented here.

This Article is organized as follows: In the first section, we demonstrate a benchmark study of three reference complexes (hereafter called spin crossover model complexes) ([Fe(H2O)6]2+, [Fe(NH3)6]2+, and [Fe(NCH)6]2+ (Figure ), where we compare DLPNO-CCSD(T) results with canonical CCSD(T) ones to establish the performance of the DLPNO method in a direct comparison to its parent. This necessitates one to choose the investigated systems such that a treatment with canonical CCSD(T) including all physically meaningful effects is still feasible. Also, we opted for compounds that are well established as spin crossover model systems to be able to draw a comparison to previous results by other groups. Furthermore, the ligand field splittings of the respective complexes are of sufficiently different sizes. As the ligand field splittings are directly connected to the spin-state splittings, this is expected to introduce some variety into the results which in turn should reduce the danger of “cherrypicking” in the study of these compounds. We hope that this approach allows us to draw more generally valid conclusions.

Figure 1

Optimized structures of the compounds investigated in this work: [Fe(H2O)6]2+, [Fe(NH3)6]2+, and [Fe(NCH)6]2+ (from left to right).

In addition to validating the DLPNO-CCSD(T) method by comparison to its canonical counterpart, we make a dedicated effort to reach the basis set limit through extrapolation techniques and, hence, to come as close to the true (experimentally unknown) spin-state splitting energy as possible. Finally, we compare our results to those obtained by DFT with a series of popular density functional approximations (DFAs) and to literature studies.

Computational Details

All calculations in this work were conducted with the ORCA program suite.[26,27] For geometry optimizations and DFT single point calculations, we used the publicly available version 4.1.2, whereas for all coupled cluster calculations the current development version was employed. All calculations described herein make use of the second-order Douglas–Kroll–Hess Hamiltonian (DKH2)[28,29] to include scalar relativistic effects, unless otherwise stated. The effects of scalar relativity are fairly limited in these systems as shown in the Results and Discussion section.

All complexes studied here were optimized separately in the HS and LS states with DFT using the BP86[30,31] functional and the DKH-def2-TZVPP[32] basis set. Furthermore, a fine grid (grid6 in ORCA nomenclature) and tight optimization criteria (tightopt) were used. To speed up convergence, the RI[33,34] approximation was made use of with the SARC/J fitting basis set which is the decontracted version of the general purpose def2/J[35] basis set. To verify that a minimum was reached, we conducted frequency calculations which showed no imaginary frequencies in any case. The optimized coordinates are listed in the Supporting Information.

DFT single point calculations were conducted with numerous density functional approximations of different types. These were the (meta-)GGA functionals BP86, TPSS,[36] OPBE,[37,38] OLYP,[39] and M06L,[40] the (meta-)hybrids TPSSh,[41] PBE0,[42,43] PW6B95,[44] B3LYP(*),[6] and M06,[45] the range-separated ωB97X-D3,[46,47] and the two double-hybrid functionals B2PLYP[48] and DSD-BLYP.[49] All of these calculations were done with the cc-pwCVQZ-DK/cc-pVQZ-DK[50,51] basis sets on Fe and the rest of the molecule, respectively, with a tight integration grid (grid6), very tight convergence criteria (verytightscf), and without RI approximations.

As reference for the CC calculations, ROHF and RHF determinants were used for the HS and LS states, respectively. All SCF calculations were performed in the absence of any RI approximations in order to exclude any possible additional error and were converged tightly down to energy changes of less than 10–9 Hartree (verytightscf). For the benchmark study (see below), the cc-pwCVTZ-DK and cc-pVDZ-DK basis sets were used on Fe and the rest of the molecule, respectively. For the basis set convergence study, cc-pwCVnZ-DK/cc-pVnZ-DK was used with n = T,Q,5.

The obtained orbitals were fed into the separate canonical and local CC codes resulting in ROHF-UCC and RHF-RCC for the HS and LS states, respectively, in keeping with the literature notation. For the DLPNO calculations, the 3s and 3p outer core orbitals are included in the correlation treatment, and the 1s and 2s inner core orbitals are kept frozen.[52] Also, the large automatically generated Autoaux fitting basis set was used, and unless stated otherwise, the tightPNO thresholds and full iterative triples ((T1) in ORCA nomenclature) were employed.

Basis set extrapolation was carried out employing the formula for the correlation energy suggested by Helgaker and coworkers[17]with the extrapolated energy, E∞corr, the two basis set hierarchies, X and Y, and the correlation energy at basis set X and Y, E(X)corr and E(Y)corr. For the SCF part, a similar expression was usedwith α = 3.9

The energy difference between the high and low spin states, ΔEHL, used throughout this work is defined as

Results and Discussion

Geometries

The geometry optimizations were carried out without symmetry restraints and, as expected, result in an approximately octahedral arrangement of the ligands (Figure ). It is noteworthy that, due to the orbital degeneracy of the 5T2g state, the HS structures feature a small Jahn–Teller effect resulting in axially compressed/elongated metal–ligand bond lengths in the cases of [Fe(H2O)6]2+ and [Fe(NCH)6]2+/[Fe(NH3)6]2+. The metal–ligand bond lengths are listed in Table .

Table 1

Fe-X (X = O,N) Bond Lengths in Å for the Three Complexes Studied Here

	[Fe(H2O)6]2+		[Fe(NH3)6]2+		[Fe(NCH)6]2+
	HS	LS	HS	LS	HS	LS
Fe-X1	2.118	2.009	2.299	2.060	2.123	1.858
Fe-X2	2.118	2.009	2.299	2.060	2.123	1.858
Fe-X3	2.155	2.009	2.278	2.061	2.143	1.858
Fe-X4	2.154	2.009	2.287	2.061	2.143	1.858
Fe-X5	2.158	2.009	2.273	2.061	2.143	1.858
Fe-X6	2.160	2.008	2.269	2.061	2.143	1.857

Benchmarks against Canonical CCSD(T)

In the first step of this work, we establish the performance of the local DLPNO-CCSD(T) method in comparison to canonical CCSD(T). To this end we chose the smallest, yet physically sound, setup to include all significant effects but still keep computing times within reasonable bounds. Since it has previously been shown that correlation of the semicore electrons and scalar relativistic effects contribute in a non-negligible fashion to spin-state splittings,[53] we chose to include both effects and adopt the necessary basis sets, i.e., cc-pVnZ-DK/cc-pwCVnZ-DK, which have been recontracted for use with the DKH2 Hamiltonian. Note that for this benchmark study a cc-pVDZ-DK/cc-pwCVTZ-DK basis set combination was used to keep computing times as manageable as possible, whereas in the basis set study (see below) only combinations with matching basis set hierarchies were used.

A crucial step in the setup of all coupled cluster calculations is the choice of a suitable set of reference orbitals. Most often these are obtained from HF, but KS or ones resulting from multireference calculations such as CASSCF may also be used. It has been shown that the choice of reference can have a substantial impact on the results of the CC calculation so extra care must be taken.[54]

The molecules investigated here feature an 1A and 5T state for the low and high spin species, respectively. Whereas the 1A state is nondegenerate and can safely be described by an RHF wave function, the 5T state is orbitally triply degenerate. This degeneracy is lifted by the Jahn–Teller effect as is also reflected in the optimized coordinates. However, the three spatial components of 5T2g parentage are still close in energy and complicate the treatment. Clearly, a UHF wave function cannot produce a proper average of the involved microstates which renders the resulting orbitals unreliable. It has been shown that restricted open-shell orbitals from ROHF or ROKS are feasible alternatives. The most rigorous approach, however, may be a state-averaged CASSCF wave function with a small active space containing all t2g orbitals which would produce properly averaged d-orbitals. These can then be used as a basis for a subsequent coupled cluster calculation of choice. Unfortunately, we observed in the present study that CASSCF reference orbitals frequently led to severe CC convergence issues which we did not manage to overcome. Consequently, in this study, we opted for the use of ROHF orbitals for a single high-spin open-shell state.

The results of the comparison of DLPNO and canonical results for the three investigated complexes are given in Tables –4.

Table 2

Correlation Energy Contributions and Spin-State Splitting of HS and LS States of [Fe(H2O)6]2+a

	HS				LS
	Ecan	EDLPNO	Δc/mEh	Ecan/EDLPNO	Ecan	EDLPNO	Δc/mEh	Ecan/EDLPNO	ΔEcan/ kcal mol–1	ΔEDLPNO/ kcal mol–1	ΔΔEd/ kcal mol–1
SD	–1.89151	–1.88793	–3.6	0.998	–1.95155	–1.94925	–2.3	0.999	37.7	38.5	–0.8
(T0)	–	–0.03493	–3.6	0.907	–	–0.04215	–5.4	0.886	–	4.5	1.2e
(T1)	–0.03853	–0.03658	–2.0	0.949	–0.04758	–0.04595	–1.6	0.966	5.7	5.9	–0.2
Ecorrb	–1.93004	–1.92451	–5.5	0.997	–1.99913	–1.99520	–3.9	0.998	43.4	44.4	–1.0

All energies in Ha unless stated otherwise. TightPNO thresholds were used for all DLPNO calculations. Ecan and EDLPNO are energies computed with CCSD(T) and DLPNO-CCSD(T), respectively.

Calculated as the sum of the SD and (T1) values.

Calculated as Ecan – EDLPNO.

Calculated as ΔEcan – ΔEDLPNO.

Calculated as (T1)(ΔEcan) – (T0)(ΔEDLPNO).

Table 4

Correlation Energy Contributions and Spin-State Splitting of HS and LS States of [Fe(NCH)6]2+a

	HS				LS
	Ecan	EDLPNO	Δc/mEh	Ecan/EDLPNO	Ecan	EDLPNO	Δc/mEh	Ecan/EDLPNO	ΔEcan/ kcal mol–1	ΔEDLPNO/ kcal mol–1	ΔΔEd/ kcal mol–1
SD	–2.39489	–2.39063	–4.3	0.998	–2.51667	–2.51381	–2.9	0.999	76.4	77.3	–0.9
(T0)	–	–0.08424	–9.2	0.902	–	–0.10116	–15.4	0.868	–	10.6	3.9d
(T1)	–0.09339	–0.09038	–3.0	0.968	–0.11654	–0.11271	–3.8	0.967	14.5	14.0	0.5
Ecorrb	–2.48828	–2.48101	–7.3	0.997	–2.63321	–2.62652	–6.7	0.997	91.0	91.3	–0.4

All energies in Ha unless stated otherwise. TightPNO thresholds were used for all DLPNO calculations. Ecan and EDLPNO are energies computed with CCSD(T) and DLPNO-CCSD(T), respectively.

Calculated as the sum of the SD and (T1) values.

Calculated as Ecan – EDLPNO.

Calculated as ΔEcan – ΔEDLPNO.

Calculated as (T1)(ΔEcan) – (T0)(ΔEDLPNO).

As is apparent, the DLPNO error (denoted Δ) in the singles and doubles (SD) is very small resulting in greater than 99.8% of the canonical CCSD correlation energy being recovered. This translates to DLPNO errors in spin-state splittings (ΔΔE) of less than 1 kcal/mol which is regarded as chemically accurate. However, upon inspection of the results of the perturbative triples, larger errors are revealed. Until recently, the standard formulation of the triples for open- and closed-shell systems in local CC codes mostly amounted to the so-called semicanonical triples (usually denoted (T0)) which neglects the coupling between different triples via the off-diagonal Fock matrix elements and can be computed very efficiently. This usually leads to relative energies with similar accuracy as the canonically computed counterparts. However, we have recently recognized that this does not always apply and that in certain cases large deviations of (T0) from canonical triples can be observed, in particular, for open-shell systems. The presented results show that this is indeed the case. The semicanonical triples are only able to recover at most 91.4% of the canonical triples correlation energy with an absolute deviation of up to 15 mEh which corresponds to only 86.8% of the canonical triples. This translates to deviations of up to 3.9 kcal/mol in the triples contribution to the spin-state splitting which is unacceptably high. It also appears that the deviation from canonical triples increases with the amount of total correlation energy as can be seen from comparing results for [Fe(H2O)6]2+ and [Fe(NCH)6]2+. It thus stands to reason that the (T0) results will deteriorate with growing system size. Finally, it is noteworthy that the relative amount of triples correlation energies recovered is consistently lower for the LS complexes by 2%–3%.

The iterative triples recently implemented in ORCA DLPNO-(T1) are able to alleviate this situation, producing results that are a more accurate approximation to the canonical (T) correction. As shown in Tables –4, between 95% and 97% of the canonical triples correlation energy is recovered. As the (T1) approach behaves consistently for both the HS and LS states, roughly the same relative amount of correlation energy is recovered, yielding accurate results for the relative spin-state energies—which are the goal of this endeavor—and resulting in relative energies of HS and LS states that differ by at most 1 kcal/mol from the canonical reference.

The main approximation in (T0) is to neglect off-diagonal internal Fock matrix elements in the local basis (Fi,j). As mentioned earlier, this becomes a poor approximation if there are low-lying electronic states that eventually lead to denominators that are (relatively) small. In this case, the diagonal elements Fi,i are not good approximations to actual orbital energies εi, and consequently, the errors become large and somewhat irregular. The same does not happen for (T1) that uses an iterative algorithm to solve for the actual triples correction. In this case, only the PNO error remains.

Generally, while the canonical CCSD correlation energy is typically reproduced with an accuracy of about 99.8%–99.9%, the same has never been the case for the triples correction in any variant. Here, one typically reaches not more than about 98% of the canonical triples correction. Since the (T) contribution to the correlation energy is much smaller than the CCSD correlation energy, typically below 5%, the error made by missing about 2% of the (T) correction is usually tolerable. In the case of the DLPNO-CCSD energy, most of the PNO error can be recovered from a perturbation–theoretical MP2 correction. Unfortunately, the same is not possible for (T), since by its very nature it is already based on perturbation theory, and there does not appear to be a straightforward way to estimate the PNO (or in the case of the triples the TNO) error.

Relativistic Effects and PNO Thresholds

As mentioned above, several factors influencing the results were included to arrive at the results given here, namely, scalar relativistic effects and the truncation thresholds of the DLPNO method. These were investigated for the case of [Fe(NH3)6]2+ at fixed geometries, and the results are given in Table .

Table 5

Comparison of Spin-State Splittings of [Fe(NH3)6]2+ with Relaxed PNO Settings and No Scalar Relativistic Correction at Fixed Geometries

	DKH2 + tightPNO	No rel + tightPNO	DKH2 + normalPNO	DKH2 + tighterPNOa
SCF	–74.27	–76.77	–74.27	–74.27
SD	51.41	51.51	48.66	51.35
(T1)	8.09	8.12	8.00	8.15
Total	–14.77	–17.14	–17.61	–14.77

TcutPNO = 10–8.

As can be seen, either neglecting scalar relativistic effects or relaxing the PNO thresholds introduces non-negligible deviations of 2–3 kcal/mol in the spin-state splitting. We emphasize that the deviation observed in the nonrelativistic calculation solely stems from the SCF part. This is to some extent surprising, as electron correlation and relativistic effects are known to not be additive,[29] although this has been previously observed by other groups.[50,56] However, at least for the spin-state energetics investigated here, this seems not to be the case, and consequently, it would, in principle, be possible to estimate the scalar relativistic effects at the SCF level.

It also becomes apparent that tightening the TCutPNO threshold by a factor of 10 does not further improve the results. While the absolute amount of correlation energy recovered is increased, this does not affect the relative energies in a significant way. This finding suggests that the TightPNO default thresholds constitute a near optimal balance between accuracy and computational expense.

Reaching the Basis Set Limit

In the second step of this work, we were concerned with approaching the basis set limit of the DLPNO-CCSD(T) calculations. This will also be important for any practical applications of the methodology proposed in this Article. To this end, calculations at triple, quadruple, and even quintuple zeta levels were conducted together with TZ/QZ and QZ/5Z basis set extrapolation. The results are shown in Table –8 and Figures –4; the absolute energies are given in the SI. We assume here that the quality of the DLPNO-CCSD(T) method with respect to canonical CCSD(T) does not deteriorate with increasing basis set size. Since demonstrating the correctness of this assumption would require comparing both methods at the respective basis set sizes, a direct proof for the given systems is computationally prohibitive. However, we note that several studies in the literature did not find the results of DLPNO-CCSD(T) to be dependent on the basis set size.[25,55] Although the cited papers do not deal with spin-state energetics, we do not anticipate the work presented here to pose a notable exception.

Table 6

Contributions to ΔEHL of [Fe(H2O)6]2+ Depending on the Basis Set Hierarchya

	TZ	QZ	5Z	TZ/QZ	QZ/5Z
SCF	–80.8	–81.1	–81.2	–81.3	–81.2
SD	36.7	39.2	40.3	41.1	41.4
(T1)	5.8	6.2	6.3	6.5	6.4
Ecorr	42.5	45.5	46.6	47.6	47.8
Total	–38.2	–35.6	–34.5	–33.6	–33.3

All energies are in kcal/mol.

Table 8

Contributions to ΔEHL of [Fe(NCH)6]2+ Depending on the Basis Set Hierarchya

	TZ	QZ	5Z	TZ/QZ	QZ/5Z
SCF	–106.8	–107.3	–107.2	–107.5	–107.1
SD	79.3	82.2	82.6	84.3	83.0
(T1)	14.5	15.2	15.3	15.7	15.4
Ecorr	93.8	97.4	97.8	100.0	98.3
Total	–13.0	–9.9	–9.3	–7.5	–8.8

All energies are in kcal/mol.

Figure 2

SCF, SD, and SD(T) contributions to ΔEHL for [Fe(H2O)6]2+ relative to the QZ/5Z results depending on the basis set hierarchy.

Figure 4

SCF, SD, and SD(T) contributions to ΔEHL for [Fe(NCH)6]2+ relative to the QZ/5Z results depending on the basis set hierarchy.

SCF, SD, and SD(T) contributions to ΔEHL for [Fe(NH3)6]2+ relative to the QZ/5Z results depending on the basis set hierarchy.

As can be seen from the results, the correlation energy shows the expected slow convergence with respect to the basis set limit. Nonetheless, in all cases, the correlation energy part of the spin-state splitting at the 5Z level is within 1.5 kcal/mol of the respective QZ/5Z extrapolated value, thus giving us some confidence that the QZ/5Z extrapolation results are sufficiently converged. We also note that the quality of the TZ/QZ extrapolation is of inferior quality if the QZ and 5Z results are already very close energetically, while the TZ and QZ are not. In such a case, the TZ/QZ value will significantly overshoot the true value as is apparent from the results of [Fe(NCH)6]2+.

From the results, it appears that, in general, it is indeed necessary to employ basis sets of 5Z quality for the correlation energy to attain results that are sufficiently converged to the basis set limit. While this is apparently not necessary for all cases, whether or not this applies cannot be determined beforehand.

As expected, the SCF results show a fast convergence to the basis set limit. In fact, the results at all basis set levels including extrapolation are within 1 kcal/mol of each other, most even within 0.5 kcal/mol. Consequently, basis set extrapolation of the SCF energy is not strictly necessary, and one could use the TZ or QZ results for this study. However, the orbitals of all basis set hierarchies are needed as input for the subsequent CC treatment. We note that in the case of [Fe(NH3)6]2+ and [Fe(NCH)6]2+ the SCF energy differences do not behave monotonously with the basis set level, thus leading to deviations in the extrapolations. However, these deviations are small enough to be of little consequence. Due to the Jahn–Teller splitting of the HS state, the 5T2g state is split, and the degeneracy of the former t2g set of orbitals is lifted. As this leads to energetically very close-lying occupied orbitals, it is important to ensure that all SCF calculations converge to the same state. To this end, we examined the highest doubly occupied orbitals of each calculation and verified that they are indeed of the same nature (see respective figures in the SI).

Comparison with Literature

Spin-state energetics have been intensely studied in the literature.[64−67] In this section, we compare the results presented here with the available literature. Naturally, only comparisons with publications containing the same compounds are meaningful. Also, many publications are concerned with vertical spin-state splittings, whereas we only investigated adiabatic ones which rules out these results for direct comparison. An overview of the available literature is given in Table . Judging the accuracy of the results is a difficult task given that the experimental singlet–quintet splittings are not known. Even if they were known, they likely would be significantly influenced by environment effects that were not included in the modeling. In fact, environment effects are known to be a major contributor to actual spin transition phenomena, and a lot of work has gone into trying to properly theoretically model these challenging collective effects.[68−70] Since definitive (Full-CI) benchmark data on electronic energies are not available, the best that appears to be presently possible is to rely on the commonly accepted assumption that CCSD(T) close to the basis set limit provides electronic energies within chemical accuracy of the unobtainable Full-CI/complete basis set results. At this point, this is an admittedly bold assumption that awaits further support from experiments and/or even higher levels of theory.

Table 9

Best Estimates for ΔEHL for All Compounds Considered Here as Taken from the Literature in Comparison to Present Resultsa

Compound	ΔEHL	Method	Basis set	Note	ref
[Fe(H2O)6]2+	–33.3	DLPNO-CCSD(T1)	cc-pwCVnZ-DK, cc-pVnZ-DK, (n = T,Q,5)	DKH2, tightPNO	This work
	–33.4	CCSD(T)	ANO:	DKH2	(57)
			Fe (7s,6p,5d,3f,2g,1h)
			O (4s,3p,2d,1f)
			H (2s)
	–35.7	CASPT2/SORCI	6-31G**/QZVP,TZVPP	Empirically corrected	(58)
	–43.3	CASPT2	ANO:	DKH2	(59)
			Fe (7s,6p,5d,3f,2g,1h)
			O (4s,3p,1d)
			H (1s)
	–46.6	CASPT2	ANO-RCC:	DKH2	(60)
			Fe (7s,6p,5d,3f,2g,1h)
			O (4s3p2s1f)
			H (3s1p)
	–60.0	DMC	Plain waves, Dirac–Fock pseudopotentials		(61)
	–41.0	DMC	cc-pVnT (n = D,T,Q)		(62)
[Fe(NH3)6]2+	–11.3	DLPNO-CCSD(T1)	cc-pwCVnZ-DK, cc-pVnZ-DK, (n = T,Q,5)	DKH2, tightPNO	This work
	–15.2	CCSD(T)	ANO:	DKH2	(57)
			Fe (7s,6p,5d,3f,2g,1h)
			N (4s,3p,2d,1f)
			H (2s)
	–28.6	CASPT2/SORCI	6-31G**/QZVP,TZVPP	Empirically corrected	(58)
	–22.6	CASPT2	ANO:	DKH2	(59)
			Fe (7s,6p,5d,3f,2g,1h)
			N 4s,3p,1d)
			H (1s)
	–21.3	CASPT2	ANO-RCC:	DKH2	(60)
			Fe (7s,6p,5d,3f,2g,1h)
			N (4s3p2s1f)
			H (3s1p)
	–35.7	DMC	Plain waves, Dirac–Fock pseudopotentials		(61)
	–28.4	DMC	cc-pVnT (n = D,T,Q)		(62)
[Fe(NCH)6]2+	–8.8	DLPNO-CCSD(T1)	cc-pwCVnZ-DK,cc-pVnZ-DK, (n = T,Q,5)	DKH2, tightPNO	This work
	–4.3	CCSD(T)	ANO:	DKH2	(57)
			Fe (7s,6p,5d,3f,2g,1h)
			C,N (4s,3p,2d,1f)
			H (2s)
	–2.0	CCSD(T)	Fe cc-pwCVnZ (n = T,Q,5)	Indirect DKH2	(56)
			C, N cc-pVTZ
			H cc-pVDZ
	–7.3	CASPT2	ANO:	DKH2	(59)
			Fe (7s,6p,5d,3f,2g,1h)
			N (4s,3p,1d)
			C (3s,2p,1d)
			H (1s)
	–31.6	DMC	Plain waves, Dirac–Fock pseudopotentials		(61)
	–19.6– −21.9	DMC	BFD/pVQZ (HF, DFT)		(63)
			ANO-RCC (CASSCF):
			Fe (8s,7p,6d,4f,3g,2h)
			N (4s,3p,2d,1f)
			C (3s,2p,1d)
			H (2s)
	–27.0	DMC	cc-pVnT (n = D,T,Q)		(62)

All energies are in kcal/mol.

In 2005, Fouqueau et al. published a study of [Fe(H2O)6]2+ and [Fe(NH3)6]2+ in which they intended to ascertain the performance of DFT for the adiabatic HS–LS energy splitting and used CASPT2, SORCI,[71] and ligand field theoretical calculations to provide benchmark results.[58] Using these methods they arrived at an estimate for ΔEHL of ca. −36 ± 1.5 and −29 ± 3 kcal/mol for [Fe(H2O)6]2+ and [Fe(NH3)6]2+, respectively. This differs from our own results by about 4 and 17 kcal/mol, respectively. This suggests that it is not equally challenging for both compounds to compute accurate spin-state splittings which is in all probability due to the different strength of the ligand field. It should, however, be noted that in this study neither relativistic effects nor correlation of the outer core electrons was included which is known to produce significant errors in the results. Also, the CASPT2 results were obtained with a 6-31G** Pople basis set which is by far too small to produce accurate results. Even the basis set for the SORCI calculations (TZVPP for ligand atoms, QZVP for Fe) must be considered too small in the light of the results in this work. As evident from Figures –4, triple-ζ basis sets are far off of the basis set limit, and even a quadruple-ζ quality basis on all atoms is not always sufficient. Similar arguments hold for all the literature studies discussed here.

One year later, Pierloot and Vancoillie also computed the adiabatic HS–LS energy splitting of [Fe(H2O)6]2+ and [Fe(NH3)6]2+ with CASPT2.[60] With a CAS(12,10) reference wave function, they arrived at an estimate of −46.6 and −31.4 kcal/mol for ΔEHL for the aquo and ammine complexes, respectively. Furthermore, a combination of large atomic natural orbital (ANO) basis sets containing up to h-functions was used in conjunction with scalar relativistic effects. The differences from our results amount to approximately 15 and 11 kcal/mol, respectively, and can probably most consistently be explained with the known tendency of CASPT2 to overstabilize the HS state.[53] Some part of the discrepancies may also be attributed to the basis set used (as mentioned above) and different methods used for geometry optimizations as well as the employed symmetry restrictions. We also note that the discrepancy to our results is of a similar magnitude for both compounds, in contrast to the study of Fouqueau et al.[58]

The de Graaf group published two studies in 2009 and 2010, in which all the model complexes considered in this work were investigated in terms of their bonding situation.[57,59] Furthermore, the adiabatic spin-state splittings were calculated with both CASPT2 and canonical CCSD(T) with a methodology similar to that of Pierloot and Vancoillie.[60] Notably, the CCSD(T) results match our own rather well, considering that the basis set cannot be saturated in the canonical calculations. The largest discrepancy amounts to ca. 5 kcal/mol in the case of [Fe(NCH)6]2+ which can be attributed to the basis set incompleteness error. The CASPT2 numbers from this study, however, show basically the same results as the ones of Pierloot and Vancoillie.[60] They demonstrate the tendency of CASPT2 to overstabilize the HS state, which is expected since important higher-order correlation effects are missing in second-order perturbation treatments. On the basis of our experience and a recent study by Pierloot et al., NEVPT2 is expected to show similar shortcomings.[53] In this study, CASPT2 is recommended as the method of choice for spin-state splittings since it showed a lower deviation from CCSD(T) results. However, the comparison of these studies with past and present results indicate that neither multiconfigurational second-order perturbation treatment is able to reliably produce accurate spin-state energetics as in many cases the errors (relative to CCDS(T)) are quite large, particularly since the goal is chemical accuracy. It was found that this is mainly due to a poor description of the correlation of the 3s and 3p outer core orbitals. In a further study by the same group, it was pointed out that a balanced basis set is necessary for the metal center and the first coordination sphere to arrive at reliable results.[72] It was thus concluded that a more accurate method than CASPT2 or NEVPT2 must be used with a sufficiently large basis set for this challenge.

In 2012, Lawson Daku et al. conducted a detailed study of [Fe(NCH)6]2+ with CCSD(T) to generate reliable reference data for a spin crossover model complex.[56] They obtained the best guess of −2.04 kcal/mol for ΔEHL which still differs significantly from our value of −9.15 kcal/mol. This is due to a number of subtle differences which nonetheless add up to a considerable amount. First, the geometry was optimized in the Oh and D2h symmetries, in contrast to this work. Also, different optimization methods were used. For the electronic energies, the ROHF-CCSD(T) method was used in conjunction with the correlation-consistent cc-pVnZ/cc-pwCVnZ basis set. Contrary to our approach, cc-pVDZ was used for H atoms, cc-pVTZ for N and C, and cc-pwCVnZ for Fe with n = T,Q,5. The energies were extrapolated, and an elaborate correction scheme was employed for the fact that the extrapolation did not include the ligand atoms. Furthermore, the scalar relativistic correction was computed separately with a UHF reference determinant which was shown to introduce little or no error. It is difficult to quantify the amount of discrepancy with regard to our calculations coming from each methodical difference, but we surmise that, in sum, they add up to the difference of ca. 7 kcal/mol that was found. We speculate that the largest part of the difference to our work stems from basis set incompleteness errors which were not canceled by the employed correction scheme.

There are several studies using diffusion Monte Carlo simulations in an attempt to provide accurate reference data for spin-state splittings on model systems. Droghetti et al. arrived at best guess values of ca. −60, −36, and −32 kcal/mol for [Fe(H2O)6]2+, [Fe(NH3)6]2+, and [Fe(NCH)6]2+, respectively, with the uncertainty of ca. 3.5 kcal/mol due to methodical variations.[61] Fumanal and coworkers provided a value of 20.7 ± 1.2 kcal/mol for [Fe(NCH)6]2+ a few years later.[63] Recently, Song et al. used DMC values as a reference in a study of density-corrected DFT on small Fe(II) complexes.[62] These provided values of −41.0, −28.4, and −27.0 kcal/mol for [Fe(H2O)6]2+, [Fe(NH3)6]2+, and [Fe(NCH)6]2+, respectively. We note that all of the DMC results for a given complex shown here do not only deviate significantly from the results presented in this Article but also from each other. The differences far exceed the quoted presumed error bars, showing that the latter cannot be valid. We conclude that the spin-state splittings calculated with DMC presented in the literature do not provide a consistent picture. Hence, it seems doubtful that the goal of providing reliable reference values was achieved.

Radón et al. conducted extensive and valuable benchmarks on iron aquo complexes but only computed vertical excitation energies which, unfortunately, prevents a direct comparison with our results.[68,73]

A comparison to the results of Feldt and coworkers is, strictly speaking, not possible since the systems in question are quite different in nature.[54] However, we would like to make a few remarks. In this work, only ground state structures are used, whereas Feldt and coworkers investigated the C–H activation of methane by an iron-oxo complex and studied the appropriate potential energy surface. Also, we would like to point out that the skeptical verdict on DLPNO-CCSD(T) of the aforementioned publication is mostly due to the shortcomings of the (T0) approximation, which was investigated here in detail. Our results demonstrate that these technical shortcomings have been remedied and that, with sufficient care, results within 1 kcal/mol of the canonical results can be obtained with DLPNO-CCSD(T1).

More generally, it should also be noted that the amount of static correlation in the ground state structures examined here is well within the capabilities of (DLPNO-)CCSD(T) which does not necessarily hold true for all transition metal complexes or transition states. Also, as has been pointed out,[54] great care must be taken with the reference orbitals, otherwise the computed results may well become meaningless. Thus, caution is always advisible.

In summary of this part, the introduction of elaborate correction schemes for basis set incompleteness appears to have a larger influence on the final results than expected. Hence, in our opinion, it is preferable to apply the same basis set hierarchy to all atoms in the system (or at least the first coordination sphere), provided that this is computationally feasible. In this respect, the DLPNO methodology opens the door for much larger systems than previously possible. This can then be followed up by standard basis set extrapolation techniques. Furthermore, the DLPNO-CCSD(T) method appears to be more successful than CASPT2 in reproducing canonical CCSD(T) results for the systems investigated here.

Comparison with Density Functional Theory

As a comparison to a much used method, we compared the results of DLPNO-CCSD(T) with several popular density functional approximations which are often chosen for problems of the kind studied in this work. We chose the [Fe(NH3)6]2+ compound for this endeavor. The results are shown in Figure . Evidently, the DFT results display a large scattering with some of the GGA functionals even predicting the wrong sign of the spin-state splitting. Moreover, no clear tendency can be seen regarding the eligibility of any particular DFA hierarchy either. However, we note that the M06L, PW6B95, and ωB97X-D3 functionals yield results almost identical to DLPNO-CCSD(T). For M06L, this can be ascribed to fortuitous error cancellation, whereas PW6B95 also yields good results for the other two compounds investigated in this paper (cf. Tables S16 and S17). Whether this, too, is due to error cancellation cannot reliably be judged on the basis of three compounds alone.

Figure 5

Comparison of ΔEHL values for [Fe(NH3)6]2+ computed with various density functional approximations with DLPNO-CCSD(T). The latter is represented by the solid line. The dotted lines show deviations of 1 kcal/mol from DLPNO-CCSD(T).

It also becomes apparent that even double hybrid functionals which include a percentage of MP2-type correlation energy are not necessarily an improvement over other, computationally less demanding ones.

All of these findings, of course, had to be expected since it has long been known that GGA functionals tend to overstabilize the LS state due to overlocalizing of electrons, whereas hybrid functionals give results which are roughly linearly dependent on the amount of Hartree–Fock exchange included. This has been shown by Reiher and coworkers who used this to reparametrize known functionals to match experimental findings.[6] MP2, on the other hand, and by extension double hybrid functionals, are known to be vulnerable to substantial amounts of static correlation which leads to inaccurate results for many transition metal compounds such as the ones studied here.

We note that over the years many authors performed sophisticated ab initio studies on SCO compounds or model complexes in order to use the results as benchmarks for DFT. While this is certainly a sensible approach for many areas of chemistry, in our opinion, this has significant shortcomings in the case spin-state splittings. It has been demonstrated many times that there is no single density functional approximation capable of describing a range of transition metals in different oxidation states and coordination environments correctly. This is true even for, say, all of the Fe(II) compounds, let alone all SCO compounds, which would include Co(III), Fe(III), or others. The optimal amount of Hartree–Fock exchange will vary significantly for all these situations, making it nearly impossible to arrive at a reliable DFT-based protocol. Considering the broad field of spin-state energetics in the field of bioinorganic chemistry, for example, in the context of multistate reactivity, reliable wave function-based protocols appear to offer new opportunities not previously available.

Conclusion

In summary, we were able to show that the results of DLPNO-CCSD(T1) are within the limits of chemical accuracy compared to canonical CCSD(T), at least for the compounds investigated here.

Furthermore, it was determined that both scalar relativistic effects in the form of the second-order DKH Hamiltonian and the use of the TightPNO thresholds as implemented in ORCA contribute non-negligibly (ca. 3 kcal/mol each) to the final results and thus should always be included. Also, as has been shown previously, and as is the default in ORCA, it is necessary to include the correlation brought in by the outer core orbital of the metal center.

It appears from our results that the basis sets are reasonably converged to the complete basis limit when performing QZ/5Z extrapolation which was only feasible with the DLPNO approximation. TZ/QZ extrapolation was shown to not give satisfactory results in all cases. Since this constitutes a serious challenge for an application to real systems, the employment of composite basis sets on different parts of the molecule at hand or explicit correlation methods may prove worthwhile in the future.

For the SCF part, the dependence on the basis set was found, as expected, to be less severe. However, given that the SCF calculations are hardly the computational bottleneck in a given study, we recommend using the values of the highest available basis set hierarchy in favor of basis set extrapolation.

It should be noted that applying (T1) to all problems of interest henceforth is not necessary as it has been demonstrated previously that the semicanonical triples (T0) yield accurate results for relative energies in most cases, e.g., for chemical reactions.[23,25] Nonetheless, in selected cases, such as in systems with small HOMO–LUMO gaps, the use of (T1) may become mandatory. It should be kept in mind that (T1) is not only more accurate but also more expensive than (T0) with the impact of the additional computing time decreasing with basis set size. Consequently, if very high accuracy is desired and large basis sets are used, (T1) is barely more expensive than (T0). In general, in case of doubt, we recommend testing the difference between both approaches in a controlled manner.

A thorough comparison with literature results once more demonstrates the difficult nature of the problem at hand as it became apparent that the methods used in the past do not produce sufficiently accurate results. This is due to either systematic method errors, as is the case for, for example, CASPT2 (and NEVPT2), or the fact that convergence to the basis set limit is not feasible due to the steep scaling with the system size (which is the case for CCSD(T)). We are pleased to note that the DLPNO results presented here agree reasonably well with the literature on canonical CCSD(T) considering several differences in basis sets, methodology, and geometries.

A comparison with a range of popular density functionals shows a large scattering of the results of almost 30 kcal/mol, depending on the functionals used, and underlines that spin-state splittings are an exceedingly difficult problem for DFT.

With the present work, it was shown that the calculation of accurate electronic energy contributions to spin-state splittings with respect to canonical CCSD(T) results under consideration of all significant physical effects including extrapolation to the basis set limit is indeed possible for the systems investigated here. While it is clear that more work is necessary to arrive at streamlined computational protocols that are usable in a blackbox fashion, we hope to inspire new confidence that the accurate calculation of spin-state splittings on the basis of correlated wave function methods is indeed feasible.

Table 3

Correlation Energy Contributions and Spin-State Splitting of HS and LS States of [Fe(NH3)6]2+a

	HS				LS
	Ecan	EDLPNO	Δc/mEh	Ecan/EDLPNO	Ecan	EDLPNO	Δc/mEh	Ecan/EDLPNO	ΔEcan/ kcal mol–1	ΔEDLPNO/ kcal mol–1	ΔΔEd/ kcal mol–1
SD	–1.85684	–1.85283	–4.0	0.998	–1.93843	–1.93476	–3.7	0.998	51.2	51.4	–0.2
(T0)	–	–0.04207	–3.9	0.914	–	–0.05208	–7.0	0.881	–	6.3	1.9d
(T1)	–0.04601	–0.04394	–2.1	0.955	–0.05909	–0.05683	–2.3	0.962	8.2	8.1	0.1
Ecorrb	–1.90285	–1.89677	–6.1	0.997	–1.99752	–1.99159	–5.9	0.997	59.4	59.5	–0.1

All energies in Ha unless stated otherwise. TightPNO thresholds were used for all DLPNO calculations. Ecan and EDLPNO are energies computed with CCSD(T) and DLPNO-CCSD(T), respectively.

Calculated as the sum of the SD and (T1) values.

Calculated as Ecan – EDLPNO.

Calculated as ΔEcan – ΔEDLPNO.

Calculated as (T1)(ΔEcan) – (T0)(ΔEDLPNO).

Table 7

Contributions to ΔEHL of [Fe(NH3)6]2+ Depending on the Basis Set Hierarchya

	TZ	QZ	5Z	TZ/QZ	QZ/5Z
SCF	–75.4	–75.8	–75.7	–76.0	–75.7
SD	51.2	53.5	54.4	55.1	55.4
(T1)	8.2	8.7	8.8	9.0	9.0
Ecorr	59.4	62.1	63.2	64.2	64.4
Total	–16.0	–13.6	–12.5	–11.8	–11.3

All energies are in kcal/mol.