The thermochemistry of london dispersion-driven transition metal reactions: getting the 'right answer for the right reason'.

Reliable thermochemical measurements and theoretical predictions for reactions involving large transition metal complexes in which long-range intramolecular London dispersion interactions contribute significantly to their stabilization are still a challenge, particularly for reactions in solution. As an illustrative and chemically important example, two reactions are investigated where a large dipalladium complex is quenched by bulky phosphane ligands (triphenylphosphane and tricyclohexylphosphane). Reaction enthalpies and Gibbs free energies were measured by isotherm titration calorimetry (ITC) and theoretically 'back-corrected' to yield 0 K gas-phase reaction energies (ΔE). It is shown that the Gibbs free solvation energy calculated with continuum models represents the largest source of error in theoretical thermochemistry protocols. The ('back-corrected') experimental reaction energies were used to benchmark (dispersion-corrected) density functional and wave function theory methods. Particularly, we investigated whether the atom-pairwise D3 dispersion correction is also accurate for transition metal chemistry, and how accurately recently developed local coupled-cluster methods describe the important long-range electron correlation contributions. Both, modern dispersion-corrected density functions (e.g., PW6B95-D3(BJ) or B3LYP-NL), as well as the now possible DLPNO-CCSD(T) calculations, are within the 'experimental' gas phase reference value. The remaining uncertainties of 2-3 kcal mol(-1) can be essentially attributed to the solvation models. Hence, the future for accurate theoretical thermochemistry of large transition metal reactions in solution is very promising.

Introduction

Noncovalent interactions (NCIs) are of utmost significance in chemistry and biology,1–3 but their accurate experimental4–6 and theoretical7–10 description remains a challenging task. Aside from H-bonding,11 London dispersion interactions,12, 13 specifically π–π and CH–π contacts are essential for the stereospecificity of natural chemical processes14–16 and molecular recognition.17, 18 Furthermore, NCIs play an increasingly important role in the rational synthetic design of coordination polymers15, 19 and supramolecular assemblies,20 thus reinforcing efforts to gain fundamental insight into their role in chemistry and reaction thermochemistry in particular.

For a long time, NCIs between closed-shell metals were considered to matter only significantly for heavy elements, mostly late transition metals1–25 for which these so-called metallophilic interactions result from a combination of short-range electron correlation, dispersion interactions, and relativistic effects of the metal center.6–30 However, this metal-centered view is incomplete, particularly for bimetallic complexes of 3 d and 4 d transition metals, where attractive dispersion interactions with surrounding ligands become more significant and often support deficient donor–acceptor bonding between the two metals.31–33 It has been shown that in reactions involving ligand–metal bond formation, as for example, encountered in catalytic processes, larger organic ligands contribute to the ligand–metal association energy by means of a significant dispersion contribution.4–40 This work addresses the issue of how to make reliable thermochemical predictions for reactions involving larger transition metal complexes where intramolecular London dispersion interactions are suspected to be essential to their stabilization. Accurate measurements and calculations of these attractive interactions are still a challenge particularly for reactions in solution where intra- and intermolecular NCIs are partly quenched.41, 42

With modern wave function theory (WFT)43–45 and density functional theory (DFT)46, 47 -based quantum chemical (QC) methods combined with new hardware and software developments, various properties of molecules and their reactions can be calculated with high accuracy. Structural DFT predictions for isolated molecules,48 and even for molecular crystals,49–51 are quite accurate and often complementarily used to interpret experimental data. The account of the long-range London dispersion forces in DFT has also advanced in recent years,7, 10, 9, 52 such that on average for medium-sized systems (<50 atoms), the best contemporary dispersion-corrected DFT methods yield relative energies with errors of only 1–2 kcal mol−1 (<5 %).53, 54, 47 These methods are also efficient enough to investigate reactions involving larger systems (100–300 atoms), although the accuracy diminishes on an absolute scale.55 ‘Chemical accuracy’ for energetic properties (errors <1 kcal mol−1, <1–2 % relative error) can be obtained routinely by canonical coupled cluster (CC)-based methods,56 but their applicability is limited to small- and medium-sized molecules (<20–30 atoms) due to their high computational cost. However, recent local CC implementations, employing for example, pair-natural orbitals and domain-based techniques (DLPNO),57, 58 make CC calculations for molecules with up to 150 atoms and reasonable triple-ζ atomic orbital (AO) basis sets feasible. In our opinion, these methods need yet to be confronted with ‘real-life’ chemical problems in order to ensure that they can be used as reference method in general. Hence, it is necessary to provide more accurate thermochemical reference data for larger complexes of chemical interest. Particularly for transition metals, only a few extensive benchmark studies exist,59, 60 and accurate thermodynamic data for ligand–metal binding in larger complexes that could serve as reference for evaluating QC methods are scarce. In this work, experimental reference data are provided by isothermal titration calorimetry (ITC;61 see below), a technique that allows the accurate measurements of reaction enthalpies as low as 1 kcal mol−1. The thermodynamic data collected by ITC will be used to assess the performance of modern and frequently used QC methods.

As an illustrative example, the dimer breaking reaction shown in Scheme 1 was investigated. It constitutes a typical case of dispersion-driven transition metal reactions (see below), even more pronounced than in a related reaction recently studied by Sieffert and Bühl.62 The μ-chlorido-bridged palladacyclic dimer (1) derived from 2-phenylpyridine (Figure 1 a) has a singlet ground state and C2 symmetry. The dimer is essentially stabilized by the bridging chlorido ligands but at least some (multi-center) binding seems to exist as indicated by the Pd−Pd Wiberg bond index63 of 0.1. The Pd−Cl−Pd bridge is quenched by phosphane ligands (PR3, where R is phenyl (Ph) or cyclohexyl (Cy)) yielding the corresponding monopalladium species (2 a and 2 b; Figure 1 b,c). If the stoichiometric factors are taken into account, the larger system contains on each side of the reaction about 180 atoms, which is considered as sizeable for any accurate first-principles treatment. Moreover, the disruption of labile μ-chlorido-bridged palladacycles is an important elementary ligand coordination reaction, particularly central in insertion reaction processes involving unsaturated organic substrates such as carbon monoxide, alkynes and alkenes, which play the role of the incoming donor ligand that can displace the labile Cl−Pd bond.64 The latter generally binds the metal center quite stereospecifically anti to the donor ligand of the metallacycle. This elementary reaction was often used for the resolution of racemic chiral phosphanes from the related diastereomeric adducts obtained from enantiopure μ-chlorido-bridged palladacycles.65

Scheme 1

The investigated transition metal reactions (denoted throughout as ‘R=Ph’ and ‘R=Cy’). All involved molecules are uncharged and have a singlet (closed-shell) ground state.

Figure 1

