Performance of density functional theory and orbital-optimised second-order perturbation theory methods for geometries and singlet–triplet state splittings of aryl-carbenes

Abstract

Carbenes are challenging molecular species for quantum chemistry because of the energetic proximity of their singlet and triplet spin states and the sensitive dependence of spin-state energetics on the geometry of the carbene site. Here we use an extended set of aryl-carbenes to evaluate the performance of density functional theory (DFT) approximations as well as of wave function based perturbation theory approaches (orbital-optimised perturbation theory methods OO-MP2 and OO-SCS-MP2) against reference coupled cluster calculations with singles, doubles and perturbative triples conducted with the aid of the domain-based local pair natural orbitals approach, DLPNO-CCSD(T). In addition to the expected functional dependence, our results document a remarkable discordance in the performance of DFT methods in the sense that the functionals that yield the best geometries do not coincide with those that provide the best spin-state energetics. Analysis of the results allows us to propose a series of methods that are expected to perform reliably within certain confidence limits for the title systems. Additionally, methodological issues regarding the reference singlet–triplet gaps obtained by the DLPNO-CCSD(T) approach are discussed.

GRAPHICAL ABSTRACT

KEYWORDS:

• Carbenes
• spin states
• coupled cluster theory
• density functional theory
• aryl-carbenes

1. Introduction

Carbenes are highly reactive organic molecules where a neutral carbon atom has two electrons less than an octet structure [1,2]. In the language of valence bond theory, carbene centres adopt sp² hybridisation and hence are bent. Their electronic structure and spin state multiplicity can be described in a simplified manner by assuming two electrons to be distributed in two nearly degenerate p and sp² orbitals (Figure 1). Formation of parallel spins in this configuration results in a high-spin triplet state (S = 1), while pairing of the orbitals in the sp² orbital leads to a singlet state (S = 0). The triplet and singlet carbenes are significantly different species both structurally and electronically. Most carbenes adopt a triplet ground state, but singlet ground state carbenes can also form, or spin-state equilibria established, through interaction with proper substituents and steric restrictions [3–6]. The significant difference in electron distribution between the two spin states results in distinctly different chemical properties and reactivity. For instance, the often higher-lying singlet state participates in reactions rather than the lower-lying triplet state. Therefore, knowledge of accurate singlet–triplet energy gaps are of high importance in carbene chemistry.

Figure 1. Schematic distribution of electrons in the triplet (left) and singlet (right) spin states of the simplest carbene, methylene.

In recent studies, it was shown that localised coupled cluster methods, specifically the domain-based pair natural orbital (DLPNO) coupled-cluster theory [7–9], can be used to assess the energetics of systems as large as a small protein and even larger. This approach can therefore bring the accuracy of the ‘gold standard’ CCSD(T) calculations to much larger systems than those accessible with conventional CCSD(T) calculations. Carbenes are one of the chemically interesting systems where canonical CCSD(T) calculations are absolutely required [10] for achieving chemical accuracy in singlet–triplet gaps (ca. 0.4 kcal/mol in case of CH₂ [11,12]). DLPNO-CCSD(T) can be used to approximate closely the results of canonical coupled-cluster when the latter is not applicable. For instance, we have recently shown that DLPNO-CCSD(T) reproduces the canonical singlet–triplet splitting for a number of aryl-carbenes (the AC-12 set) with a root-mean-squared (RMS) error of only 0.4 kcal/mol [10]. Furthermore, DLPNO calculations can be used to decompose the reference and correlation energy of bonded and non-bonded carbenes [13]. This is achieved by the application of the local energy decomposition scheme (LED) for restricted and unrestricted DLPNO-CCSD [14,15]. However, lower-cost wave function methods or even more practically DFT calculations, despite being characterised by inconsistent behaviour [10,16], should in principle remain a practical tool for geometry optimizations and at least for exploratory calculations of singlet–triplet gaps, provided the performance of different functionals is carefully assessed and evaluated.

In this work we take advantage of the DLPNO approach in order to gauge the performance of a number of widely used density functionals and of orbital-optimised second-order Møller–Plesset perturbation theory (OO-MP2) and its spin-component-scaled variant (OO-SCS-MP2) [17–22]. For this purpose, we have set up a benchmark to examine the performance of the tested methods both in predicting geometries and singlet–triplet gaps. 18 different aryl-carbenes were selected (AC-18 set, Figure 2) slightly extending our previous AC-12 benchmark set with larger molecules. This set is a mixture of singlet and triplet ground state aryl- and diaryl-carbenes with different electron-donating and electron-withdrawing substituents stabilising different spin states. Using the DLPNO-based coupled cluster calculations, we first analyse the quality of geometry optimisation provided by the methods evaluated in this article. Although geometry optimizations at the DLPNO-CCSD(T) using analytic gradients are not yet feasible, the total DLPNO-CCSD(T) energy at DFT and MP2 optimised geometries is an adequate measure of how closely the DFT optimisation may approach the DLPNO-CCSD(T) minimum—provided of course that the same local minimum is targeted. Hence, we take a lower total DLPNO-CCSD(T) energy as a criterion for assuming that a given DFT geometry is closer to the (unknown) DLPNO-CCSD(T) minimum than DFT derived structures that lead to higher DLPNO-CCSD(T) energies. Furthermore, critical carbene geometrical features such as carbene bond, angle and dihedral angles are analysed both chemically, in terms of substituent dependency, and computationally, in terms of method dependency. In a second step, we focus on the singlet–triplet gaps. First, a discussion about the approximations used in DLPNO coupled cluster is provided in view of the recent development of the iterative perturbative triples excitation correction that is more accurate than the previously employed semi-canonical approximation in the case of small HOMO–LUMO gaps. Finally, the accuracy of singlet–triplet gap prediction by DFT functionals is evaluated both at the structures optimised by each DFT functional itself (relaxed) and for specific reference geometries (unrelaxed) .

