Steven M.
Maley‡
a,
Doo-Hyun
Kwon‡
a,
Nick
Rollins
a,
Johnathan C.
Stanley
a,
Orson L.
Sydora
b,
Steven M.
Bischof
*b and
Daniel H.
Ess
*a
aDepartment of Chemistry and Biochemistry, Brigham Young University, Provo, Utah 84602, USA. E-mail: dhe@chem.byu.edu
bResearch and Technology, Chevron Phillips Chemical Company LP, 1862, Kingwood Drive, Kingwood, Texas 77339, USA. E-mail: bischs@cpchem.com
First published on 21st August 2020
The use of data science tools to provide the emergence of non-trivial chemical features for catalyst design is an important goal in catalysis science. Additionally, there is currently no general strategy for computational homogeneous, molecular catalyst design. Here, we report the unique combination of an experimentally verified DFT-transition-state model with a random forest machine learning model in a campaign to design new molecular Cr phosphine imine (Cr(P,N)) catalysts for selective ethylene oligomerization, specifically to increase 1-octene selectivity. This involved the calculation of 1-hexene:1-octene transition-state selectivity for 105 (P,N) ligands and the harvesting of 14 descriptors, which were then used to build a random forest regression model. This model showed the emergence of several key design features, such as Cr–N distance, Cr–α distance, and Cr distance out of pocket, which were then used to rapidly design a new generation of Cr(P,N) catalyst ligands that are predicted to give >95% selectivity for 1-octene.
We recently reported the development and use of a density functional theory (DFT) transition-state model that provided quantitative prediction of molecular Cr catalysts for controllable selective ethylene trimerization and tetramerization (Scheme 1).30 This selective catalyst design effort is important because the ubiquity of polyethylene resulting from robust Phillips31 and Ziegler–Natta catalysts32,33 has led to an increase in global need for linear α-olefin (LAO) polymerization comonomers.34 Additionally, these short-chain LAOs, specifically 1-hexene and 1-octene, are used in the manufacture of plasticizers, lubricants, detergents, and plastomers/elastomers. Prior to our work, Sydora and co-workers reported a series of aryl and benzyl substituted phosphine imine (P,N) catalysts that experimentally produced ∼30% 1-octene.35 In our computational catalyst design effort, we used the cationic high-spin transition-states TS1 and TS2 shown in Scheme 1b to develop a linear correlation model between DFT computed values and experimental 1-hexene:1-octene ratios (see later Scheme 2 for complete catalytic cycle). This allowed using transition-state calculations to computationally design a new general class of phosphine monocyclic imine Cr(P,N) catalysts where changes in the ligand structure control 1-hexene versus 1-octene selectivity (Scheme 1b). Experimental ligand and catalyst synthesis, and reaction testing, quantitatively confirmed our transition-state predictions.
Scheme 2 Outline of Cr-catalyzed metallacycle mechanism for selective ethylene oligomerization to 1-hexene and 1-octene. |
While our DFT transition-state model is practical, accurate, and successfully identified new ligands that were experimentally validated, it did not offer general or specific catalyst design guidance to further enhance 1-octene selectivity. Moreover, the interpretation of singular controlling transition-state features is unobvious since the energy difference between 1-hexene and 1-octene selectivity is relatively small. Therefore, we decided to combine our transition-state model with data science methods that can potentially provide the emergence of chemical features to enhance 1-octene selectivity. While machine learning and related multi-dimensional methods are beginning to be used for molecular catalyst design,36 there are currently no examples of an experimentally verified quantum-mechanical transition state model merged with machine learning methods for catalyst selectivity. Contemporary to our work, Balcells and Azpuru-Guzik very recently reported combining machine learning and transition states for analysis of stoichiometric H2 oxidative addition reactivity by Vaska-type complexes.37 Currently, our work reported here and the work of Balcells and Azpuru-Guzik represent the only reports of uniquely combining transition states and machine learning. Generally, data science approaches to molecular catalyst design emphasize ground-state properties of either pre-catalysts or ligands without metal centers. For example, Fey and Pringle developed databases of ground-state ligand properties38–43 that enabled the prediction and experimental verification of new fluorophosphine ligands for hydroformylation and hydrocyanation.44 Recently, Denmark reported a workflow where a subset of a library of synthetically accessible catalysts are selected and tested. The data obtained was then used to train statistical learning models to optimize a chiral catalyst for thiol addition to N-acylimines.45 Machine learning methods have also been used to predict reaction barrier heights in heterogeneous catalyst applications.46–48
Outlined in Scheme 1c, our approach reported here involves using DFT-computed transition-state features and selectivities for machine learning analysis. The analysis of >100 Cr(P,N) catalysts and 14 molecular descriptors through machine learning regression algorithms with multifold cross validation resulted in a low root mean square error (RMSE) and emergence of three critical design elements to enhance 1-octene selectivity. The utility of these machine-learning identified selectivity features was demonstrated by the design and calculation of several new ligands that are predicted to give >90% 1-octene selectivity.
As discussed in the introduction, in our recent computational campaign, we developed a DFT transition-state model that provided quantitative prediction of molecular Cr(P,N) catalysts (Scheme 1).30 This resulted in the computational design and experimental verification of a new general class of phosphine monocyclic imine Cr(P,N) catalysts where slightly more than 50% 1-octene was formed. Importantly, this DFT transition-state model is accurate and relatively fast to use. New ligands can be virtually screened by calculating the energies of transition-state conformation ensembles for TS1 and TS2 (Scheme 1b),52–63 which can then be used in our linear correlation model to provide a predicted 1-hexene:1-octene ratio for given set of reaction conditions. Our transition-state selectivity model was based on our DFT mechanistic calculations, and studies by Britovsek and McGuinness, Cheon, Liu and others.52,53,56–63 Briefly, pre-catalyst activation results in a bis-ethylene CrI species A and oxidative C–C bond formation then results in chromacyclopentane B (Scheme 2). A third ethylene coordination gives C and migratory insertion leads to the chromacycloheptane intermediate D where 1-hexene and 1-octene pathways diverge. Through β-hydrogen transfer TS1 there is formation of 1-hexene and Cr catalyst reduction. The 1-octene pathway requires a fourth ethylene coordination and intermediate E followed by migratory insertion TS2 to give the chromacyclononane intermediate F that leads to 1-octene.
While our DFT transition-state model is extremely useful, because it is a correlation model with small energy differences, it did not reveal general catalyst design principles that would easily lead to the design of extremely high 1-octene selectivity. Because our transition-state selectivity model is accurate, this provided a platform for combining it with machine learning data science methods that can provide the emergence of general chemical features to enhance 1-octene selectivity for Cr(P,N) catalysts. To our knowledge, there has not been a previous report of combining an experimentally verified quantum mechanical transition-state model for selectivity with machine learning analysis as a workflow to design new catalysts.
While machine learning is beginning to find significant application as a tool to aid organic synthesis,64–74 there are relatively few examples of machine learning applied to inorganic or organometallic reactions, especially heterogeneous46,75–87 and homogeneous catalysis.88,89 Kulik trained an artificial neural network to predict the high-spin to low-spin splitting energies of ∼2700 transition metal complexes.90 They also demonstrated the usage of a kernel ridge regression model for predicting spin-splitting, bond lengths, and redox potentials for a relatively large collection of transition metal complexes.91 Related, Corminbouef trained machine learning models to screen over 18000 potential homogeneous catalysts for the Suzuki–Miyaura C–C cross-coupling reaction.92 Sunoj used a combination of a neural network and random forest model to identify the regioselectivity of catalytic difluorination of alkenes.93 Brgoch screened over 100000 compounds using a support vector machine regression to identify novel highly compressible metal materials,94 and Buehler used convolutional neural networks to search for new composite metal materials.95 Xin used artificial neural networks to identify heterogeneous metallic catalysts for CO capture and reduction.96 Most directly related to our work, a recent and important contribution by Balcells and Azpuru-Guzik demonstrated the use of machine learning for reactivity to screen ligands of Vaska-type complexes to identify optimal ligands for H2 oxidative addition. This work showed that machine learning identified atom size and electronegativity as key features impacting H2 oxidative addition reactivity.37
With the success of previous machine learning studies predicting spin-splitting energies and redox potentials,90,91 we were relatively confident that one or more machine learning algorithms would perform well for our workflow. Similar to our previous computational studies,30 and computational assessments by McGuinness suggesting the good performance of M06L,52,53 we used the unrestricted M06L density functional97 for describing the electronic structure of Cr(P,N) catalysts. The UM06L/Def2-TZVPP//UM06L/6-31G**[LANL2DZ] level of theory was combined with the SMD continuum model98 for cyclohexane to estimate the free energies of TS1 and TS2 (see ESI† for computational details). In this transition state model, the relative free energies of TS1 and TS2 provide selectivity under the assumption of Curtin–Hammett type conditions. All transition-state structures were optimized, and vibrational frequencies were computed to verify the stationary points as first-order saddle points. Normal rigid-rotor harmonic oscillator approximations were applied with free energies at 1 atm and 298 K. Because the transition-state model is a linear correlation scheme, no temperature or pressure corrections were applied. All DFT calculations were performed using Gaussian 09.99 Machine learning analyses were performed using scipy,100 numpy,101 pandas102 and scikit-learn103 in Python 2.7. A detailed description of the machine learning analysis along with the source code are available in the ESI.†
We used 105 unique (P,N) ligands in our transition-state training data set, which included 14 experimentally measured ligands. Scheme 2 outlines these (P,N) ligands, which have a variety of different functional groups, but retain the phosphine and imine, or imine-like, ligand coordination for which the DFT transition-state model was developed. A major motivation in the selection of these 105 ligands was to further optimize the five-membered imine ring system we previously designed and to stay within the bounds of predictability for our correlation model. As shown in Scheme 3, this set includes a variety of substituted heterocycles such as pyrroles, imidazoles, and oxazoles. In each of these cases we also examined their combination with alkyl, fluoroalkyl, aryl, and amido phosphines as well as phospholanes.
Fourteen atomic and molecular descriptors/features were extracted from TS1 and TS2 for each of the 105 ligands shown in Scheme 3. Described in Fig. 1, extracted features included geometric parameters such as bond lengths, angles, dihedrals, percent volume buried, and Cr metal center distance out pocket. Several electronic features were also harvested, such as electrostatic-based atomic charges. Percent volume buried describes the extent to which the first coordination sphere of the Cr metal center is occupied by a (P,N) ligand.104 The distance out of pocket describes the how far the Cr metal is situated from the (P,N) ligand.
Fig. 1 Descriptors extracted for machine learning analysis. (A) Geometric descriptors and electrostatic charges. (B) Definition of percent volume buried. (C) Definition of distance out of pocket. |
The Scikit-Learn python library was used to set up and train regressors on this transition-state data set, which was split into 25% training and 75% testing sets. Seven regression algorithms were tested including: random forest, Gaussian process regression, LASSO, elastic-net, ridge regression and support vector regression with both a linear and radial basis function kernel. Multifold cross validation was performed to protect against model overfitting common in small datasets. This random sampling was performed 10 times and 20-fold cross validation was used at each iteration to determine regression accuracy. A brief test of different amounts of random sampling and cross validation was also performed (see ESI†). The RMSE of each model determined using cross validation averaged across iterations is shown in Fig. 2.
The machine learning regression algorithms were used to evaluate the use of the 14 atomic and molecular features to quantitatively correlate with the DFT calculated energy differences between TS1 and TS2. The RMSE of the regression algorithms ranged from 0.344 to 0.568 (Fig. 2). The best performing model was random forest (RMSE = 0.344) and the poorest performing model was elastic-net (RMSE = 0.568). Unsurprisingly the LASSO and ridge algorithms, which are related to elastic-net, also performed poorly with RMSEs of 0.565 and 0.494, respectively. The performance of support vector regression improved by almost 10% when changing from a linear (RMSE = 0.491) to a radial basis function (RMSE = 0.408) kernel. Gaussian process regression performed comparable to SVR-rbf.
Related to the top performance of random forest for our correlation of 1-hexene:1-octene ratios, Doyle reported that random forest outperformed linear models, SVR, k-nearest neighbor, and artificial neural networks at predicting yields of Pd-catalyzed Buchwald–Hartwig cross-coupling of aryl halides with 4-methylaniline.70 Also, for spin-splitting energies and bond lengths of transition metal complexes, Kulik achieved exceptionally low mean unsigned errors using random forest.91 Palmer et al. showed that random forest outperforms SVR and artificial neural networks in predicting the aqueous solubility of organic molecules.105
Because the random forest algorithm performed well in our case, and the applications mentioned above, this algorithm was chosen for further hyperparameter optimization using the GridSearch CV method from SciKit-Learn. Different permutations of hyperparameters and five-fold cross validation were tested in order to determine the set of hyperparameters that maximized the performance of the model. The number of trees in each forest was varied from 20 to 210 and the trees were split from 5 to 125 times. Both mean signed error and mean absolute error were considered when determining the quality of each split and between three and 14 features were examined when considering the best split. The optimized random forest model was then re-fit to the training data to validate the hyperparameter optimization. The RMSE of the random forest model decreased from 0.344 to 0.272 after optimization. The RMSE of the 1-hexene to 1-octene selectivities are 0.275 and 0.269, respectively.
The middle of Fig. 3 plots the selectivities determined from the optimized random model against those determined from the DFT selectivity model. In our data set, overall 1-hexene selective (i.e. >50% 1-hexene vs. 1-octene) is labeled as a negative value and overall 1-octene selective is labeled as a positive value (i.e. >50% 1-octene vs. 1-hexene). The random forest model correctly predicted the overall 1-hexene versus 1-octene selectivity for 83 ligands and incorrectly predicted the overall selectivity for 22 ligands. This incorrect assignment occurs in cases where the DFT computed 1-hexene selectivity of a ligand is less than 1%. The random forest model tends to perform best for ligands ranging from 20:80 to 50:50 1-hexene:1-octene selectivity (see ESI-Fig. 1†). Importantly, the majority of ligands in our data set are overall 1-octene selective, which is useful for our goal of improving the percentage of 1-octene, but it is likely that the RMSE of the random forest model would be reduced if the data set were more evenly distributed between overall 1-hexene and overall 1-octene selective.
At the top and the bottom of Fig. 3, the 1-octene selectivities calculated using the DFT selectivity model and the optimized random forest model are plotted against the experimentally determined selectivities. The DFT selectivity model agreed with experiment very well (R2 = 0.91, mean absolute deviation = 4.4%). The DFT model underestimates the experimental selectivity, however, this is overall advantageous with the goal to increase 1-octene production. In contrast to the DFT calculated values, the random forest model, not unexpectedly, is unable to quantitatively reproduce the experimental selectivities with a high degree of linear correlation. The lack of very high quantitative correlation between random forest and experiment values is likely due to the relatively small sample size of experimentally studied ligands. Despite this, the random forest can be used to determine critical chemical features that are responsible for enhancing 1-octene selectivity. The relative importance of the 14 features included in our dataset can be determined by replacing data with random values and observing the impact on the RMSE value. If replacing data of a feature with random values results in a small change to the RMSE then it has a low degree of importance. Conversely, if there is a large change in the RMSE then this feature has a large importance. Fig. 4 displays this feature importance analysis using the optimized random forest model.
Fig. 4 Normalized feature importance determined from random forest model with 95% confidence intervals (red bars). |
Inspection of Fig. 4 shows that the Cr–N distance, Cr–α distance, and distance out of pocket were identified as being most important in enhancing 1-octene selectivity. The Cr–N–C2 angle, which is related to the Cr–α distance, was also identified as an important 1-octene enhancing feature. Interestingly, despite the proposed importance of the ligand bite angle, especially for phosphine catalysts,106 we found that the P–Cr–N1 ligand bite angle is among several lesser important features, which is consistent with Sydora previously showing that for Cr-phosphinoamidine catalyst ligands with similar bite angles resulted in significantly different 1-hexene:1-octene selectivities.35 After we completed this machine learning analysis and identified the importance of the Cr-distance out of pocket for the 105 ligands examined here, Liu reported a DFT analysis of ethylene oligomerization by a Cr(2,2′-dipicolyamine) catalyst and also found that enhanced tetramerization was likely due to the Cr-distance out of pocket.
With the emergence of chemical features by the random forest machine learning model, we then turned to the final step in our workflow, which is using this information to virtually identify new catalyst ligands (Scheme 1c). Based on the important features, we rapidly designed ligands L9–L15 shown in Fig. 5. The machine learning features directed us to change the 4-membered (P,N) ligand scaffold found in previously reported ligands of generations 1 (L1–L5) and 2 (L6–L8) to a 3-membered (P,N) ligand, which would potentially alter the Cr–N distance, Cr–α distance, and Cr distance out of ligand pocket. Briefly, generation 1 ligands were reported by Sydora35 and based on experimental trial and error with only a single example of 1-octene production in the vicinity of 30%. Generation 2 ligands were computational designed by us using only transition-state energies, and then experimentally verified to show our DFT calculations were extremely accurate. However, this only led to the improvement of 1-octene selectivity to about 50%. Moreover, these transition states did not provide an obvious way to design a new, more selective generation of catalysts. Our machine-learning driven modification led to the proposal of ligands L9–L11, and with our transition-state model they are predicted to be 97–99% 1-octene selective. With the rapid success of this new generation 3 type of (P,N) ligands we further decreased the ligand to have direct phosphine-nitrogen connection, which naturally led to the proposal of 2-membered (P,N) ligands L12–L14. Based on our transition state model, ligands L12–L14 have predicted selectivities of >95% for 1-octene. Importantly, as plotted at the bottom of Fig. 5, the use of our transition-state model combined with translation of machine learning features to new catalyst ligands, resulted in increasing the prediction of 1-octene from between <35% for generation 1 and ∼50% for generation 2 to >95% for generation 3.
With the rapid success of designing ligands L9–L15 it is clear that there are a number of new candidates that are now available for experimental testing, and several more ligands can now be virtually and rapidly designed. As one experimental confirmation of our results, without prior knowledge, and subsequent to our design of L9–L15, a literature search of all reported Cr-phosphine catalysts for ethylene oligomerization revealed that Yang reported that L12 is indeed highly 1-octene selective.107 In this experimental test, which is signficantly different than the Chevron Phillips Chemical reaction conditions that our transition-state model was developed for, in situ catalyst formation with the combination of Cr(acac)3, MAO, and L12 oligomerized under 50 bar of ethylene gave a 1-hexene:1-octene ratio of 28:70 with trace production of polyethylene.
Footnotes |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sc03552a |
‡ Authors contributed equally to this work. |
This journal is © The Royal Society of Chemistry 2020 |