Matthew D.
Wodrich
a,
Michael
Busch
b and
Clémence
Corminboeuf
*a
aLaboratory for Computational Molecular Design, Institute of Chemical Sciences and Engineering, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland. E-mail: clemence.corminboeuf@epfl.ch
bLaboratory for Computational Molecular Design (LCMD) and National Center for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, Station 12, 1015 Lausanne, Switzerland
First published on 3rd June 2016
Volcano plots are frequently used as aids in the search for new heterogeneous and electrochemical catalysts. These tools successfully predict catalytic processes based solely on thermodynamic descriptions, which also capably describe many aspects of the catalytic cycles of homogeneous species. However, homogeneous catalysts also frequently depend upon the kinetic influences brought about by steric interactions to promote or prevent specific chemical reactions. Here, a prototypical transformation facilitated by a homogeneous catalysis, the hydroformylation of an olefin using CO and H2, is examined to establish the viability of creating kinetic volcano plots and to determine their ability to ascertain the influences steric bulk plays on catalytic cycle energetics. Similar to their thermodynamic counterparts, kinetic volcanoes successfully reproduce many experimentally known facets of the hydroformylation reaction. In contrast to thermodynamic volcanoes, kinetic volcanoes emphasize changes in the height of the different activation barriers brought about by steric interactions. This crucial information, however, comes with considerable computational cost, since the transition states of catalysts bearing large bulky ligands must be identified and characterized. To overcome this drawback, a procedure is proposed that relates a simple steric parameter, the Tolman cone angle, with the descriptors used to create the kinetic volcano plots. In this way, the activation barriers of bulky catalysts can be estimated without requiring expensive transition state computations. These newly derived structure–activity relationship volcano plots represent useful tools for identifying new homogeneous catalysts.
Scheme 1 Current state of methods to computationally assess the thermodynamic and kinetic profiles of homogeneous catalysts. |
Faced with challenges of this type, the heterogeneous catalysis and electrochemistry communities developed an intuitive tool, volcano plots, that facilitate the rapid analysis of large data sets with the aim of identifying potential catalysts with appealing thermodynamic profiles. These graphical depictions, named for their abstract resemblance to the geological formation, represent catalytic activity with respect to the magnitude of catalyst/intermediate interactions.1,2 The hallmark volcano shape identifies thermodynamically appealing catalysts, which appear near the top of the plot, and are considered to be “most ideal” because they closely match Sabatier's concept of an ideal catalyst.3 Sabatier's principle states that interactions between a catalyst and a substrate should be neither too weak nor too strong. In other words, placing reactants onto a catalyst and dissociating the resulting products should be equally facile. Deviations from Sabatier's ideal behaviour result in catalysts appearing progressively lower on the volcano plot. In this way, volcano plots allow for the rapid, simultaneous visualization of the key energetic aspects of the catalytic cycle of many catalysts. Indeed, when used as a tool for screening catalysts, volcano plots can be tremendously valuable for experimental/computational collaborations. Data derived from first principles computations serves to narrow the pool of potential catalysts to only those having the best free energy profiles, which can then be tested experimentally. Importantly the behaviour of heterogeneous catalysts has been successfully described using only thermodynamics, a situation that greatly simplifies computational investigations of different potential catalysts since the activation barriers can be ignored. As a result of their ease of use and remarkable ability to assist in rapidly ascertaining data concerning numerous catalysts, volcano plots are now commonplace in the heterogeneous catalysis and electrochemistry literature.4–6
An essential step toward constructing volcano plots is establishing linear free energy scaling relationships (LFESRs)7 that map the properties of different catalytic cycle intermediates onto a single variable (e.g., the binding energy of a single intermediate).8,9 These relationships are based on the notion that the stabilities of the different intermediates depend upon one another and are restricted from varying independently. Thus, if the stability of the single variable descriptor is known, it becomes possible to estimate the stability of all other catalytic cycle intermediates using the corresponding LFESRs. Indeed, the power of such relationships have long been used to better understand physical organic phenomena,10–12 and are still regularly employed to garner information regarding the behaviour of heterogeneous13 and homogeneous14,15 catalysts.
Determining the existence of LFESRs is a relatively straightforward process. Assuming one possesses knowledge of the reaction mechanism for the chemical transformation of interest, key intermediates can be identified from which it is possible to compute (or determine experimentally) the free energy profile for an individual catalyst (Fig. 1a). During the screening process, multiple free energy profiles will be created for catalysts having different metal/ligand combinations. The data from these multiple free energy profiles can then be analysed for the existence of any LFESRs (Fig. 1b). The mathematical equations defining the scaling relationships can then be plotted relative to one another, which leads to a volcano plot (Fig. 1c).16 Individual catalysts are then placed onto the final plot based on the free energies associated with key transformations, thereby identifying species possessing appealing thermodynamic profiles.
While volcano plots are regularly used for identifying potential new catalysts in the heterogeneous and electrochemistry domains,17 applications to homogeneous catalysts are virtually unknown.18 In 2015, we first explored the viability of volcano plots for homogeneous catalysis19 by examining a prototypical example, Suzuki cross-coupling. The resulting homogeneous volcano plots capably reproduced a variety of experimentally known trends and successfully served as a proof-of-principle example that confirmed the earlier abstractly proposed20 framework. Despite this success, our study considered only thermodynamic aspects of the reaction and ignored the important effects that activation barriers have on overall catalytic performance (Scheme 1, bottom row). This point raises a logical question: can the predictive power of volcano plots be increased through the inclusion of kinetic data? It is well established that the success or failure of homogeneous catalysts often lies in the delicate interplay between thermodynamic and kinetic factors. This is particularly true for regio- or enantioselective reactions, where kinetics fundamentally alter the final isomer ratio. The purpose of this contribution is to explore the viability of constructing “kinetic volcanoes” (Scheme 1, 2nd and 3rd rows) and interpreting the resulting data by examining the hydroformylation of an alkene as a prototypical kinetically driven reaction.
(1) |
The hydroformylation of alkenes via the addition of CO and H2 (eqn (1)) is an industrially important reaction21–27 that employs transition metal containing homogeneous catalysts. While the earliest transformations relied upon Co catalysts,28 Wilkinson's 1965 discovery29,30 of a Rh-based alternative that utilized milder conditions paved the way for the patented low-pressure oxo process of Union-Carbide, which today is used to produce enormous quantities of different aldehydes. To varying degrees this reaction can also be catalysed by other metals,31 but Rh is the most active.21 The proposed reaction mechanism for Rh-catalysed hydroformylation is depicted in Fig. 2. CO dissociation from 1 yields 2 and serves as the entry point into the catalytic cycle. From 2 olefin coordination leads to the π-complex 3, which rearranges to the σ-complex 4. Addition of CO yields 5 and then 6 following CO insertion into the metal–carbon bond. The addition of H2 then produces 7, which undergoes reductive elimination to give 2, completing the catalytic cycle. For the purposes of this study, ethene was selected as the olefin. This choice is motivated by two factors: first, using ethene eliminates conformational problems that arise from using longer alkenes, second, due to the equivalent nature of the carbon atoms, the hydroformylation of ethene results in only one product (propanal), thereby removing competing pathways that lead to the linear and branched regioisomers for asymmetric alkenes. While issues of regioselectivity are certainly of interest, for the purposes of this proof-of-principle study we chose to minimize the occurrence of such complicating factors and focus on assessing the viability and value of creating kinetic volcano plots.
Since the goal of this work is to probe the ways in which different metal/ligand combinations influence the thermodynamic and kinetic aspects of olefin hydroformylation through the use of density functional theory computations, we examined group 8 (Fe, Ru, Os), 9 (Co, Rh, Ir), and 10 (Ni, Pt)32 metals coupled with a series of monodentate phosphine ligands of varying sizes (PH3, PMe3, PPh3, PCy3). The different metal/ligand combinations produced 32 catalysts, from which linear free energy scaling relationships were built and, ultimately, volcano plots representing both the thermodynamics and kinetics constructed. For the sake of simplicity and to more closely follow the Fig. 2 mechanism, the oxidation states of the catalysts were adjusted to comply with the known 16e−/18e− nature of the reaction mechanism. Thus, each of the group 8 and 10 catalysts was represented as a charged system, isoelectronic with the group 9 metal.
(2) |
(3) |
(4) |
(5) |
(6) |
Fig. 3 Linear free energy scaling relationships amongst intermediates for the thermodynamic steps of olefin hydroformylation without consideration of different ligand size. The free energies are relative to a reference state, defined in eqn (2)–(6). The mathematical equations that define the linear scaling relationships are used to create the simulated reaction profile and, ultimately, the final volcano plot (vide infra). The zero slope observed between intermediates 5 and 4 indicates that the binding energy is independent of the choice of catalyst. |
While linear free energy scaling relationships ably describe the thermodynamic interplay amongst intermediates for some homogeneous catalytic cycles, the existence of similar relationships that delineate kinetic aspects is less clear. Certainly amongst the most celebrated works along this line is the Bell–Evans–Polanyi (BEP) model11,12 that describe linear relationships between activation energies and a descriptor, thereby permitting determination of kinetic parameters. For instance, BEP relationships are useful in describing some organic reactions51,52 and have shed light on processes occurring in heterogeneous catalysis.53,54 On the other hand, these relationships have shortcomings of their own, including problems describing certain aspects of organometallic complexes.55 Thus, it is necessary to determine if kinetic aspects of the hydroformylation catalytic cycle, namely the free energy needed to overcome the activation barriers between intermediates, correlate with a purely thermodynamic descriptor (e.g., ΔGRRS of an intermediate). Fortunately, Fig. 4 shows that correlations do exist between the thermodynamic descriptor, ΔGRRS(5), and the relative free energies associated with each of the transition states (TS3,4, TS5,6, TS6,7, TS7,2). Owing to this fact, it is possible to construct volcano plots that describe aspects of both reaction thermodynamics and kinetics, which will be shown in the proceeding section.
Fig. 5 (a) Overview of thermodynamics of olefin hydroformylation without consideration of ligand size, as derived from the linear free energy scaling relationships that define the catalytic cycle depicted in Fig. 2. (b) Volcano plot illustrating the thermodynamic suitability of catalysts derived from (a). Reactions that define the potential determining step are: 7 → 2 for region I, 3 → 4 for region II, 6 → 7 for region III, and 4 → 5 for region IV. Lines defining the volcano are obtained by taking the lowest −ΔG(pds) value amongst all the reaction for each ΔGRRS(5) value. Note that the deviations of the points representing individual catalysts from the free energies predicted through linear scaling relationships (depicted as lines) are expected and acceptable. Significant deviations from the scaling relationships [such as seen in Ni(PMe3)2 and Ni(PPh3)2] may be indicative of unusual behaviour by a specific catalyst. |
Extracting the most relevant information from Fig. 5a leads to the volcano plot shown in Fig. 5b in which only the relationships for the most difficult reaction step to complete, the potential determining step, are depicted. The potential determining step is ascertained from eqn (7).
ΔG(pds) = max[ΔGRxn(2 → 3), ΔGRxn(3 → 4), ΔGRxn(4 → 5), ΔGRxn(5 → 6), ΔGRxn(6 → 7), ΔGRxn(7 → 2)] | (7) |
Put more simply, eqn (7) states that the reaction between catalytic cycle intermediates having the most negative −ΔG(pds) value (y-axis) for any descriptor value (x-axis) represents the potential determining step. Thus, the nature of the potential determining step directly relates to how strongly a catalyst binds a substrate [i.e., the x-axis value, ΔGRRS(5) delineates which ΔGRxn is the most energetically difficult]. When the catalyst binds substrates very strongly [i.e., very exergonic binding energies for ΔGRRS(5)] the reaction moving between 7 → 2 is potential determining (green line). If a catalyst still overbinds substrates, but to a lesser extent then 3 → 4 is potential determining (red line). When a catalyst weakly underbinds substrates 6 → 7 is potential determining (purple line) and when substrates are strongly underbound 4 → 5 is potential determining (pink line). These reactions divide the Fig. 5b volcano plot into four regions (I–IV), which correspond to two unique situations. The two regions located to the left of the volcano peak (I and II), are commonly known as the volcano's “strong binding” slope. As their name implies, catalysts in this region bind substrates too strongly. Thus, the release of the final product (7 → 2) or hydrogen atom release by the metal associated with moving from a π- to a σ-complex (3 → 4) are thermodynamically difficult and serve as potential determining steps. Conversely, the right slope (regions III and IV), commonly referred to as “weak binding”, characterizes situations where the substrates bind too weakly to the metal centre. Catalysts falling in those regions have problems completing mechanistic steps requiring the addition of new ligands, including H2 (6 → 7) and CO (4 → 5) coordination.
Once the slopes of the volcano plot representing the potential determining steps have been identified, visualizing the points representing catalysts only involves plotting Cartesian coordinates [ΔGRRS(5), −ΔG(pds)] obtained from the previously computed free energy profiles.57 Examining the placements of the individual catalysts in Fig. 5b already confirms known experimental trends, even from this solely thermodynamic picture. For example, Rh catalysts (indicated by “▼” in Fig. 5b) appear near the volcano peak, a placement corresponding to a thermodynamic free energy plot for these catalysts in which all steps are exergonic or roughly thermoneutral, including the potential determining step.58 Indeed, the energetic assessments indicated in Fig. 5b align well with the experimentally established enhanced activity of Rh-based catalysts. Other catalyst appearing high on the volcano, such as those containing Ir (+) and Co (▲) metals likewise show good thermodynamic free energy profiles, which echoes their experimental capabilities.59–61 Note that all of the group 9 metals tend to cluster around the top of the volcano, indicating that catalysts incorporating these metals are more ideal in terms of Sabatier's principle, as the intermediates are neither bound too strongly nor too weakly.
While group 9 metals closely follow Sabatier's principle and, thus, appear high on the volcano plot, metals belonging to groups 8 and 10 uniformly appear in the lower regions. The placement of catalysts into these regions means the free energies required to move between catalytic cycle intermediates are energetically uneven (i.e., a combination of more highly endergonic and exergonic steps) for these species.54 Correspondingly, the potential determining step is more energetically costly and these catalysts have overall less favourable thermodynamics. Group 8 metals uniformly appear on the strong binding (left) side of the volcano, meaning that reducing the number of coordinated ligands, such as releasing products, is energetically very costly. Steps of this nature are particularly problematic for Fe and Os, with Ru appearing to suffer less. This slight shift of Ru-catalysts toward the volcano peak corroborates their known uses for hydroformylation,31,62,63 while Fe (ref. 64) and Os (ref. 65–67) are less commonly employed. Group 10 metals suffer from the opposing problem of their group 8 counterparts: rather then binding intermediates too strongly; adding the necessary intermediates (H2 and CO) to complete the catalytic process is energetically costly. As a result, group 10 metals uniformly appear on the weak binding (right) side of the volcano plot. Indeed, Ni,68 Pd,69,70 and Pt (ref. 71) are all capable of catalysing hydroformylation reactions, but generally require co-catalysts72,73 or proceed through alternative mechanisms than that studied here.31
While the influence of the metal centre is quite clear, the effect of ligand size is more difficult to discern. In general, adding bulkier ligands destabilizes the intermediates, which is seen by a shift to weaker binding (i.e., a left to right movement on the volcano plot as steric bulk is increased). From a purely thermodynamic perspective, group 8 catalysts appearing on the strong-binding side of the volcano (left side) will benefit from larger ligands, since bulkier ligands destabilize the intermediates (causing a rightward shift on the volcano plots) that would push these species closer to the volcano peak and, correspondingly, reduce the endergonicity of the potential determining step. On the other hand, group 10 metals already bind too weakly. Further destabilization of the intermediates via the inclusion of bulkier ligands would only exacerbate the problem, leading to even less favourable thermodynamic energy profiles (i.e., more endergonic free energies for the potential determining step).
Looking at the influences of ligand size for particular metal species also reveals key pieces of information regarding how volcano plots can best be used for screening new catalysts. For Rh, the two intermediate-sized ligands (PMe3 and PPh3) have better thermodynamics than the smallest (PH3) and bulkiest (PCy3) ligands, while for Ir and Co, thermodynamic factors appear less influenced by ligand bulk, as the potential determining steps for each species are predicted to be roughly equivalent. Of course, the lack of any discernable trends regarding ligand choice is unsurprising, given that differences in ligand bulkiness should minimally influence thermodynamic aspects of the catalytic cycle for the case studied here. However, some aspects of ligand bulkiness, particularly when favourable non-bonded interactions involving large ligands are used to promote reactions or form specific products,74–78 as is known to occur in organocatalysis,79,80 may be observed in thermodynamic assessments of catalytic behaviour. Fortunately, checking the behaviour of bulkier catalysts from a thermodynamic perspective can be accomplished relatively easily, since only minima on the potential energy surface need to be computed. By plotting the LFESRs for a series of catalysts containing a combination of ligands with greater and lesser amounts of steric bulk any steric factors that arise in the thermodynamics would quickly be visible as outlier points in the LFESRs.
Importantly, the generally small impact of ligand bulk indicates that in future studies aiming to better understand the thermodynamics of different catalysts, smaller, less bulky ligands should be used (e.g., PMe3 to represent phosphines, OMe2 to represent ethers, etc.) since their bulkier alternatives often give similar thermodynamic profiles but require considerably more computational time to obtain. However, the same situation should not be observed in kinetic volcanoes, which should show significant influences that arise from ligand bulk. To determine if this is true, we turn toward constructing kinetic volcano plots to analyse the different catalysts.
ΔG‡(kds) = max[ΔG‡(3 → TS3,4), ΔG‡(5 → TS5,6), ΔG‡(6 → TS6,7), ΔG‡(7 → TS7,2) | (8) |
Fig. 6 (a) Overview of the olefin hydroformylation kinetics without considering ligand size, as derived from the linear free energy scaling relationships that define the catalytic cycle depicted in Fig. 2. (b) Volcano plot illustrating the kinetic suitability of different catalysts derived from (a). Reactions that define the kinetic determining step are 3 → TS3,4 for region I and 6 → TS6,7 for region II. Lines defining the volcano are obtained by taking the lowest −ΔG‡(kds) value amongst all the reaction for each ΔGRRS(5) value. |
As can be seen in the comprehensive kinetic overview provided in Fig. 6a, the LFESRs defining the free energies of 3 → TS3,4, 5 → TS5,6, and 7 → TS7,2 run parallel to one another and are also very close energetically. This means that the barriers to overcome the transition states associated with converting the olefin from a π-bound to a σ-bonded alkyl group (3 → TS3,4), inserting CO into the metal–alkyl bond (5 → TS5,6), and inserting an H atom into the metal CO bond (7 → TS7,2), are roughly equivalent. Since the largest reaction barriers are associated with moving from 3 to TS3,4 (i.e., it is the lowest line on the left slope of Fig. 6a), this reaction step is used to define the kinetic determining step for region I in Fig. 6b.81 Indeed, this matches the strong-binding nature of the left side of the volcano plot, since moving from 3 to 4viaTS3,4 should be influenced by how strongly the metal centre binds the H atom. Likewise, 6 → TS6,7 defines the slope of the kinetic volcano in region II (Fig. 6b), which aligns with the weak binding nature of the H2 addition that characterizes moving from 6 to 7viaTS6,7. Using these newly created volcano plots, it is possible to assess the different catalytic systems based on a kinetic picture and to compare and contrast trends predicted by thermodynamic and kinetic pictures.
Many of the same trends identified in the thermodynamic volcano (Fig. 5b) are also seen in its kinetic counterpart (Fig. 6b). For instance, group 9 metals tend to lie higher on the volcano than those of group 8 and 10. The peaks of both the kinetic and thermodynamic volcanoes each indicate the ideal strength of the catalyst/substrate interaction (x-axis) is the same [i.e., ΔGRRS(5) value of ∼−10 kcal mol−1]. Therefore, individual catalysts appear in the same areas of both volcano plots (i.e., group 8 – strong binding, group 10 – weak binding). However, the relative energetic ordering of the metals has changed: thermodynamically Rh-based species had the best energetic profile, yet catalysts containing Co centres appear slightly better when kinetics are considered (i.e., they lie higher on the volcano plot). As an example, the kinetic determining step for RhPCy3 (6 → TS6,7) has a barrier height for the potential determining step that is ∼8 kcal mol−1 higher than the CoPCy3 species. Ir-based catalysts show less appealing kinetic profiles and uniformly appear in the lower regions of the kinetic volcano.
Disappointingly, the effects of ligand bulkiness are only somewhat clear from the kinetic volcano. For Co and Rh the bulkiest ligand (PCy3) have a more favourable kinetic determining step than for less bulky ligands, despite being shifted away from the top of the volcano. Indeed, this trend is relatively robust, also appearing for Ru (group 8) and Ni (group 10), but not for Ir, which shows an opposing trend. Of course, using LFESRs that combine the energetic descriptions of the kinetic aspects of all phosphine ligands without consideration of their size may lead to a somewhat muddled picture, in which the influences of ligand bulk cannot be properly ascertained. Considering the LFESRs and volcano plots in a different way may unravel the effects that sterically bulky ligands play on reaction kinetics, thus, it is necessary to take a closer look at these relationships.
Fig. 7 Revised linear free energy scaling relationships amongst intermediates for the kinetic steps of olefin hydroformylation separated by ligand type. |
Also note that these newly derived relationships benefit from greatly improved correlations [higher R2 values and smaller standard deviations (2σ), see ESI Fig. S2†] that arise when the ligands are classified individually. Because of this, volcano plots created from ligand-specific LFESRs will provide a more accurate description of the catalytic cycle. The importance of this point cannot be understated, as having LFESRs capable of accurately predicting the relative free energies of catalytic cycle intermediates and transition states for catalysts with unknown free energy profiles is paramount for the rapid, yet accurate, screening of new catalysts. Having distinguished that different ligand sizes produce distinct LFESRs, we now construct ligand-specific kinetic volcano plots.
Fig. 8 shows these new ligand-specific kinetic volcanoes, which appear dramatically different from one another. Notably, the height of each volcano increases moving from the least to the most bulky ligand. This situation illustrates that the barrier height of the kinetic determining step decreases when bulkier ligands are used. While the shifts are rather minor amongst PH3, PMe3, and PPh3, the change is vivid for the bulkiest PCy3 ligand, corresponding to a transition state barrier of around 12 kcal mol−1, as opposed to 18 kcal mol−1 for the smaller PH3 ligand. Importantly, this finding shows that increasing the bulkiness of the phosphine ligand results in a more appealing catalysts, which is only identified when examining kinetic aspects of the catalytic cycle. Fig. 8b illustrates that the location of actual catalysts on the volcano plots closely correspond to the predicted energies based on LFESRs, indicating that kinetic volcanoes would also be useful for predicting new catalysts. Having established that both general thermodynamic volcano plots and ligand-specific kinetic volcano plots identify two important yet distinct aspects of the catalytic cycle energetics, we proceed to demonstrate how these type of plots can be used to identify and predict the behaviour of new catalysts.
Fig. 8 (a) Overview of the hydroformylation kinetics separated by ligand type, as derived from the linear free energy scaling relationships that define the catalytic cycle depicted in Fig. 2. (b) Volcano plot illustrating the kinetic suitability of catalysts derived from (a). |
While thermodynamic studies will often provide the necessary amount of data to move directly to the laboratory, it can be envisioned that in some instances more information regarding the role steric interactions play in dictating activation barrier heights will be desired. In these cases it becomes necessary to construct kinetic volcano plots. While we have already shown that kinetic volcanoes distinguish subtle aspects related to catalytic cycle energetics and reproduce known experimental trends, their principle drawback is their considerable computational cost. As such, an ideal and cheap tool for screening homogeneous catalysts would facilitate access to the same type of information obtained from kinetic volcano plots but without the need to undertake the time-intensive task of computing all the necessary data for numerous catalysts.
One attractive solution to more quickly access kinetic data would be establishing a structure–activity relationship (SAR) between a parameter capable of describing ligand sterics (e.g., Tolman's cone angle,82 Verloop's Sterimol,83etc.) and the descriptor used to create the kinetic volcano plot. Indeed, chemists often use steric parameters to better understand, describe, and even predict different aspects of homogeneous catalysis.77,84–88 A hypothetical procedure to create kinetic volcanoes based on structure–activity relationships (SAR volcanoes) is outlined in Fig. 9. First, the free energy profiles of a series of catalysts bearing the same ligand must be computed (Fig. 9a). The number of individual catalysts tested bearing this specific ligand need only be sufficient to establish reliable linear free energy scaling relationships (LFESRs) in the manner previously outlined (Fig. 9b, black lines). This procedure needs to be completed for several small to medium sized ligands (absolute minimum of two but including more should give more accurate predictions, vide infra), which leads to different sets of LFESRs for each ligand type (Fig. 9b, red versus black lines). Using these ligand-dependent LFESRs hypothetical volcanoes for each ligand can be created following the earlier described procedures (Fig. 9c). Each ligand-specific volcano has a different peak point (black and red squares, Fig. 9c). Additionally, each ligand type has a specific steric descriptor, such as the Tolman cone angle. Using the peak point from each ligand-specific volcano along with the ligand-specific steric descriptor, a mathematical relationship can be established that relates the Cartesian coordinates of the peak position [in our case representative of ΔGRRS(5) and −ΔG‡(kds)] to the steric parameter (in our case the Tolman cone angle), as seen in Fig. 9d. The steric parameters for other ligands, for which no data has been computed, can easily be obtained from a database or determined as needed. By combining the steric parameter values of these untested ligands with the mathematical relationships relating the peak position coordinates to the steric parameters the new peak location for the hypothetical volcano of any ligand can be established (Fig. 9e). Finally, the same slopes obtained from the volcano plots of computed ligands can be used as approximations to construct a hypothetical SAR kinetic volcano (Fig. 9f).
As an example of the procedure to create SAR kinetic volcanoes, the LFESRs for the PH3 and PPh3 ligands were obtained by computing the free energy profiles (including transition states) of just eight catalysts [Fe, Ru, Os, Co, Rh, Ir, Fe, Pt] for each ligand type. The resulting LFESRs describing those ligands are depicted in Fig. 7 (PH3: black lines, PPh3: purple lines). Using only those LFESRs, the corresponding volcano plots can be created, which appear in Fig. 8a (PH3: black, PPh3: purple). Cartesian coordinates defining the peak positions [i.e., ΔGRRS(5) and −ΔG‡(kds)] of each of the Fig. 8a simulated kinetic volcanoes can then be related to the Tolman cone angle using a linear fit (analogous to Fig. 9d). Knowing the Tolman cone angle for other phosphine ligands (e.g., see database in ref. 82), the peak points for the new ligand-specific kinetic volcano plots can be derived. As a final step, the sides of the volcanoes can be approximated by using the slopes from computed volcanoes of ligands with the closest Tolman cone angle (for this example, slopes from either PH3 or PPh3). Following this procedure produces the SAR kinetic volcanoes for the PMe3 and PCy3 ligands shown in purple in Fig. 10. As can be seen, the agreement between the black volcanoes, which were derived directly from the computed data for the PMe3 and PCy3 ligands, and their estimated counterparts is quite satisfactory. While some deterioration is seen for the largest ligand, particularly in the prediction of the x-coordinate, the projected value of kinetic determining steps (y-axis) differ only slightly. Despite the observed reduction in accuracy, it should be stressed that these estimations were obtained with a considerable reduction of computational time since no data for either the PMe3 or the PCy3 ligands was needed to derive the purple Fig. 10 plots.
To demonstrate the full power of the above-described methodology, we created SAR kinetic volcanoes for several phosphine ligands for which we computed no data (purple lines, Fig. 11). The three representative ligands include one having very little bulk, P(OMe)3, one extremely bulky, P(Mes)3, and one of intermediary bulk, P(NMe2)3. Using these SAR volcanoes, the predicted kinetic determining step for any catalyst bearing a phosphine ligand can be estimated solely by computing the descriptor value [ΔGRRS(5)]. Such a tool is remarkably appealing, as the approximate value of the highest activation barrier encountered in the catalytic cycle for any ligand can be quickly identified. Thus, accessing the kinetic information of synthetically appealing catalysts identified from thermodynamic volcanoes can be accomplished with very little computational expense.
Fig. 11 Structure–activity relationships (SAR) volcanoes for selected phosphine ligands computed using structure–activity relationships fit from four (green) and two (purple) data points. |
One potential strategy for improving the accuracy of predictions made from SAR volcanoes is including additional ligands into the fitting procedure (e.g., complete the Fig. 9a–d steps with more ligand types). Indeed, using a 4-point (green lines, Fig. 11), as opposed to the earlier employed 2-point fit (purple lines, Fig. 11), changes the peak location of the volcano that corresponds to the free energy associated with the kinetic determining step (y-axis) and the ideal descriptor value (x-axis). The quality of the estimates are less affected for P(OMe)3, while more significant deviations exist for P(Mes)3. It must again be stressed, however, that the improvements made by using a fit from more ligand points come at a significant increase in computational time. Users are encouraged to undertake studies that maximize the predictive ability of SAR kinetic volcanoes while simultaneously minimizing computational expense.
Footnote |
† Electronic supplementary information (ESI) available: Details on deriving the linear scaling relationships, creating the volcano plots, and building the structure–activity relationships volcano plots, as well as computed values of the catalytic cycles using the PBE0-dDsC and M06 functionals are included. Cartesian coordinates of relevant compounds are included (.xyz). See DOI: 10.1039/c6sc01660j |
This journal is © The Royal Society of Chemistry 2016 |