Figure 2. The AC-18 test set of aryl-carbenes used in this study.

2. Computational methods

Geometries were optimised using a set of DFT functionals (listed in Table 1), orbital-optimised second-order Møller–Plesset perturbation theory (OO-MP2), derived by minimisation of Hylleraas functional, and its spin-component-scaled implementation OO-SCS-MP2 [17–19]. Previous studies suggested that conventional UHF MP2 is unreliable for carbenes due to spin contamination of the triplet state [10], therefore only the orbital-optimised variant is discussed here. Both OO-MP2 and OO-SCS-MP2 methods are employed with the resolution of the identity approximation (RI-MP2). SCS-MP2 is a variant of MP2 where the correlation energy is partitioned into same-spin and opposite-spin components that are scaled separately to achieve better performance [23]. DFT calculations benefitting from the D3 dispersion correction [24] with Becke–Johnson damping [25] are denoted by ‘D3BJ’ [25,26]. All geometry optimizations have used the def2-TZVPP [27] basis set and the resolution of identity approximation (RI) for Coulomb integrals together with the def2/J [28] auxiliary basis sets. The chain-of-spheres approximation (COSX) to exact exchange was used in geometry optimizations with hybrid functionals [29]. All integration grids were increased substantially with respect to the default values (Grid6 and GridX8 in ORCA convention). The auxiliary basis sets used set in MP2 calculations was the def2-TZVPP/C [30]. Final energies for DFT methods were computed using the def2-QZVPP basis sets without the RI approximation in order to eliminate any numerical errors. SCF convergence was set to ‘VeryTight’ for all calculations reported in the present work.

Table 1. DFT functionals used in this work.

Type	Functional	Ref.
GGA	BP86	[33,34]
	BLYP	[33,35]
	PBE	[36]
	REVPBE	[36–38]
	RPBE	[39]
Hybrid	B3LYP	[35,40]
	O3LYP	[35,41]
	X3LYP	[42]
	PBE0	[43]
	ωB97-X	[44]
Meta GGA	TPSS	[45]
	M06-L	[46]
Meta hybrid	TPSSh	[47]
	M06-2X	[48]
Double hybrid	B2PLYP	[49]
	B2GP-PLYP	[50]
	mPW2PLYP	[51,52]

DLPNO coupled-cluster calculations for singlet and triplet states [7–9] were performed either using tight or normal PNO thresholds. In the initial evaluation of geometries (Section 3.1) the relevant DLPNO calculations employed normal (default) PNO thresholds and made use of the (T0) triples corrections. Final evaluation of energies (Section 3.2) and all reference values for the singlet–triplet gaps reported in this work used tight PNO thresholds and (T1) triples. The correlation consistent basis sets cc-pVTZ and cc-pVQZ were used for the DLPNO coupled-cluster calculations [31], in combination with their respective auxiliary basis sets. In DLPNO calculations the keyword ‘Randomize 0’ was used in order to ensure the consistency of localised orbitals in different calculations. We define the S-T gap as ΔE_ST= E_T − E_S, therefore negative values of splitting, when mentioned, mean a triplet ground state. All calculations were performed with ORCA [32].

3. Results and discussion

3.1. Evaluation of optimised geometries

The quality of geometries optimised using OO-MP2, OO-SCS-MP2, and various DFT functionals are evaluated by calculating the DLPNO-CCSD(T0)/cc-pVTZ energies of the optimised geometries. This approach is based on the assumption that the optimum structure will have the lowest DLPNO-CCSD(T0) energy. Considering the high number of optimised geometries with various methods, one can assume that one of the optimised geometries will be sufficiently close to the optimal DLPNO-CCSD(T0)/cc-pVTZ structure to give a good qualitative overall picture. Therefore, those structures are selected as the reference structures for subsequent evaluations. Each spin state is treated separately in terms of which method performs best, and hence which method provides the reference structure. It is noted that for the initial evaluation of geometries the semi-canonical triples were used in the coupled-cluster calculations. As shown in Table 2, TPSS and TPSSh provide five of the reference structures, while B2GP-PLYP-D3BJ, X3LYP and B2PLYP-D3 each provide one reference structure. The rest of the lowest-energy structures are all provided by OO-SCS-MP2.