Optimized structures (PBE-D3(BJ)/def2-TZVP) of a) the investigated dipalladium complex (1), b) the adduct with the PPh3 ligand (2 a), and c) the product with the PCy3 ligand (2 b). Color code: green Cl, blue N, emerald Pd, grey C, white H.

The reaction enthalpies in chlorobenzene (PhCl) solution were determined using ITC measurements. This dimer breaking reaction is known to work quite selectively.66 The ligand addition is kinetically favored in apolar solvents leading to an adduct that can possibly undergo a second ligand addition. However, the activation energy for this second addition is high in apolar solvents and hence the cationic double addition product is not observed under normal conditions (298.15 K, atmospheric pressure) in PhCl. Since the measured reaction enthalpies are not directly comparable to calculated gas-phase energies of the molecules in question, we have decided to ‘back-correct’ the experimental enthalpies to pure electronic reaction energies in the gas phase at 0 K. To this end, we have used a scheme that was already successfully applied to calculate binding affinities of supramolecular complexes.67, 68 It is based on a modified rigid-rotor-harmonic-oscillator approximation (RRHO)67 and Gibbs free solvation energy corrections calculated by continuum solvation models. It will be shown that particularly the latter introduces a sizeable error in such a scheme, and hence obtaining reliable reference energies for reactions in solution remains quite challenging. The ‘back-corrected’ experimental reaction energies will then be used to benchmark various DFT and WFT methods. The investigated molecules and similar systems that change their size and/or spatial shape in the course of a reaction are challenging to both DFT and local correlation methods since NCIs have a sizeable long-range contribution.34, 62 We will investigate whether the atom-pairwise D3 dispersion correction scheme is applicable in such a case and how dependent the results are on the chosen density functional. Furthermore, it will be discussed whether the recently published DLPNO methods57, 58 can be used as reference of benchmark quality for subsequent large-scale computational studies or screening applications with DFT methods. In this context, it will be investigated whether an MP2-based error estimate is problematic for the DLPNO methods if MP2 itself yields large errors. Moreover, the question of basis set convergence will be addressed. Having an accurate experimental reference at hand in combination with careful investigation of various error sources, it is possible to discriminate between fortuitous error cancellation and correct physical behavior of approximate methods that are feasible for larger systems of chemical interest. Consequently, it is possible to judge on the predictive power of the corresponding methods, particularly for transition metal reactions in solution with significant long-range dispersion contribution.

After a brief description of the experimental and computational details, we will discuss the role of dispersion and the challenge to obtain reliable reference values in detail followed by a comprehensive discussion of the reaction energies obtained with various DFT and WFT methods. The results will substantiate our general hypothesis that both modern DFT and WFT methods are able to predict the thermochemistry of larger transition metal complexes correctly provided that long-range London dispersion interactions are properly treated. Nonetheless, it will become obvious from our results that the term ‘chemical accuracy’ has to be redefined for computational studies of larger molecules.

Experimental and Computational Details

ITC measurements of the reaction enthalpy

Typical ITC measurements can be carried out with high accuracy from −10 °C up to +50 °C. They consist of fractioned injections of small portions (a few μL of solution) of the reagent into a cell containing the counterpart, each injection generating a thermographic response.61 This work was performed using a Waters-SAS nano-ITC61 device (two stainless steel hastealloy cells of 1 mL volume each) housed in a clean inert argon atmosphere glove box equipped with an ultrasonic degasser and an ambient atmosphere conditioning device. The reference cell was filled with pure PhCl (freshly distilled and subsequently degassed prior to usage), which is a particularly well-suited solvent since its binding capabilities to metal centers are weak. This is a crucial prerequisite since ITC measurements require clean reactions that produce no complex figures in the resulting ITC thermographs as shown in Figure 2 for the investigated reactions. They display a significant exothermic response for both reactions where 1 mm solution of the dipalladium complex 1 in PhCl was titrated at least in three independent experiments by concentrated solutions of the phosphane ligands, triphenylphosphane (PPh3) or tricyclohexylphosphane (PCy3), in the same solvent. Applying the independent substrate reaction thermodynamic model A+B⇌C (multiple set of independent binding sites69) to the measured ITC thermographic data provides an accurate fitted reaction enthalpy (ΔH) and, in a slightly less accurate way, the constant of association (Kasso), Gibbs free reaction energy (ΔG) and the reaction entropy (ΔS); for results, see Table 1. The fitted reaction enthalpies can be cross-checked against the values extracted from the raw thermogram (Figure 2).

Figure 2

Thermographic traces: a) heat released in the titration of 1 (c=1 mm) dissolved in PhCl with PPh3 (c=26.3 mm) yielding adduct 2 a: ΔH298.15 =−23.6±1.1 kcal mol−1; b) heat released in the titration of 1 (c=1 mm) dissolved in PhCl with PCy3 (c=23.3 mm) yielding adduct 2 b: ΔH298.15 =−25.3±0.3 kcal mol−1. In both cases: T=298.15 K, 600 s interval between injections of the respective ligand into 1, 90 injections of 2 μL each.

Table 1

Results of the isotherm titration calorimetry (ITC) measurements and the subsequent fitting procedure (T=298.15 K).

	ΔH [kcal mol−1]	Kasso	ΔG [kcal mol−1]	ΔS [kcal mol−1]
R=Ph	−23.6±1.1	(6.4±0.6)*106	−18.6±0.1	(−1.7±0.4)*10−2
R=Cy	−25.3±0.3	(5.0±1.0)*107	−20.9±0.4	(−1.5±0.2)*10−2

Computational setup

For the palladium atoms, an effective core potential (ECP) (SD(28,MWB)70) (and matching ECP basis sets) was used throughout this work, which will not be mentioned further. The geometries of the investigated molecules were optimized at the PBE71-D3(BJ)72/def2-TZVP73, 74 level of theory using the Turbomole suite of programs (version 6.4).75 For comparison, the dipalladium complex was also optimized without dispersion correction (see below). The Cartesian coordinates of the minimum structures (confirmed by analytical frequency calculations at the same level of theory) are listed in the Supporting Information. The thermostatistical corrections necessary for determining gas-phase reaction energies from measured reaction enthalpies were calculated with the RRHO approximation67 (see below) and the conductor like screening model for real solvents, COSMO-RS,76–78 (BP86/TZVP) or solvation model density, SMD, [79] (BP86/def2-TZVP) continuum solvation models. The latter calculations were carried out with the May 1st, 2013 (R1) release of General Atomic and Molecular Electronic Structure System (GAMESS)80, 81 employing the standard parameters for solvation in PhCl. Most single-point energy (SPE) calculations were performed with the ORCA82 program (version 3.0.083) using either the def2-TZVP73, 74 or the def2-QZVP74 basis sets (except for HF-3c84) and corresponding auxiliary basis sets (def2-TZVP/C,73 def2-QZVP/C73, 85) if applicable. Hartree–Fock (HF), RI-MP286 and RI-SCS-MP287 calculations were conducted with Turbomole 6.4 to save computation time by exploiting molecular symmetry. Furthermore, RI-MP2/def2-TZVPD88 and RI-MP2/def2-QZVPD88 calculations were carried out with Turbomole 6.4 using very tight integral accuracy criteria (scftol=10−16) for the latter to avoid self-consistent field (SCF) convergence problems. Additional HF and RI-MP2/aug-cc-pVDZ89, 90 (aug-cc-pVDZ/C91) calculations were performed with ORCA using the def2-SVP74 ECP basis set for the palladium atoms. Extrapolation to the complete basis set limit (CBS) was carried out according to Helgakers two-point extrapolation scheme separately for Kohn–Sham (KS)-DFT/HF energies92 and MP2 correlation energies.93

