Hydrogen Activation by Silica-Supported Metal Ion Catalysts: Catalytic

Dec 5, 2018 - Materials Science Division, Argonne National Laboratory , 9700 South ... functionals have the lowest mean absolute deviation (MAD) value...
0 downloads 0 Views 2MB Size
Subscriber access provided by Gothenburg University Library

A: Spectroscopy, Molecular Structure, and Quantum Chemistry

Hydrogen Activation by Silica Supported Metal Ion Catalysts: Catalytic Properties of Metals and Performance of DFT Functionals Cesar Plascencia, Larry A Curtiss, and Cong Liu J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b08340 • Publication Date (Web): 05 Dec 2018 Downloaded from http://pubs.acs.org on December 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Hydrogen Activation by Silica Supported Metal Ion Catalysts: Catalytic Properties of Metals and Performance of DFT Functionals Cesar Plascenciaa,b, Larry A. Curtissa*, and Cong Liuc* a) Materials Science Division, Argonne National Laboratory, 9700 South Cass Ave., Lemont IL 60439 USA b) Department of Chemistry, Michigan State University, 220 Trowbridge Rd, East Lansing, MI 48824 USA c) Chemical Sciences and Engineering Division, Argonne National Laboratory, 9700 South Cass Ave., Lemont IL 60439 USA