Table 2. Methods that provide the reference structure for each spin state of the carbenes considered in the present study, as evaluated using DLPNO-CCSD(T0) energies.

Carbene	Singlet state	Triplet state
1	OO-SCS-MP2	OO-SCS-MP2
2	B2PLYP-D3	TPSSh-D3BJ
3	OO-SCS-MP2	OO-SCS-MP2
4	OO-SCS-MP2	OO-SCS-MP2
5	OO-SCS-MP2	OO-SCS-MP2
6	TPSSh-D3BJ	OO-SCS-MP2
7	X3LYP	TPSSh-D3BJ
8	OO-SCS-MP2	OO-SCS-MP2
9	OO-SCS-MP2	OO-SCS-MP2
10	OO-SCS-MP2	TPSS-D3BJ
11	TPSS-D3BJ	OO-SCS-MP2
12	OO-SCS-MP2	OO-SCS-MP2
13	OO-SCS-MP2	OO-SCS-MP2
14	OO-SCS-MP2	OO-SCS-MP2
15	OO-SCS-MP2	OO-SCS-MP2
16	OO-SCS-MP2	OO-SCS-MP2
17	B2GP-PLYP-D3BJ	OO-SCS-MP2
18	OO-SCS-MP2	OO-SCS-MP2

Table 3 shows the absolute average energy difference between different structures to the minimum energy found among all available structures for a specific spin state of a particular species, along with the maximum deviation observed and the species for which it was observed. We note that the RMS deviation in going from NormalPNO to TightPNO thresholds is about 0.13 kcal/mol for the AC-12 benchmark subset that was reported previously [10].

Table 3. Average and maximum absolute errors (kcal/mol) for each evaluated method, defined as DLPNO-CCSD(T0)/cc-pVTZ energy differences between the AC-18 carbene geometries optimized by a given method and the reference geometry that yields the lowest DLPNO-CCSD(T0) energy for each carbene in each spin state.

		Singlet			Triplet
	Average	Max	Max No.	Average	Max	Max No.
OO-SCS-MP2	0.22	3.47	17	0.03	0.51	2
OO-MP2	0.53	2.03	17	0.24	0.99	2
mPW2PLYP	0.39	0.96	16	0.49	1.02	16
B2GP-PLYP	0.37	0.86	16	0.75	1.32	14
B2GP-PLYP-D3BJ	0.39	0.89	16	0.79	1.42	14
B2PLYP-D3	0.23	0.72	16	0.32	0.76	16
B3LYP	0.37	1.22	14	0.33	0.92	14
B3LYP-D3BJ	0.34	1.02	14	0.34	0.87	14
X3LYP	0.37	1.20	14	0.34	0.93	14
O3LYP	1.37	3.05	17	1.06	2.04	14
TPSSh-D3BJ	0.27	1.12	17	0.12	0.40	16
wB97X-D3	0.46	0.90	14	0.48	1.11	14
B3PW91-D3BJ	0.53	1.18	14	0.47	0.99	16
PBE0-D3BJ	0.60	1.20	14	0.57	1.05	16
M06-2X	0.70	2.03	17	0.64	1.57	15
PBE-D3BJ	1.07	2.58	17	0.35	0.95	16
RPBE-D3BJ	0.80	2.19	2	0.31	0.72	18
REVPBE-D3BJ	1.00	2.36	2	0.44	0.89	18
TPSS-D3BJ	1.03	6.37	2	0.12	0.42	16
BP86-D3BJ	0.88	2.27	2	0.34	0.89	16
M06-L	0.89	1.99	16	0.67	1.39	16
BLYP-D3BJ	0.97	2.46	14	0.54	1.08	14

Inspection of Table 3 reveals that the average errors are mostly below 1 kcal/mol. Maximum errors are, with very few exceptions, only slightly higher (ca. 2–3 kcal/mol). O3LYP produces structures with the largest deviation from the optimal minimum-energy DLPNO-CCSD(T0) structure for both singlet and triplet states. Non-hybrid functionals typically produce worse results for the singlet carbene structure compared to the triplet. However, in the case of the triplet state, the results vary significantly. Interestingly, both TPSS and TPSSh produce very good triplet geometries as judged from these energy differences, however, for the singlet geometries, TPSSh clearly surpasses the non-hybrid TPSS. Most hybrid functionals provide similar accuracies, with the exception of O3LYP as mentioned above. TPSSh provides the second-best results overall, which is important because it surpasses the performance of double-hybrid functionals at lower computational cost and better scaling. As for TPSSh, the accuracy of the geometries provided for both spin states are similar for most hybrid functionals. This can be a significant advantage because it can result in beneficial error cancelation in the calculation of S-T gaps. With regards to double-hybrid functionals, the original B2PLYP outperforms the other double hybrids included in the present work and provides very balanced results for both spin states with respect to DLPNO-CCSD(T0) energies. Despite that, double-hybrid functionals do not specifically improve upon the conventional hybrid functionals on an average sense, even though the maximum errors are reduced. Taking note of this fact and considering the higher cost of double hybrid calculations, one may conclude that application of double hybrid functionals for geometry optimizations of carbenes may not offer any advantages when it comes to creating reference geometries for the most reliable subsequent evaluation of accurate energetics.