Various KS-DFT methods were applied in this work, namely the semi-local (generalized gradient approximation (GGA) PBE,71 the meta-GGA TPSS,94 the hybrid functionals PBE095 and B3LYP,96, 97 as well as the hybrid-meta functionals PW6B95,98 M06,99 M06-2X,99 and the double hybrid functional B2PLYP.100 The M06 and M06-2X calculations were carried out with the largest possible grid (ORCA grid7) since the results for the M06 family of functionals are known to be rather strongly grid-dependent.101, 54 All other DFT calculations were performed with grid4 (final grid: grid5). The total SCF energy was converged to 10−8 E, and the one-electron energy to 10−5 E (ORCA: tightscf; Turbomole: scfconv 7). The atom-pairwise dispersion correction (D3102, 7) with Becke–Johnson damping (D3(BJ)72) was applied for all functionals except M06 and M06-2X, where the D3 correction with ‘zero’-damping (D3(0)102) was used since medium-range NCIs are already captured at the electronic level by these functionals. The D3M+ calculations (see below) were performed with a modified standalone DFT-D3 program.103 The nonlocal (NL) van der Waals functional approach by Vydrov and van Voorhis (VV10,104 DFT-NL105) was employed non-self-consistently (i.e., calculating the NL energy once in a post-SCF manner) for PBE, PBE0 and B3LYP using grid4.

Aside from the already mentioned wave-function-based methods HF-3c, HF, MP2 and SCS-MP2, the recently developed local correlation methods DLPNO-CEPA/157 and DLPNO-CCSD(T)58 were applied. Two different values for the electron pair cut-off (T) have been used for comparison, the default value 10−4 E and a tighter value of 10−5 E. The DLPNO energies were extrapolated to the basis set limit using an additive scheme based on focal-point analysis106, 107 (denoted as δCBS). All correlation calculations were done with the frozen-core approximation (chemical core, that is, only the valence electrons were correlated).

Results and Discussion

Back-correction of the experimental reference values

The ITC measurements in solution (room temperature, PhCl) revealed that the reaction for R=Cy (ΔH=−25.3 kcal mol−1) is slightly more exothermic than that for R=Ph (ΔH=−23.6 kcal mol−1). This is in agreement with chemical intuition since PCy3 is a stronger Lewis base than PPh3. Moreover, there is no strong steric cluttering around the palladium center and thus it is easily accessible for both ligands. The measured reaction enthalpies are, however, not directly comparable to results from a standard calculation since the SPEs of the molecules in question represent pure electronic energies in the gas phase at 0 K. Hence, solvent and ro-vibrational effects as well as their temperature dependence have to be properly taken into account. Calculating reaction enthalpies from reaction energies, which on first sight seems to be the most straightforward way of comparing the experimental and calculated values, is not advisable since continuum solvation models are usually parametrized for the description of Gibbs free energies and not for enthalpies. Hence, we have used a different approach that is based on ‘back-correction’ of the Gibbs free reaction energies to pure electronic energies [Equation (1)], which then can be directly compared to calculated reaction energies. In this scheme, the thermostatistical (ro-vibrational) corrections from energy to Gibbs free energy (ΔG) are calculated from analytic PBE-D3(BJ)/def2-TZVP frequencies within a modified RRHO approximation to avoid large errors for low-frequency modes.67 The Gibbs free solvation energy contributions (ΔδG(X)) are determined using two different continuum solvation models (COSMO-RS76–78 and SMD79) for the solvent X=PhCl at T=298.15 K.

ΔE(exptl.,gas) denotes the zero-point-vibrational exclusive ‘experimental’ reaction energy for the isolated molecules. Since this ‘back-correction’ scheme is based on Gibbs free energies, the measured reaction enthalpies had to be converted to Gibbs free energies by applying a fitting procedure (see above). The theoretical corrections—as the SPEs themselves—do of course also depend on how well the optimized structures represent the actual geometries of the molecules in solution. Note that dimer 1 in the solvent cavity, as well as in the gas phase, is not planar but shows a slightly folded C2 symmetric structure. The geometry optimizations were performed at the PBE-D3(BJ)/def2-TZVP level and for comparison also with plain PBE and the same basis set. While the structures of the ligands and monopalladium complexes 2 a and 2 b do not change significantly due to the dispersion correction, the Pd−Pd distance in dimer 1 is shortened due to dispersion by about 5 % from 3.32 Å to 3.14 Å. Furthermore, the ligand moieties are folded stronger in the PBE-D3(BJ) structure (the closest C−N distance between the two sides shrinks from 6.08 Å to 5.72 Å) thus indicating the importance of interligand dispersion interactions. Additional consideration of COSMO108 in the optimizations did not change the structures significantly. Hence, we have used the gas phase PBE-D3(BJ)/def2-TZVP geometries for all further calculations (see Figure 1). The effect of different geometries on the calculated reaction energies was estimated to about 1 kcal mol−1 by comparing four methods (TPSS-D3(BJ), PW6B95-D3(BJ), B3LYP-NL, DLPNO-CCSD(T)) using the two different geometries for the reaction in which R=Ph (see the Supporting Information).

The ro-vibrational corrections to the reactions have sizeable values of 16.2 kcal mol−1 and 17.7 kcal mol−1 for R=Ph and R=Cy, respectively, with an estimated uncertainty value67 of ±1 kcal mol−1 (see the Supporting Information). However, a larger source of error stems from the Gibbs free solvation energy correction. Although the absolute value of ΔδG298.15 K(PhCl) is much smaller (−4.0 kcal mol−1 and −2.5 kcal mol−1 for R=Ph and R=Cy, respectively, with COSMO-RS) than the ro-vibrational correction, a comparison with the respective SMD values revealed significant differences. The SMD solvation model yields only −2.3 kcal mol−1 and −0.7 kcal mol−1 for R=Ph and R=Cy, respectively, from which we estimate an error of ±2 kcal mol−1 in the Gibbs free solvation energy contribution to the reactions.

In order to check for consistency, we also tried an analogous ‘back-correction’ using COSMO-RS solvation enthalpies instead of Gibbs free solvation energies. These values are estimated indirectly from the temperature dependence of calculated δG values in a 100 K window around room temperature. This very approximate procedure yields a decent ‘back-corrected’ reaction energy of −35.1 kcal mol−1 for R=Cy but the value for R=Ph of −27.0 kcal mol−1 deviates by about 4 kcal mol−1 from the result obtained with Gibbs free solvation energies. This confirms that the empirical continuum solvation models in their present form should only be used to calculate Gibbs free solvation energies. This conclusion is further supported by the results of the ab initio DLPNO-CCSD(T) calculations as discussed below.

The best estimates of the experimental reference values for both ligands are given in Table 2. We arrive at a relatively large uncertainty of ±3 kcal mol−1 (see the Supporting Information), which also reflects the size and complexity of the involved chemical species. Hence, the term ‘chemical accuracy’, which for small molecules is typically defined as 1 kcal mol−1 error has to be adjusted for the large systems that are the focus of this work.

Table 2

‘Back-corrected’ experimental reference values (in kcal mol−1).

	[R=Ph]	[R=Cy]
ΔG(exptl.,PhCl) (fitted ITC values)	−18.6	−20.9
ΔE(exptl.,gas) (E→G (COSMO-RS, PhCl))	−30.8	−36.1
ΔE(exptl.,gas) (E→G (SMD, PhCl))	−32.5	−37.9
ΔE(exptl. reference) (best estimate)	−32±3	−36±3

The special role of London dispersion interactions

Significance of the long-range NCI contribution

The HF method predicts both reactions energies qualitatively wrong (endothermic; see below), which clearly shows that their driving force is electron correlation. However, it is not the relatively large medium-range contribution but the significant long-range London part to the NCI energy, which renders this reaction particularly difficult for all QC methods that are feasible for molecules of this size. The dispersion contributions to the respective reaction energies are visualized in Figure 3 as a function of the interatomic pair distance. This information is accessible from the atom-pairwise D3 scheme by simple grouping of the dispersion energies to several distance regimes considering the stoichiometric factors of the reaction (i.e., adding up the data from all reaction partners). Note that the term ‘dispersion energy’ is only well defined in the region where two charge distributions do not overlap while at distances, where this overlap is crucial for the shape of the potential energy surface, the dispersion energy is only defined by extrapolation. Nonetheless, it is obvious that attractive correlation effects for distances larger than 4 Å have a sizeable stabilizing influence on the product side. In fact, inter-ligand, cross ligand-metal and inter-metal dispersion interactions are all important for the stabilization of these palladium complexes. Due to the larger size of the PCy3 ligand, the total D3 dispersion contribution (Edisp) to the reaction energy calculated with for example the TPSS functional is slightly larger (−25.4 kcal mol−1) for the reaction with R=Cy than for R=Ph (−24.6 kcal mol−1). A qualitative difference is, however, not recognizable from Figure 3, although their individual π–π and CH–π contributions are expected to be considerably different.

Figure 3

London dispersion contribution to the reaction energy from TPSS-D3(BJ)/TZVP124 calculations as a function of the interatomic pair distance. The contributions for distances shorter than 2 Å are summed to a single value.

In the following, dispersion-corrected DFT methods (specifically the D3 correction with the original and with modified C6 coefficients; see below) and local correlation methods (specifically DLPNO-CEPA/1 and DLPNO-CCSD(T) with different pair cut-off values) are assessed with respect to their ability to capture these long-range London dispersion interactions in a physically correct manner.

Molecular environment dependence of the C6 coefficients for (organo)metallic compounds

In many organometallic compounds the metal atoms carry a positive partial charge, for example, for the reaction investigated in this work the CM5 charges109 (PBE/TZVP) of the Pd atoms are 0.37 and 0.40 in dimer 1 and monomer 2 a, respectively. The dispersion forces associated with neutral metal atoms and metal atoms with a positive (partial) charge are quite different. Extreme examples are alkaline and earth alkaline atoms for which this oxidation state induced change of the C6 dispersion coefficients is up to three orders of magnitude.110

In the atom-pairwise D3 dispersion correction,102 the C6 coefficients are determined for different coordination numbers of an element that is using the free atoms and different neutral hydrides as model systems. In the case of metals, the electron density around the metal atom in these model systems will be larger than it is the case for typical organometallic compounds. This leads to too large C6 coefficients and to a systematic overestimation of the related dispersion energy. However, it is not well known how large this physical effect in typical organometallic compounds is, or how well it is described by for example density-based dispersion correction schemes like VV10.

Due to the importance of dispersion forces in the Pd complexes discussed in this work, we calculated new C6 coefficients based on model systems with a net charge of +1. A comprehensive investigation of analogously recomputed C6 coefficients for all metals will be published elsewhere. This modified version of the D3 dispersion correction will be called D3M+ dispersion correction in the further context. To show the difference to the standard D3 dispersion correction, the resulting C6 coefficients are shown in Table 3 along with the standard C6 coefficients.

Table 3

Coordination number (CN)-dependent C6 coefficients (in atomic units) for homogeneous atom pairs.[a]

	CN	C6[Pd−Pd]	CN	C6[Pd−Pd]	CN	C6[Pd−Pd]
D3	0	608.5	1.863	287.0	5.710	265.9
D3M+	0	67.5	1.618	133.1	2.945	99.7

[a] The C6(Pd−Pd) coefficients refer to London dispersion interactions between two Pd atoms with the same coordination number. The difference between D3 and D3M+ concerning the coordination numbers results from the use of different hydride compounds.

For the free atom (CN=0), the neutral and ionic values of the C6(Pd−Pd) differ by about an order of magnitude and thus already much less than for alkaline metals. For more realistic coordinations typically found in molecules, the neutral and ionic coefficients differ even less (only by a factor of two to three; c.f. 5th column in Table 3). It should be pointed out that only the C6 coefficients of atom pairs containing Pd are modified here; that is, the dispersion forces between for example two carbon atoms remain identical in D3 and D3M+ and also the same BJ damping scheme is applied. However, from the results given in Table 3, we expect a sizeable effect on computed reaction energies when the overall change of the dispersion energy is large and the metal atom(s) contribute significantly to the latter.

Calculated reaction energies

The results for the tested DFT and WFT methods are shown in Table 4 and Table 5, respectively (additional results can be found in the Supporting Information).

Table 4

Reaction energies obtained by density functional theory (DFT) methods.[a]

Method/Basis	ΔEreac [kcal mol−1]
	[R=Ph]	[R=Cy]
Exptl. Ref. (E→G)	−32±3	−36±3
TPSS/def2-TZVP	−18.0	−19.7
TPSS/def2-QZVP	−17.7	−19.0
TPSS/CBS	−17.7	−18.8
TPSS-D3(BJ)/CBS	−42.3	−44.2
TPSS-D3 M+(BJ)/CBS	−39.9	−41.7
PBE/def2-TZVP	−22.5	−24.6
PBE-D3(BJ)/def2-TZVP	−41.0	−44.0
PBE-D3 M+(BJ)/def2-TZVP	−39.3	−42.2
PBE-NL/def2-TZVP	−39.4	−43.6
PBE0/def2-TZVP	−22.9	−25.0
PBE0-D3(BJ)/def2-TZVP	−41.5	−44.3
PBE0-D3 M+(BJ)/def2-TZVP	−39.8	−42.5
PBE0-NL/def2-TZVP	−38.3	−42.2
B3LYP/CBS	−11.8	−12.7
B3LYP-D3(BJ)/CBS	−41.4	−43.2
B3LYP-D3 M+(BJ)/CBS	−38.5	−40.1
B3LYP-NL/CBS	−34.7	−38.3
M06-2X/CBS	−20.2	−23.6
M06-2X-D3(0)/CBS	−23.1	−27.3
M06-2X-D3 M+(0)/CBS	−23.1	−27.2
M06/def2-TZVP	−30.3	−32.8
M06/def2-QZVP	−28.8	−31.4
M06/CBS	−28.5	−31.0
M06-D3(0)/CBS	−34.8	−38.4
M06-D3 M+(0)/CBS	−34.6	−38.1
PW6B95/CBS	−21.1	−22.1
PW6B95-D3(BJ)/CBS	−34.9	−36.2
PW6B95-D3 M+(BJ)/CBS	−33.7	−35.0
B2PLYP/CBS	−27.3	−27.7
B2PLYP-D3(BJ)/CBS	−42.7	−43.4
B2PLYP-D3 M+(BJ)/CBS	−41.2	−41.9

[a] Complete basis set limit (CBS) extrapolation with def2-TZVP and def2-QZVP basis sets).

