Accurate Description of Catalytic Selectivity: Challenges and Opportunities for the Development of Density Functional Approximations

Abstract

Open AccessCCS ChemistryCOMMUNICATION1 Nov 2021Accurate Description of Catalytic Selectivity: Challenges and Opportunities for the Development of Density Functional Approximations Zhe-Ning Chen†, Tonghao Shen†, Yizhen Wang and Igor Ying Zhang Zhe-Ning Chen† Collaborative Innovation Center of Chemistry for Energy Materials, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, MOE Key Laboratory of Computational Physical Sciences, Shanghai Key Laboratory of Bioactive Small Molecules, Department of Chemistry, Fudan University, Shanghai 200433 State Key Laboratory of Structural Chemistry, Fujian Institute of Research on the Structure of Matter, Chinese Academy of Sciences, Fuzhou 350002 , Tonghao Shen† Collaborative Innovation Center of Chemistry for Energy Materials, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, MOE Key Laboratory of Computational Physical Sciences, Shanghai Key Laboratory of Bioactive Small Molecules, Department of Chemistry, Fudan University, Shanghai 200433 , Yizhen Wang Collaborative Innovation Center of Chemistry for Energy Materials, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, MOE Key Laboratory of Computational Physical Sciences, Shanghai Key Laboratory of Bioactive Small Molecules, Department of Chemistry, Fudan University, Shanghai 200433 and Igor Ying Zhang *Corresponding author: E-mail Address: [email protected] Collaborative Innovation Center of Chemistry for Energy Materials, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, MOE Key Laboratory of Computational Physical Sciences, Shanghai Key Laboratory of Bioactive Small Molecules, Department of Chemistry, Fudan University, Shanghai 200433 https://doi.org/10.31635/ccschem.020.202000635 SectionsSupplemental MaterialAboutAbstractPDF ToolsAdd to favoritesDownload CitationsTrack Citations ShareFacebookTwitterLinked InEmail An accurate description of catalytic selectivity poses an enormous challenge for theoretical simulations. Due to the absence of a well-defined benchmark set on the catalytic selectivity, the performance of even the widely used density functional approximations (DFAs) is yet to be validated. This work reports a test set based on the selective hydrogenation of α,β-unsaturated aldehydes catalyzed by ruthenium (Ru) hydride complexes. Special attention is paid to benchmark the regioselectivity of aldehydes to either unsaturated alcohols or saturated aldehydes. Accurate reference data were calculated by the massive parallel implementation of the coupled-cluster single, double, and perturbative triple excitations [CCSD(T)] approach, based completely on set limits. Furthermore, we performed the microkinetic simulation based on the CCSD(T) energy profiles, serving the most direct criteria for the performance of DFAs on catalytic selectivity of hydrogenation reactions. Using this test set, we uncovered the intrinsic difficulty of semilocal and hybrid functionals for such a purpose. In the context of the XYG3-type double hybrid (xDH) framework, we showed that this particular challenge could be addressed by the top rung functionals only when the many-body nondynamic correlation effect was accounted for adequately. A recently proposed xDH, namely scsRPA, showed unprecedented accuracy with only 5% error on average. The predicted kinetic selectivity by scsRPA is in close agreement with the reference value, revealing a unique versatility of the top rung DFAs for a reliable description of catalytic selectivity of hydrogenation reactions. Download figure Download PowerPoint Introduction Activity and selectivity are the most key concepts at the heart of catalysis.1–5 Rational designs, along with a good understanding of catalysts with the desired activity and selectivity, put forward a high demand on the theoretical multiscale kinetic study.6,7 To prepare the catalytic reaction network as the first-principle-based inputs for the followed microkinetic simulations, density functional approximations (DFAs) in the framework of (generalized) Kohn–Sham density functional theory (DFT) are the current methods of choice.8–11 However, the accuracy of the microkinetic study is known to be highly sensitive to the input data.12,13 Uncertainty remains for the reliability of the existing DFAs for treating chemical bonds, in particular, the transition metal–ligand (M–L) bonds,14–27 which play a prominent role in the activity and/or the selectivity of catalysts.28,29 To date, the majority of the transition-metal involved benchmarks focus only on the thermochemical and thermodynamic properties such as the M–L bonding energies, reaction energies, and/or reaction barriers. These are the relative energies, ΔE, characterizing each elementary reaction.14–27 However, these elementary reactions belong to a complicated catalytic reaction network, while the selectivity in catalysis is more of a matter of two or more competitive reaction pathways. To the best of our knowledge, there is no test set reported that explicitly aims to differentiate the competitive reaction pathways, which are the differences of the aforementioned properties, ΔΔE. The lack of accurate energetic and kinetic reference data from theory has prevented an unbiased benchmark of the intrinsic errors of DFAs for a reliable description of catalytic selectivity. Results and Discussion In this work, we have established a test set with accurate reference data to capture the fundamental bonding character in the selective hydrogenation of α,β-unsaturated aldehydes, which is not only of great industrial importance but also of broad academic interest.30–34 The schematic diagram of the corresponding catalytic reaction network is presented in Figure 1, for which a detailed description is provided, along with the listing of all reactions in Supporting Information Figure S1. Briefly, there exist essentially competing organic events among Reactions 1a–1c. The difficulty in achieving a high selectivity of C=O hydrogenation to unsaturated alcohols (UOL; Reaction 1b) can be inherently traced back to the thermodynamically more favorable side reactions yielding the corresponding saturated aldehydes (SAL; Reaction 1a) or fully saturated alcohols (SOL; Reaction 1c). The catalytic processes involve the elementary steps, as catalyzed by the transition-metal hydride species. By using four α,β-unsaturated aldehydes of acrolein, crotonaldehyde, cinnamaldehyde, and citral, in conjunction with the models of common metal hydride catalysts in hydrogenation, that is, [RuHCl(PH3)3(H2O)] and [RuH2(PH3)4],31,35–39 we have compiled a test set that includes 12 net hydrogenation reaction energies in the main-group chemistry ( MG12, ΔE associated with Reactions 1a–1c) and 50 Ru-mediated reaction energies ( TM50, ΔE associate with Reactions 2–9). Figure 1 | Typical processes for selective hydrogenation of α,β-unsaturated aldehyde catalyzed by transition-metal complexes. The upper part involves the main-group chemistry in the MG12 subset, while the lower part involves the transition-metal chemistry in the TM50 subset with 12 and 50 data points (ΔE), respectively. The corresponding SMG12 and STM16 subsets contain 12 and 16 pairwise energy differences (ΔΔE), respectively, to describe the regioselectivity. The model reactions have been elaborated in Supporting Information Figure S1, which are indexed by the same labels shown here. Codes for net reactions are UOL for C=O hydrogenation to UOL, SOL for full hydrogenation to SOL, and SAL for C=C hydrogenation to SAL. Download figure Download PowerPoint To set up a benchmark by specifically targeting catalytic selectivity, we then extracted 12 data points from the MG12 subset as pairwise differences ΔΔE of the hydrogenation energies among Reactions 1a–1c to depict the intrinsic regioselectivity of the problem ( SMG12). Similarly, the influence of transition-metal Ru on the regioselectivity could be quantified by the other 16 data points ( STM16). Furthermore, we performed microkinetic simulations on the relevant model reaction network (see Supporting Information for more details about the kinetic simulations). The resulting hydrogenation selectivity of cinnamaldehyde to the corresponding UOL and SAL provided the most direct criteria for the performance of DFAs on the catalytic selectivity of the target systems. To produce accurate reference data, in particular, for the energy differences representative to the intricate regioselectivity ( SMG12, STM16), and the predicted kinetic selectivity, we performed the calculations using the coupled-cluster single, double, and perturbative triple excitations [CCSD(T)] approach and extrapolated to the complete basis set (CBS) limit, which has long been recognized as the gold standard in quantum chemistry and has been used widely as the benchmark of DFAs for varying chemical and physical properties.40,41 Note that the average system size of the test set is ∼25 atoms per molecule, with the largest system consisting of 42 atoms, for which the canonical CCSD(T) calculations were extremely time-consuming. We employed the massive-parallel implementation of CCSD(T) in the Fritz Haber Institute ab initio molecular simulations (FHI-aims) package,40 to take advantage of the modern supercomputer powers without resorting to any local approximations to the correlation.41 We have used the FHI-aims software package for the single-point calculations in addition to molecular geometries for all DFT methods as well as the CCSD(T) reference. Thus, a strict comparison could be achieved for the intrinsic errors of DFAs and immunes from the influences due to the numerical inconsistency and the DFT errors in the geometrical optimization. However, it should be noted that the good performance of DFAs on energetic properties might not always (despite often) guarantee their accuracy for the geometry prediction. It is worthy of a careful benchmark on the influence of the geometry optimization for the complicated catalytic selectivity in the future. The computational details are found in the Supporting Information, while all calculation results (see Supporting Information Figures S1-S3, and Tables S1-S3, including the Cartesian coordinates for all systems are provided in the Supporting Information. Here, we first examined six DFAs, namely, PBE,39 SCAN,42 PBE0,43 M06L,44 M06,45 and M06-2X.45 These are widely used DFAs in chemistry and/or materials science in the hierarchy of DFAs46 from the second rung of generalized gradient approximation (GGA) through the third rung of meta-GGA to the fourth rung of hybrid (meta-)GGA. The former three DFAs corresponded to the nonempirical functionals, constructed by Perdew and co-workers42 to satisfy as many exact constraints as possible on each rung, and believed to be crucial for achieving a generally good description of diverse bonding systems. The latter three are empirical DFAs developed by Truhlar's group.45,46 They often contained several tens of empirical parameters, which were optimized directly with respect to some accurate experimental and/or theoretical thermochemical/thermodynamic properties for bonding systems of varying natures. To put the errors for thermochemistry and regioselectivity of both the main-group and the transition-metal-catalyzed systems on the same scale and ensure the scalability of their benchmark results from model studies to real complex systems, we took the mean absolute percentage error (MAPE) hereafter to compare the performance of the DFAs against the CCSD(T)/CBS reference. (See Table 1 for MAPEs of the six DFAs mentioned above. Supporting Information Table S3 presents more results with other DFAs, as well as those obtained using some other statistical descriptors.) The detailed computing reaction energies as well as the reference values can be found in Supporting Information Table S1 and S2. Table 1 | Percent MAPEa of 90 Calculated Results for 11 DFAs and CCSD against the CCSD(T)/CBS Reference Data for Each Subset and Total Test Set PBE SCAN PBE0 M06L M06 M06-2X MG12 9.9 3.1 7.1 9.2 4.6 3.9 SMG12 5.9 4.2 7.3 12.4 6.0 5.7 TM50 38.4 17.8 64.4 63.9 52.6 101.5 STM16 72.0 55.1 41.7 99.9 62.6 25.1 Overall 36.2 20.7 45.1 56.1 41.7 62.2 CCSD XYG3 XYGJ-OS xDH-PBE0 ZRPS scs-RPA MG12 5.9 2.2 1.4 2.5 2.6 2.5 SMG12 3.0 3.7 3.2 6.5 4.9 4.4 TM50 74.5 24.2 17.4 17.3 7.0 5.5 STM16 18.9 46.0 31.3 23.2 36.9 5.7 Overall 45.9 22.4 15.9 14.9 11.4 5.0 aMAPEs are given by 1 N ∑ i = 1 N | A i DF A − A i Ref | / | A i Ref | × 100 % with | A i DF A − A i Ref | denoting the absolute deviation of DFA results A i DF A with respect to the CCSD(T)/CBS reference A i Ref . To avoid undesired numerical noise that decreases the statistical value of MAPEs, we have replaced the denominator of | A i Ref | by 1 kcal/mol (chemical accuracy) if | A i Ref | <1 kcal/mol. The best three performers are indicated in bold. As shown in Table 1, the widely used DFAs are competent for main-group chemistry. Their MAPEs for MG12 are all below 10%, with the best performance provided by the state-of-the-art meta-GGA SCAN. The nonempirical family of DFAs showed better scalability for the regioselectivity with the main-group elements, providing similar or even better MAPEs for the SMG12 subset. For comparison, empirical DFAs, including M06L, M06, and M06-2X, all delivered degraded performances on SMG12, compared with MG12. The statistics presented in Table 1 show that both families of these lower rung DFAs underachieved for transition-metal chemistry. Again, the best performance on TM50 was provided by SCAN, resulting in a MAPE of 17.8%. However, the error of SCAN is more than tripled for STM16 when benchmarking on the transition-metal involved regioselectivity. The comparison on TM50 and STM16 subsets revealed that none of the commonly used DFAs on the lower rungs (no matter how nonempirical or highly parametrized) could provide a consistently well (or balanced) description for both TM50 and STM16 subsets. Carefully inspection of the benchmark results on TM50 suggested that the errors correlated positively with a fraction of the exact exchange, revealed by a notable increase of the errors from PBE to PBE0 and from M06 to M06-2X. However, the trend was reversed entirely with STM16. The best performance on STM16 (and the worst for TM50) was given by M06-2X, which contained a doubled amount of the exact exchange, compared with that of M06. Although there should still be plenty of room for improvement on the lower rung DFAs,47 a balanced and satisfactorily accurate description of the catalytic selectivity undoubtedly imposed a great challenge to the (hybrid) meta-GGAs, which was challenging to overcome by merely tuning the existing parameters in the functional formula. This observation was further confirmed by an extended benchmarking on other 15 (dispersion-corrected) DFAs on rungs 2–4 (see Supporting Information Table S3). An empirical dispersion correction improved the performance of PBE (but not for PBE0) on MG12, SMG12, and TM50, while showing less influence on STM16 for transition-metal involved regioselectivity. To study this unexpected challenge, we depicted a detailed distribution of the percentage errors for 28 regioselectivity relevant data in SMG12 and STM16 (Figure 2). Interestingly, none of the six lower rungs DFAs provided well-behaved error distributions. We did not see clear decay error frequencies in the regions with larger errors (see Figures 2a and 2b). M06 gave the best error distribution with 12 out of 28 test cases (12/28) falling within the −5% to 5% interval while having a great number of outliers (7/28) with percentage errors larger or lower than 25% or −25%. M06-2X shows a worsened error distribution. Considering that the best MAPEs for SMG12 and STM16 are given by M06-2X, it further confirmed the intrinsic difficulty to improve the lower rung DFAs for the reliable description of catalytic selectivity. Figure 2 | (a–d) Histogram of percentage errors of various methods for 28 reaction energy differences in SMG12 and STM16 (see Supporting Information Figure S2 for the error distributions for MG12 and TM50). Each vertical bar represents errors in the 10% range. The cases with errors larger or lower than 25% or −25% are grouped as outliers. Download figure Download PowerPoint Next, we examined a series of the XYG3-type double hybrids (xDHs) on the top-fifth rung of DFAs. The results from the wave-function method CCSD are also presented for comparison in Table 1. XYG3,48 XYGJ-OS,49 and xDH-PBE050 used the nonlocal correlation model in the form of the second-order perturbation theory (PT2) or the opposite-spin part only (osPT2) quantum state. The top rung correlation model designed for ZRPS explicitly explicitly introduces the electron-pair nondynamic correlations based on the osPT2 model51; while a recently proposed xDH, namely scsRPA, further included many-body nondynamic correlations in the context of spin-pair distinct adiabatic-connection fluctuation-dissipation theorem.52 The xDHs featured the use of self-consistent orbitals and density from a lower rung DFA with a full portion of the semi-local exchange correlation energy functional [i.e., B3LYP for XYG3 and XYGJ-OS; PBE0 for xDH-PBE0, ZRPS, and scsRPA]. They contained at most four empirical parameters in the functional formula and had the same (or even lower) computational scaling as that of the canonical PT2 method (N5 with N denoting the system size), which is formally one order higher than that of hybrid DFAs (N4), but one order lower than CCSD (N6). Top rung xDHs have been repeatedly proved to be unprecedentedly accurate in the description of the individual reactions in the main-group chemistry.48–50 As seen in Table 1, all five xDHs yield MAPEs lower than 3% for MG12. Their accuracy for the main-group chemistry involved the relative energy differences for pairs of competitive reactions on SMG12, which is also satisfactory. Except for xDH-PBE0 (6.5%), the resulting MAPEs were all <5%. It has been well recognized that the deficiency of the (os)PT2-based top rung functionals, including XYG3, XYGJ-OS, and xDH-PBE0, for transition-metal systems is closely related to the large nondynamic correlation error in the (os)PT2 approximation (see Supporting Information Table S3 for the MP2 results).53 Indeed, as shown in Table 1, MAPEs associated with these functionals were notably larger for TM50 than those for MG12, while their performances were even more degraded for STM16. Notably, the advanced wave-function method of CCSD yielded an even worse performance for TM50, whose MAPE was as high as 74.5%. However, the errors associated with CCSD were quite systematic, resulting in a much smaller, although still unsatisfactory MAPE of 18.9% for STM16, highlighting the difficulty toward a balanced description of catalytic selectivity. As shown in Table 1, ZRPS significantly reduced the MAPE for TM50 to only 6.6%, by virtue of including nondynamic correlation effect explicitly at the two-electron level, which, unfortunately, showed little help for transition-metal involved regioselectivity. The MAPE of ZRPS remained ∼40% for STM16. This error could not be effectively reduced until the many-body nondynamic correlation had been properly included in the context of the xDH framework, resulting in a well-balanced and outstanding performance of scsRPA on all four subsets, as shown in Table 1. Moreover, the well-behaved error distribution (Figures 2c and 2d) added weight to distinguish scsRPA as an accurate and robust method to address our challenge. There were 18/28 test cases falling in the −5% to 5% interval with no outlier! Our last benchmark was the kinetic selectivity of cinnamaldehyde toward the SAL and UOL products. The percentages of products are shown in Figure 3. The CCSD(T) energy profile predicted that the ratio between the SAL and UOL products was 7∶3. PBE, PBE0, SCAN, M06L, M06-2X, and XYGJ-OS were erroneously overestimating the selectivity of the UOL product. By contrast, M06, XYG3, xDH-PBE0, and ZRPS corrected the ordering but, notably, overshot the selectivity of the SAL product. Figure 3 revealed that the best performer was scsRPA, resulting in an almost identical ratio to that of CCSD(T). Figure 3 | Hydrogenation selectivity of cinnamaldehyde toward the SAL and UOL products. The model reaction network contains two competitive hydrogenation processes, as shown in Supporting Information Figure S3. Kinetic simulations were conducted based on the energy profiles at 300 K (see Supporting Information Section 2 for more details to account for the steric effect and the rate equations used in the kinetic modeling). The initial concentrations of cinnamaldehyde and catalysts were set to be 1.0 and 0.1 mol/L, respectively. Numeric time integral continued until all cinnamaldehyde reactants were consumed. The dotted line is to guide the eyes, showing a 1∶1 selectivity to SAL and UOL. The blue ovals indicate that the product ratio from scsRPA is almost identical to that of CCSD(T). Download figure Download PowerPoint Conclusion We have developed the first test set to study the capability of DFAs for the description of catalytic selectivity quantitatively. This goal was achieved based on the Ru-catalyzed selective hydrogenation of α,β-unsaturated aldehydes using the calculation results from our massive parallel implementation of the CCSD(T)/CBS method as references. We found that semilocal and hybrid DFAs did not provide a satisfactory description for the test set, revealing the intrinsic difficulty encountered in improving lower rung DFAs for a reliable description of the catalytic selectivity. The excellent performance of scsRPA demonstrated the great potential of top rung DFAs toward a universally higher accuracy and broader applicability. These characteristics were achieved when the many-body nondynamic correlation effects were properly taken into account to go beyond the commonly used DH functionals for transition-metal involved chemistry. We expect that the test set and the relevant benchmark study in this work would serve as the first step toward quantitatively understanding the performance of DFAs for the catalytic selectivity of hydrogenation reactions and could motivate more constructive test sets for other kinds of catalytic selectivity. Supporting Information Supporting Information is available and includes a detailed description of the test set, computational methods, kinetic simulations, calculation results of CCSD, CCSD(T), along with various DFAs, Cartesian coordinates of the systems in the test set, and the histogram of percentage errors for MG12 and TM50. Conflict of Interest The authors declare no conflicts of interest. Funding Information This work was supported by the National Natural Science Foundation of China (nos. 21688102, 21973015, 21973094, and 91427301), the Science Challenge Project (no. TZ2018004), an of local in and a key of the of Shanghai (no. and the Natural Science Foundation of Fujian (no. The authors to for and of this 1. in for Catalysis of Catalytic 3. for and of the of for 1, of in to Activity and the Computational of 1, Computational in between and Computational and Computational to and A to and Computational for of and of and Uncertainty in Kinetic of Catalysis and An of DFT and for the of to of Density for and An of from Density Functional to and DFT of the Structure and of in the of Density Functional for Reaction as of the of Density Functional for with Set of Reactions for the of Density of Density for of the of Density Functional to with Catalytic Activity toward by We Reactions by of Density Functional for by a Set of Accurate DFT for for for Density Functional a Density Functional The to Hydrogenation of Hydrogenation of and and Catalytic In The of of A Reaction for of Wang of A in Catalysis Hydrogenation of by of on the Molecular of and on the of the Catalytic Hydrogenation of on the of the Hydrogenation of by of the C=C Hydrogenation of by in the in Hydrogenation of A Molecular with Numeric Zhang of the in the Numeric for Molecular Perdew and Density Density Functional without The PBE0 Density Functional for and M06 of Density for and and of and 12 Perdew of Density Functional Approximations for the Institute of of Density Functional in Computational An and of Density Zhang Density Functional for Accurate of and Zhang Density Functional to Using a Zhang Density Functional xDH-PBE0 from a Zhang Perdew Density for and Zhang of and in Density Functional A Zhang the Top of of Density Functional the of and Information Chinese functional hybrid authors to for and of this

Topics & Concepts