At this point, it is useful to discuss the effect of the dispersion correction term in geometry optimisation. We focus on two functionals, B2GP-PLYP and B3LYP, with which we have optimised all structures in the AC-18 set with and without the D3BJ dispersion correction. The difference in both cases is insignificant and the small changes effected by the empirical dispersion correction do not always improve the structures. Although one may not necessarily expect a considerable advantage from including dispersion effects for the carbene molecules of this study, the above results suggest that empirical corrections such as the D3BJ approach employed here are not beneficial for medium-sized carbene systems on an average sense.

Having discussed the performance of DFT functionals, we now turn to the orbital-optimised OO-MP2 and OO-SCS-MP2 methods. The orbital-optimised variant of the MP2 method is selected as it is remedying the spin contamination of the reference wavefunction [10]. Spin component scaling empirically counters the bias toward same spin electrons of the SCF calculations by adjusting the contribution of same-spin and opposite-spin electron correlation. The overestimation of non-bonding interactions by MP2 is also suggested to be remedied by the SCS modification [24]. For the present comparison against DLPNO-CCSD(T0) energies, OO-MP2 provides relatively good geometries especially for the triplet states, however OO-SCS-MP2 improves on the OO-MP2 structures in all cases but one, thus providing the optimal reference structures. Therefore, between the two options we can recommend uniformly the application of the OO-SCS-MP2 variant. It has been noted that OO-MP2 may tend to overestimate correlation effects with concomitant effects on bond lengths [17,53], an effect that could be damped in the SCS variant consistent with the results obtained here. Although not tested in the present work, it is also worth mentioning the regularised OO-MP2 approach of Head-Gordon and co-workers [54,55], which can successfully counter divergences that may occur during the optimisation procedure for OO-MP2 when the frontier orbitals approach degeneracy.

Despite the above findings, considering the relatively high cost of orbital-optimised MP2 compared to any DFT alternative, it is recognised that depending on the type of application it may be questionable whether the higher computational cost always justifies the increase in accuracy, particularly given the good performance of standard hybrid functionals. Therefore, as an alternative to OO-SCS-MP2, the recommended DFT approach for obtaining reliable structures is TPSSh. Figure 3 provides a closer look at the quality of the structures provided by TPSSh per species of the AC-18 set. The most problematic cases in terms of unfavourable DLPNO-CCSD(T0) energies seem to be species 14, 17, and 18. Interestingly, these are all species which feature a distinct deviation from planarity between the two carbene substituents. Therefore, this issue is likely related to the correlated prediction of correct carbene bond distance, carbene angle, and dihedral angle between the carbene substituents. In the following, we take a closer look individually at these structural parameters.

Figure 3. Energy errors (kcal/mol) of TPSSh-optimized geometries for the carbenes of the AC-18 set. The energy error is defined as the DLPNO-CCSD(T0) energy difference between the TPSSh geometry and the optimized geometry that provides the lowest DLPNO-CCSD(T0) energy for each carbene in each spin state.

3.1.1. Carbene bond length

The variation in the length of the C–C bond between the carbene centre and the phenyl group in substituted carbenes may be qualitatively connected to Bent’s rule [56], according to which the σ bond length decreases as s orbital character increases. This would imply that the C–C bonds involving the carbene centre might be shorter in the singlet state because this features a smaller carbene angle and, hence, larger s character in the carbene sp² orbital. However, the opposite is observed, with the triplet state having shorter C–C bonds. The longer bond length in the singlet state is presumably due to the additional electron repulsion between the σ-bond substituents and the doubly occupied σ-sp² orbital [57]. It is also noted that chemical substitution of stabilising groups for each spin state reduces the bond length (electron donors for singlet and electron withdrawal groups for the triplet carbenes).

Deviations in the optimised bond lengths of each method in comparison to the overall optimal reference structure are shown in Figure 4. All DFT functionals and wave function-based methods predict bond lengths within 0.03 Å of the optimal, the only exception being the TPSS-optimised carbene–carbon bond in species 2 and 11. This seems to be an anomaly as the TPSS singlet wavefunction is spin-contaminated (specifically, <S²> = 0.84 and 0.79 for carbenes 2 and 11, respectively) because it converges to a broken-symmetry state (‘open shell’ singlet). This issue occurs also with TPSSh but to a much smaller extent. No other functionals showed such behaviour for these two or any other carbenes.

Figure 4. Average and maximum errors of tested methods in carbene–carbon bond lengths (in Å) against the minimum DLPNO-CCSD(T0) energy reference structures.