Table 5

Reaction energies obtained by wave function theory (WFT) methods.[a]

Method/Basis	ΔEreac [kcal mol−1]
	[R=Ph]	[R=Cy]
Exptl. Ref. (E→G)	−32±3	−36±3
HF-3c	−42.9	−39.0
HF/def2-TZVP	3.6	3.4
HF/def2-QZVP	4.1	4.2
HF/CBS	4.2	4.4
HF-D3(BJ)/CBS	−35.6	−40.4
HF-D3 M+(BJ)/CBS	−32.4	−37.5
MP2/def2-TZVP	−53.5	−50.6
MP2/def2-QZVP	−53.9	−51.1
MP2/CBS	−54.4	−54.4
SCS-MP2/def2-TZVP	−40.6	−37.8
SCS-MP2/def2-QZVP	−40.8	−38.1
SCS-MP2/CBS	−41.1	−41.1
TCutPairs=10−4 Eh
DLPNO-CEPA-1/def2-TZVP	−33.5	−31.8
DLPNO-CEPA-1/δCBS	−34.4	−35.6
DLPNO-CCSD/def2-TZVP	−32.9	−31.0
DLPNO-CCSD/δCBS	−33.8	−34.8
DLPNO-CCSD(T)/def2-TZVP	−35.4	−32.9
DLPNO-CCSD(T)/δCBS	−36.3	−36.7
TCutPairs=10−5 Eh
DLPNO-CEPA-1/def2-TZVP	−27.2	−27.4
DLPNO-CEPA-1/δCBS	−28.1	−31.2
DLPNO-CCSD/def2-TZVP	−27.3	−26.9
DLPNO-CCSD/δCBS	−28.1	−30.7
DLPNO-CCSD(T)/def2-TZVP	−31.2	−30.0
DLPNO-CCSD(T)/δCBS	−32.1	−33.8

