High-Throughput Experimentation, Theoretical Modeling, and Human Intuition: Lessons Learned in Metal–Organic-Framework-Supported Catalyst Design

We have screened an array of 23 metals deposited onto the metal–organic framework (MOF) NU-1000 for propyne dimerization to hexadienes. By a first-of-its-kind study utilizing data-driven algorithms and high-throughput experimentation (HTE) in MOF catalysis, yields on Cu-deposited NU-1000 were improved from 0.4 to 24.4%. Characterization of the best-performing catalysts reveal conversion to hexadiene to be due to the formation of large Cu nanoparticles, which is further supported by reaction mechanisms calculated with density functional theory (DFT). Our results demonstrate both the strengths and weaknesses of the HTE approach. As a strength, HTE excels at being able to find interesting and novel catalytic activity; any a priori theoretical approach would be hard-pressed to find success, as high-performing catalysts required highly specific operating conditions difficult to model theoretically, and initial simple single-atom models of the active site did not prove representative of the nanoparticle catalysts responsible for conversion to hexadiene. As a weakness, our results show how the HTE approach must be designed and monitored carefully to find success; in our initial campaign, only minor catalytic performances (up to 4.2% yield) were achieved, which were only improved following a complete overhaul of our HTE approach and questioning our initial assumptions.


■ INTRODUCTION
Metal−organic frameworks (MOFs) constitute a class of porous crystalline materials that can act as finely tuned heterogeneous catalyst supports. 1−5 They can enhance catalysis through confinement or containment effects, 25,6 support effects, high surface area, and framework modulation of molecular transport. 7,8 These all represent phenomena that are not present in traditional inorganic metal oxide supports. Through the deposition of catalytically active metal sites onto MOF nodes via atomic layer deposition in MOFs (AIM) 9−12 or surface organometallic chemistry (SOMC), 13 the porous MOF NU-1000 14 has proven to be an effective support for catalytic sites capable of ethylene hydrogenation 10,15 and oligomerization, 16−18 propane oxidative dehydrogenation, 19 alkene epoxidation, 20 and propyne isomerization. 9 imprNonetheless, the design of novel MOF catalysts is quite challenging: deposited active sites in MOFs are by their nature quite difficult to characterize, 21 and the active site and catalytic activity are known to change depending on deposition technique. 19 As such, conventional computational approaches may struggle to convincingly describe the underlying chemistry, although some success has been demonstrated. 10,18,22 In particular, MOF-deposited nanoscale catalysts, 23 especially for highly fluxional metals such as Cu, 11,21,24−27 provide a very difficult challenge for theory despite their versatile reactivity. On top of this, the activity of MOF catalysts can vary heavily due to synthetic or operating parameters, such as metal weight loading, reduction temperature, reaction temperature, space velocity, and operating pressures�all difficult parameters to model via standard theoretical means (i.e., density functional theory (DFT)). Thus, high-throughput experimentation (HTE) appears well-positioned to play an important role in MOF catalyst design, as it can rapidly explore chemical space with minimal (i.e., data-driven or machinelearned) theoretical input. Despite this, high-throughput approaches are quite rare in experimental studies of metal− organic frameworks 28−31 and in catalytic applications of MOFs. 9 In this work, we present the first application of this highthroughput approach to developing MOF-supported catalysts. By empirically learning the interaction between synthetic and operating conditions over 23 different metals, we use datadriven algorithms to identify and optimize interesting catalytic activity prior to in-depth characterization and theoretical modeling. Focusing on the reaction of propyne dimerization, important for the removal of alkyne compounds from polymerization feedstocks, 33,34 we are able to improve hexadiene yields in flow reactors from 0.4 to 24.4% on the MOF NU-1000 ( Figure 1). Characterization of the final optimized catalysts reveals the necessity of Cu nanoparticle formation in the MOF for achieving significant catalytic activity, which is supported by theoretical results.
However, these results were not obtained without significant intellectual effort�initial assumptions regarding the nature of metal active sites and the range of relevant experimental variables proved incorrect, resulting in over half the experiments conducted giving no significant hexadiene yield. Improvement was only obtained through questioning our initial assumptions and a redesign of our high-throughput campaign, with the full story only revealing itself through careful posthoc characterization and modeling of the successful catalyst. Thus, our results show how HTE is not a fully automated approach�although useful for finding unexpected catalytic activity, high-throughput campaigns often must be monitored and redesigned to find success.