In general, the deviations among methods are not pronounced. Most methods slightly underestimate the C–C bond length compared to the reference geometry. In the singlet state, the results present a wider scatter between different GGA, hybrid, and double hybrid methods compared to the triplet state. Most hybrid and GGA functionals provide bond lengths within 0.01 Å for the singlet state, with the exception of O3LYP that yields the shortest bond lengths among all methods for both spin states. In the triplet state, most of the hybrid and GGA functionals show larger deviations from the reference geometry. However, TPSSh provides on average results that are balanced for both spin states and close in comparison to the more expensive wave function-based and double hybrid methods. The Minnesota hybrid M06-2X is the opposite extreme compared to O3LYP, as it provides the longest bond lengths for the singlet, in this respect better-resembling OO-SCS-MP2. Double hybrid functionals and M06-2X also yield in general longer bond lengths than the GGA and hybrid functionals in the triplet state, which is closer to the MP2-based methods.

3.1.2. Carbene angle

The carbene angle is the most critical structural parameter of carbenes. If a carbene was linear then both π and σ orbitals would be degenerate, resulting in a triplet ground state. As the carbene centre bends, the σ orbital gains s character and relaxes. A competition between electron–electron repulsion and orbital relaxation then determines the carbene angle. Therefore, singlet states have a more bent structure to maximise orbital relaxation and reduce degeneracy. Prediction of carbene angle presents a challenging problem for most methods in specific aryl-carbenes. In the singlet state, most species bonded with a single phenyl ring provide acceptable angles with a maximum error below 5 degrees. Addition of a second phenyl ring slightly decreases the accuracy. Considerable errors are observed for species 17, where errors up to 10 degrees are observed for the typically most accurate method OO-SCS-MP2. Table 4 lists deviations from the structure providing the lowest energy at coupled-cluster level for species 17 (in this case, B2GP-PLYP-D3BJ). As can be seen, orbital-optimised MP2 methods show significant divergence from the optimum angle as judged by minimal DLPNO-CCSD(T0) energies. Double hybrid methods on the other hand provide uniformly very good angles.

Table 4. Absolute errors in optimized carbene angles (in degrees) for species 17 in its singlet state and species 14 in its triplet state.

Method	17 singlet	14 triplet
OO-SCS-MP2	10.47	0.00
OO-MP2	12.51	0.74
mPW2PLYP	1.60	0.52
B2GP-PLYP	0.67	0.17
B2GP-PLYP-D3BJ	0.00	0.92
B2PLYP-D3	1.61	0.73
B3LYP	4.86	6.91
B3LYP-D3BJ	2.82	3.64
X3LYP	4.32	6.23
O3LYP	8.33	8.42
TPSSh-D3BJ	4.31	3.29
wB97X-D3	0.19	3.29
B3PW-D3BJ	2.99	2.77
PBE0-D3BJ	2.62	2.85
M06-2X	3.71	0.40
PBE-D3BJ	8.63	5.30
RPBE-D3BJ	5.70	2.97
REVPBE-D3BJ	6.96	4.04
TPSS-D3BJ	6.48	3.97
BP86-D3BJ	7.30	4.79
M06-L	3.39	4.44
BLYP-D3BJ	6.88	6.06

Calculation of the optimal angle of the triplet state appears to be even more challenging for most methods. Overall maximum errors are higher in the triplet state in comparison to the singlet state. The most problematic species in the triplet state is 14, where the maximum error can be up to 8 degrees as shown in Table 4. In this case, the only straightforward solution in achieving the optimum angle is the use of a double hybrid functional or orbital-optimised MP2. The results of these methods are approximated very well by M06-2X. It is noted that the addition of exact exchange does not necessarily improve the performance of GGA functionals, as judged by comparison of related pairs of hybrid and non-hybrid functionals.

Figure 5 shows the absolute average error of all methods used in this study. OO-SCS-MP2 provides the best overall carbene angles although it faces large errors in the specific case of the singlet state of species 17. As was evident from the detailed analysis of the worst-case scenarios of both spin states, the addition of exact exchange does not clearly improve the prediction of carbene angle. However, all methods benefiting from perturbation theory provide highly accurate carbene angles. Despite the occasional large errors in specific cases, absolute average errors can be considered small.

Figure 5. Absolute average and maximum errors in carbene angles (in degrees) of tested methods against the minimum DLPNO-CCSD(T0) energy reference structures.

We note that the carbene angle is not dependent on the empty p orbital, therefore, it is not considerably affected by the substitution of electron donor or acceptor groups. This is not the case for the out-of-plane distortions that will be discussed in the following.

3.1.3. Carbene dihedral angle

As mentioned above, dihedral angles involving the carbene site seem to be among the trickiest geometrical parameters to predict correctly. The out-of-plane conformation of the groups attached to the carbene centre is affected by steric and electronic effects. Principally the dihedral angle is determined by the steric demands of aryl groups; addition of electron donors and acceptors to the aromatic ring introduces further modifications. The phenyl rings in diaryl-carbenes are always forced out of a plane in both spin states. A few substituted aryl-carbene species also contribute in tilting the substituent at the carbene centre out of the plane of the other phenyl ring, such as with species 18. These twists are spin-state specific and may be suppressed or entirely absent in one of the two spin states. Specifically, for substituent effects, the dihedral angle appears to correlate with change in the effective electron donation or withdrawal to the phenyl ring (resonance adjustment) and when this happens the dihedral angle can be significantly different between the two spin states. The more in-plane conformations increase the magnitude of the effect of the ring substituent and the opposite occurs for more out-of-plane conformations [58].