[a] Complete basis set limit (CBS) extrapolation with def2-TZVP and def2-QZVP basis sets; δCBS calculated according to [Eq. (2)].

First, basis-set-related issues will be discussed. The investigated systems are difficult in this respect due to their size, which makes it impossible to apply huge basis sets very close to complete basis set (CBS) limit, at least for the correlated WFT methods. Furthermore, very large basis sets, although much less prone to basis set superposition errors (BSSE), quickly become linear dependent, and hence SCF convergence problems can occur. On the other hand, smaller basis sets have to be checked for the remaining basis set incompleteness error (BSIE) and BSSE. While the BSIE can be estimated with extrapolation schemes, the intramolecular BSSE is not directly accessible for the investigated reactions. However, the latter can be significant as shown for a different transition metal triphenylphosphane complex.62 Due to the exponential basis set convergence of the total energy in HF and DFT methods, the well balanced def2-TZVP73, 74 basis set (denoted ‘TZ’ in the following) is usually already close enough to the basis set limit. However, it is still advisable to check for the remaining BSIE in cases where this is computationally feasible. For the investigated reactions, the estimated BSIE (as obtained with extrapolation according to Halkier et al.92 using the def2-TZVP and def2-QZVP74 (denoted ‘QZ’) basis sets) for the TPSS functional amounts to only 0.3 kcal mol−1 and 0.9 kcal mol−1 for R=Ph and R=Cy, respectively. B3LYP, PW6B95 and also the HF method show essentially the same BSIE as TPSS. While the reaction energy for R=Ph is practically converged at the quadruple-ζ level, the reaction in which R=Cy shows a more pronounced basis set dependence. This difference is difficult to predict by inspection of the chemical structures beforehand, which once again underlines that computational studies aiming at high accuracy should always include a careful investigation of basis set effects. The highly parametrized M06 hybrid-meta functional, however, shows a larger BSIE of 1.8 kcal mol−1 for both reactions and the reaction energy is not converged at the QZ level. Similar observations were recently reported111 for the M06 L and M06-2X functionals that yielded an artificially large BSSE even with very extended AO basis sets thus indicating some numerical instability of the M06 family of DFT functionals. Surprisingly, our M06-2X results do not show any basis set issue and the reaction energies were practically converged at the QZ level. The double hybrid B2PLYP was expected to show a larger BSIE due to the MP2 correlation part but this was not be observed for the investigated reactions (0.3 kcal mol−1 for both using the frozen (chemical) core approximation). Although double hybrid calculations are usually performed with correlated core electrons in the MP2 part this is not advisable for transition metals since standard basis sets are not optimized for accurate correlation calculations of core levels. The B2BLYP reaction energies without frozen core differ by −1.2 kcal mol−1 and −1.3 kcal mol−1 for R=Ph and R=Cy, respectively, from the frozen core results and also show a larger BSIE of 0.8 kcal mol−1 for both reactions. All DFT methods compute the energetic ordering of two reactions correctly already at triple-ζ level thus indicating that BSIE and BSSE effects are not crucial in this case.