■ RESULTS AND DISCUSSION
Our initial search space consisted of 22 different metals, due to their known propensity to catalyze alkyne functionalization and upgrading (Ag, Al, Ce, Cu, Fe, Ga, In, Ir, La, Mo, Ni, Pd, Pt, Re, Rh, Ru, Sm, Sr, Ta, Ti, V, and Zn), deposited onto NU-1000 at five different weight loadings (0.5−5 wt %) and tested in a 16-reactor parallel plug flow reactor (Flowrence from Avantium). In addition to the metal identity and loading, five different reaction variables were also identified. These included two pretreatment variables related to the reducibility of the organometallic complexes to the active metal site, hydrogen reduction temperature (150−225°C) and steamed air treatment temperature (125−225°C) used to remove residual organic material, and operating conditions that influence conversion and selectivity including total flow (2.5−20 mL/ min), which in turn will affect the space velocity, reaction temperature (100−250°C), and the amount of H 2 cofed during the reaction (0− 5 vol %). These initial parameters were chosen based on domain knowledge 9,16,35−38 and in consideration of the thermal and chemical stability of NU-1000. 15,39 Given this search space, our initial selection algorithm employed Bayesian optimization 40−45 to maximize the total hexadiene yield (defined as the ratio of moles of hexadiene formed to the number of moles of propyne that have been consumed) for each metal individually, prioritizing In and Cu over other metals due to literature precedent. 9 A summary of this approach is given in Figure S1. However, after several iterations of this approach totaling 721 experiments (taking roughly 6 months of work), the results showed very little activity or selectivity for any metal in the parameter space outlined, obtaining a maximum hexadiene yield of 4.2% over 4 wt % Cu/NU-1000 (H 2 reduction T = 150°C, T rxn = 250°C, space velocity = 2000 mL/min/g cat , 5 vol % H 2 cofed) and little improvement of predictions with number of iterations. Furthermore, Cu was far-and-away the best-performing catalyst, with the next-highest yields obtained with Ga (1.0%) and In (0.4%). Thus, a decision was made to abandon the initial approach.
Prior studies have shown that an increase in H 2 partial pressure in the inlet results in a higher selectivity to propylene, and that oligomers are preferentially produced at lower H 2 :C 3 H 4 ratios in rather large concentrations, [35][36][37]46 and yet, this trend was not observed in our initial results. Instead, the highest hexadiene yields were obtained at the maximum 5 vol % H 2 concentrations. Based on this unexpected trend, we made the decision to screen all metals (slightly biased toward Cu, based on its larger activities) at a higher concentration of H 2 (40 vol %, SI Section 1.2). To our surprise, this immediately boosted hexadiene yields on Cu from 4.2 to 15.4%. This key result allowed us to redesign our highthroughput campaign to focus on Cu activity and higher inlet H 2 concentrations (0−80%) while eliminating pretreatment variables (air temperature and H 2 reduction temperature), which did not seem to have an impact on hexadiene yields. Figure 2 summarizes the approach used in our redesigned campaign. Instead of optimizing conditions for all metals symmetrically, a dual-model scheme was employed, in which model 1 selects operating parameters (loading, reaction temperature, inlet H 2 concentration, and space velocity) for Cu and model 2 selects other metals and metal loadings to test at these operating conditions. This approach enabled both (i) a more efficient use of the high-throughput reactor bed and (ii) the possibility of learning trends across different metals with model 2. To achieve this latter goal, the adsorption of propyne and H 2 onto different single-atom catalysts deposited onto NU-1000 was computed with DFT. Thermochemical features of these adsorptions (e.g., ΔE ads , ΔH ads , ΔG ads ) were combined with atomic data 47,48 of different metals to interpolate between different metals. Full details of these models are available in the Supporting Information.
Using this approach, yields on Cu were increased from the 4.2% found in our initial campaign to 24.4% over the course of 227 experimental trials at 79 unique operating conditions (with several conditions tested multiple times to ensure reproducibility). Furthermore, significant yields of 7.3 and 7.0% were obtained for Pd and Ni, respectively, showing the ability of our second approach to discover interesting reactivity on new metals while simultaneously optimizing Cu. Interestingly, the highest hexadiene yields for Pd and Ni were found at considerably different conditions. Pd had the highest yield at 5 wt %, 200°C, 460 mL/min/g cat , and 13.5 vol % H 2 cofed, while Ni had the highest yield at 3 wt %, 175°C, 162 mL/ min/g cat , and 38.5 vol % H 2 cofed. These in turn also differ dramatically from Cu, which exhibited high yields with much higher space velocities, between 1000−2000 mL/min/g cat , and lower metal loading. Figure 3 shows different Voroni tesselations (K-nearestneighbor models with K = 1) of the maximum hexadiene yield (green), conversion (red), and selectivity (blue) achieved in two-dimensional cuts of experimental parameter space (metal loading, H 2 vol %, space velocity, and reaction temperature) for the 79 trial-averaged Cu experiments carried out in our second campaign. Each set of plots demonstrates one trend found in the final data set. The first set of plots, showing yield as a function of loading and temperature ( Figure 3a) and H 2 vol % and space velocity (Figure 3b), emphasizes the sparsity of the catalytic activity observed in the experimental parameter space. High catalytic activity for Cu requires a highly specific reaction temperature and loading as well as space velocity and H 2 vol %, and small changes in any of these parameters can drastically affect results. Interestingly, a trend in catalytic activity appears when H 2 vol % and space velocity are increased in turn with each other (i.e., increasing one or the other alone results in diminished yields). Figure 3c,d demonstrates the trade-off observed between conversion and selectivity as loading and temperature are changed. Higher conversions are noted at the intersection of high temperature and high metal loading. However, higher selectivity for hexadienes is observed at higher dispersion of metals. Thus, the optimal yields occur at the combination of high temperatures (200°C) and medium metal loadings (2 wt %). This demonstrates the importance of considering conversion and selectivity in tandem for the optimization of yield, as the conditions that promote formation of one product may not be thermodynamically favorable. Interestingly, this trade-off appears to occur differently for Ni and Pd, which is discussed in the Supporting Information ( Figure S7).
Finally, Figure 3e shows the highly dependent behavior of the catalytic system (M + MOF) on both loading and H 2 vol %. These dependencies point strongly toward nanoparticle-like catalytic behavior, in which one needs a large enough nanoparticle so that two propyne molecules can react on its surface but small enough such that excess H 2 is not adsorbed to over hydrogenate the reactants. If one has in mind separated metal centers in the MOF, which carry out catalysis independently from one another, then yield should increase roughly linearly with loading, up to a point at which the increased loading results in sintering or other destructive processes. Instead, a sharply peaked activity at a loading of 2 wt % was observed. This implies the existence of a strongly optimal loading, which is difficult to explain without including interactions between Cu sites in one model, most likely through nanoparticle formation.
Indeed, postreaction microscopy analysis of the Cu/NU-1000 catalysts shows the formation of nanoparticles, which is unique to the Cu/NU-1000 catalysts and not exhibited for the other deposited metals. SOMC produces highly dispersed metal species on supports, and XRD analysis of Cu/NU-1000 prior to reaction or preactivation indicates that the Cu species are below the limit of detection for XRD for a catalyst containing 1.52 wt % Cu by inductively coupled plasma-optical emission spectroscopy (ICP-OES, nominal loading of 1.0 wt %). Migration of isolated Cu sites has been previously reported on NU-1000 8,11,21 where the size of the nanoparticle formed is constrained by the size of the hexagonal (3 nm) and triangular (1 nm) channels of the pore when installed through AIM. Here, the Cu particle size appears to be manipulated by controlling the duration of the time, during which the catalysts are reduced under H 2 .  (Figure 4e). During the reaction, the particle size continues to increase due to cofed H 2 , increasing to an average of 13.9 nm (Figure 4d). This nanoparticle size greatly exceeds the pore size of the MOF, implying the formation of Cu nanostructures spanning several MOF pores. Both the effect of particle size and inlet H 2 concentration on turnover frequency (TOF) at 200°C were explored (Figure 4f). In all cases, larger Cu nanoparticles and higher H 2 vol % cofeed resulted in a more active catalyst, up to 40%, at which point the TOF was observed to decrease slightly. No other higher oligomers were observed as evidenced by the H balance close to 100% ( Figure  S11). Yield for hexadiene is greater than that of Cu nanoparticles deposited onto SiO 2 via SOMC ( Figure S13), despite the presence of coke on the Cu/NU-1000 ( Figure  S15). This higher yield appears to be the result of NU-1000 stabilizing larger nanoparticles�a high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) image of Cu/SiO 2 taken after the reaction shows smaller particle size ( Figure S14). Further studies are being conducted to determine the growth kinetics of these nanoparticles and whether nanoparticles may exist within the pores of the MOF or decorate the exterior.
Finally, to rationalize the high performance of Cu nanoparticles in comparison to single-atom (low% loading) Cu catalysts, reaction mechanisms for propyne dimerization and hydrogenation were computed using DFT ( Figure 5). Specifically, dimerization vs hydrogenation reaction pathways on a single-atom catalyst (SAC) model of Cu/NU-1000 ( Figure 5, left) were compared to a nanoparticle (NP) model ( Figure 5, right, estimated as a slab). These two models are  Figures S1 and  S2), hydrogenation is favored by 48 kJ mol −1 on the SAC model and by only 35 kJ mol −1 on the NP model. This suggests that the mechanistic explanation for the higher hexadiene yield observed for larger Cu loadings and H 2 cofeed % is that higher H 2 cofeed assists in formation of NPs in NU-1000, which in turn have increased selectivity toward hexadienes compared to SAC Cu/NU-1000, thereby increasing the overall hexadiene yield. Beyond this, coverage and confinement effects 8 of the MOF may also play important roles in engendering hexadiene selectivity.