*Corresponding authors: [email protected] (CL), [email protected] (LAC) Abstract: Single-site heterogeneous catalysts (SSHC) have received increasing attention due to their well-defined active sites and potentially high specific activity. Detailed computational studies were carried out on a set of potential SSHC’s i.e., silica supported metal ions, to investigate the reactivity of these catalysts with H2 as well as to evaluate the performance of density functional theory (DFT) methods in conjunction with triple-ζ quality basis sets (i.e., cc-pVTZ) on reaction energetics. The ions considered include 4d and 5d metals as well as several post-transition metal ions. A representative cluster model of silica is used to calculate reaction free energies of the metal hydride formation that results from the heterolytic cleavage of H2 on the M-O bond. The hydride formation free energy is previously shown to be strongly correlated with the catalytic activity of such catalysts for alkene hydrogenation. ONIOM calculations (CCSD(T)//MP2) are used to assess the accuracy and reliability of the MP2 results and it is found that MP2 is a suitable level of theory for gauging the performance of DFT functionals. The performance of various DFT functionals is assessed relative to MP2 results and it is found that the wB97xd and PBE0 functionals have the lowest standard deviation (STD) value while the MN12SX and PBE functionals have the lowest mean absolute deviation (MAD) values. The B3LYP functional is shown to have the similar MAD and STD values as the top performing functionals. Potential active SSHC’s for exergonic hydrogen activation predicted in this study include mostly late and post transition metal ions, i.e., Au3+, Pd2+, Pt4+, Pd4+, Ir4+, Hg2+, Rh3+, Pb4+, Tl3+, In3+, Ir3+, Os4+, Cd2+, Ru2+, and Ga3+. This study provides important guidance to future computational studies of such catalyst systems.

Introduction Single-site heterogeneous catalysts (SSHCs) that consist of a metal atom, ion, or small cluster of atoms, held by surface ligands to a rigid framework, have received increasing attention.1 The most salient features of such catalysts2–4, in addition to potentially high specific activity due to the isolated active site and higher recyclability relative to homogeneous catalysts, are the spatially well-defined nature of its active sites. These characteristic features allow for comprehensive structure-function investigations that can reveal kinetic and mechanistic information about the active site; information that can be used to derive rational design principles for either improving current SSHC’s or designing neoteric ones.1–8 The rational design of SSHC’s based on its active site properties is supported by studies which show that it is typically the few atoms composing the SSHC’s active site and its immediate support structure, often low-coordinated and unsaturated sites, that directly contribute to the catalytic activity of the SSHC.9–15 These active sites are typically composed of metals or organometallic complexes that are distributed across high-area surfaces and can exist either as clusters, as a coordination 1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 28

compound, or as single atoms/ions.16–23 Previous studies have shown that supported singlesite metal atoms/ions heterogeneous catalysts (SAHC) could possess favorable properties relative to SSHC’s with supported clusters or nanoparticles as the active sites.6,24 Among the choices for the support materials for SSHC’s are silica16,17,20–23 and zeolite25,26 supports as well as metal organic frameworks (MOF)24,27–29 and 2D structures30. While MOF’s and 2D systems allow for a greater control of the homogeneity of the SSHC’s active site due to their crystalline-type or regularly ordered structures,31–35 silica remains attractive due to its high surface area, thermal stability and structural flexibility. Additionally, silica is earth-abundant and readily binds36 with a wide selection of elements. The focus of this computational study are SAHC’s that use silica as the supporting material. Specifically, the SAHC’s are composed of bare metal ions on a silica support20,22 that are synthesized through the use of strong electrostatic adsorption (SEA) methods37,38 or surface organometallic chemistry (SOMC)5 techniques and are catalytically active for the hydrogenation of alkenes. Modeling the amorphous nature of the silica support requires the careful application of computational chemistry modeling techniques. The morphology and topology of a silica support is highly dependent on the manner in which it is prepared,15,36 and consequently, silica can have siloxane rings of varying sizes. The modeling techniques available to describe such an amorphous structure have been discussed in previous studies.15,39 In one procedure, the slabs are generated either from quartz40 or melted glass structures41 with a similar density of surface species as the target silica. In another procedure, the amorphous silica structure is obtained through the optimization of starting silica structure with siloxane rings of varying size.42 The first procedure does not readily reproduce the amorphous nature of silica, whereas the second procedure quickly becomes computationally intractable due to the substantial number of conformations available to the starting structure with varying siloxane ring sizes. In addition, the modeling technique chosen by Martynowycz et. al39 involves modeling the silica support as a 2D film; however, due to the regular ordering of a 2D film, this technique is incapable of reproducing silica’s amorphous structure. An alternative modeling technique is to use cluster models.43 Cluster models represent the local structure relevant to the property being investigated through the truncation of a larger, periodic structure. Truncation of the larger, periodic structure necessarily removes the presence of long-range interactions. However, several studies on SSHC that use a silica support show that cluster models of such systems recover similar energetics as slab models,15,42,44–46 indicating that the structural and electronic effects of these silica supported SSHC’s are mostly confined within a local boundary. For this study, the cluster approach is adopted. Silsesquioxane cages are chosen to model the silica support as silsesquioxane cages have been previously shown to be good cluster representations of amorphous silica.47 Accurate and reliable quantum chemistry methods are required to arrive at results that ultimately can be used to derive rational design principles for SSHC’s or SAHC’s. While wavefunction (WF) methods based on the Hartree-Fock approximation offer such accuracy and reliability, their computational cost scales exponentially with respect to the basis set size and thus quickly becomes computationally intractable for larger systems.48 For this reason, density functional theory (DFT) is usually chosen since its computational cost is lower and can reproduce qualitative trends. While there has been an attempt to create a systematic series of DFT functionals with increasing accuracy (“Jacob’s ladder”)49, choosing an appropriate functional is still not black box. Instead, choosing the appropriate functional for a target species relies on either the similarity of the target species with the training set used to derive the functional or the existence of benchmark studies that prove that the functional performs well with the target species.50 While there have been previous benchmark studies for the electronic and geometric properties of silica51, transition metals (TM)52,53, or SSHC’s54, to the knowledge of the authors, no study exists that measures the performance of DFT on SAHC’s with silica 2 ACS Paragon Plus Environment

Page 3 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

supported metal ions. As such, the performance of DFT functionals on the target systems of the present study must be explored. Previously, we have explored the catalytic alkene hydrogenation of five ions (In3+, Ga3+, Zn2+, Mn2+, Ti4+) supported by silica using a combined computational and experimental study.55 Propylene hydrogenation was used as a probe reaction to explore the catalytic properties of these SAHC systems. Since olefin hydrogenation is the microscopic reverse reaction of alkane dehydrogenation, the study of the catalytic hydrogenation reactions can provide useful insights into reaction intermediates and structure reactivity relationships while having favorable thermodynamics under mild conditions.21 In addition, similar to propylene hydrogenation, many hydrogenation reactions (e.g., CO2 hydrogenation) are initiated by the activation of H2 by the catalyst.56,57 Thus, the study of hydrogen activation could provide insights into the catalytic properties of the supported metal ions for hydrogenation reactions of other molecules. In the aforementioned combined computational and experimental study on SAHC’s,55 it was discovered that for such catalysts, the hydrogen activation undergoes a heterolytic cleavage pathway to form a metal hydride intermediate. The study also found that the formation of this hydride intermediate is critical for the activity of alkene hydrogenation. An activity-descriptor relationship was established where experimentally-derived reaction rates were correlated with computationally-derived reaction free energies of metal hydride formation at 200 °C, hereafter referred to as ΔGM-H,200 °C; The ΔGM-H,200 °C values are calculated at 200 °C to mirror experimental conditions. Essentially, the more negative the ΔGM-H,200 °C values, the more active the silica supported metal ions are for the alkene hydrogenation (See Figure 1). This activity-descriptor relationship can be used to predict the catalytic activity of any silica supported single-site metal ions given that ΔGM-H,200 °C is known. Therefore, in this work, the ΔGM-H,200 °C of a broad range of metal ions was investigated; including the five ions previously studied (In3+, Ga3+, Zn2+, Mn2+, and Ti4+) along with Pb, Tl, and 4d and 5d TM ions. The ΔGM-H,200 °C is calculated using twelve DFT functionals that are ultimately calibrated against ΔGM-H,200 °C calculated using ab initio WF methods. From these results, the performance of DFT functionals is measured, the model chemistry that is used to derive the activity-descriptor relationships is evaluated, and novel SAHC’s using bare transition metal ions are proposed.

Figure 1. Activity-descriptor relationship of propylene hydrogenation on M/SiO2 systems at 200 °C from study by Liu et al.55 The curve is a fitted correlation between the calculated ∆GM-H and the log of the experimental turn-over frequencies (TOFs). Figure reproduced using data from ref[55].

Methodology First principles calculations based on DFT and WF methods have been performed on a series of silica cluster models to validate and predict their catalytic activity for olefin hydrogenation. The structure of the silica cluster is taken from a recent study by Das, U. et al.,15 in which the authors 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 28

explored the effect of siloxane ring strain on the catalytic properties of single-site catalyst. The study demonstrates that a cluster model composed of six siloxane rings, in fact the cluster model used in this study (Figure 2), is sufficient to reproduce results from a periodic slab model; and, validates that the electronic and steric effects of silica are local and can be adequately represented by a cluster model. Furthermore, a previous study47 shows that silica clusters composed of different size siloxane rings (silsesquioxane) are a good approximation of an amorphous silica surface. The series of model catalysts considered in this work includes the five ions previously studied (In3+, Ga3+, Zn2+, Mn2+, Ti4+) along with the 4d and 5d metals as well as Tl and Pb. All calculations are performed using the Gaussian 09 software.58

Figure 2. Example of silsesquioxane cage used in the study, before the addition of the metal ion active site. The silsesquioxane cage used for this study has six siloxane rings. The oxygens on the hydroxyl groups can form bonds with the metal ion via surface organometallic chemistry.

The silica clusters are designed to model SAHC’s surfaces where the catalyst is a bare TM ion. The stability of possible surface species/sites for these type of silica clusters have been previously investigated by our group55, and it was found that under catalytic conditions low-coordinated surface sites where the metal is exposed is either more thermodynamically favorable or is in equilibrium with higher-coordinated surface sites that have absorbed water or hydroxyl molecules. For this reason, the silica clusters used in this study all feature the exposed metal site as the active site. Consequently, fourcoordinated active sites are the maximum coordination number that this model can achieve. In this study, only 3-coordinated and 4-coordinated bare-ion active site models are considered. The silica active site cluster models were first screened at the B3LYP59–62/CEP-31G63 level of theory. At this level of theory, all the active site cluster models were optimized at various possible oxidation states, coordination numbers, and multiplicities. The most common oxidation states for each ion were considered, based on literature research of reported coordination compounds. Only the lowest energy structures were used for higher-level calculations. The geometries of the lowest energy structures were then further optimized at the B3LYP/cc-pVTZ(-PP)64–69 level of theory; 4d and 5d metal atoms are modeled with the pseudopotential based aug-cc-pVTZ-PP basis set while all other atoms are modeled with the all-electron aug-cc-pVTZ basis set. The optimized structures at both the B3LYP/CEP-31G and B3LYP/cc-pVTZ(-PP) are characterized through a vibrational analysis to ensure the optimized structure corresponds to global minima. The B3LYP/cc-pVTZ(PP) optimized active site clusters were then used to obtain the hydrogenated cluster models. This is done by adding a hydrogen atom to the metal site as 4 ACS Paragon Plus Environment

Page 5 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

well as to an adjacent oxygen. These hydrogenated cluster models were then reoptimized at the B3LYP/cc-PVTZ(PP) level and characterized through a vibrational analysis to ensure the optimized structure corresponds to global minima. Note, that the coordination number increases by one upon hydrogenation and the multiplicity of the active site structure is assumed to be the same for the hydrogenated cluster model. The free energy model used in this work is the standard statistical mechanical model employed by Gaussian 09. For both the active site and hydrogenated clusters models, corrections to the zero-point energy (ZPE) and for Gibbs free energy are calculated at the B3LYP/ccpVTZ(-PP) level of theory. These corrections were added to the calculated total energies using different functionals and methods. Thus, the errors of each method only result from the calculations of total energies. The optimized structure at the B3LYP/cc-pVTZ(-PP) level of theory for both the active site cluster model and the hydrogenated active site cluster model is what is used for all subsequent singlepoint calculations using other DFT functionals and WF methods. The level of theory chosen for the optimizations is supported by a recent benchmark study51 on siloxane and siloxane derivatives that demonstrated that B3LYP with a triple-ζ basis set lead to good structural and thermochemical results. Note that similar to the original study by Liu et.al55, The bottom half of the cluster is frozen during optimization (6 H, 11 O and 6 Si atoms). Single point calculations are used to calculate ΔGM-H,200 °C, i.e. reaction free energies of metal hydride formation at 200 °C, based on the following equation: (1) In Eq (1), the metal hydride results from the heterolytic cleavage of H2 on the M-O bond.23,55,70–73 Our previous experimental and theoretical studies have explored both heterolytic and homolytic addition of H2 for all 3d metals, for selected 4d metals, and selected post-TM ions.55 It was found that that all of these SSHC’s go through heterolytic H2 cleavage on the M-O bond. Therefore, in this paper we assume that heterolytic addition mechanism is dominant for the metal ions presented. Values for ΔGM-H,200 °C are calculated using a series of DFT functionals and WF based methods. An effort was made to benchmark at least one functional in each rung of “Jacob’s ladder”; The series of functionals include: the local spin density approximation (LSDA) functional SVWN561,74,75; the generalized gradient approximation (GGA) functionals BLYP60,76, PBE77,78 along with their hybrid functional counterparts B3LYP59–62 and PBE079; the meta-GGA functional TPSS80 along with the corrected version for solid-state calculations revTPSS81,82; the meta-GGA functional M06-L83 along with its hybrid functional counterpart M0684; The empirical dispersion corrected GGA functional B97D85; The range-separated hybrid functional MN12-SX86; and finally, the long-range corrected functional with empirical dispersion corrections wB97xd87. The aim is to calibrate the performance of the DFT functionals using the WF methods as the standard. The WF methods used are the second-order Møller-Plesset perturbation theory (MP2)88, and coupled-cluster theory with inclusion of single and double excitations and a perturbative treatment of the triple excitations{CCSD(T)}.89,90 The CCSD(T) method was used in a two-layer ONIOM (“our own nlayer integrated molecular orbital and molecular mechanics”)91–94 model to assess the necessity of highlevel electron correlation methods and to validate the accuracy of the MP2 energies. The two-layer ONIOM method allows for the partitioning of the chemical system into two layers (the inner and outer layer), each layer calculated at a different level of theory. Generally, the inner layer is treated with a higher level of theory relative to the outer layer. The energy is then “extrapolated” according to the following equation: (2)

𝐸𝑂𝑁𝐼𝑂𝑀 = 𝐸𝑚𝑜𝑑𝑒𝑙ℎ𝑖𝑔ℎ +𝐸𝑟𝑒𝑎𝑙𝑙𝑜𝑤 −𝐸𝑚𝑜𝑑𝑒𝑙𝑙𝑜𝑤

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 28

The “real” system refers to the entire cluster while “model” refers to the inner layer. In this study, the “model” system is capped using hydrogen as the link atom. The manner in which to partition the system has been previously studied and guidelines exist for general use.95,96 The partitioning has been studied for silica clusters in particular and it has been demonstrated that link atom distances are not crucial for obtaining correct results97, therefore default G09 values are used for link atom distances for all calculations.

Figure 3. Example ONIOM 2-layer model for clusters where the metal ion is three-coordinated. Atoms in ball-and-stick model are inner layer. Atoms in tube model are outer layer. The cluster on the left is the active site model while the cluster on the right is the model for the hydrogenated active site.

Figure 4. Example ONIOM 2-layer model for clusters where the metal ion is four-coordinated. Atoms in ball-and-stick model are inner layer. Atoms in tube model are outer layer. The cluster on the left is the active site model while the cluster on the right is the model for the hydrogenated active site.

In this study, the inner layer for the silica clusters in which the TM is three-coordinated includes the TM and any oxygen and hydroxyl groups that are directly coordinated to the TM (See Figure 3). Partitioning schemes that use a larger inner layer in silica systems have been previously studied by Roggero, et al.,98 where it was demonstrated that there are no significant improvements to the energy when using a larger inner layer than that used for the current study. A different layering scheme is used for the four-coordinated silica models, these models contain a strained 4-membered ring which includes the TM. Thus, a larger inner layer was used that included the strained ring to avoid breaking the bonds involved in the strained ring (See Figure 4). Several inner layer models were tested for one of these four-coordinated silica models and the model performance was assessed via the S-value test99 (values available in supplementary material in 6 ACS Paragon Plus Environment

Page 7 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table S1). The S-value measures the adequacy of the partitioning scheme used to divide the inner from the outer layer. (3)

𝑆𝑙𝑒𝑣𝑒𝑙 = 𝐸𝑟𝑒𝑎𝑙𝑙𝑒𝑣𝑒𝑙 −𝐸𝑚𝑜𝑑𝑒𝑙𝑙𝑒𝑣𝑒𝑙

Due to the way the ONIOM method is derived, the error (ΔS) in the overall ONIOM model can be measured and is related to the S-values previously mentioned. (4)

𝛥𝑆 = 𝐸𝑂𝑁𝐼𝑂𝑀 −𝐸𝑟𝑒𝑎𝑙ℎ𝑖𝑔ℎ = (𝐸𝑟𝑒𝑎𝑙ℎ𝑖𝑔ℎ −𝐸𝑚𝑜𝑑𝑒𝑙ℎ𝑖𝑔ℎ )−(𝐸𝑟𝑒𝑎𝑙𝑙𝑜𝑤 −𝐸𝑚𝑜𝑑𝑒𝑙𝑙𝑜𝑤 ) = 𝑆ℎ𝑖𝑔ℎ −𝑆𝑙𝑜𝑤

The inner layer chosen for the four-coordinated silica model is based on the S-value results, ΔS results, and chemical intuition. Values for ΔGM-H,200 °C are calculated with two ONIOM schemes for each active site. In one scheme, CCSD(T) is used for the inner layer while MP2 is used for the outer layer. In the other, MP2 is used for the inner layer while B3LYP is used for the outer layer. Since the geometries as well as the initial ΔGM-H,200 °C were obtained using the B3LYP functional, comparisons of the two ONIOM schemes and the MP2 results allow for the calibration of the initial study by Liu et al.55 as well as for the assessment for the need of high level electron correlation methods for the type of silica clusters included in this study.

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 28

Results and Discussion The reaction free energies of metal hydride formation at 200 °C , i.e. ΔGM-H,200 °C, were calculated using Eqn. 1 for a series of TM and post-TM containing SAHC cluster models using both DFT and WF methods. ΔGM-H,200 °C is an important activity descriptor for alkene hydrogenation, identified in our previous study on propylene hydrogenation using SAHC.55 The activity-descriptor relationship showed a strong correlation between the calculated ΔGM-H,200 °C and the experimental hydrogenation activity, i.e., the lower the ΔGM-H,200 °C, the higher activity of the catalyst (See Figure 1).55 Therefore, the catalytic activity of novel catalysts for alkene hydrogenation can be predicted by calculating the ΔGM-H,200 °C on the catalyst active site. In the present study, we not only make predictions of the catalytic properties of different metal ions on silica for alkene hydrogenation, but also examine the performance of ONIOM models, DFT functionals, and the MP2 method on trends and accuracy of ΔGMH,200 °C. The investigation proceeds as follows: First, to assess the necessity of high-level electron correlation wave-function methods above the MP2 level, the performance of two 2-layer ONIOM models is compared against MP2 results. Second, the DFT predicted ΔGM-H,200 °C trends are examined and potential active SAHC’s are identified with various functionals being tested. Third, to gauge the validity of ΔGM-H,200 °C trends predicted by the various functionals, the MP2 ΔGM-H,200 °C trends are examined and compared to the DFT ΔGM-H,200 °C trends. Finally, the results from MP2 calculations are used to calibrate the performance of various functionals for the prediction of ΔGM-H,200 °C. Because all the 3d TM ions were studied in our previous work,55 in this study we consider mostly 4d and 5d TM ions, as well as selected 3d and post-TM ions. 1. Performance of ONIOM Models against MP2 Results. Two ONIOM models are tested in this study. The partitioning of the inner and outer layers for the two models is as follows: [MP2/cc-pVTZ(-PP):B3LYP/cc-pVTZ(-PP)] which is shortened to {ONIOM(MP2)} and [CCSD(T)/cc-pVTZ(-PP):MP2/cc-pVTZ(-PP)] which is shortened to {ONIOM(CC)}. Due to the costly nature of coupled-cluster calculations, only a subset of cluster models is calculated using the {ONIOM(CC)} model. The difference between MP2 and {ONIOM(MP2)} values correspond to the measure of the {ONIOM(MP2)} model error. This relationship is described by equation 4 which is reproduced below: (4) 𝛥𝑆 = 𝐸𝑂𝑁𝐼𝑂𝑀 −𝐸𝑟𝑒𝑎𝑙ℎ𝑖𝑔ℎ = (𝐸𝑟𝑒𝑎𝑙ℎ𝑖𝑔ℎ −𝐸𝑚𝑜𝑑𝑒𝑙ℎ𝑖𝑔ℎ )−(𝐸𝑟𝑒𝑎𝑙𝑙𝑜𝑤 −𝐸𝑚𝑜𝑑𝑒𝑙𝑙𝑜𝑤 ) = 𝑆ℎ𝑖𝑔ℎ −𝑆𝑙𝑜𝑤 Since the lower level for {ONIOM(MP2)} is B3LYP, ΔS for {ONIOM(MP2)} also represents the additional correlation energy recovered by using MP2 on the inner layer. There is no straightforward definition for the difference between MP2 and {ONIOM(CC)}; However, since MP2 is the lower level for {ONIOM(CC)} then the difference between MP2 and {ONIOM(CC)} represents additional correlation recovered by using CCSD(T) on the inner layer.

The ΔGM-H,200 °C values calculated using the ONIOM models and MP2 are shown as a bar graph in Figure 5. The bar graph shows that the trends between ΔGM-H,200 °C calculated using MP2 and {ONIOM(CC)} are similar. There are some changes in trend order between MP2 and {ONIOM(CC)} namely for 3c-Ag+, 4c-Ir4+, and 3c-Pd2+. These metal ions contain half-filled and full d-shell configurations that most likely require post-HF methods due to the existence of differential correlation effects.100 The discussion on differential correlation effects is expanded in a later section, but for now it is sufficient to note that the differential correlation effects might be responsible for the change in trend order. The change in trend order for these systems indicates that CCSD(T) is recovering additional correlation effects that is not being recovered by MP2. However, this change in trend order is not 8 ACS Paragon Plus Environment

Page 9 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

associated with a large shift in ΔGM-H,200 °C magnitude. For example, 3c-Ag+ is out of order for MP2 with a ΔGM-H,200 °C magnitude of 7.4 kcal mol-1, however its nearest neighbors have a ΔGM-H,200 °C magnitude of 9.7 kcal mol-1 and 8.7 kcal mol-1 so the change in order does not significantly change the predicted catalytic activity. On average, the use of the {ONIOM(CC)} method affects the ΔGM-H,200 °C values by ~ ±2 kcal mol-1 relative to the MP2 values, however these differences can get as high as 8 kcal mol-1. Thus, in general, using CCSD(T) for the inner layer in {ONIOM(CC)} has small effects on ΔGM-H,200 °C and does not significantly change the order of the MP2 trend. The previous discussion is evidence that MP2 is sufficient to recover the necessary correlation effects present in the SAHC cluster models tested.

Figure 5. Trends for ΔGM-H,200 °C calculated with {ONIOM(CC)}, MP2, and {ONIOM(MP2)} methods for a subset of cluster models. The data is ordered according the results from the {ONIOM (CC)} model. Units in kcal/mol.

Further analysis of Figure 5 results in a comparison between MP2 and {ONIOM(MP2} results. There are some changes in trend order between MP2 and {ONIOM(MP2)} namely 3c-Ag+, 4c-Rh4+, and 3c-Os2+. Again, it is observed that the metal ions that have a change in trend order contain half-filled and full d-shell configurations. The difference between MP2 and {ONIOM(MP2)}, i.e. ΔS, is further analyzed in Figure 6 which is a Letter-Value (LV) plot of the ΔS values. The data of Figure 6 is given in Table 1. Letter-value plots are similar to box-plots, but with improved description of the tail region for small data sets and the inclusion of more percentile ranges.101 A notable feature of the LV plot in Figure 6 is that it contains only negative values. This indicates that the {ONIOM(MP2)} model calculations tend to predict more exergonic values relative to MP2. Since the only difference between the two is the use of the B3LYP functional, the overprediction is directly attributable to the B3LYP functional and is most likely caused by known issues with metals for the LYP portion of the functional.102 According to Figure 6, 50% of the errors lie between -8 kcal mol-1 and -11.5 kcal mol-1 while 75% of the errors lie between -7 kcal mol-1 and -12.5 kcal mol-1. While the average ΔS value is ~ -10 kcal mol-1. These average values represent both the average value for the failure of {ONIOM(MP2)} to reproduce MP2 results as well as the average value of additional correlation energy recovered when using MP2. However, it is important to note that the ΔS value can get as low as ~-24 kcal mol-1. The large ΔS value along with the changes in the trend order for {ONIOM(MP2)} indicate that at least the MP2 level of theory is necessary to get correct magnitudes and trends for ΔGM-H,200 °C. There are even more changes in the trend order when comparing {ONIOM(MP2)} with {ONIOM(CC)}, however, these changes also correspond to species with half-filled or full d-shell configurations. From all the evidence presented so far, including the similarities among the ONIOM and MP2 trends, the similarity in the magnitude of ΔGM-H,200 °C produced by {ONIOM(CC)} and MP2, and the stark improvements made when going from {ONIOM(MP2)} to MP2; it is reasonable to conclude that 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 28

the MP2 level is sufficient to recover the relevant effects present in the SSHC cluster models, and consequently that it can be used as reference method against which all other DFT results can be calibrated. Table 1. ONIOM Model Error, ΔΔGM-H,200 °C, for ONIOM model {ONIOM(MP2)}. The ONIOM model error gauges the effectiveness of the ONIOM model to reproduce results using only the high level method, i.e. MP2. Units in kcal/mol. Statistical information is visualized in Figure 6.

structure

ΔΔGM-H,200 °C(kcal mol−1)

structure

ΔΔGM-H,200 °C(kcal mol−1)

3c-Ru2+

1.6

4c-Nb4+

10.2

4c-Mn4+

2.6

4c-Mo4+

10.3

3c-Zn2+

4.8

3c-In3+

10.3

3c-Tc3+

6.7

3c-Mo2+

10.4

4c-Ti4+

6.9

3c-Rh3+

10.4

3c-Mn3+

6.9

4c-Hf4+

10.8

3c-Re3+

7.5

3c-Pd2+

11.0

3c-Os2+

7.7

3c-Tl3+

11.0

3c-Mo3+

7.8

4c-Tc4+

11.2

3c-Ru3+

7.8

4c-Tl3+

11.6 12.0

3c-Ta3+

7.9

4c-Os4+

4c-Zr4+

8.3

4c-Pd4+

12.0

3c-Mn2+

8.5

3c-Ag1+

12.2

3c-Pb2+

8.6

3c-Au3+

12.3

3c-Ir3+

8.9

4c-Pt4+

12.7

3c-Cd2+

9.1

4c-Ta4+

13.1

4c-Ru4+

9.4

4c-Ir4+

13.1

3c-Ga3+

9.7

4c-Re4+

13.3

3c-Hg2+

9.8

4c-W4+

18.6

3c-Nb3+

10.0

4c-Rh4+

24.1

4c-Pb4+

10.1

2. DFT Trends The predicted ΔGM-H,200°C trends produced by various functionals are shown in Figure 7. The SVWN5 functional generally results in the lowest ΔGM-H,200 °C for each SAHC model i.e. for each active site, SVWN5 predicts a stronger M-H bond relative to the rest of the functionals. However, this result is to be expected as the SVWN5 functional has the well-known issue of over-binding.103 LSDA-type functionals, such as SVWN5, are considered the first rung in the so called “Jacob’s Ladder” of functionals.80 Sophisticated functionals that contain corrections to the LSDA functional are the subsequent rungs; where each successive rung offers a systematic improvement on the previous rung. Despite this, the distribution of DFT calculated ΔGM-H,200 °C for any particular SAHC model does not follow the “Jacob’s Ladder” scheme. There is, however, a general symmetry in the ordering. For example, for any particular metal ion, the revTPSS calculated ΔGM-H,200 °C regularly is the highest value compared to those calculated with other functionals. Similarly, the MNS12SX functional regularly produces a ΔGM-H,200 °C value that is close or at the median value for the distribution of DFT calculated 10 ACS Paragon Plus Environment

Page 11 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ΔGM-H,200 °C for any metal ion. The ordering for the rest of the functionals is not exact. However, there are again general regularities. For example, the M06 and M06L functionals are usually near the bottom of the distribution of DFT calculated ΔGM-H,200 °C since the resulting ΔGM-H,200 °C are relatively more exergonic. In contrast, the TPSS and B97D functional are usually near the top of the distribution of DFT calculated ΔGM-H,200 °C since the resulting ΔGM-H,200 °C are relatively more endergonic.

Figure 6. Letter-Value plot of the ONIOM Model Error, ΔΔGM-H,200 °C, calculated for ONIOM model {ONIOM(MP2)}. The ONIOM model error gauges the effectiveness of the ONIOM model to reproduce results using only the high level method, i.e. MP2. Units in kcal/mol. The values for ΔGM-H,200 °C calculated with DFT range from ~ -45 kcal mol-1 to ~ 55 kcal mol-1. According to the activity-descriptor relationship described by Liu et. al.55, ΔGM-H,200 °C is strongly

correlated with the catalyst activity for propylene hydrogenation. Essentially, a lower i.e. more exergonic ΔGM-H,200 °C value corresponds to a higher catalytic activity (See Figure 1). Thus, the relative hydrogenation activities among different metal ions are highly dependent on the trend of the calculated ΔGM-H,200 °C. Generally speaking, all the DFT functionals predicted similar trends of ΔGM-H,200 °C as a function of metal ion. Thus, the DFT predicted ΔGM-H,200 °C values are partitioned by metal type—early TM vs late TM vs post-TM—in Table 2 in an effort to isolate interesting trends (No notable trends observed based on formal oxidation state or coordination number, see supplementary information Tables S7 and S8). From Table 2, it is observed that all the early TM ions (Groups 3-7) tend to have highly endergonic ΔGM-H,200 °C values, regardless of which DFT functional was used. This suggests that these silica supported early TM ions have very little to zero catalytic activity. In contrast, most of the post-TM ions (Groups 13 & 14) are predicted to have an exergonic ΔGM-H,200 °C value, regardless of the DFT functional used. This indicates that post-TM ions could be good candidates for propylene hydrogenation. Indeed, the work by Liu et. al.55 predicted and confirmed the favorable catalytic activity of the novel In3+/SiO2 catalyst for propylene hydrogenation. The only post-TM SAHC that is predicted by all DFT functionals to have poor catalytic activity is the 3c-Pb2+ SAHC. In contrast, the 4c-Pb4+ model is predicted to have high catalytic activity. A literature search indicates that Pb(II) has been previously adsorbed onto various type of silica104–106 while there is no evidence for Pb(IV) silica. This indicates that Pb(II) could have greater stability on silica compared to Pb(IV). The more stable Pb(II) on silica, however, would engender a higher ΔGM-H,200 °C, simply because ΔGM-H,200 °C is calculated relative to the free energy of the catalyst active site (Eqn. 1). While there are no synthesized or predicted structures of Pb(IV) on a surface reported in the literature, there are reports of Pb(IV) metallocanes107. The presence of weak transannular interactions stabilize the Pb(IV) metallocanes albeit to a highly distorted structure. In essence, the structure and stability of the Pb4+ metallocane is very sensitive to donor-acceptor type transannular interactions. The siloxane rings present in the silica could provide similar donor-acceptor stabilizing effects.108,109 However, the feasibility of the synthesis of a Pb4+/SiO2 SAHC is currently not accessible to the authors. Although different DFT functionals show similar trends of ΔGM-H,200 °C as a function of metal ion, it is difficult to determine potential SAHC candidates based solely on the results from the DFT functionals. For example, in Table 2, the cluster models containing late TM ions (Groups 8-12), have 11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 28

DFT predictions ranging from -44.2 kcal mol-1 to 44.7 kcal mol-1. Furthermore, the number of exergonic and endergonic ΔGM-H,200 °C results is different for each DFT functional. For example, SVWN5 predicts that all but three cluster models containing late TM’s will have an exergonic ΔGM-H,200 °C value. In contrast, revTPSS predicts that only five cluster models containing late TM’s will have an exergonic ΔGM-H,200 °C value. The consistent trends predicted by all the DFT functionals for early TM and post-TM ions are not present for the late TM ions. Therefore, a reference method is needed to calibrate the DFT results in order to determine which of the late TM ions are likely to possess favorable catalytic activity. The calibration will also help confirm the previous discussion and analysis of the DFT results for the early and post-TM ions. The reference method chosen is MP2 due to its affordability and accuracy. The MP2 calculated ΔGM-H,200 °C values will help determine which cluster models possess exergonic or very low endergonic ΔGM-H,200 °C values that in turn helps identify potential catalytically active SAHC’s.

Figure 7. DFT predicted ΔGM-H,200 °C trends. A MP2 ΔGM-H,200 °C trendline is overlaid on the DFT points. Any deviation from the MP2 trend can then be visually and directly compared. Units in kcal/mol. Table 2. DFT predicted ΔGM-H,200 °C trends. DFT Trends partitioned by metal type. Units in kcal/mol. NOTE: ΔGM-H,200 °C = G in this table. Early – Group 3-7, Late – Groups 8-12, Post-TM – Groups 13,14. type

structure

G(B3 LYP)

G(B97 d)

G(BL YP)

G(M06 )

G(M06 L)

G(MN12 SX)

G(PBE )

G(PB E0)

G(SV WN5)

G(TPSS)

G(revTP SS)

G(wB97x d)

Early TM

3c-Tc3+

14.7

13.9

19.5

4.8

6.6

10.8

14.2

10.0

5.4

18.5

19.8

8.5

3c-Mn3+

16.7

16.1

21.2

6.6

11.5

14.1

16.1

12.1

8.3

21.4

23.0

12.1

2+

3c-Mo

18.9

18.8

21.8

13.7

14.1

18.7

18.1

15.7

11.8

22.7

24.8

15.9

4c-Mn4+

23.0

25.3

27.1

17.1

20.9

21.7

23.6

19.8

15.4

28.2

30.4

19.5

3c-Re3+

29.4

31.6

33.2

25.5

25.4

25.4

29.8

26.3

22.8

33.4

36.3

26.9

4c-Re

29.5

32.1

33.0

27.7

26.6

29.9

29.5

25.9

22.7

33.1

35.8

26.6

3c-Mn2+

30.1

30.5

33.1

25.8

27.1

31.4

30.1

27.6

25.3

34.1

35.9

28.3

4+

12 ACS Paragon Plus Environment

Page 13 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry 3c-Mo3+

31.2

32.1

34.6

29.6

28.6

29.4

31.9

28.8

25.2

35.1

37.3

28.2

4c-Tc

31.7

33.2

35.0

27.7

29.4

31.4

31.4

28.6

22.6

36.4

38.9

27.6

3c-Y 3+

33.3

35.4

35.7

32.3

32.4

34.9

33.9

31.7

29.6

37.5

38.7

31.4

4c-Mo4+

35.1

36.6

37.9

30.5

32.2

35.0

35.3

33.0

27.6

39.3

41.3

31.8

4c-Zr4+

35.8

38.7

38.4

34.1

34.1

37.0

36.7

34.4

30.6

39.6

41.1

34.0

4c-Nb4+

37.4

39.6

40.3

34.9

35.0

36.9

37.6

35.1

30.6

40.8

42.7

34.8

3c-Nb3+

37.6

37.6

40.7

34.0

33.5

35.2

37.5

34.8

31.0

41.0

42.9

34.6

4c-Hf4+

37.7

41.4

40.3

35.1

34.8

38.5

38.5

36.3

32.4

41.2

42.5

36.4

4c-W4+

39.2

41.4

42.4

34.7

35.7

38.4

40.0

37.1

32.8

43.4

45.6

36.2

4c-Ti4+

39.5

42.4

41.9

36.7

36.8

39.0

39.9

37.9

33.7

43.3

45.0

37.9

4c-Ta4+

40.0

43.0

42.8

37.0

36.4

39.4

40.7

38.2

34.4

43.5

45.4

38.1

3c-Ta3+

47.7

49.7

50.7

44.1

43.9

46.5

48.2

45.7

42.2

51.0

52.8

44.7

3c-Au3+

-29.2

-26.0

-21.9

-35.2

-29.9

-32.7

-28.8

-35.5

-44.2

-23.1

-20.3

-36.3

4c-Pt4+

-22.2

-19.9

-17.1

-27.7

-24.7

-24.9

-22.3

-26.6

-34.5

-16.5

-13.8

-27.4

3c-Pd2+

-18.1

-18.7

-15.6

-26.0

-22.6

-15.1

-20.9

-21.9

-34.4

-13.3

-10.7

-21.9

4c-Pd4+

-14.5

-12.6

-9.1

-23.1

-18.1

-16.4

-15.1

-19.8

-29.5

-8.6

-5.8

-20.5

3c-Hg2+

-10.1

-8.3

-7.9

-10.7

-10.2

-10.8

-10.7

-12.3

-16.8

-4.5

-1.2

-11.1

3c-Rh3+

-5.5

-4.9

-0.3

-12.0

-8.4

-9.1

-5.9

-10.7

-15.3

0.0

2.4

-11.6

4c-Ir4+

-5.3

-3.8

-0.9

-10.6

-8.8

-7.9

-5.9

-9.7

-17.1

-0.3

2.2

-10.4

3c-Ag1+

2.3

3.3

2.0

4.2

5.5

5.2

-0.6

1.1

-10.5

6.0

9.3

3.4

3c-Cd2+

2.9

5.2

4.8

0.3

1.7

2.7

2.0

1.0

-5.8

7.9

10.4

1.1

3c-Ir3+

4.9

4.4

9.6

-1.5

-0.3

-0.2

4.1

0.1

-6.5

9.0

11.1

-0.7

3c-Ru2+

6.0

5.1

8.1

0.2

-1.6

0.6

3.8

2.9

-4.3

9.1

11.5

3.2

3c-Zn2+

7.9

9.3

9.2

2.9

3.9

6.6

6.5

6.2

-0.1

12.0

14.1

6.0

4c-Os4+

8.0

9.9

11.7

3.9

4.7

6.3

7.4

4.3

-2.1

12.2

14.6

4.1

4c-Rh4+

9.8

11.4

14.7

1.9

5.8

6.1

10.3

5.9

-1.0

15.8

18.3

4.6

3c-Ru3+

13.0

14.9

17.4

6.9

9.4

11.6

13.3

9.2

4.5

18.0

20.3

8.8

4+

4c-Ru

13.6

15.2

17.5

7.7

10.4

12.3

13.4

10.1

3.8

18.8

21.3

9.5

3c-Os2+

29.0

44.7

31.6

26.8

26.5

30.3

29.6

26.9

29.6

34.2

37.4

29.3

3c-Tl3+

-14.6

-10.5

-9.3

-18.8

-15.6

-18.3

-12.9

-17.8

-23.0

-7.1

-4.4

-18.8

4c-Pb

-13.1

-7.7

-8.3

-14.0

-10.7

-14.4

-9.8

-14.8

-15.3

-4.2

-0.3

-15.1

4c-Tl3+

-8.9

-3.0

-5.5

-4.2

-2.0

-6.2

-6.2

-10.0

-9.0

-0.7

3.8

-7.7

3c-In3+

-4.1

-1.4

-0.5

-8.8

-6.4

-5.3

-3.2

-6.1

-11.9

1.9

4.1

-7.3

3c-Ga

5.2

6.4

7.6

0.0

1.2

4.2

4.9

3.3

-3.0

9.3

11.2

1.9

25.3

28.4

27.7

24.6

24.9

23.4

24.5

22.3

19.2

28.6

31.1

24.4

4+

Late TM

PostTM

4+

3+

3c-Pb2+

3. Comparison of MP2 and DFT Trends The results for ΔGM-H,200 °C calculated with MP2 is shown in Figure 8 where the MP2 ΔGM-H,200 °C values are partitioned by metal type (No notable trends are observed based on formal oxidation state or 13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 28

coordination number, see supplementary information Figure S1). Some general features of the MP2 trends are the same as those of the DFT trends. Namely, that all early TM ions are predicted to have an endergonic ΔGM-H,200 °C while all except one of the post-TM containing cluster models are predicted to have favorable catalytic activity (i.e., exergonic ΔGM-H,200 °C). The MP2 ΔGM-H,200 °C’s predicted for late TM ions have a mixture of endergonic and exergonic values similar fashion as the DFT predicted values. The reproduction of the general features of the DFT trends gives further confirmation that MP2 is a reliable choice for the reference method. More importantly; together, the MP2 and DFT results indicate that late TM and early p-block (i.e. post-TM) metal ions have favorable catalytic properties over early TM ions.

14 ACS Paragon Plus Environment

Page 15 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 8. MP2 predicted ΔGM-H,200 °C trends. MP2 Trends partitioned by metal type. Units in kcal/mol. Early – Group 3-7, Late – Groups 8-12, Post-TM – Groups 13,14.

A more detailed comparison of the trends is given in Figures 9 and 10 that shows a heatmap and letter-value plot for both the MP2 and DFT trends (See supplementary information for Table representation, Table S2). The heatmap representation (Figure 9) allows for a straightforward comparison of the MP2 and DFT trends; herein, the following description is provided to help with interpreting the heatmap. A row in the heatmap contains the ΔGM-H,200 °C values for a particular metal ion 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 28

calculated with various methods. The order for the rows is determined by the trend predicted by MP2. Since the columns represent the trends predicted by a particular method, analysis of the color gradient of a column results in a direct comparison of the DFT trend represented by that column with the MP2 trend. The main reasons for possible differences in the column color gradient are two-fold. First, the spread of values can be different between methods. Therefore, even if a functional were to reproduce the exact trend predicted by MP2, the color gradient could be different due to the different spread of values. Second, an abrupt change in color in the color gradient of a column represents a difference in the ordering between the MP2 trend and the trend for the functional represented by the column. It is apparent from the heatmap in Figure 9 that the trends predicted by the MP2 method are generally reproduced by all the functionals tested except for SVWN5 and revTPSS. The following discussion excludes the SVWN5 and revTPSS functionals; the discussion for the SVWN5 and revTPSS outliers is presented shortly in a later section. In general, any difference in the trend order between a functional trend and MP2 trend are small. For example, 4c-W4+ is out of order for most of the functionals tested. However, the magnitude of ΔGM-H,200 °C for 4c-W4+ is very close to the magnitude of ΔGM-H,200 °C for the next two nearest metal ions. Thus, the change in ordering probably does not result in a meaningful change in the activity. There are some noticeable exceptions to this general observation. The most egregious example is the 45 kcal/mol ΔGM-H,200 °C value predicted by the B97D functional for the three-coordinated Os2+ (3c-Os2+), the magnitude of ΔGM-H,200 °C for the next two nearest metal ions are 32 kcal/mol and 30 kcal/mol respectively which is clearly seen in Figure 9. Other notable exceptions, include the 4c-Rh4+, 3c-Zn2+, and 3c-Ag+ cluster models. The change in trend order for these metal ions is also large. For example, the ΔGM-H,200 °C for the next two nearest neighbors of 3c-Zn2+ calculated with M06L are 11 kcal/mol and 6 kcal/mol where the value for the 3c-Zn2+ cluster model is 3.9 kcal/mol. The common thread among these clusters with substantial changes in the trend order are that, given their formal oxidation state and hence d-count, they possess either a half-filled or full d shell. The violin plots in Figure 10 shows the distribution of ΔGM-H,200 °C values for each method. Violin plots are essentially box plots with a rotated kernel density plot on each side. Thus, the shape of the violin plot represents all possible results while the thickness indicates how common that result is. The three lines in the violin plot, from top to bottom, represent the 75% quartile, 50% quartile i.e., the median, and the 25% quartile. The vertical shifts of the various DFT violin plots relative to the MP2 violin plot shows the difference in the magnitude ΔGM-H,200 °C. While there are vertical shifts for all the functionals, the violin plots show that the overall distribution shape in addition to the statistically important values, i.e., the maximum, minimum, and interquartile ranges, are similar across all methods, indicating that the trends for all the methods are very similar. The previous discussion demonstrated that the MP2 predicted trends are mostly conserved within all the functionals tested. The discussion that follows concentrates on the accuracy of ΔGM-H,200 °C for the functionals tested through a comparison with ΔGM-H,200 °C calculated with MP2.

16 ACS Paragon Plus Environment

Page 17 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 9. Heatmap of ΔGM-H,200 °C calculated with MP2 and various DFT functionals. Units in kcal/mol.

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 28

Figure 10. Violin Plot of ΔGM-H,200 °C calculated with MP2 and various DFT functionals. Units in kcal/mol.

4. Evaluation of the accuracy of DFT functionals. The mean absolute deviation (MAD) and standard deviation (STD) of ΔGM-H,200 °C calculated using DFT from ΔGM-H,200 °C calculated with MP2 are shown in Figure 11, partitioned by the functionals tested. Figure 11 reveals that the wB97xd functional has the lowest STD values followed closely by the PBE0 functional while MN12SX functional has the lowest MAD values followed closely by the PBE functional. While STD and MAD values have similar meaning, STD values are more affected by large deviations than are MAD values. Hence, the fact that MN12SX and PBE have the lowest MAD but not STD values, indicate that there are large deviations present in the MN12SX and PBE data sets that are not present in the wB97xd and PBE0 data sets. The low deviations resulting from using PBE and PBE0 functionals is to be expected; previous benchmark calculations on silica structures51 show that PBE and PBE0 result in low deviations. Additionally, benchmark studies on reactions involving late transition metals53 and on metal containing catalysts110 show that PBE and PBE0 lead to good results. The low deviation from the MN12SX and wB97xd functionals is most likely a consequence of the rangeseparated and long-range corrected feature, respectively, embedded in the design of these functionals.

Figure 11. The mean absolute deviation and standard deviation from ΔGM-H,200 °C calculated with MP2 is given by functional tested. Units in kcal/mol.

18 ACS Paragon Plus Environment

Page 19 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The range-separated and range-corrected functionals allow for the capture of long-range interactions. The many hyperconjugation interactions present in the silica structure108,109 are long-range interactions that benefit from the use of a range-separated or range-corrected functionals. Conversely, the revTPSS functional has the highest STD and MAD values followed by the SVWN5 functional. The failure of the SVWN5 functional is expected since it represents the simplest functional and is known to fail for systems where the electron densities deviate from the density of a uniform electron gas i.e. most chemical systems. However, the substandard performance of the revTPSS functional warrants further discussion. The revTPSS functional is a modification of the TPSS functional that accounts for solid-state effects not present in molecular environments through parametrization.81 However, while there are longrange effects in the SSHC through the presence of hyperconjugation in the silica support, these effects are molecular in nature. Further, the properties of the active site of the SAHC are highly local in nature. Thus, the failure of the revTPSS functional is most likely due to the explicit parametrization of the revTPSS method to reproduce solid-state parameters such as lattice energy. The solid-state effects are missing in these supported metal ions and therefore the application of the revTPSS functional leads to spurious corrections and hence inconsistent results. The original TPSS functional also shows a high STD value, however the MAD value is much lower, indicating the presence of some large deviations that are skewing the STD results. Finally, B3LYP which is used in the original study by Liu et. al.55 does just as well as the top performing functionals. The results for the B3LYP functional have the 5th lowest STD value and 3rd lowest MAD value and are close in magnitude to the lowest STD and MAD values, thus validating the choice in model chemistry used in the original study by Liu et. al55 that developed the activity-descriptor relationship used throughout the present study.

Figure 12. The signed deviation from ΔGM-H,200 °C calculated with MP2 is given for each functional tested as well as a box plot of the statistics. The deviations are partitioned by the cluster models tested. Units in kcal/mol.

The signed deviations of ΔGM-H,200 °C calculated using various DFT functionals from ΔGM-H,200 °C calculated with MP2 i.e. ΔΔGM-H,200 °C values are shown in Figure 12, partitioned by the metal ions tested. 19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 28

From Figure 12, it is observed that the spread in the distribution of ΔΔGM-H,200 °C for each cluster model tested varies significantly. Early TM ions such as 4c-Zr4+ display the smallest spread in the distribution of ΔΔGM-H,200 °C, indicating a high precision in the calculation of ΔGM-H,200 °C across all functionals for cluster models that contain an early TM ion. The accuracy of the DFT functionals in the calculation of ΔGM-H,200 °C for early TM containing cluster-models is good as well, as the entire distribution of ΔΔGM-H,200 °C values produced by the functionals tested usually lie within ±5 of the median value, where the median value is usually around zero. The accuracy and precision begin to suffer as the d-count for the metal ion is raised, i.e. the accuracy and precision is generally worse for late TM ions. For example, the 3c-Au3+ and the 3cPd3+ cluster models have large errors and a large spread of ΔGM-H,200 °C values. The increasing error with increasing d-count indicates that the DFT functionals are not properly account for the differential correlations present in TM. TM atoms and TM containing species require post-HF methods due to the existence of differential correlation effects.100 Most notable of these differential correlation effects are the valence-valence and core-valence correlation effects that increase as the d-count increases. Consequently, the large deviations and large spread of values resulting from DFT calculations on late-TM systems are likely due to the failure of the DFT functionals in describing these differential correlation effects that are appropriately accounted for in the MP2 method. The failure of the DFT functionals in describing the differential correlation effects is also most likely responsible for the substantial changes from the MP2 trend order observed earlier for metal ions that have half-filled and full d-shell configurations. However, it is important to note that DFT and MP2 capture correlation effects in distinct ways, thus the conjecture presented in this paper requires further validation with more careful studies comparing correlation effects recovered through DFT and MP2. The mean signed deviations (MSD) of ΔGM-H,200 °C calculated using various functionals relative to the MP2 ΔGM-H,200 °C are shown in Table 3, partitioned by various TM properties. Regardless of TM property, SVWN5 consistently gave the most negative deviation (reflecting the tendency of SVWN5 to over-bind), revTPSS consistently gave the most positive deviation, while the MN12SX functional consistently produces MSD values close to zero. These are the same trends described for Figure 7, i.e., the general DFT trends. There is no clear trend among the oxidation and coordination number categories as each functional usually produces similar MSD values for all oxidation states and coordination numbers (1+ state category is excluded since there is only one cluster model that is in a 1+ state.). The MSD for the functionals categorized by the type of metal gives more interesting results. Notably, the functionals that display the smallest MAD and STD results all have MSD values close to zero for late TM ions; namely the wB97xd and PBE0 functionals followed closely the M06L, M06, MN12SX and B3LYP functionals. The superior performance of these functionals on late TM containing species is most likely the reason why they produced the lowest MAD and STD values overall. Table 3. Signed deviation from MP2, ΔΔGM-H,200 °C, given by functional tested. The signed deviations are partitioned by oxidation state, coordination number, and type of metal. The data is a further breakdown of the data presented in Figure 12. Units in kcal/mol. NOTE: ΔΔGMH,200 °C = G in this table. property

Oxidation State

Coordination Number

Type of Metal

G(SVWN5)

G(M06)

G(wB97xd)

G(PBE0)

G(M06L)

G(MN12SX)

G(B3LYP)

G(PBE)

G(B97d)

G(BLYP)

G(TPSS)

G(revTPSS)

1+

-18.0

-3.3

-4.0

-6.3

-1.9

-2.3

-5.1

-8.1

-4.1

-5.4

-1.5

1.9

2+

-6.7

-3.0

-1.1

-1.8

-2.4

0.3

0.8

-0.2

3.3

3.1

5.1

7.6

3+

-6.3

-3.1

-2.4

-1.9

-1.4

-0.7

1.4

1.8

3.0

5.5

6.3

8.6

4+

-7.2

-3.2

-2.4

-1.9

-1.6

0.0

1.0

1.3

3.4

4.7

5.7

8.0

3c

-7.0

-3.2

-1.9

-1.9

-1.9

-0.3

1.2

0.8

3.0

4.4

5.7

8.0

4c

-7.1

-3.1

-2.5

-2.2

-1.5

-0.2

0.6

1.1

3.2

4.3

5.6

8.0

Early

-7.0

-4.5

-3.4

-3.1

-3.8

-1.3

-0.5

-0.3

1.1

2.7

3.4

5.4

20 ACS Paragon Plus Environment

Page 21 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Late

-6.7

-1.3

0.0

0.0

0.7

2.0

3.6

2.9

5.8

7.3

8.6

11.2

PostTM

-7.6

-4.0

-4.3

-4.4

-1.9

-3.3

-2.2

-1.0

1.5

1.4

4.1

7.1

Conclusion The reactivity of silica supported metal ions with H2 i.e. reaction free energies of metal hydride formation at 200 °C, ΔGM-H,200 °C was studied using both wavefunction and DFT methods. The considered metal ions include mostly 4d and 5d TM ions, as well as selected 3d and post-TM ions. MP2 is shown to capture sufficient electron correlation effects to produce reliable ΔGM-H,200 °C values. Through comparison with MP2, it is found that the wB97xd and PBE0 functionals produce the lowest STD value while the MN12SX and PBE functionals produce the lowest MAD values. In addition, the model chemistry used to derive the activity-descriptor relationship—the B3LYP functional in conjunction with a triple-ζ basis set—is shown to have the similar MAD and STD values as the top two performing functionals, validating the original choice in model chemistry. Based on these results, it is recommended to use a triple-ζ basis set in conjunction with one of the following functionals: wB97xd, MN12SX, PBE, PBE0 or B3LYP for SAHC’s that contain a bare metal ion as the active site. Previous studies by Liu et. al55 showed that the calculated ΔGM-H,200 °C (an activity descriptor) is strongly correlated with the experimental activity of propylene hydrogenation for such catalyst systems, i.e., lower ΔGM-H,200 °C corresponds to higher catalytic activity (See Figure 1). Using the same activitydescriptor relationship, our calculations suggest that early TM ions would possess poor catalytic activity for propylene hydrogenation while late TM ions have a wide range of predicted catalytic activity. The best candidates are observed to be post-TM (p-block semimetals) ions (i.e. Tl3+, Pb4+, In3+, and Ga3+) where all but one metal ion (i.e. Pb2+) tested are predicted to be catalytically active. Since Ga3+ is known to be catalytically active,55 it is used as a natural cutoff point for assigning SAHC’s with possible active catalytic activity. Based on this natural cutoff point, the following ions could be potentially catalytically active towards the hydrogenation of olefins, ordered from most catalytically active to least: Au3+, Pd2+, Pt4+, Pd4+, Tl3+, Ir4+, Hg2+, Rh3+, Pb4+, In3+, Ir3+, Os4+, Cd2+, Ru2+, and finally Ga3+. An important caveat to these predictions is the stability of the bare ion as an active site. For example, the Pt4+/SiO2 SAHC is difficult to synthesize experimentally since the Pt4+ active sites quickly agglomerate to form nanoparticles instead of highly dispersed active sites. Thus, the accuracy of these predictions is contingent on the stability of the bare ion as an active site. ASSOCIATED CONTENT Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org. Details of the ONIOM Models tested, S-Value results for 4-coordinated SAHC ONIOM models, raw values for single-point calculations used to make Figures, pivot tables for DFT results, additional figures that visualize pivot table for MP2 results and ONIOM model error, as well as the coordinates for both the active site and hydrogenated cluster models in xyz format. AUTHOR INFORMATION Corresponding Author C. Liu ([email protected]), L. A. Curtiss ([email protected]) 21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 28

Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. ACKNOWLEDGEMENT This work was financially supported by the U.S. Department of Energy (DOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, under Contract DEAC0206CH11357 (Argonne National Laboratory). Cesar Plascencia would like to acknowledge the financial support from the Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education (ORISE) for the DOE. ORISE is managed by ORAU under contract number DE-SC0014664. The calculations were performed using the computing resources provided by the Laboratory Computing Resource Center (LCRC) at Argonne National Laboratory and the National Energy Research Scientific Computing (NERSC) Center at Lawrence Berkeley National Laboratory.

References (1) (2) (3) (4) (5)

(6) (7) (8) (9) (10) (11) (12)

Pelletier, J. D. A.; Basset, J. M. Catalysis by Design: Well-Defined Single-Site Heterogeneous Catalysts. Acc. Chem. Res. 2016, 49, 664–677. Hlatky, G. G. Heterogeneous Single-Site Catalysts for Olefin Polymerization. Chem. Rev. 2000, 100, 1347–1376. Thomas, J. M.; Raja, R.; Lewis, D. W. Single-Site Heterogeneous Catalysts. Angew. Chemie - Int. Ed. 2005, 44, 6456–6482. Thomas, J. M.; Raja, R. Mono-, Bi- and Multifunctional Single-Sites: Exploring the Interface Between Heterogeneous and Homogeneous Catalysis. Top. Catal. 2010, 53, 848–858. Copéret, C.; Comas-Vives, A.; Conley, M. P.; Estes, D. P.; Fedorov, A.; Mougel, V.; Nagae, H.; Núnez-Zarur, F.; Zhizhko, P. A. Surface Organometallic and Coordination Chemistry Toward Single-Site Heterogeneous Catalysts: Strategies, Methods, Structures, and Activities. Chem. Rev. 2016, 116, 323–421. Yang, X. F.; Wang, A.; Qiao, B.; Li, J.; Liu, J.; Zhang, T. Single-Atom Catalysts: A New Frontier in Heterogeneous Catalysis. Acc. Chem. Res. 2013, 46, 1740–1748. Serna, P.; Gates, B. C. Molecular Metal Catalysts on Supports: Organometallic Chemistry Meets Surface Science. Acc. Chem. Res. 2014, 47, 2612–2620. Dal Santo, V.; Liguori, F.; Pirovano, C.; Guidotti, M. Design and Use of Nanostructured SingleSite Heterogeneous Catalysts for the Selective Transformation of Fine Chemicals. Molecules 2010, 15, 3829–3856. Judai, K.; Abbet, S.; Wörz, A. S.; Heiz, U.; Henry, C. R. Low-Temperature Cluster Catalysis. J. Am. Chem. Soc. 2004, 126, 2732–2737. Lei, Y.; Mehmood, F.; Lee, S.; Greeley, B.; Siefert, S.; Winans, R. E.; Elam, J. W.; Meyer, R. J.; Redfern, P. C.; Teschner, D.; et al. Increased Silver Activity for Direct Propylene Epoxidation via Subnanometer Size Effects. Science. 2010, 328, 224-228. Turner, M.; Golovko, V. B.; Vaughan, O. P. H.; Abdulkin, P.; Berenguer-Murcia, A.; Tikhov, M. S.; Johnson, B. F. G.; Lambert, R. M. Selective Oxidation with Dioxygen by Gold Nanoparticle Catalysts Derived from 55-Atom Clusters. Nature 2008, 454, 981–983. Herzing, A. A.; Kiely, C. J.; Carley, A. F.; Landon, P.; Hutchings, G. J. Identification of Active Gold Nanoclusters on Iron Oxide Supports for CO Oxidation. Science. 2008, 321, 1331–1335. 22 ACS Paragon Plus Environment

Page 23 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(13)

(14) (15)

(16)

(17) (18) (19) (20)

(21) (22)

(23) (24) (25)

(26)

(27)

The Journal of Physical Chemistry

Vajda, S.; Pellin, M. J.; Greeley, J. P.; Marshall, C. L.; Curtiss, L. A.; Ballentine, G. A.; Elam, J. W.; Catillon-Mucherie, S.; Redfern, P. C.; Mehmood, F.; et al. Subnanometre Platinum Clusters as Highly Active and Selective Catalysts for the Oxidative Dehydrogenation of Propane. Nat. Mater. 2009, 8, 213–216. Remediakis, I. N.; Lopez, N.; Nørskov, J. K. CO Oxidation on Rutile-Supported Au Nanoparticles. Angew. Chemie - Int. Ed. 2005, 44, 1824–1826. Das, U.; Zhang, G.; Hu, B.; Hock, A. S.; Redfern, P. C.; Miller, J. T.; Curtiss, L. A. Effect of Siloxane Ring Strain and Cation Charge Density on the Formation of Coordinately Unsaturated Metal Sites on Silica: Insights from Density Functional Theory (DFT) Studies. ACS Catal. 2015, 5, 7177–7185. Sohn, H.; Camacho-Bunquin, J.; Langeslay, R. R.; Ignacio-de Leon, P. A.; Niklas, J.; Poluektov, O. G.; Liu, C.; Connell, J. G.; Yang, D.; Kropf, J.; et al. Isolated, Well-Defined Organovanadium( iii ) on Silica: Single-Site Catalyst for Hydrogenation of Alkenes and Alkynes. Chem. Commun. 2017, 137, 6770–6780. Schweitzer, N. M.; Hu, B.; Das, U.; Kim, H.; Greeley, J.; Curtiss, L. A.; Stair, P. C.; Miller, J. T.; Hock, A. Propylene Hydrogenation and Propane Dehydrogenation by a Single-Site Zn2+ on Silica Catalyst. ACS Catal. 2014, 4, 1091–1098. Camacho-Bunquin, J.; Ferrandon, M.; Das, U.; Dogan, F.; Liu, C.; Larsen, C.; Platero-Prats, A. E.; Curtiss, L. A.; Hock, A. S.; Miller, J. T.; et al. Supported Aluminum Catalysts for Olefin Hydrogenation. ACS Catal. 2017, 7, 689–694. Camacho-Bunquin, J.; Siladke, N. A.; Zhang, G.; Niklas, J.; Poluektov, O. G.; Nguyen, S. T.; Miller, J. T.; Hock, A. S. Synthesis and Catalytic Hydrogenation Reactivity of a Chromium Catecholate Porous Organic Polymer. Organometallics 2015, 34, 947–952. Camacho-Bunquin, J.; Aich, P.; Ferrandon, M.; Getsoian, A.; Das, U.; Dogan, F.; Curtiss, L. A.; Miller, J. T.; Marshall, C. L.; Hock, A. S.; et al. Single-Site Zinc on Silica Catalysts for Propylene Hydrogenation and Propane Dehydrogenation: Synthesis and Reactivity Evaluation Using an Integrated Atomic Layer Deposition-Catalysis Instrument. J. Catal. 2017, 345, 170–182. Getsoian, A. G. B.; Hu, B.; Miller, J. T.; Hock, A. S. Silica-Supported, Single-Site Sc and Y Alkyls for Catalytic Hydrogenation of Propylene. Organometallics 2017, 36, 3677–3685. Getsoian, A.; Das, U.; Camacho-Bunquin, J.; Zhang, G.; Gallagher, J. R.; Hu, B.; Cheah, S.; Schaidle, J. A.; Ruddy, D. A.; Hensley, J. E.; et al. Organometallic Model Complexes Elucidate the Active Gallium Species in Alkane Dehydrogenation Catalysts Based on Ligand Effects in Ga K-Edge XANES. Catal. Sci. Technol. 2016, 101, 12–18. Hu, B.; Getsoian, A.; Schweitzer, N. M.; Das, U.; Kim, H.; Niklas, J.; Poluektov, O.; Curtiss, L. A.; Stair, P. C.; Miller, J. T.; et al. Selective Propane Dehydrogenation with Single-Site CoII on SiO2 by a Non-Redox Mechanism. J. Catal. 2015, 322, 24–37. Ranocchiari, M.; Lothschutz, C.; Grolimund, D.; van Bokhoven, J. A. Single-Atom Active Sites on Metal-Organic Frameworks. Proc. R. Soc. A Math. Phys. Eng. Sci. 2012, 468, 1985–1999. Wu, L.; Navrotsky, A. Synthesis and Thermodynamic Study of Transition Metal Ion (Mn2+ , Co2+ , Cu2+ , and Zn2+ ) Exchanged Zeolites A and Y. Phys. Chem. Chem. Phys. 2016, 18, 10116– 10122. Bernales, V.; Yang, D.; Yu, J.; Gümüşlu, G.; Cramer, C. J.; Gates, B. C.; Gagliardi, L. Molecular Rhodium Complexes Supported on the Metal-Oxide-like Nodes of Metal Organic Frameworks and on Zeolite HY: Catalysts for Ethylene Hydrogenation and Dimerization. ACS Appl. Mater. Interfaces 2017, 9, 33511–33520. Li, Z.; Schweitzer, N. M.; League, A. B.; Bernales, V.; Peters, A. W.; Getsoian, A. B.; Wang, T. 23 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(28)

(29) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39) (40) (41) (42) (43) (44) (45) (46) (47)

Page 24 of 28

C.; Miller, J. T.; Vjunov, A.; Fulton, J. L.; et al. Sintering-Resistant Single-Site Nickel Catalyst Supported by Metal-Organic Framework. J. Am. Chem. Soc. 2016, 138, 1977–1982. Bernales, V.; League, A. B.; Li, Z.; Schweitzer, N. M.; Peters, A. W.; Carlson, R. K.; Hupp, J. T.; Cramer, C. J.; Farha, O. K.; Gagliardi, L. Computationally Guided Discovery of a Catalytic Cobalt-Decorated Metal-Organic Framework for Ethylene Dimerization. J. Phys. Chem. C 2016, 120, 23576–23583. Comito, R. J.; Fritzsching, K. J.; Sundell, B. J.; Schmidt-Rohr, K.; DincǍ, M. Single-Site Heterogeneous Catalysts for Olefin Polymerization Enabled by Cation Exchange in a MetalOrganic Framework. J. Am. Chem. Soc. 2016, 138, 10232–10237. Zhao, J.; Liu, X.; Chen, Z. Frustrated Lewis Pair Catalysts in Two Dimensions: B/Al-Doped Phosphorenes as Promising Catalysts for Hydrogenation of Small Unsaturated Molecules. ACS Catal. 2017, 7, 766–771. Mas-Ballesté, R.; Gómez-Navarro, C.; Gómez-Herrero, J.; Zamora, F. 2D Materials: To Graphene and Beyond. Nanoscale 2011, 3, 20–30. Coudert, F. X.; Fuchs, A. H. Computational Characterization and Prediction of Metal-Organic Framework Properties. Coord. Chem. Rev. 2015, 307, 211–236. Bosch, M.; Yuan, S.; Rutledge, W.; Zhou, H. C. Stepwise Synthesis of Metal-Organic Frameworks. Acc. Chem. Res. 2017, 50, 857–865. Rao, C. N. R.; Maitra, U. Inorganic Graphene Analogs. Annu. Rev. Mater. Res. 2015, 45, 29–62. Stock, N.; Biswas, S. Synthesis of Metal-Organic Frameworks (MOFs): Routes to Various MOF Topologies, Morphologies, and Composites. Chem. Rev. 2012, 112, 933–969. Gajan, D.; Copéret, C. Silica-Supported Single-Site Catalysts: To Be or Not to Be? A Conjecture on Silica Surfaces. New J. Chem. 2011, 35, 2403-2408. Jiao, L.; Regalbuto, J. R. The Synthesis of Highly Dispersed Noble and Base Metals on Silica via Strong Electrostatic Adsorption: II. Mesoporous Silica SBA-15. J. Catal. 2008, 260, 342–350. Regalbuto, J. R.; Navada, A.; Shadid, S.; Bricker, M. L.; Chen, Q. An Experimental Verification of the Physical Nature of Pt Adsorption onto Alumina. J. Catal. 1999, 184, 335–348. Martynowycz, M. W.; Hu, B.; Kuzmenko, I.; Bu, W.; Hock, A.; Gidalevitz, D. Monomolecular Siloxane Film as a Model of Single Site Catalysts. J. Am. Chem. Soc. 2016, 138, 12432–12439. Guesmi, H.; Tielens, F. Chromium Oxide Species Supported on Silica: A Representative Periodic DFT Model. J. Phys. Chem. C 2012, 116, 994–1001. Tielens, F.; Gervais, C.; Lambert, J. F.; Mauri, F.; Costa, D. Ab Initio Study of the Hydroxylated Surface of Amorphous Silica: A Representative Model. Chem. Mater. 2008, 20, 3336–3344. Handzlik, J.; Grybos, R.; Tielens, F. Structure of Monomeric Chromium(VI) Oxide Species Supported on Silica: Periodic and Cluster DFT Studies. J. Phys. Chem. C 2013, 117, 8138–8149. Evarestov, R. A.; Bredow, T.; Jug, K. Connection between Slab and Cluster Models for Crystalline Surfaces. Phys. Solid State 2001, 43, 1774–1782. Guesmi, H.; Grybos, R.; Handzlik, J.; Tielens, F. Characterization of Tungsten Monomeric Oxide Species Supported on Hydroxylated Silica: a DFT Study. RSC Adv. 2016, 6, 39424–39432. Rimola, A.; Costa, D.; Sodupe, M.; Ugliengo, P. Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments Chem. Rev. 2013, 113, 4216–4313. Islam, M. M.; Costa, D.; Calatayud, M.; Tielens, F. Characterization of Supported Vanadium Oxide Species on Silica: A Periodic DFT Investigation. J. Phys. Chem. C 2009, 113, 10740– 10746. Feher, F. J.; Newman, D. A.; Walzer, J. F. Silsesquioxanes as Models for Silica Surfaces. J. Am. Chem. Soc. 1989, 111, 1741–1748. 24 ACS Paragon Plus Environment

Page 25 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(48) (49) (50) (51) (52) (53) (54)

(55)

(56) (57) (58) (59) (60) (61) (62) (63) (64) (65)

The Journal of Physical Chemistry

Bartlett, R. J.; Musial, M. Coupled-Cluster Theory in Quantum Chemistry. Rev. Mod. Phys. 2007, 79, 291–352. Perdew, J. P.; Ruzsinszky, A.; Tao, J.; Staroverov, V. N.; Scuseria, G. E.; Csonka, G. I. Prescription for the Design and Selection of Density Functional Approximations: More Constraint Satisfaction with Fewer Fits. J. Chem. Phys. 2005, 123 (6). Becke, A. D. Perspective: Fifty Years of Density-Functional Theory in Chemical Physics. J. Chem. Phys. 2014, 140, 18A301. Cypryk, M.; Gostyński, B. Computational Benchmark for Calculation of Silane and Siloxane Thermochemistry. J. Mol. Model. 2016, 22, 1–19. Harvey, J. N. On the Accuracy of Density Functional Theory in Transition Metal Chemistry. Annu. Reports Sect. “C” (Physical Chem.) 2006, 102, 203-226. Quintal, M. M.; Karton, A.; Iron, M. A.; Boese, A. D.; Martin, J. M. L. Benchmark Study of DFT Functionals for Late-Transition-Metal Reactions. J. Phys. Chem. A 2006, 110, 709–716. Rohmann, K.; Hölscher, M.; Leitner, W. Can Contemporary Density Functional Theory Predict Energy Spans in Molecular Catalysis Accurately Enough to Be Applicable for in Silico Catalyst Design? A Computational/Experimental Case Study for the Ruthenium-Catalyzed Hydrogenation of Olefins. J. Am. Chem. Soc. 2016, 138, 433–443. Liu, C.; Camacho-Bunquin, J.; Ferrandon, M.; Savara, A.; Sohn, H.; Yang, D.; Kaphan, D. M.; Langeslay, R. R.; Ignacio-de Leon, P. A.; Liu, S.; et al. Development of Activity–Descriptor Relationships for Supported Metal Ion Hydrogenation Catalysts on Silica. Polyhedron 2018, 152, 73-83. Graciani, J.; Mudiyanselage, K.; Xu, F.; Baber, A. E.; Evans, J.; Senanayake, S. D.; Stacchiola, D. J.; Liu, P.; Hrbek, J.; Sanz, J. F.; et al. Highly Active Copper-Ceria and Copper-Ceria-Titania Catalysts for Methanol Synthesis from CO2. Science 2014, 345, 546–550. Liu, C.; Yang, B.; Tyo, E.; Seifert, S.; DeBartolo, J.; von Issendorff, B.; Zapol, P.; Vajda, S.; Curtiss, L. A. Carbon Dioxide Conversion to Methanol over Size-Selected Cu4 Clusters at Low Pressures. J. Am. Chem. Soc. 2015, 137, 8676–8679. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision D.01. Gaussian Inc.: Wallingford, CT 2009. Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648-5652. Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789. Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58, 1200– 1211. Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623–11627. Stevens, W. J.; Basch, H.; Krauss, M. Compact Effective Potentials and Efficient SharedExponent Basis Sets for the First- and Second-Row Atoms. J. Chem. Phys. 1984, 81, 6026–6033. Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007-1023. Dunning, T. H.; Peterson, K. A.; Wilson, A. K. Gaussian Basis Sets for Use in Correlated Molecular Calculations. X. The Atoms Aluminum through Argon Revisited. J. Chem. Phys. 2001, 114, 9244-9253. 25 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(66)

(67) (68) (69) (70) (71) (72) (73) (74) (75) (76) (77) (78) (79) (80) (81) (82) (83) (84)

Page 26 of 28

Hill, J. G.; Peterson, K. A. Explicitly Correlated Coupled Cluster Calculations for Molecules Containing Group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) Elements: Optimized Complementary Auxiliary Basis Sets for Valence and Core-Valence Basis Sets. J. Chem. Theory Comput. 2012, 8, 518–526. Balabanov, N. B.; Peterson, K. A. Systematically Convergent Basis Sets for Transition Metals. I. All-Electron Correlation Consistent Basis Sets for the 3d Elements Sc-Zn. J. Chem. Phys. 2005, 123, 064107. Peterson, K. A.; Figgen, D.; Dolg, M.; Stoll, H. Energy-Consistent Relativistic Pseudopotentials and Correlation Consistent Basis Sets for the 4d Elements Y-Pd. J. Chem. Phys. 2007, 126, 124101. Figgen, D.; Peterson, K. A.; Dolg, M.; Stoll, H. Energy-Consistent Pseudopotentials and Correlation Consistent Basis Sets for the 5d Elements Hf-Pt. J. Chem. Phys. 2009, 130, 164108. Zhang, L.; Wang, A.; Miller, J. T.; Liu, X.; Yang, X.; Wang, W.; Li, L.; Huang, Y.; Mou, C. Y.; Zhang, T. Efficient and Durable Au Alloyed Pd Single-Atom Catalyst for the Ullmann Reaction of Aryl Chlorides in Water. ACS Catal. 2014, 4, 1546–1553. Misumi, Y.; Seino, H.; Mizobe, Y. Heterolytic Cleavage of Hydrogen Molecule by Rhodium Thiolate Complexes That Catalyze Chemoselective Hydrogenation of Imines under Ambient Conditions. J. Am. Chem. Soc. 2009, 131, 14636–14637. Hulley, E. B.; Welch, K. D.; Appel, A. M.; Dubois, D. L.; Bullock, R. M. Rapid, Reversible Heterolytic Cleavage of Bound H2. J. Am. Chem. Soc. 2013, 135, 11736–11739. Zhang, S.; Appel, A. M.; Bullock, R. M. Reversible Heterolytic Cleavage of the H-H Bond by Molybdenum Complexes: Controlling the Dynamics of Exchange Between Proton and Hydride. J. Am. Chem. Soc. 2017, 139, 7376–7387. Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. B 1964, 136, 863-871. Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140,1133. Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098–3100. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made SimpleERRATA. Phys. Rev. Lett. 1997, 78, 1396. Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158–6170. Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Non-Empirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401. Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. Workhorse Semilocal Density Functional for Condensed Matter Physics and Quantum Chemistry. Phys. Rev. Lett. 2009, 103, 026403. Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. Erratum: Workhorse Semilocal Density Functional for Condensed Matter Physics and Quantum Chemistry (Physical Review Letters (2009) 103 (026403)). Phys. Rev. Lett. 2011, 106, 179902. Zhao, Y.; Truhlar, D. G. A New Local Density Functional for Main-Group Thermochemistry, Transition Metal Bonding, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Phys. 2006, 125, 194101. Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group 26 ACS Paragon Plus Environment

Page 27 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215–241. (85) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787–1799. (86) Peverati, R.; Truhlar, D. G. Screened-Exchange Density Functionals with Broad Accuracy for Chemistry and Solid-State Physics. Phys. Chem. Chem. Phys. 2012, 14, 16187–16191. (87) Chai, J.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Improved Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615–6620. (88) Møller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618–622. (89) Purvis, G. D.; Bartlett, R. J. A Full Coupled‐ Cluster Singles and Doubles Model: The Inclusion of Disconnected Triples. J. Chem. Phys. 1982, 76, 1910–1918. (90) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A Fifth-Order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479–483. (91) Maseras, F.; Morokuma, K. IMOMM: A New Integrated Ab Initio + Molecular Mechanics Geometry Optimization Scheme of Equilibrium Structures and Transition States. J. Comp. Chem. 1995, 1170–1179. (92) Humbel, S.; Sieber, S.; Morokuma, K. The IMOMO Method: Integration of Different Levels of Molecular Orbital Approximations for Geometry Optimization of Large Systems: Test for NButane Conformation and SN2 Reaction: RCl+Cl−. J. Chem. Phys. 1996, 105, 1959-1967. (93) Svensson, M.; Humbel, S.; Froese, R. D. J.; Matsubara, T.; Sieber, S.; Morokuma, K. ONIOM: A Multilayered Integrated MO + MM Method for Geometry Optimizations and Single Point Energy Predictions. A Test for Diels-Alder Reactions and Pt(P(t-Bu)3)2 + H2 Oxidative Addition. J. Phys. Chem. 1996, 100, 19357–19363. (94) Dapprich, S.; Komáromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. A New ONIOM Implementation in Gaussian98. Part I. The Calculation of Energies, Gradients, Vibrational Frequencies and Electric Field Derivatives. J. Mol. Struct. THEOCHEM 1999, 461–462, 1–21. (95) Clemente, F. R.; Vreven, T.; Frisch, M. J. (2010) Getting the Most Out of ONIOM: Guidelines and Pitfalls In Quantum Biochemistry; Wiley & Co; pp 61–83. (96) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; et al. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678–5796. (97) Handzlik, J. Application of the ONIOM (QM/QM) Method in the Study of Molybdena–Silica System Active in Olefin Metathesis. Int. J. Quantum Chem. 2007, 107, 2111-2119. (98) Roggero, I.; Civalleri, B.; Ugliengo, P. Modeling Physisorption with the ONIOM Method: The Case of NH3 at the Isolated Hydroxyl Group of the Silica Surface. Chem. Phys. Lett. 2001, 341, 625–632. (99) Morokuma, K.; Musaev, D. G.; Vreven, T.; Basch, H.; Torrent, M.; Khoroshun, D. V. Model Studies of the Structures, Reactivities, and Reaction Mechanisms of Metalloenzymes. IBM J. Res. Dev. 2001, 45, 367–395. (100) Harrison, J. F. Electronic Structure of Diatomic Molecules Composed of a First-Row Transition Metal and Main-Group Element (H-F). Chem. Rev. 2000, 100, 679–716. (101) Hofmann, H.; Wickham, H.; Kafadar, K. Letter-Value Plots: Boxplots for Large Data. J. Comput. Graph. Stat. 2017, 26, 469–477. (102) Solomon, E. I.; Scott, R. A.; King, R. B. Computational Inorganic and Bioinorganic Chemistry; Wiley, 2009. 27 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 28

(103) Becke, A. D. Density‐ functional Thermochemistry. I. The Effect of the Exchange‐ only Gradient Correction. J. Chem. Phys. 1992, 96, 2155–2160. (104) Chiron, N.; Guilet, R.; Deydier, E. Adsorption of Cu(II) and Pb(II) onto a Grafted Silica: Isotherms and Kinetic Models. Water Res. 2003, 37, 3079–3086. (105) MacHida, M.; Fotoohi, B.; Amamo, Y.; Mercier, L. Cadmium(II) and Lead(II) Adsorption onto Hetero-Atom Functional Mesoporous Silica and Activated Carbon. Appl. Surf. Sci. 2012, 258, 7389–7394. (106) Tran, H. H.; Roddick, F. A.; O’Donnell, J. A. Comparison of Chromatography and Desiccant Silica Gels for the Adsorption of Metal ions—I. Adsorption and Kinetics. Water Res. 1999, 33, 2992–3000. (107) Cea-Olivares, R.; García-Montalvo, V.; Moya-Cabrera, M. M. The Importance of the Transannular Secondary Bonding Strength in the Molecular Structures of Metallocanes of Type [X(CH2CH2Y)2MRR′] and [X(CH2CH2Y)2M′R] (M = Ge(IV), Sn(IV), Pb(IV), M′ = As(III), Sb(III) and Bi(III); X = NR″, O, S; Y = O, S). Coord. Chem. Rev. 2005, 249 (7–8), 859–872. (108) Weinhold, F.; West, R. Hyperconjugative Interactions in Permethylated Siloxanes and Ethers: The Nature of the SiO Bond. J. Am. Chem. Soc. 2013, 135, 5762–5767. (109) Moraru, I. T.; Petrar, P. M.; Nemeş, G. Bridging a Knowledge Gap from Siloxanes to Germoxanes and Stannoxanes. A Theoretical Natural Bond Orbital Study. J. Phys. Chem. A 2017, 121, 2515–2522. (110) Steinmetz, M.; Grimme, S. Benchmark Study of the Performance of Density Functional Theory for Bond Activations with (Ni,Pd)-Based Transition-Metal Catalysts. ChemistryOpen 2013, 2, 115–124.

TOC Graphic

28 ACS Paragon Plus Environment