For correlated WFT methods, however, the results suffer much more from BSIE and BSSE effects due to the slow convergence of the correlation energy with respect to the size of the basis set. Note that for all correlated WFT calculations, the frozen (chemical) core approximation was employed, that is also the 4 s and 4 p orbitals of the palladium atoms were excluded from the correlation treatment since the applied basis sets are also not optimal for semi-core orbitals. For MP2 and SCS-MP2, the respective HF and correlation energies at TZ and QZ level were extrapolated to the basis set limit according to the scheme proposed by Halkier et al.93 Since this extrapolation scheme93 was originally developed for the correlation consistent Dunning basis sets112 we have cross-checked against extrapolation with optimized exponents113 for def2-TZVP/def2-QZVP basis sets but the results were virtually identical (deviations <0.1 kcal mol−1). Opposed to the plain DFT and HF methods, the BSIE for correlated WFT methods is expected to be negative since a larger basis also implies more electron correlation and thus stronger binding of the dipalladium complex 1. Note that in the context of larger molecules the CBS limit should be rather understood as large basis set limit since the true MP2 CBS limit cannot be reached for molecules of this size due to computational limitations. Using density fitting (RI) techniques86 it is, however, possible to perform MP2 calculations with 4000 basis functions routinely. Hence, it was also possible to check the effect of diffuse functions by employing the def2-TZVPD and def2-QZVPD88 basis sets (denoted TZD and QZD, respectively) for comparison (see the Supporting Information). It should be noted that for the latter, the SCF convergence of the palladium complexes is rather slow. This is even more pronounced for aug-cc basis sets89, 90 that contain additional diffuse functions and hence an aug-cc-pVQZ basis could not be applied. A comparative extrapolation of the MP2 energy based on the aug-cc-pVDZ and aug-cc-pVTZ basis sets was discarded due to the large BSSE of the aug-cc-pVDZ results as evident from the large error of the MP2/aug-cc-pVDZ reaction energy for R=Ph (see the Supporting Information). As expected, the BSIE for MP2 is significantly larger than for DFT methods (−0.9 kcal mol−1 and −3.8 kcal mol−1 for R=Ph and R=Cy, respectively) and at the MP2/QZ level, the reactions energies are definitely not converged. Particularly when R is Cy, the large difference between the QZ and CBS result of −3.3 kcal mol−1 indicates that enlarging the basis set merely corrects for the intramolecular BSSE, which again indicates that the reaction in which R is Cy is more prone to BSSE issues. The CBS extrapolation for SCS-MP2 yielded the same findings although the BSIE is slightly smaller in this case (−0.5 kcal mol−1 and −3.3 kcal mol−1 for R=Ph and R=Cy, respectively). The semi-diffuse basis functions included in the QZ basis set are able to capture the long-range dispersion interactions of the investigated molecules, which is reflected in the convergence of the reaction energies towards the experimental reference value for both high-level WFT and modern DFT methods. There are, however, cases for which semi-diffuse basis functions are not sufficient as for example, shown in a recent study of a frustrated Lewis pair.114 Since the largest system of the present study has 3737 basis functions (8246 auxiliary basis functions) at the QZ level, DLPNO-CCSD(T)/QZ with tight pair cut-off would be very computational demanding. However, DLPNO-CCSD(T)/TZ calculations are efficient enough to be applied routinely for molecules of this size (SPE calculations of the largest system take only about a day even for tight thresholds). Furthermore, the TZ basis set benefits from error compensation between BSIE and BSSE.68 Hence, the BSIE of the DLPNO-CCSD(T) results (and analogously for DLPNO-CEPA-1 and DLPNO-CCSD) was estimated (denoted δCBS) using focal point analysis106, 107 for the electronic energy E in which a correction from a TZ calculation is added to the respective MP2/CBS results.

Such an additivity scheme is only valid if the BSIE is similar for MP2 and CCSD(T)115, 116 and the MP2 methods itself yields reasonable results for the system in question. Although the MP2 method overestimates the reaction energies of both reactions, we are confident that this additivity scheme can be still applied since the BSIE correction repaired the wrong energetic ordering at MP2/TZ level and the less overestimated SCS-MP2 results showed the same trend. This assumption could be confirmed by the DLPNO-CCSD(T)/δCBS results that predict the small energetic difference of these two very similar reactions correctly while at the TZ level the reaction with PPh3 is wrongly predicted to be more exothermic.

Since the investigated reactions differ only by the phosphane ligand used to quench the dimer 1, the difference of their reaction energies (≈10 %) is rather small. Hence, it poses a challenging test case for computational thermochemistry as consistency for varying ligands can be extremely important in practice. Investigation of two similar reactions also provides an internal consistency check of our theoretical procedures.

Both reaction energies are predicted qualitatively wrong at the HF level (endothermic; see Table 5), which is expected since the reactions are driven by electron correlation and in particular, long-range London dispersion interactions as discussed above. Hence, not surprisingly, DFT methods without dispersion correction are also not able to reproduce the experimental reference values. All tested uncorrected DFT functionals underestimate the reaction energies although the ligand effect is predicted correctly. Comparing the (meta-)GGA functionals TPSS and PBE with the hybrid-(meta) functionals PBE0, PW6B95 and B3LYP there is no systematic improvement visible and the latter even shows the largest deviations of 17.2 kcal mol−1 and 20.3 kcal mol−1 for R=Ph and R=Cy, respectively. Although the medium-range NCI energy is already included in the highly parametrized M06-2X and M06 functionals by construction, they still underbind due to the missing long-range part. The better performance of the latter agrees with the recommendation by Zhao and Truhlar that M06 is fitted for systems involving both NCIs and transition metals while the M06-2X functional was proposed to be accurate for aromatic stacking interactions of main group chemistry.117 A similar observation is made for the double hybrid B2PLYP, which includes NL correlation through its scaled MP2 contribution. Since the latter is strongly overestimated for the investigated reactions the underbinding of B2PLYP is less pronounced. The long-range London dispersion interaction is, however, missing in all functionals and fortuitous agreement with the experimental reference values would be just a matter of uncontrollable error cancellation.