As seen from the results in Figure 6, the absolute average error for all methods in the value of the dihedral angle is below 6 degrees. Local functionals are characterised by large deviations in the singlet state. This seems to be the primary reason for the considerably poor performance of GGA functionals for the singlet state. Surprisingly, double hybrid methods produce worse triplet dihedral angles than hybrid and GGA functionals. This might be related to MP2 problems with spin contamination. Orbital-optimised MP2 methods have rather large average errors in the singlet state, however this is mostly due to a large error in a single species (17), where the dihedral angle is about 25 degrees off. Otherwise, the orbital-optimised MP2-based methods provide accurate dihedral angles for all the other carbenes of the test set.

Figure 6. Absolute average and maximum errors in carbene dihedral angles (in degrees) of tested methods against the minimum DLPNO-CCSD(T0) energy reference structures.

3.1.4. Qualitative analysis of dihedral change by substituent

It is informative to take a closer look at specific cases for a qualitative discussion of variation of dihedral angle by the presence of electron donor and acceptor groups on the aryl rings. Focusing first on species 2, the triplet state is stabilised by the electron-withdrawing nitro group. The nitro group helps to restore the resonance of the phenyl ring that is reduced in the triplet state by delocalisation of the carbene p electron on the aryl. However, this effect is preferentially reduced in the singlet state so the phenyl ring can instead donate electron density to the empty p orbital on the carbene centre. Therefore, in the singlet state, the nitro group is slightly twisted out of plane to reduce the electron-withdrawing effect. Prediction of the twisted singlet state varies significantly among different methods and a twist of up to 20–30 degrees is predicted by most GGA methods. However, the reference structure (B2PLYP-D3) predicts a more modest out-of-plane angle of ca. 5 degrees. In the case of the para-substitution with an amine group (species 5) both amine hydrogens are bent out of the aryl ring plane in the same direction for both spin states. The bending is more pronounced in the singlet state. This can be rationalised by the different extent of conjugation of the nonbonding electron pair of the nitrogen with the aromatic system of the phenyl, depending on the occupation of the carbene p orbital and its coupling through the π-network of the benzene ring.

The methoxy group bound at the carbene centre in species 8 forms a planar structure in the singlet state but the methyl bends out of plane in the triplet state. This can be attributed to the interaction between the oxygen lone pairs with the unpaired electrons in the triplet carbene, whereas the singlet state is in-plane to maximise the π–π interaction between the oxygen lone pair and the empty p orbital of carbene. Prediction of this twist is overall quite accurate, with errors within a few degrees for most methods.

Substitution of electron donor and acceptor into the diphenylcarbene results to small effects, with a single exception of O^-. All cases follow the same trend mentioned before, the electron-donating group O^- significantly increases the dihedral angle in the triplet state, while in the singlet state the dihedral angle is significantly reduced to increase electron donation. Electron withdrawing groups CN and NH₂ have a much smaller structural effect. In the case of species 17 and 18, where the carbene centre is adjacent to a carbon atom bonded to one or two oxygen atoms, the triplet and singlet states have significantly different structures. In the singlet state, the structure is bent so the empty p orbital can interact with the oxygen electron pair. However, the triplet state adopts a planar structure.

3.2. Singlet–triplet gaps

3.2.1. Approximations relevant to the reference calculations

In DLPNO coupled-cluster calculations, the reference energy value (restricted HF for the singlet, quasi-restricted for the triplet) is the same as in the canonical counterpart. However, there is a deviation in the value of the correlation energy because the virtual space is localised in order to enable the lower scaling that characterises the DLPNO approach. We call this deviation ‘localisation error’. Table 5 lists deviations for the different excitation components of coupled-cluster theory and additionally compares deviations due to distinct approaches to the triples corrections. As can be seen, the error introduced by the DLPNO approximation is very small at the CCSD level for most cases, with the maximum deviation observed for species 17, the only case exceeding 0.5 kcal/mol. In the case of triple excitations, since DLPNO calculations use quasi-restricted orbitals, in canonical coupled-cluster theory one needs to either iteratively calculate the perturbative triples amplitudes or diagonalise the internal and virtual spaces separately. The latter case is not possible in case of localised orbitals. This means one either ignores the off-diagonal elements, i.e. the (T0) approximation, or iteratively calculates the perturbative triple excitation amplitudes, i.e. (T1) [59]. From Table 5, it is clear that the (T0) approximation may have a non-negligible effect on singlet–triplet gaps, however (T1) only slightly suffers from the localization error, to an extent usually less than the associated CCSD error. The reference values for this work were all obtained with perturbative (T1) triples.

Table 5. Deviations in singlet–triplet gaps (kcal/mol) for various comparisons of approaches relevant to the DLPNO-CCSD(T) calculations.