■ CONCLUSIONS
In this work, we have screened an array of 23 metals deposited onto NU-1000 for the catalytic dimerization of propyne to hexadienes under multiple reaction conditions. Using datadriven algorithms and high-throughput experimentation, hexadiene yields on Cu/NU-1000 were improved in series from 0.4 to 15.4 to 24.4%. Characterization of the catalyst reveals that selectivity toward hexadiene is likely due to large Cu nanoparticles, embedded in the MOF framework, which grow in size when exposed to H 2 under reductive conditions. This activity is supported by DFT calculations comparing single-atom models of the Cu active site to nanoparticle models, which predict nanoparticles to be more selective for dimerization to hexadienes.
To the best of our knowledge, this project serves as the first application of the HTE approach to post-modification catalysis in MOFs and provides key takeaways for further research in this direction. It is clear that the improved activity obtained in this synergistic study could not have been discovered a priori; the high-throughput approach was critical to finding the precise experimental and synthetic conditions at which the catalysts showed meaningful selectivity and conversion to hexadiene. The discovered trend of large amounts of cofed H 2 drastically increasing hexadiene yields goes against both literature precedent [35][36][37]46 and chemical intuition, which dictates that lower H 2 :C 3 H 4 ratios should facilitate hexadiene formation. This trend could only be explained via post-hoc analysis with microscopy and DFT, which shows that high H 2 % cofeed encourages nanoparticle formation, which is more selective for propyne dimerization than single-atom copper sites. Furthermore, our initial single-atom models of the active site proved at odds with the final nanoparticle catalysts, and even if our DFT models had been more in line with this final characterization, it is unlikely that any ab initio model could satisfyingly incorporate the effects of reaction temperature, metal loading, H 2 cofeed%, and space velocity, all of which proved vital for optimizing hexadiene yield. It is worth noting that features of these single-atom models were utilized in the machine learning scheme developed to select new experiments (Figure 2), which likely hindered their contribution to improving catalyst yield.
On the other hand, this does not mean the HTE approach is void of complications. Designing and initiating a HTE campaign requires committing to a parameter space ahead of time which may or may not contain significant results. Although one attempts to design their campaign wisely to maximize the chances of success, using literature precedent and chemical intuition to select the choices and ranges of the experimental variables, initial assumptions can be proven wrong. For example, literature precedent dictated that a small H 2 vol % would encourage hexadiene selectivity, which led us to select an initial parameter range of 0−5 vol % H 2 . This unfortunate initial choice led to about six months of work (over half the experiments conducted in this project) resulting in very low hexadiene yields. It required us questioning our initial assumptions and following the data to drastically expand the range of hydrogen concentrations being tested in order to find significant catalytic activity. This overhaul of our HTE approach also provided an opportunity to drop insignificant variables and redesign our experimental selection algorithm, which greatly increased the efficiency of our second campaign. Furthermore, HTE alone cannot be used to explain the catalytic activity it finds�the data must be analyzed, characterization must be done, and the active site must be modeled to provide a full and satisfying picture. We suggest that high-throughput campaigns in general may benefit from a series of low-commitment exploratory searches followed by increasingly targeted and refined searches once interesting behavior has been found.
Thus, our story shows that although data-driven HTE can successfully be used as a tool to discover novel and highperforming materials, it is by no means a fully automated approach. Much care and thought must be put into the design and operation of the campaign throughout its lifetime, and the campaign must be monitored and redesigned if it is not providing satisfactory results. Of course, the parameter space of any HTE campaign can in principle be expanded ad inf initum to include all of chemical space (at exponentially increasing cost), but in reality, some (potentially ill-founded) assumptions must be made, and chemical intuition must be applied to design an initial search space and monitor its success. Human intuition and input will continue to be vital to the success of improving catalytic yields. Finally, we hope that this journey shows that it is important to share both successful and unsuccessful results to advance knowledge as a community. ■ METHODOLOGY Metal Deposition. NU-1000 was synthesized as previously described 9,32 and used as prepared. An automated synthesis platform (CM3 Core Module deck, Unchained Laboratories Inc.) housed in a custom-built N 2 -filled glovebox (MB 200B, MBruan) was used for catalyst synthesis. The CM3 performed both solid and liquid dispensing within a 0.5% tolerance. A surface organometallic chemistry (SOMC) process was used for grafting metals onto NU-1000. 13 This synthetic approach was used over ALD 9 and other cation exchange techniques 49 due to its generalization to many different precursors and feasibility in a high-throughput environment. In detail, 20 mg of NU-1000 was first dispensed into 4 mL vials. Various organometallic precursors (trimethylaluminum, gallium(II) acetylacetonate, trimethylindium, titantium tetraisopropoxide, t e t r a k i s ( e t h y l m e t h y l a m i n o ) v a n a d i u m ( I V ) , b i s -(ethylcyclopentadienyl) manganese, bis(N,N′-di-t-butylacetamidinato) nickel(II), bis(dimethylamino-2-propoxy) copper-(II), diethylzinc, bis(t-butylimido) bis(dimethylamino) molybdenum(VI), methylcyclopentadienyl)(1,5-cylooctadiene) iridium(I), (trimethyl)methylcyclopentadienyl platinum(IV), bis(ethylcyclopentadienyl) ruthenium(II), rhodium(II) acetylacetonate, palladium(II) hexafluoroacetylacetonate, trimethylphosphine(hexafluoroacetylacetonato) silver(I), tris(2,2,6,6-tetramethyl-3,5-heptanedionato) lanthanum(III), tris(2,2,6,6-tetramethyl-3,5-heptanedionato) samarium(III), pentakis(dimethylamino) tantalum(V), or methyltrioxorhenium(VII)) were then dispensed into each vial correlating to either 0.5, 1.0, 2.0, 4.0, or 5.0 wt % metal, after which 4 mL of dried toluene was added. The vials were then placed manually onto a heated unit atop a shaker plate and left to metalate for 72 h at 60°C and 400 rpm. After cooling, the samples were washed in toluene and centrifuged (Speedvac Concentrator, SPD121P, ThermoElectron) five times. Solvent exchange using pentane was done in the last washing step. After removal of the supernatant fluid, the vials were then removed from the glovebox and vacuum-dried at 60°C overnight.
High-Throughput Screening. Dimerization reactions were performed in a 16-reactor fixed bed flow system (Flowrence, Avantium). Typically, 5 mg of as-prepared catalyst was loaded into a quartz reactor (inner diameter = 2 mm, outer diameter = 3 mm, length = 300 mm). Reactions were performed between 100 and 250°C at a heating rate of 10°C/ min at 1 atm. All gases were purchased from Airgas. Propyne (2 vol %) in Ar was used with pure H 2 in various ratios for hydrogenation reactions. All catalysts were reduced in 100 vol % H 2 for 2 h at 200°C. For each data point collected, the catalysts were first exposed to the propyne/H 2 /Ar mixture for 1 h before measurements, after which two or three measurements were taken per condition tested, and reported values are the average of these data points. N 2 (UHP) was added to sweep gas to the gas chromatograph (GC) after reaction. He (UHP) was used as an internal standard. The effluent of each reactor was collected in a unique vial for each reactor to avoid cross contamination of effluent streams. The effluent was then analyzed with GC (7890B, Agilent Technologies), equipped with a thermal conductivity detector (TCD) and two flame ionization detectors (FIDs). The gas products identified were methane, ethane, ethylene, propane, propylene, propadiene, acetylene, 2,3-dimethyl-1,3-butadiene, t,t-2,4-hexadiene, c,c-2,4-hexadiene, c,t-2,4-hexadiene, t-2-methyl-1,3-pentadiene, and 2,4-hexadiyne.
Microscopy. Nanoparticle imaging was conducted using a FEI Talos F200X high-angle annular dark-field scanning transmission electron microscope scanning transmission electron microscope (HAADF-STEM) operated at 200 kV. This microscope was used in coordination with the Center for Nanoscale Materials (CNM) at Argonne National Laboratory. In preparation for analysis by electron microscopy, powdered samples (∼20 mg) were sonicated in ethanol (10 mL) for 15 min. The resulting suspension was dropcast onto a lacey carbon TEM grid (Ted Pella, Inc., UC-A on Lacey 400 mesh Cu). Particle size was measured using ImageJ software.
Single-Atom Catalyst Model. A cluster model representing a single Zr 6 O 8 (OH) 4 (H 2 O) 4 node of NU-1000 was extracted from a previously optimized unit cell of NU-1000. The hydroxy and aqua ligands on one face of the node were used to deposit different catalytic transition metals, varying in their oxidation state from +1 to +4. The transition metals were ligated to the appropriate number of hydroxy ligands to maintain the overall neutrality of the MOF. 51 The linkers connecting this node of NU-1000 to the adjacent nodes were truncated as formate groups. The carbon atoms of the formate groups were spatially frozen to mimic the structural rigidity of NU-1000. These cluster models representing single-atom catalytic metals supported on NU-1000 (M-NU-1000) were then employed to compute the adsorption energies of the reactants (propyne, H 2 ) on the catalytic metal centers and to study the reaction mechanism for propyne dimerization occurring on a single-atom Cu-NU-1000.
The above cluster models were optimized at the density functional theory level in the Gaussian 16 software 52 program using the M06-L 53 local exchange-correlation density functional. The def2SVP basis set was used for O, C, and H, while the def2TZVPP basis set was used for the catalytic metals and the support metal (Zr). 54 The SDD effective-core pseudopotential was used for metals present beyond the fourth row of the periodic table to enhance the computational efficiency. Vibrational frequencies were computed in the rigid rotor harmonic oscillator (RRHO) approximation for the optimized structures to determine the nature of the stationary point and to compute the entropic contribution to the free energies. No imaginary frequencies were determined for the intermediates (minimum on the potential energy surface), while only one imaginary frequency along the reaction coordinate was determined for the transition states. All frequencies below 50 cm −1 were corrected to 50 cm −1 to avoid errors due to anharmonicity. 55,56 Different possible spin states corresponding to the oxidation state of each transition metal were considered, and only the spin states lowest in energy are reported here.
Nanoparticle Model. Periodic density functional theory calculations to study the reaction mechanism of propyne dimerization on Cu NP model were performed in the Vienna Ab Initio Simulation software (VASP 5.4.4). 57,58 The Cu nanoparticles formed in NU-1000 upon pretreatment in H 2 were modeled using a 5 × 5 × 3 Cu(111) periodic slab (75 Cu atoms) with a lattice constant of a = 2.522 Å. Such periodic models have previously been used to model large Cu nanoparticles (NPs) in NU-1000. 50 An additional vacuum of at least 18 Å was used along the z-direction. The Perdew− Becke−Ernzerhof (PBE) 59 exchange-correlation density functional was used along with Grimme's D3 correction with zero damping for dispersion energy correction. 60 A plane-wave basis set with a cutoff energy of 400 eV was employed. All structures were optimized using an energy and force convergence criteria of 10 −4 eV and 0.08 eV/Å respectively. A 3 × 3 × 3 Γ-centered k-points grid was used for the Brillouin zone sampling. The climbing image nudged elastic band (CI-NEB) method 61 with at least eight images between the optimized reactant and product along the reaction coordinate was used for locating the transition state. All calculations were performed without any spin polarization. The electronic energies of the optimized structures were used to represent the potential energy interface along the reaction cycle. The relative energies were calculated using eq 1. The energies of the reference gas phase molecules were obtained by placing a single gas phase molecule at the center of a simulation box having the same dimensions as the Cu NP model. The absolute and relative electronic energies are reported in Table S5