The atom-pairwise D3 dispersion correction7 (with Becke–Johnson damping72) corrects for a large part of the error. It yields a sizeable contribution to the reaction energies and for TPSS and B3LYP the D3 contribution is even larger than that of the functional itself (e.g. −30.5 kcal mol−1 for R=Cy with B3LYP). The TPSS-D3(BJ), PBE-D3(BJ), PBE0-D3(BJ), B3LYP-D3(BJ), and B2PLYP-D3(BJ) reaction energies differ much less from each other than for the plain functional, which is a very positive result. The values are all too large (about 5–7 kcal mol−1) but the reaction in which R=Ph is still correctly predicted as being less exothermic. Besides the intrinsic errors of the functionals, the overbinding on the product side (which involves the larger dispersion energy) also indicates a slight overshooting of the D3 correction. The PW6B95-D3(BJ) results, however, are already within the error bars of the experimental reference values thus indicating that this functional shows clearly smaller errors for the investigated molecules and that the overestimation of the D3 dispersion energy is only a small part of the remaining error for the other functionals. The M06 and M06-2X energies can be improved by adding the missing long-range part of the dispersion energy via a D3 correction with ‘zero’-damping to avoid double counting of medium-range contributions (D3(0)102) and only after correcting for this missing physics, the reaction energies are approaching the experimental reference value for the right reason, at least with the M06 functional. As discussed above the slight overestimation of the D3 correction can be attributed to the too large C6 coefficients of the Pd atoms. The new C6 coefficients (denoted D3M+) based on positively charged model systems are assumed to be more appropriate to describe the partially positive charged metal in the complexes. This is proven by the systematic improvement of all DFT results with D3M+ correction compared to the D3 results. The overestimation of the reaction energies is reduced by up to 10 % (3.1 kcal mol−1 for R=Cy with B3LYP) and consequently all tested DFT functionals combined with D3 M+ deviate by less than 5 kcal mol−1 from the experimental reaction energies. Moreover, the PW6B95-D3 M+(BJ) results perfectly agree with the experimental reference (similar for M06-D3 M+(0)). Only the B2PLYP-D3 M+(BJ) reaction energy for R=Ph is still too large by ≈6 kcal mol−1 but this can be attributed to the strong overestimation of the MP2 correlation energy for molecules with significant π–π interactions and potentially also to problems associated with perturbation theory applied to metal containing systems. Although scaled down in double hybrids, these MP2 related errors directly show up in the B2PLYP results and deteriorate the usually superior performance of double hybrid functionals in this case. As a sanity check for D3M+, the NL van der Waals functional approach by Vydrov and van Voorhis (VV10,104 DFT-NL105) was employed for PBE, PBE0, and B3LYP. The results obtained with the density dependent dispersion correction are in good agreement with the D3M+ results thus indicating that the new C6 coefficients for the palladium atoms are physically sound. This argument can be reversed by stating that the local dipole polarizability model in VV10 seems to reflect the change in the electronic environment in the reaction quite well. Moreover, the B3LYP-NL reaction energies are within the error bound of the experimental reference.

Adding dispersion contributions to the plain HF reaction energies via the D3M+ method leads to remarkably accurate results, which are within the error bars of the experimental reference values and only the energetic difference between the two reactions is predicted slightly too large. This again emphasizes that the investigated reactions are driven by London dispersion interactions, which can be accurately estimated by the atom-pairwise correction scheme. Notably, the D3M+ contributions for the HF method (−36.6 kcal mol−1 and −41.9 kcal mol−1 for R=Ph and R=Cy, respectively) are even larger than the reaction energies itself. The minimal basis set method HF-3c84 which incorporates a D3 and BSSE correction yields surprisingly accurate reaction energies for R=Cy, while the overbinding of about 8 kcal mol−1 for R=Ph is comparable to the B2PLYP-D3(BJ) result. Although the energetic ordering is predicted wrongly, HF-3c offers a very fast possibility (the entire calculations take only a few minutes) for a preliminary investigation of large but electronically simple complexes where pure semi-empirical methods are not applicable due to missing or improper parametrization.

Besides the already discussed basis set and metal related issues, the MP2 method and to a much lesser extent, SCS-MP2, suffer from strong overestimation of π–π and CH–π dispersion interactions118 which is most evident for monopalladium complex 2 a. Consequently, the reaction energy for R=Ph is spuriously identical to that of the reaction in which R=Cy and about 20 kcal mol−1 too exothermic. Notably, SCS-MP2 reduces the overbinding by 13 kcal mol−1 compared to standard MP2 but the difference between the two ligands is still not described correctly. Obviously, correlated methods beyond MP2 such as the ‘gold standard’ CCSD(T) are mandatory in this case to get accurate results with WFT methods but their unfavorable scaling of the computational demands with system size necessitates further approximations. One promising approach in this direction is the recently published DLPNO scheme, which allows CCSD(T) calculations of larger molecules with only minor loss of accuracy.58 The DLPNO-CCSD(T) reaction energies are slightly too large if the default value for the pair cut-off (T=10−4 E) is applied. Tightening this threshold by one order of magnitude, however, yields reaction energies that are in excellent agreement with the experimental reference values. A closer inspection reveals that the difference of the reaction energy obtained with the two threshold values is significantly larger for R=Ph (4.2 kcal mol−1) than for R=Cy (2.9 kcal mol−1) which furthermore underlines that the reaction with the PPh3 ligand is more problematic for MP2. The calculations with the tighter pair cut-off are approximately twice as expensive in the latter case. About 50 % of the electron pairs enter the post-MP2 correlation treatment while for the default value only about 25 % of the pairs are kept, which in turn increases the MP2 error estimate contribution to the total correlation energy from ≈0.5 % to ≈1 %. In order to determine the importance of connected triple excitations in CC calculations, the DLPNO-CCSD(T) results were compared to DLPNO-CCSD and DLPNO-CEPA/1. Previous studies on smaller molecules have shown that CEPA/1 clearly outperforms CCSD for thermochemical applications,119 which was ascribed to an implicit inclusion of triple excitation contributions. However, as can be seen from the results in Table 5, this is not the case for the reaction investigated in this work and DLPNO-CEPA/1 and DLPNO-CCSD with tight threshold values yield almost identical reaction energies that are about 1–2 kcal mol−1 below the lower error bound of the reference values. DLPNO-CEPA/1 and DLPNO-CCSD with default threshold values give the right answer for the wrong reason since they merely benefit from error cancellation due to overestimation of the MP2 correction for neglected pairs and the missing triples contribution. Comparing the DLPNO-CCSD and DLPNO-CCSD(T) results with the T value of 10−5 E reveals significant contributions from the connected triple excitations, which are almost twice as large for R=Ph (−4.0 kcal mol−1) than for R=Cy (−2.6 kcal mol−1). Hence, only if connected triple excitations are included at least perturbatively and a sufficient number of electron pairs are kept thus treating also the long-range London dispersion interactions accurately, the ‘right answer for the right reason’ is obtained provided that the BSSE and BSIE issues could be ruled out by for example, proper CBS extrapolation.