Carbene	DLPNO-CCSD vs canonical CCSD	(T0) vs canonical (T)	(T1) vs canonical (T)	DLPNO-CCSD(T1) vs canonical CCSD(T)
1	0.08	0.63	0.06	−0.14
2	0.07	0.79	0.02	0.10
3	0.02	0.18	0.06	−0.08
4	0.02	0.08	0.05	−0.07
5	0.04	0.23	0.09	−0.13
6	0.01	0.21	0.05	−0.04
7	0.15	0.45	0.12	−0.27
8	0.02	0.64	0.01	−0.04
9	0.03	0.62	0.02	−0.05
10	0.08	0.73	0.07	−0.15
11	0.08	0.72	0.05	−0.13
12	0.33	1.35	0.20	−0.53
13	0.30	0.32	0.02	0.29
14	0.42	1.58	0.21	−0.63
15	0.07	0.44	0.21	−0.14
16	0.11	0.33	0.16	−0.26
17	0.56	1.55	0.51	−1.07
18	0.20	0.68	0.19	−0.39

Notes: The present comparisons were made using the cc-pVDZ basis set, which allowed canonical CCSD(T) calculations to be performed for all carbenes of the test set. TightPNO setting were applied for the DLPNO calculations and a quasi-restricted orbital reference was employed in the canonical calculations to ensure comparability of results.

3.2.2. Singlet–Triplet gaps by DFT methods

To compare the accuracy of DFT methods for spin-state energetics, we have considered two different approaches that represent likely scenarios for practical applications: (a) singlet–triplet gaps based on structures optimised with each specific method, and (b) singlet–triplet gaps computed with single-point calculations using a common set of reference structures. In the present case, these common structures were obtained with OO-SCS-MP2. The reference singlet–triplet gaps in both cases are the same and are obtained using DLPNO-CCSD(T1) energies with cc-pVTZ/cc-pVQZ basis set extrapolation on the OO-SCS-MP2 optimised geometries. In the following we will focus on the overall average errors of the methods without going into details for individual carbenes of the test set. The results are presented in Table 6. Detailed results are provided as Supporting Material.

Table 6. Evaluation of density functional theory methods for singlet–triplet gaps using either geometries optimized individually with each functional, or reference OO-SCS-MP2 geometries: absolute average errors (AAE), root-mean-squared errors (RMS), and maximum errors (Max) in kcal/mol against reference cc-pV(T/Q)Z basis-set-extrapolated DLPNO-CCSD(T1) calculations.

	Optimized geometries			Reference geometries
Functional	AAE	RMS	Max	AAE	RMS	Max
mPW2PLYP	1.10	1.38	3.24	1.05	1.26	2.45
B2GP-PLYP	2.41	2.97	6.27	2.66	3.22	7.15
B2GP-PLYP-D3BJ	2.50	3.11	6.56	2.74	3.35	7.45
B2PLYP-D3	1.53	1.98	4.58	1.51	1.89	3.94
B3LYP	2.20	2.45	4.46	2.26	2.57	4.26
B3LYP-D3BJ	1.97	2.32	4.66	1.86	2.26	4.63
X3LYP	2.37	2.62	4.56	2.50	2.75	4.41
O3LYP	3.81	4.10	6.06	3.98	4.26	6.19
TPSSh-D3BJ	5.24	5.49	7.71	5.55	5.74	7.99
ωB97X-D3	3.13	3.34	4.34	3.60	3.68	4.60
B3PW-D3BJ	5.06	5.36	7.57	5.32	5.51	7.59
PBE0-D3BJ	6.08	6.29	8.25	6.36	6.47	8.35
M06-2X	4.04	4.14	5.35	4.34	4.44	6.72
PBE-D3BJ	3.54	3.95	6.65	3.65	4.00	6.53
RPBE-D3BJ	3.31	3.77	6.94	3.40	3.86	6.89
RevPBE-D3BJ	3.68	4.10	7.32	3.79	4.19	7.23
TPSS-D3BJ	4.55	4.92	8.10	4.47	4.80	7.20
BP86-D3BJ	3.51	3.93	6.86	3.60	4.01	6.75
M06-L	5.02	5.29	8.53	5.18	5.40	8.49
BLYP-D3BJ	1.67	2.19	5.32	1.70	2.11	4.82

Note: All DFT calculations were performed with the def2-QZVPP basis sets.

A crucial observation we would like to highlight from the outset is that the functionals which provide good structures as discussed in the preceding sections of this paper, do not necessarily provide good energetics according to the results of Table 6. Most prominently, TPSSh provides some of the worst singlet–triplet gaps whilst providing very good geometries across the AC-18 set. Therefore, it matters a great deal for practical applications whether the target is a high-quality geometry that can be the basis for further investigations with higher-level methods, or whether the target is a reliable screening for singlet–triplet gaps.

