Coupled Cluster Benchmark of New Density Functionals and of Domain Pair Natural Orbital Methods: Mechanisms of Hydroarylation and Oxidative Coupling Catalyzed by Ru(II) Chloride Carbonyls Irena Efremenko1 and Jan M.L. Martin1, a) 1 Department of Organic Chemistry, Weizmann Institute of Science, 76100 Reḥovot, Israel a) Corresponding author: gershom@weizmann.ac.il Abstract. In the present work we tested the performance of several new functionals for studying the mechanisms of concurrent reaction of hydroarylation and oxidative coupling catalyzed by Ru(II) chloride carbonyls. We find that DLPNOCCSD(T) is an acceptable substitute for full canonical CCSD(T) calculations; that the recent ωB97X-V and ωB97M-V functionals exhibit superior performance to commonly used DFT functionals; and that the revised DSD-PBEP86 double hybrid represents an improvement over the original, even though transition metals were not involved in its parametrization. INTRODUCTION Synthetic methods allowing one-step C-C bond formation through homogenous-catalyst-mediated transformation of C-H bonds have become increasingly important in both industry and academia.1,2 Potential routes for selective CH bond activation and subsequent C-C bond formation in alkenes and aromatic compounds include hydroarylation by addition of aromatic C–H bonds across an unsaturated C=C bond, or an oxidative coupling that preserves the double bond. The latter reaction, while a highly desirable industrial goal, is challenging from the synthetic point of view. Pioneering examples of Ru-catalyzed coupling of aromatic carbon-hydrogen bonds with olefins were reported two decades ago.3,4 Since then (for a review, see Ref. 2), an increasing number of examples catalyzed by Rh, Ru and Pd have been published, but the mechanistic aspects of the reactions have only been addressed by experimental methods. Motivated by the experimental results of Milstein and coworkers,4 over a decade ago we started to explore the mechanisms of the concurrent reactions of oxidative coupling and hydroarylation of methyl acrylate (MA) catalyzed by Ru carbonyl complexes (Scheme 1) using hybrid and double hybrid DFT families. These calculations showed that either proton elimination by chloride ion, or hydrogen transfer to coordinated olefin, can serve as the initial step of aromatic C-H bond activation, while oxidative addition mechanism could be excluded. Proton elimination proceeds via transition state TS1 and initiates COOMe oxidative coupling with olefin (TS2) COOMe + + COOMe according to Scheme 2a. Next, methyl oxidative cinnamate coupling interaction of Ru hydride intermediate E isomer RuCl2, CO (1) + Int3 with MA and HCl regenerates the COOMe 180oC methyl acrylate initial RuCl2 carbonyl (TS4a). COOMe hydroCOOMe (MA) arylation Alternatively, the catalytic cycle could or methyl be closed by interaction of Int3 with hydrocinnamate (2) (3) linear branched MA and benzene, yielding Int1 (TS4b). Scheme 1. Possible interactions between MA and benzene in presence of Hydrogen transfer to MA (TS1a and TS1b) causes coexistence of phenyl and Ru(II) chloride carbonyl complexes in benzene solution: oxidative one of the two isomeric alkyl ligands in coupling (observed) and hydroarylation (not observed). the Ru coordination sphere (Int4 and Int5). This mainly leads to the hydroarylation products via TS3 and TS3a (Scheme 2,b), but coordination of the second olefin molecule followed by oxidative coupling is also possible. Scheme 2. Mechanisms of MA interactions with benzene in presence of Ru(II) chloride carbonyl complexes: oxidative coupling (a) and hydroarylation (b). We found that the activation barriers, the relative energies of the key intermediates, and the overall RuCl2 RuCl2 RuCl2 direction of the catalytic reaction strongly depend on CO CO CO R R R the composition of the Ru coordination sphere. Ru Cl complexes that could form in the reaction mixture, and CO serve as initial species of catalytic cycles ,are shown in RuCl2 RuCl2 RuCl Cl2Ru RuCl2 CO OC OC CO Scheme 3. In all the complexes, chloride anions are OC CO CO CO CO strongly bound to the metal atom; the only exception is Ru(CO)2Cl2 RuCl2(CO)4(benzene) (rightmost complex in the 2nd Ru(CO)Cl2 Ru(CO)4Cl2 Ru(CO)3Cl2 row of Scheme 3) showing CO insertion into one of the Ru-Cl bonds. Moreover, the energetic results, O O O R particularly calculated activation barriers, differ O O O RuCl2 RuCl2 RuCl2 Ru between computational approaches and were hard to CO OC Cl2 OC CO CO CO reconcile with experimental observations. Therefore, for this specific case, we will attempt to Scheme 3. Calculated RuCl2 complexes with benzene, CO obtain rigorous first-principles results by means of and MA. Arrows indicate that complex [RuCl2(benzene)] is coupled cluster theory near the complete basis set limit, used as a reference. and use these results to assess the performance of more affordable computational methods. We believe our results have broader relevance for modeling mechanisms of catalytic reactions mediated by transition metal complexes. The CCSD(T) method5 is considered the “gold standard” of quantum chemistry. However, the computational cost of canonical CCSD(T) calculations scales as O(N7) and becomes prohibitively high for mechanistic studies of practical transition metal catalysis problems. Recently developed domain pair natural orbital methods, such as DLPNOCCSD(T) of Neese and coworkers6 and PNO-LCCSD(T) of Werner and coworkers,7 scale almost linearly with system size (at least for closed-shell cases) and, in the main group, provide similar accuracy to the corresponding canonical calculation. Recently, benchmark studies of the performance of density functionals for transition metal problems, using DLPNO-CCSD(T) for calibration, have started appearing for reaction energies8 and barrier heights.9 Since the Ru complexes shown in Schemes 2 and 3 are still barely tractable by canonical methods, this enables us to assess the “domain error” for real-size transition metal complexes. In this paper, we assess both DFT and PNO methods against canonical CCSD(T) for the hydroarylation and oxidative coupling of benzene and methyl acrylate (MA) catalyzed by RuCl2-carbonyl complexes, as a representative example for the complex mechanisms of homogeneous catalytic reactions. Having in mind specific computational problems usually addressed in such mechanistic studies, the calculations were divided into four groups: (i) overall reaction energies (Scheme 1) that do not involve transition metals; (ii) relative energies of stable RuCl2 complexes with CO, benzene and MA (Scheme 3) that should reproduce breaking and formation of the metal-ligand bonds and deep alterations of the coordination sphere and electronic structure of the metal atom; (ii) energies of key intermediates along reaction pathways catalyzed by different Ru complexes (Scheme 2) and (iv) barrier heights along these reaction pathways. It would be natural to calculate relative energies of the carbonyl complexes relative to RuCl2, however, our calculations revealed that it has a triplet ground state, and that the singlet is essentially purely biradical (which is also reflected in the pathological D1 diagnostic10 value of 0.435). Discussion of the open-shell calculations are beyond the scope of this preliminary report; therefore, we used RuCl2(benzene) as a reference, as RuCl2 will anyway have no independent existence under the experimental conditions (benzene solvent). At all levels, 1:1 exchange of benzene with CO or MA is energetically unfavorable; all other complexes are exothermic with respect to the reference. The relative energies of key intermediates and transition states along each reaction path were calculated relative to the initial form of the catalyst, (C6H6)(CO)nRuCl2 (n=0-4). COMPUTATIONAL METHODS The Weigend-Ahlrichs basis set family11 def2-TZVP, def2-TZVPP, and def2-QZVPP was used throughout. Reference geometries were optimized at the PBE0-D3BJ/def2-TZVP level12,13 using Gaussian 09;14 identities of transition states were verified by frequency and intrinsic reaction coordinate calculations. At the final geometries, canonical CCSD(T)/def2-TZVPP single-point energy calculations were performed using MOLPRO 2018,15 both using default frozen cores and including Ru(4s,4p) subvalence orbitals (which are correlated by default in ORCA. We found in this work that the mean absolute effects of Ru(4s,4p) subvalence correlation on carbonyl ligand energies and transition states are both a nontrivial 1.0 kcal/mol.) DLPNO-CCSD(T)16 and the version with improved iterative triples, DLPNO-CCSD(T1),17 were calculated with the def-TZVPP basis set using ORCA,18 likewise DLPNO-CCSD(T)/def2-QZVPP calculations were done for basis set extrapolation using the simple L–3 formula.19 TightPNO cutoffs20 were used to reduce domain discretization error; we found in the present work that DefaultPNO causes errors up to 3.5 kcal/mol in energy differences, and hence do not recommend its use. In addition, single-point DFT calculations with a number of DFT functionals were carried out using ORCA. Aside from PBE0 already mentioned, these include: (a) the Berkeley “combinatorially optimized” 21 B97M-V, ωB97X-V and ωB97M-V; (b) the M06 family:22 M06-L, M06 and M06-2X; (c) TPSS23 and two different hybrids thereof, namely, TPSSh and TPSS0 (10% and 25% HF exchange, respectively); (d) both the original double-hybrid DSDPBEP86-D3BJ24 and its reparametrized version revDSD-PBEP86-D4;25 in the latter, D3BJ also replaced with the very recently published next-generation D4 model.26 As basis set convergence of double hybrids tends to be dominated by the MP2-like term, we carried out def2-TZVPP and def2-QZVPP calculations and applied L–3 basis set extrapolation;27 for the remaining DFT functionals we applied def2-TZVPP except for ωB97, which were accurate enough that we also tried def2-QZVPP. (Changes are on the order of 1 kcal/mol.) GRID6 was used in all DFT calculations. In all Orca calculations, the RIJCOSX approximation28 was employed, as well as the RI-MP2 approximation29 for the double hybrids, in conjunction with the respective appropriate auxiliary basis sets30 ,31 for the def2 family. RESULTS AND DISCUSSION MAD (mean absolute deviation) and RMSD (root mean square deviation) error statistics with respect to our best CCSD(T) basis set limit estimates are shown in Figure 1 for the four types of energy differences. As expected, the smallest deviations were found for the overall reaction energies, which involve only main group elements. The hybrid functionals of the M06 family (MAD=0.68 and 0.57 kcal/mol for M06 and M06-2X, respectively) and double hybrid revDSD-PBEP86 functionals (MAD=0.64) show the best performance in this group. However, overall for the four criteria, the best results were obtained using ωB97M-V and ωB97X-V range-separated hybrids as well as by the revDSD-PBEP86 double hybrid, with accuracy between DLPNO-CCSD and DLPNO-CCSD(T). On the 3rd rung (meta-GGA) of the Jacob’s Ladder, B97M-V performed best for reaction energies and carbonyl complex stabilities; however, M06-L showed similar performance, and TPSS outperformed them for barrier heights (MAD 3.69 vs. 4.54 for M06-L and 5.47 for B97M-V). The revDSD-PBEP86 functional outperforms the original DSD-PBEP86 for all four groups of calculations, whereas DSD-SCAN-D4 improved on DSD-PBEP86-D3BJ only for the carbonyl complexes. The DLNPO-CCSD(T) approach shows very close agreement (MAD=0.35 kcal/mol) with canonical CCSD(T) for the reaction energies, while for other groups there is a bit more daylight between them (MAD=1.04 for carbonyls, 1.56 for intermediates, and 1.58 kcal/mol for TSes). Using the DLNPO-CCSD(T1) approach with improved perturbative triples decreases these statistics to 0.60, 0.80, and 1.02 kcal/mol, respectively. The RMSD/MAD ratios are close to the theoretical value32 (for a normal distribution) of √(π/2)≈1.2533, for carbonyls and intermediates but much larger for TSes, indicating an outlier (RMSD=0.80, 0.86, 2.05 kcal/mol): we note that TS2-CO2 is essentially biradical (and has D1>0.4). If we exclude it, MAD and RMSD drop to quite pleasing values of 0.63 and 0.73 kcal/mol, respectively. TSes minima CO ∆E_ligand ∆Ee,reac TSes minima EM E- P2 E- CCS CC D SD (T ) 3B J TP E0 -D PB D D D 0 06 6 -2 X M M M 0D PB E CO ∆E_ligand D 0.00 D 0.00 D 3.00 B9 7M B9 -V w 7X B9 -V D 7M SD -P -V BE re P8 vD 6 SD D 3 D -P BJ SD B -S EP8 CA 6 N D 4 6.00 3.00 w 6.00 S TP S SS h SC TP A SS N 0 D 3B J 9.00 EM E- P2 C E- CS CC D SD (T ) 9.00 B9 7 w MB9 V w 7X B9 -V D 7M SD -P -V BE re P8 vD 6 SD -D 3 D -P BJ SD BE -S P8 CA 6 N D 4 12.00 06 L 12.00 3B J TP S TP S SS h SC TP A SS N 0 D 3B J 15.00 M 06 L M M 06 06 -2 X MAD relative to CCSD(T)/CBS (kcal/mol) RMSD relative to CCSD(T)/CBS (kcal/mol) 15.00 ∆Ee,reac Figure 1. RMSD and MAD relative to CCSD(T)/CBS for the four types of energetics considered. At all levels except for DLPNO-MP2, the highest RMSD and MAD values were found for the first reaction, i.e. dissociation of benzene and association of CO and MA ligands. Most density functionals tend to overbind the ligands; detailed analysis shows that for one group of functionals (PBE0-D3BJ, SCAN-D3BJ, TPSS, M06L and DSD-PBEP86 family, Fig. 2a) the error becomes greater with increasing number of CO ligands; the other group (M06, M06-2X, TPSS0, TPSSh, and B97, Fig. 2b) exhibits a less pronounced opposite trend, especially in the presence of coordinated benzene. Carbonyls: Energy Deviation relative to CCSD(T)/CBS (kcal/mol) 10 PBE0-D3BJ SCAN D3BJ DSD-PBEP86-D3BJ revDSD-PBEP86 DSD-SCAN D4 M06L a b 5 TPSS Number of CO ligands (n) Number of CO ligands (n) 0 0 2 3 4 0 1 2 3 4 0 1 2 3 0 1 2 1 ΔE, kcal/mol ΔE, kcal/mol 1 -10 -20 2 3 4 0 1 2 3 4 0 2 3 0 1 2 -5 -10 (C6H6) (CO)nRuCl2 -30 (MA) (CO)nRuCl2 -15 (CO)nRuCl2 -40 1 (C6H6) (CO)nRuCl2 (MA) (CO)nRuCl2 (C6H6)(MA) (CO)nRuCl2 -20 (CO)nRul2 M06 M06-2X B97M-V wB97M-V TPSS0 TPSSh wB97X-V (C6H6)(MA) (CO)nRuCl2 Figure 2. Energy deviation relative to CCSD(T)/CBS in different RuCl2 (CO)n complexes as a function of n. Our findings are consistent (see also our companion paper in the present volume) with the findings of Najibi and Goerigk33 and ourselves25 for the very large GMTKN55 main-group benchmark34 and of Iron and Janes9 for the MOBH35 transition metal reaction benchmark: Notably, that the range-separated hybrids ωB97X-V and ωB97M-V acquit themselves particularly well, that revDSD represents an improvement over the original DSD not just for the main group but also transition metals, and that unlike for the main group where empirical double hybrids are clearly superior, they offer no clear advantage over ωB97M-V for transition metal reactions. Unlike Iron and Janes, however, who found the new DSD-SCAN double hybrid25 to be among the best performers for MOBH35, we find it to be inferior to revDSD-PBEP86 and the ωB97n-V family for the present problem (Figure 1). ACKNOWLEDGMENTS This research was supported by the Israel Science Foundation (grant 1358/15), the Minerva Foundation, and the Helen and Martin Kimmel Center for Molecular Design (Weizmann Institute of Science). REFERENCES 1 S. Murai, F. Kakiuchi, S. Sekine, Y. Tanaka, A. Kamatani, M. Sonoda, and N. Chatani, Nature 366, 529 (1993); J. Kischel, I. Jovel, K. Mertins, A. Zapf, and M. Beller, Org. Lett., 8, 19 (2006); H. A. Zhong, J. A. Labinger, and J. E. Bercaw, J. Am. Chem. Soc. 2002, 124, 1378; C. Jia, D. Piao, J. Oyamada, W. Lu, T. Kitamura, and Y. Fujiwara, Science 287, 1992 (2000). 2 (a) From C-H to C-C Bonds: Cross-Dehydrogenative-Coupling, Chao-Jun Li, ed. Royal Society of Chemistry, 2014. (b) Catalytic Arylation Methods: From the Academic Lab to Industrial Processes, A. J. Burke, C. S. Marques, Ed.. Wiley, 2015. 3 F. Kakiuchi, T. Sato, M. Yamauchi, N. Chatani and S. Murai, Chem. Lett. 28, 19 (1999). 4 H. Weissman, X. Song, and D. Milstein, J. Am. Chem. Soc. 123, 337 (2001). D. Milstein, D. Pure Appl. Chem. 75, 445 (2003). 5 I. Shavitt and R.J. Bartlett, Many – Body Methods in Chemistry and Physics (Cambridge University Press, Cambridge, 2009). 6 C. Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, J. Chem. Phys. 139, 134101 (2013). 7 Q. Ma and H. Werner, WIRes Comput. Mol. Sci. 8, e1371 (2018) and references therein. 8 S. Dohm, A. Hansen, M. Steinmetz, S. Grimme, and M.P. Checinski, J. Chem. Theory Comput. 14, 2596 (2018). 9 M. A. Iron, and T. Janes. J. Phys. Chem. A 123, 3761-3781 (2019). 10 C. L. Janssen and I. M. B. Nielsen, Chem. Phys. Lett. 290, 423–430 (1998). 11 F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005) 12 C. Adamo and V.Barone, J. Chem. Phys. 110, 6158 (1999) and references therein. 13 S. Grimme, S. Ehrlich, and L. Goerigk, L. J. Comput. Chem. 32, 1456 (2011) and references therein 14 M. J.Frisch, et al., Gaussian09, Revision C.01wis4; Gaussian, Inc.: Wallingford, CT, 2009. Website: http://www.gaussian.com 15 H.-J. Werner et al., MOLPRO 2018.2, a package of ab initio programs. Cardiff, UK, 2018. Website: http://www.molpro.net. 16 M Saitow, U Becker, C Riplinger, E Valeev, F Neese, J. Chem. Phys. 146 164105 (2017). 17 Y. Guo, C. Riplinger, U. Becker, D.G. Liakos, Y. Minenkov, L. Cavallo, and F. Neese, J. Chem. Phys. 148, (2018). 18 F. Neese, F. WIREs Comput. Mol. Sci. 8, e1327 (2018) and references therein. 19 A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, and A.K. Wilson, Chem. Phys. Lett. 286, 243 (1998). 20 D.G. Liakos, M. Sparta, M.K. Kesharwani, J.M.L. Martin, and F. Neese, J. Chem. Theory Comput. 11, 1525 (2015). 21 N. Mardirossian and M. Head-Gordon, PCCP 16, 9904 (2014); J. Chem. Phys. 142, 074111(2015); ibid. 144, 214110 (2016). 22 Y. Zhao and D. G. Truhlar, J. Chem. Phys. 2006, 125, 194101; idem, Theor. Chem. Acc. 2008, 120, 215-241. 23 J. Tao, J. Perdew, V. Staroverov, and G. Scuseria, Phys. Rev. Lett. 91, 146401 (2003). 24 S. Kozuch and J.M.L. Martin, PCCP 13, 20104 (2011); J. Comput. Chem. 34, 2327 (2013). 25 G. Santra, N. Sylvetsky, and J. M. L. Martin, J. Phys. Chem. A, 123, ASAP (2019) (Leo Radom festschrift). DOI: 10.1021/acs.jpca.9b03157 26 E. Caldeweyher, S.Ehlert, A.Hansen, H.Neugebauer, S.Spicher, C.Bannwarth, S. Grimme, J. Chem. Phys. 150, 154122 (2019). 27 A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, and A.K. Wilson, Chem. Phys. Lett. 286, 243 (1998). 28 S. Kossmann and F. Neese, J. Chem. Theory Comput. 6, 2325–2338 (2010) and references therein 29 For a review, see: D. E. Bernholdt, Parallel Computing 26, 945 (2000) 30 F. Weigend, PCCP 8, 1057 (2006); F. Weigend, J. Comput. Chem. 29, 167 (2008). 31 C. Hättig, PCCP 7, 59 (2005); A. Hellweg, C. Hättig, S. Höfener, and W. Klopper, Theor. Chem. Acc. 117, 587 (2007). 32 R. C. Geary, Biometrika 27, 310-332 (1935) 33 A. Najibi and L. Goerigk, J. Chem. Theory Comput. 14, 5725 (2018). 34 L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi, and S. Grimme, PCCP 19, 32184 (2017).