In summary, good performers yielding reaction energies within the ±3 kcal mol−1 error bars (≈10 %) of the experimental reference for both reactions are B3LYP-NL/CBS, M06-D3 M+(0)/CBS, HF-D3 M+(BJ)/CBS, PW6B95-D3 M+(BJ)/CBS, and DLPNO-CCSD(T)/δCBS (T=10−5 E) in order of increasing accuracy. This agreement between experimental and theoretical results is very encouraging, particularly considering the size and complexity of the investigated molecules and the large errors observed for other frequently used DFT or WFT methods as for example, M06-2X or MP2. The good performance of B3LYP-NL squares with the observation made by Hujo et al.105 that there is a beneficial compensation of overrepulsive (hybrid GGA) and slightly overbinding VV10 parts. However, the PW6B95-D3 M+(BJ) method (or PW6B95-D3(BJ)) generally provides superior results at a similar level of empiricism and numerical complexity. Since a few years this method serves as our default density functional on the hybrid-meta level for SPE calculations. After incorporating long-range interactions via the atom-pairwise dispersion correction with ‘zero’-damping (D3(0) or D3 M+(0)) the highly parametrized M06 functional also yielded good results but it is numerically less stable. The good performance of HF-D3 M+(BJ) is certainly not generalizable, while the most accurate reaction energies predicted by the DLPNO-CCSD(T)/δCBS method with tight pair cut-off of 10−5 E indicates the great potential of the latter method to serve as reference in case of missing experimental data.

Conclusions

This joint experimental and theoretical investigation addressed the thermochemistry of two exemplary reactions in solution involving larger palladium complexes where long-range intramolecular London dispersion forces are essential for the exothermicity of the reaction. Such systems and other larger 3 d and 4 d transition metal complexes do play an important role for example, in catalytic processes and are in focus of current research.39, 40, 62, 120, 4–38 Our results obtained for this challenging ‘real-life’ chemical problem substantiate that significant progress has been made in recent years for both DFT and WFT methods. They are in agreement with a previous study about ligand substitution energies for transition metal complexes in solution.42 According to our experience over the last years, accurate predictions of the geometrical structure and electronic energies of molecules do not pose a fundamental problem to modern QC methods anymore at least for closed-shell molecular systems in the gas phase with up to about 200 atoms. Modern density functionals as for example, PW6B95 perform very well if properly corrected for long-range London dispersion effects via for example, the atom-pairwise D3 (or D3M+) correction. Contrary to small transition metal complexes, larger ones are significantly stronger stabilized by long-range intramolecular dispersion interactions, which is particularly important for dissociation reactions. These interactions are missing in the M06 family of functionals, which only include medium-range NCIs and hence they tend to under-bind for larger molecules,34, 62 and the good performance for small organometallic systems121, 117 is probably not fully generalizable to large complexes with significant long-range interactions. The atom-pairwise D3 dispersion correction slightly overshoots in our example because the metal has a positive (partial) charge and hence a reduced effective number of electrons for the dispersion interaction. However, with re-computed C6 coefficients derived from charged systems (DM3+), a good agreement with the non-local density-dependent VV10 (NL) energies is obtained. Whether this is a general phenomenon will be investigated in a subsequent comprehensive study of D3M+ published elsewhere.

It is encouraging to see that the recently proposed DLPNO-CCSD(T) method can serve as accurate alternative applicable for larger systems where CCSD(T) is not feasible. It is very helpful to judge on the performance of still more efficient (and less basis set dependent) DFT methods. However, more experience has to be gained with respect to the choice of the only critical threshold, the pair cut-off. We propose a tighter T value of 10−5 E in cases where MP2 does not provide accurate correlation energies in order to decrease the influence of MP2-related errors. Nonetheless, we are confident that the DLPNO-CCSD(T) method has great potential to serve as general reference method also in theoretical organometallic thermochemistry.

It should be noted that for small-molecule thermochemical benchmark sets,59, 60, 122 the issues associated with long-range London dispersion interactions are less evident and there is need for further transition metal chemistry benchmarks particularly for larger systems of chemical interest. In this context, it is of great advantage to have accurate experimental thermodynamical data (e.g., provided by ITC measurements), which is part of an ongoing collaboration of our groups. However, experimental data for larger molecules are typically obtained for reactions in solution and hence comparisons to the calculated gas-phase data require corrections for thermal and entropic effects as well as solvation contributions. The thermostatistical calculation of reaction entropies calculated with for example, the RRHO model and the Gibbs free solvation energy contribution to a reaction calculated with continuum solvation models show uncertainties of ≈1 kcal mol−1 and ≈2 kcal mol−1, respectively, for larger systems. In most cases, entropic and solvation contributions do not cancel each other, and a more accurate description of particularly the latter would be very beneficial. New developments of explicit solvation models may improve the situation in the future.

Although nowadays many experimental transition metal studies are complemented by theoretical thermochemistry, the level of theory used is not always sufficient, which can easily result in wrong conclusions.123 Cancellation of different error sources may lead to the correct result, but reliable predictions are only possible if all physical interactions are captured properly. This is particularly important if experimental reference values are missing as it is the case for most larger molecule studies especially in transition metal chemistry. The fortuitous agreement of numerical results should not lead one to believe that physically important contributions can be neglected and aiming at the ‘right answer for the right reason’ remains an important issue in quantum chemistry. Nonetheless, it is encouraging that both modern dispersion-corrected density functionals and local CC approaches are able to predict the thermochemistry of larger transition metal complexes with significant long-range London dispersion interactions correctly and, provided that all physical interactions are all described properly, the modern DFT and WFT methods converge practically to the experimental reference value. The general validity of our findings can finally only be judged through extensive application of the proposed thermochemistry protocol to various chemical reactions in combination with further ITC measurements, which is ongoing work in our groups. However, it should be reiterated that realistic error estimates of sizeable (‘back-corrected’) experimental reaction energies for large molecules do not comply with standard notion of a ‘chemical accuracy’ of 1 kcal mol−1. Instead, an accuracy of 2–3 kcal mol−1 (≈5–10 % of ΔE) seems more realistic and is indeed achievable with modern dispersion-corrected density functionals as for example, PW6B96-D3(BJ) in combination with nowadays possible large-scale DLPNO-CCSD(T) calculations, thus opening a promising future for theoretical thermochemistry.

Experimental Section

Complex 1 was prepared according to a published procedure,40 while the phosphane ligands were purchased and purified. The synthesis, spectroscopic characterization (1H NMR, 13C NMR, ESI-MS) of the adducts 2 a und 2 b and structural X-ray diffraction analysis for 2 a is described in the Supporting Information.

CCDC 1015924 contains the supplementary crystallographic data for the structural X-ray diffraction analysis for 2 a; these data can be obtained free of charge from The Cambridge Crystallographic Data Centre via http://www.ccdc.cam.ac.uk.