Most GGA functionals provide singlet–triplet gaps with average errors of more than 2 kcal/mol and maximum errors of over 6 kcal/mol. Interestingly BLYP, one of the cheapest methods assessed here, provides very reasonable average errors of less than 2 kcal/mol, with a maximum error of less than 6 kcal/mol. The functional is directly competitive to hybrid functionals such as B3LYP and X3LYP, which do not uniformly improve on BLYP over the whole set. An interesting observation is that the D3BJ dispersion corrected version of B3LYP, which is the best hybrid functional among the teste methods, performs better for singlet–triplet gaps than the uncorrected form even when using the OO-SCS-MP2 geometries. Double-hybrid methods, specifically mPW2PLYP and B2PLYP, can be recommended as the best possible choices among DFT methods. It is important to note that on average they do not improve dramatically on the singlet–triplet gaps obtained by, for example, the simple BLYP functional. However, they achieve much better error distribution. The results change little when considering singlet–triplet gaps calculated using OO-SCS-MP2 structures. In fact, the average and maximum errors are generally increased, but a gain is observed when using OO-SCS-MP2 geometries for mPW2PLYP and TPSS.

Overall, the above results suggest that BLYP is a fine choice at least for exploratory calculations on singlet–triplet gaps, whereas a conceivable composite DFT-based method that should result in both, better geometries and better energies, might combine TPSSh for geometry optimisation and mPW2PLYP for calculation of final energies.

A final comment is in order with regard to the expected accuracy of DFT methods versus wave function based approximations. OO-SCS-MP2 and OO-MP2 provide RMS errors of 3.82 and 3.49 kcal/mol for singlet–triplet gaps (AAE of 3.27 and 2.29 kcal/mol, respectively) using self-consistently optimised geometries. The average error of OO-MP2 is practically the same when it used with OO-SCS-MP2 geometries (AAE and RMS of 2.25 and 3.33 kcal/mol, respectively). The typical errors of the two methods are within 2–3 kcal/mol of the reference in most cases. A prominent exception is the challenging carbene 17, for which they both exhibit maximum absolute errors of over 10 kcal/mol compared to the reference, predicting an inverse S-T gap, i.e. a singlet ground state. A consistent difference among the two is that OO-SCS-MP2 favours the singlet state in comparison to OO-MP2, but it is likely that there is overcompensation of the triplet preference of OO-MP2, resulting in slightly worse singlet–triplet gaps on average for the SCS variant. An interesting comparison can be made with the restricted open-shell version of MP2 (RO-MP2). The performance of RO-MP2 and OO-MP2 is very similar (see Supporting Material for explicit comparisons of S-T gaps at OO-SCS-MP2 geometries), but the former improves on average over the latter with AAE and RMS errors of 1.48 and 1.60 kcal/mol, respectively. Importantly, RO-MP2 correctly predicts a triplet ground state for carbene 17, with an S-T gap within 1.5 kcal/mol of the reference. This is not the only source of improvement, however, because RO-MP2 has a slight edge over OO-MP2 even when this specific carbene is excluded from the comparison (AAE of 1.48 vs 1.69 kcal/mol). Overall, it can be concluded that MP2 methods do not match the performance of double hybrid functionals, and are surprisingly challenged even by a functional such as BLYP.

Turning to coupled-cluster approaches, it is noted that the AAE, RMS and Max values for CCSD, i.e. with the exclusion of any triples corrections, are 1.43, 1.57, and 3.90 kcal/mol using the OO-SCS-MP2 geometries. Considering the results presented in Table 6, these average metrics are comparable with, but inferior to the results of double-hybrid methods, therefore the choice of CCSD over mPW2PLYP would be questionable. A systematic improvement on singlet–triplet gaps over double-hybrid DFT can therefore only be achieved with inclusion of perturbative triples to CCSD.

4. Conclusion

A range of DFT methods and orbital-optimised versions of MP2 and SCS-MP2 were evaluated for geometric parameters and singlet–triplet gaps of aryl-carbenes against DLPNO-based coupled cluster calculations. In the case of geometry optimisation, we found that TPSSh is the most practical method that provides accurate structures at a reasonable cost. If the system size allows the application of more expensive methods, OO-SCS-MP2 provides the most accurate structures as judged against DLPNO-CCSD(T) energetics. In case of very large systems, the faster non-hybrid TPSS could be the method of choice for good geometries. Further investigation into geometrical parameters of carbenes demonstrates that the in- and out-of-plane dihedral rotations are the most challenging parameters to capture in the aryl-carbene series. Crucially, these parameters have a significant energetic consequence due to resonance stabilisation of the carbene centre. The accuracy of DFT functionals in calculation of singlet–triplet gaps is mixed. Importantly, using better or worse structures does not significantly change on average the quality of singlet–triplet gaps computed by DFT functionals. Interestingly, BLYP provides very good singlet–triplet energy gaps and we recommend this functional for exploratory calculations and the mPW2PLYP double hybrid or DLPNO-CCSD(T) when higher accuracy is desired.

Acknowledgement

Support by the Max Planck Society is gratefully acknowledged.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Max Planck Society. This work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2033–390677874 – RESOLV.