Tianbai
Huang
a,
Robert
Geitner
*b,
Alexander
Croy
*a and
Stefanie
Gräfe
*a
aInstitute for Physical Chemistry (IPC) and Abbe Center of Photonics, Friedrich Schiller University Jena, Helmholtzweg 4, 07743 Jena, Germany. E-mail: alexander.croy@uni-jena.de; s.graefe@uni-jena.de
bInstitute of Chemistry and Bioengineering, Technical University Ilmenau, Weimarer Str. 32, 98693 Ilmenau, Germany. E-mail: robert.geitner@tu-ilmenau.de
First published on 28th May 2024
Transition metal complexes have played crucial roles in various homogeneous catalytic processes due to their exceptional versatility. This adaptability stems not only from the central metal ions but also from the vast array of choices of the ligand spheres, which form an enormously large chemical space. For example, Rh complexes, with a well-designed ligand sphere, are known to be efficient in catalyzing the C–H activation process in alkanes. To investigate the structure–property relation of the Rh complex and identify the optimal ligand that minimizes the calculated reaction energy ΔE of an alkane C–H activation, we have applied a Δ-machine learning method trained on various features to study 1743 pairs of reactants (Rh(PLP)(Cl)(CO)) and intermediates (Rh(PLP)(Cl)(CO)(H)(propyl)). Our findings demonstrate that the models exhibit robust predictive performance when trained on features derived from electron density (R2 = 0.816), and SOAPs (R2 = 0.819), a set of position-based descriptors. Leveraging the model trained on xTB-SOAPs that only depend on the xTB-equilibrium structures, we propose an efficient and accurate screening procedure to explore the extensive chemical space of bisphosphine ligands. By applying this screening procedure, we identify ten newly selected reactant–intermediate pairs with an average ΔE of 33.2 kJ mol−1, remarkably lower than the average ΔE of the original data set of 68.0 kJ mol−1. This underscores the efficacy of our screening procedure in pinpointing structures with significantly lower energy levels.
Computational chemistry serves as a potent tool for ligand design for organometallic catalysts. Reaction energies and activation barriers can readily be obtained by analyzing the energy profiles of molecular configurations on the potential energy (hyper-)surface, calculated by quantum chemical (QC) methods, such as density functional theory (DFT). These values can be linked to the reaction rate by means of the transition state theory.9 Moreover, with the aid of QC methods, the key intermediates and transition states (TSs) of the catalytic reaction can be identified. Once the ligand sphere of the transition metal complex is specified, this process can also be accomplished via automated exploration of the chemical reaction network (CRN).10–16 However, the massive number of possible combinations of the building blocks of the ligands leads to a vast chemical space with varying properties.17 Due to the high computational cost, QC methods become impractical to screen thousands of key intermediates and TSs for the reaction of interest to find an optimized ligand structure.
By virtue of high computational efficiency, machine learning (ML) techniques have emerged as complements to QC methods, and have been successfully applied in drug discovery,18,19 as well as in screening the properties of metal–organic frameworks20 and transition metal complexes.21 ML techniques are also widely applied in the investigation of chemical reactions.22 For instance, Choi et al.23 predicted the activation barriers of reactions from the RMG-py12 database with a mean absolute error (MAE) of 1.95 kcal mol−1. Ye et al.24 predicted the activation barriers of a reductive elimination step of Pd-catalyzed C–N coupling reactions with a MAE of 0.48 kcal mol−1. In addition, the activation barriers in the CRN of formamide, aldol, and unimolecular decomposition of 3-hydroperoxypropanal were accurately predicted using an artificial neural network (ANN) model.25 It is noteworthy that the ML models used in these studies incorporated the thermodynamic properties of the products and reactants as descriptors for the reaction, aligning with the Bell–Evans–Polanyi principle26 and leading to high prediction accuracy. The utilization of QC-free molecular descriptors was also proven to be successful in predicting the activation barriers of glutathione adduct formation27 and dihydrogen activation.28
Despite the successes in predicting the activation energy of elementary reactions, the prediction of energy differences between the reactant and the intermediates in the entire catalytic cycle, which is important for multistep catalytic reactions, is uncommon in literature. The challenge remains in effectively modeling the reaction energy ΔE, namely, the energy difference between the reactant and the intermediate within the same elementary reaction. This challenge may stem from the large number of conformers of both reactants, intermediates, and products, making it difficult to pair the proper reactant-product or reactant–intermediate pairs. Consequently, the improper pairing can lead to inaccuracies in the training data that are fed into the ML models.
The present study aims to predict the energy difference, denoted as ΔE, between the 6-coordinated metal alkyl-hydride (Rh(PLP)(Cl)(CO)(H)(alkyl)) and the 4-coordinated precursor (Rh(PLP)(Cl)(CO)) featuring different bidentate phosphine ligands PLP. According to the theoretical and experimental mechanistic studies,29–31 the Rh complex undergoes a series of transformations including CO dissociation, C–H activation, CO recombination, C–C coupling and reductive elimination, where C–H activation is identified to be the one of the rate-limiting steps. According to the Bell–Evans–Polanyi principle, the activation energy ΔE‡ is highly correlated to the reaction energy. A complex with low ΔE for the C–H activation usually features a low ΔE‡ of this step, which benefits the entire catalytic cycle. In addition, the 6-coordinated intermediate is one of the relatively stable intermediates after the C–H activation and CO recombination, which can proceed with a subsequent C–C formation step in the carbonylation process. In this regard, lowering the ΔE could also increase the equilibrium concentration of the 6-coordinated intermediate, accelerating the subsequent C–C coupling reaction. Therefore, to minimize ΔE of the C–H activation by varying the ligand structure of the complexes is the primary objective for designing a complex with high catalytic activity. In this context, an efficient screening scheme of ΔE between the 4-coordinated reactant and the 6-coordinated intermediate is of great importance.
We demonstrate in the present study that the Δ-ML approach,32 which has been successfully applied in predicting the activation barriers and reaction enthalpies of the “breaking two bonds, forming two bonds” type reactions,33 is a possible way for the reaction energy prediction. By training on diverse sets of descriptors, the Δ-ML models can obtain a well-balanced performance between the efficiency and the accuracy of the prediction. We demonstrate that this screening approach allows us to identify ten new bisphosphine ligands, which correspond to ten reactant–intermediate pairs with an average ΔE remarkably lower than the average ΔE of the original data set. Thus, we are able to identify ten new promising Rh-bisphosphine catalysts for C–H activation.
ΔPt(Rt) ≈ Δbt(Rb) = b(Rb) + FML(σ(Rb)). | (1) |
In our study, we have employed this Δ-ML methodology to investigate reaction energies associated with the C–H activation process mediated by Rh complexes. In this context, the properties obtained at the DFT level serve as the target values to be predicted while those obtained at the computationally much more efficient semi-empirical GFN2-xTB34 level of theory (further denoted as xTB throughout this study) are used as the baseline values.
Firstly, to account for the errors that are solely introduced by the different levels of theory, namely DFT vs. xTB, we examine the upper bound of the reaction driving force, denoted as (see Fig. 1a). We denote the corresponding Δbt-model as . The prediction task can be formulated as follows,
(2) |
(3) |
Beyond the baseline value , the feature vector σ of the reactant–intermediate pair is constructed by incorporating the structural information from both reactant and intermediate. In this context, we define σ as σ(RDFT,r,RxTB,i), which helps to account for the difference between and . Referring to Fig. 1a, RDFT,r and RxTB,i represent the geometries of the DFT-equilibrium reactant and xTB-equilibrium intermediate, respectively. Topological-based descriptors, such as autocorrelation functions (ACs)35,36 or position-based descriptors, such as smooth overlap atomic positions (SOAPs),37,38 can be used to construct the feature vector σ. If electronic information is also taken into consideration, the feature vector can be further expanded as σ = σ(RDFT,r, RxTB,i, ρDFT,r, ρxTB,r, ρxTB,i), where ρ depends on the electron density information obtained at different level of theories (see AIM-AC in Section 2.2 for details).
Secondly, another Δbt-model, denoted as ΔDx-model, is trained to predict the reaction energy ΔEDFT (see Fig. 1b)
ΔEDFT = EDD,i − EDD,r ≈ ΔDx(ΔExTB, σ), | (4) |
ΔExTB = Exx,i − Exx,r. | (5) |
(6) |
(7) |
Furthermore, ACs can be evaluated on a specified atom i, which helps to describe the atomic environment. This is defined as the atom-specific ACs (a-ACs):
(8) |
(9) |
The illustrations of standard ACs, rest-ACs, l-ACs and a-ACs are shown in Fig. S1.† To offer a detailed description that focuses on the local reaction site, the a-ACs of selected atoms (Rh, P1, P2, Cl, R11, R12, R21, R22, L1, L2) were employed to construct the feature vector for the reactant, while a broader set of atoms (Rh, P1, P2, Cl, R11, R12, R21, R22, L1, L2, H, C) was used for the intermediate (see Fig. S3† for details). The chemical environment of the reaction site, i.e., Rh and Cl ions in our study, is expected to undergo significant changes before and after reacting with the C–H bond. Therefore, including the a-ACs of these specified atoms, which emphasizes the description of the atomic environment, is deemed beneficial for characterizing the reactant–intermediate pairs and the associated reaction energy. In addition, we also incorporated the global information of the Rh complex, by including the rest-ACs and l-ACs. All ACs are solely dependent on the graph structures of the complex, therefore the feature vectors in both - and ΔDx-models take the form of σ(Rr,Ri).
In the framework of AIM theory,40,41 the nuclei in a molecule are naturally separated in atomic basins Ω, the boundary of which is defined by a zero-flux surface in the gradient vector field of the electron density. The partitioning of the molecular space into atomic basins enables the assignment of the molecular properties to individual atoms. Details of the atomic properties used in this study are listed in Table S1 and explained in the ESI.† The AIM analysis, as implemented in Multiwfn (3.8),42 can be applied to both, the DFT and xTB calculations. Consequently, including the AIM-ACs of the reactant calculated at DFT level of theory provides additional information for describing the reaction energy. In this context, both feature vectors for - and ΔDx-models not only depend on geometries RDFT,r, RxTB,i and RxTB,r, but also on the electron information ρDFT,r, ρxTB,r and ρxTB,i.
Given SOAPs' ability to capture the 3D structure of the molecules, SOAPs evaluated on the xTB-equilibrium reactant are included to account for the energy difference caused by the change in geometry, when predicting the reaction energy ΔEDFT (recall Fig. 1b). Consequently, the feature depends on RDFT,r, RxTB,r and RxTB,i in the ΔDx-model, whereas in the -model, the feature vector of the reactant–intermediate pair depends only on RDFT,r and RxTB,i, given that RDFT,r and RxTB,r are equivalent.
The geometry guesses required to optimize the reactant and intermediate states were generated using a combination of RDKit (2022.03)45 and OpenBabel (3.0).46 We initiated the construction of reactants and intermediates geometries from an octahedral Rh(PLP)(Cl)4 prototype, where the linkers L, as well as the residual groups R1 and R2 were selected from a comprehensive set of 19 different linkers and 62 residual groups (see Fig. S5†). To address potential conformational variations, conformer searches using CREST47 were performed within a subset of prototypes. In addition to the conformer with the lowest energy, several other minima on the xTB potential energy (hyper)surface (PES) were identified in the conformer search processes. Each of these minima represented different prototypic molecules, and would subsequently serve as individual data points for ML analysis.
After the prototypes were optimized at GFN2-xTB (6.6.1)34 level of theory, the initial guesses of 4-coordinated reactants Rh(PLP)(Cl)(CO) were achieved by removing the excess chloride ions in the axial positions. Furthermore, one additional chloride ion in the equatorial plane is replaced by CO. Additionally, the initial geometries of 6-coordinated intermediates Rh(PLP)(Cl)(CO)(H)(propyl) were generated by substituting one chloride ion in the axial position with H, and two equatorial chloride ions with a propyl group and a CO, respectively.
Calculations at density functional of theory (DFT), as well as at xTB level of theory, were carried out sequentially for every system, to obtain the target values ΔEDFT and , as well as the baseline values ΔExTB and for our ML objectives.
Firstly, xTB optimizations were performed on the initial guesses of reactants and intermediates, yielding the xTB single point energies Exx,r and Exx,i. Subsequently, both xTB-optimized species were further optimized at the level of DFT, yielding the single-point DFT energies EDD,r and EDD,i, respectively.
The DFT calculations, employing the range-separated ωB97XD functional48 and the def2-SVP basis set, along with the respective effective core potential,49 were carried out in the Gaussian 16 software package.50 A vibrational analysis was carried out for DFT-equilibrium structures to verify that a minimum was obtained on the PES. Furthermore, the DFT single-point energies ExD,i, were obtained by performing a calculation for each xTB-optimized intermediate while the xTB single-point energies EDx,r, were obtained for each DFT-optimized intermediate. Recall that ExD,i is used to calculate the target upper bound of reaction energy (eqn (2)), and EDx,r is used to calculate the corresponding baseline value (eqn (3)). Structural assessments were performed on reactants and intermediates optimized at both DFT and xTB levels of theory, to validate that the coordination numbers of Rh ion are 4 and 6 in reactant and intermediate states, respectively. In addition, the phosphine ligands of both reactant and intermediate states of the molecules are ensured to be in cis-position.
The artificial neural network (ANN) implemented within the NeuralFastAI54 framework, was employed to train the two types of Δ-ML models on the features after dimensionality reduction. The ANN models utilize the rectified linear functions (ReLU) to capture potential nonlinear relationships between the features on the target values. For optimization, the Adam optimizer,55 a method for efficient stochastic optimization, was employed to fit the parameters of the model. Hyperparameters, including the network architectures, embedding layers dropout rate, linear layers dropout rate, and number of epochs were optimized efficiently using the Bayesian optimization56,57 in an automated fashion, as implemented in the AutoGluon package (1.0).58 Other classical ML models, such as XGBoost59 and CatBoost,60 are also available in the AutoGluon package and can be implemented conveniently in an automated fashion as well.
Furthermore, the permutation importance61,62 of the features was assessed on the best-performing model, in order to gain insights into the relationships between the descriptors and the deviation of the baseline values from the target values.
To account for the difference arising from the different calculation levels, -models were trained using AIM-ACs and ACs of varying depths as molecular descriptors. The dataset of 1743 reactant–intermediate pairs were split into training, validation and test sets with a ratio of 64:16:20 in two steps: first, 20% of the total data were randomly selected to form the test set for all models in our study. The remaining data were randomly split into training and validation set with a ratio of 8:2 (64% and 16% of all data, respectively) during the automated training procedure as implemented in the Autogluon package.58 The performance of the -models is depicted in Fig. 4. Having comprised the electron information computed at both DFT and xTB levels of theory, all models trained on AIM-ACs with varying maximum depths dmax exhibited notably good performance, with a RMSE of approximately 11.0 kJ mol−1 and a R2 exceeding 0.90. When increasing the maximum feature depth from 1 to 3, the RMSE decreased from 11.1 kJ mol−1 to 10.6 kJ mol−1, and the R2 increased from 0.905 to 0.913. The optimal performing model, achieved with dmax = 3, had 4 hidden layers with 1042, 367, 150, 121 neurons, respectively, along with an embedding layer dropout rate of 0.651, linear layer dropout rate of 0.002 and 24 epochs. The hyperparameter sets and other detailed results of the optimal models trained on different dmax are summarized in Table S2.† Besides the ANN models, the identical training set was also fed to other models, such as XGBoost59 and CatBoost.60 The results for these models are summarized in Table S3.† The performance of the ANN model is slightly better than these classical ML models. Therefore, we mainly focus on the discussion of the performance of ANN models throughout this study. Further increasing the maximum feature depth did not improve the performance of the -models. This suggests that the discrepancy between and can be attributed to the difference in description of electron information in the short range (d = 0, 1, 2, 3). Fig. 4c provides insight into the fraction of AIM-ACs with different dmax after dimensionality reduction. As dmax increased from 3 to 9, the fraction of AIM-ACs of d = 0, 1, 2, 3 remained consistently over 0.72. This observation supports the conclusion that including AIM-AC features with larger dmax (dmax > 3) does not provide additional information for addressing the difference between and .
Although -models trained on AIM-ACs have exhibited excellent performance, it is important to note that the calculations for AIM-based atomic properties are computationally demanding, rendering AIM-ACs less practical in accelerating high-throughput screening for the desired catalysts. In contrast, the original ACs based on the isolated atomic properties demand significantly lower computational costs.
In our study, a substantial number of isolated atomic properties derived from DFT calculations with a ωB97XD functional48 and the parameters used in xTB calculations were used to construct the ACs (see ESI†). However, in comparison to AIM-ACs, the -models trained on ACs generally had poorer performance in predicting the . The best performing model was obtained at dmax = 7, with a RMSE of 14.2 kJ mol−1 and a R2 of 0.845. Fig. 4d illustrates the importance of including the atomic properties with larger dmax. As dmax increased from 1 to 7, the peak of the fraction also moves towards a higher d values, with d = 5 and d = 6 contributing the largest fraction (0.23 and 0.24, respectively) to the entire feature vector. On the contrary, ACs of d = 0 and d = 1 were less important and discarded via the dimensionality reduction procedure as dmax increased.
The details of the best performing models trained on AIM-ACs (dmax = 3) and ACs (dmax = 7), along with the corresponding feature importance analysis, are presented in Fig. 5. Regarding the model trained on AIM-ACs with dmax = 3, data points in the test and training sets were distributed evenly around the y = x line in the parity plot. Compared to the prediction performance using the pure baseline value (as shown in Fig. 3a), this -model displays a notable enhancement, with the RMSE decreasing from 23.2 to 10.6 kJ mol−1. Feature importance analysis indicates that the a-AIM-ACs contributed the most to the improvement, among a total of 364 features. On the contrary, the l-AIM-ACs, which are designed to describe the ligands, play less crucial roles in improving the model. Although this optimal model was trained on the AIM-ACs with dmax = 3, features of d = 0 and 1 exhibited the strongest predictive power, with AIM-ACs of d = 2, 3 contributing only minor corrections to the model.
Compared to AIM-ACs, -model trained on ACs exhibited lower accuracy, and notable deviating points were observed in both test and training sets. This disparity could be attributed to the inherent inflexibility of ACs compared to AIM-ACs, as the distribution of ACs in is much sparser than that of AIM-ACs. Nevertheless, when compared to the pure baseline model, the performance of the -model trained on ACs has substantially improved at a considerably lower cost. In contrast to the model trained on AIM-ACs, the ACs of d = 0 or d = 1 barely provided distinctive information for different complexes. Consequently, ACs of greater depth (d = 2 ∼ 6) have higher predictive power in addition to the baseline value . It is noteworthy that the features with high importance are dependent on the xTB parameters, particularly the anisotropic XC scaling parameter fXCθA, which is parametrized to describes the quadrupole expansion of the electron density of a specific element and account for the anisotropic exchange-correlation effect.34
In addition to the topology-based features such as AIM-ACs and ACs, the -model was also trained on the position-based feature set, such as SOAPs, which can be computed very efficiently in the describe43,44 package. Compared to the model trained on ACs, the model trained on SOAPs exhibited superior prediction performance, with a RMSE of 13.0 kJ mol−1 and a R2 of 0.871. The optimal set of hyperparameters is summarized in Table S2.† In contrast to the model trained on AIM-ACs features, where the xTB charge of the Cl ion in the intermediate (chg0xTB(Cl)) emerged as the most important feature, the atomic environments around Rh (p(Rh)) played the most significant role among a total of 463 features, in enhancing the predictive performance of the model trained on SOAPs.
In summary, the -models trained on AIM-ACs exhibit strong predictive performance. However, it is important to note that the calculation of the AIM properties demands high computational power (roughly equivalent to a single point calculation), due to the heavy numerical integration of the electron density. Nevertheless, by comprising the electron information into the feature sets, these models point out that the charge in the Cl basin calculated at the xTB level of theory (chg0xTB(Cl)) is the most relevant feature accounting for the difference between and . Alternatively, ACs and SOAPs provide a much faster route to predict with relatively high accuracies. In the case of ACs, features with higher d value have stronger predictive power while in the case of SOAPs, the atomic environments of Rh play significant roles in enhancing the performance of the model.
The ΔDx-models trained on ACs generally exhibit lower accuracies than the models trained on AIM-ACs. Unexpectedly, the best performing model was obtained at dmax = 1, with a RMSE of 12.7 kJ mol−1 and a R2 of 0.723. Fig. 6d illustrates that the fraction of ACs of d = 0, 1 dropped drastically as dmax increases. However, in contrast to the prediction, including the AC features with d > 1 does not improve the accuracy of the model.
Details and the feature importance analysis of the ΔDx-model trained on AIM-ACs with dmax = 9 are shown in Fig. 7. Data points in the test and training sets were distributed uniformly around the y = x line in the parity plot, with only a few notably deviating points present in the high energy region where ΔEDFT exceeds 100 kJ mol−1. The corresponding feature importance analysis reveals that a-AIM-ACs are the most significant contributors to the high predictive performance of the model. Among these features, the electron information around Cl and Rh ions are of the highest significance among a total of 507 features, particularly the AIM-charge of Cl in the intermediate calculated at xTB level of theory (chg0xTB,i(Cl)). As discussed in Section 3.1, this feature primarily accounts for the difference arising from the different calculation levels of theory. In contrast to the -model, where electronic information from xTB calculations plays a more prominent role, DFT information from the reactant is more important in this model, in addition to the already known importance of chg0xTB,i(Cl). Note that the changes in geometries in both reactants and intermediates can introduce errors independently to the prediction in the energy difference. The structural change of the reactants can be captured by the AIM analysis on the DFT-equilibrium structure and xTB-equilibrium structure of the reactant, namely, on RDFT,r and RxTB,r, respectively. However, the change in geometries of the intermediates remains hard to describe solely through the AIM analysis of RxTB,i. Therefore, higher importance of the AIM-ACs which depend on the DFT calculations is observed. On the contrary, it is less practical to discuss the results for ΔDx-models trained on ACs due to their low prediction accuracy. Nevertheless, the parity plot and the corresponding feature importance ranking are shown in Fig. S6.†
Furthermore, the ΔDx-model was trained using SOAPs that depend on the structure of reactants optimized at both xTB and DFT levels of theory, as well as the structure of the intermediate optimized at the xTB level of theory. It is noteworthy that SOAPs can be calculated much more efficiently compared to AIM-ACs, and the performance of the model trained on these features even slightly exceeds that of the model trained on AIM-ACs, with a RMSE of 10.3 kJ mol−1 and a R2 of 0.819. Given that SOAPs are less accurate in predicting the difference in different levels of theory than the AIM-ACs, as demonstrated in the comparative study on the -models, this result suggests that SOAPs may have stronger predictive power regarding the change in geometries. Among all the SOAP features, the atomic environment of the Rh ions in DFT-equilibrium structures of the reactants and in xTB-equilibrium structures of the intermediates is the most significant factor for the accuracy of this ΔDx-model.
The SOAPs can be further simplified by excluding the information of the DFT structure, although it plays an important role according to the results of the feature ranking. In this manner, the energy difference of the reaction obtained at the DFT level of theory ΔEDFT can be predicted exclusively using the information from xTB calculations, which eliminates the need for computationally expensive structural optimizations at the DFT level of theory. This simplified feature set is denoted as xTB-SOAPs. The DFT calculation for a molecule containing 72 atoms typically requires more than two days of CPU time, while a xTB optimization for the same molecule requires less than one minute on the same machine. In addition, the generation of the SOAP features requires less than one second. In general, xTB-SOAPs can be generated very efficiently. Note that xTB-SOAPs only gain efficiency in predicting ΔEDFT but not . This is due to the fact that xTB-SOAPs and SOAPs are identical for the prediction of , since the energies of the complexes at xTB and DFT levels were evaluated for identical geometries. As expected, the performance of the ΔDx-model slightly deteriorated after excluding the information from DFT-equilibrium structures, with the RMSE increasing to 11.1 kJ mol−1, and the R2 reducing to 0.789. However, compared to the prediction from pure baseline values ΔExTB, this is already a considerable improvement with only a minor increase in computational cost. Although this result is less accurate than the predictive study of activation energies of different types of elementary reactions,23,24,28 where the errors range from 2.0 to 8.1 kJ mol−1, our study allows the prediction of energy differences between the reactant and intermediate. The structural difference between the reactant and the intermediate is greater than the difference between the reactant and the TS, because they are further apart on the 3N-6 potential energy surface. This could be the reason for the difficulties in predicting the reaction energy.
Without the information of the atomic environments of Rh obtained in the DFT-equilibrium structure, the feature ranking shows an increased importance of the environment around P1 in the reactant, which emphasizes information from the ligands. A large size and high bulkiness of a ligand usually implies a large deviation between DFT-equilibrium and xTB-equilibrium structures. Therefore, the importance of a more detailed description on the ligand structure is heightened, when the direct descriptions on the change in geometries, such as SOAPs evaluated on DFT-equilibrium structure, are absent.
In summary, the ΔDx-model trained on AIM-ACs and SOAPs exhibits good performance in predicting the energy difference ΔEDFT. However, these two feature sets may account for different aspects of the difference between the baseline value and the target value: AIM-ACs are more related to the difference in levels of theory, while SOAPs are more associated to the change in geometry. Furthermore, xTB-SOAPs features are highly recommended for efficient prediction of the energy difference, owing to the low computational costs for xTB optimization and SOAPs calculation.
Fig. 8 Relationship between the reaction energy (ΔEDFT) and activation energy ΔE‡DFT of C–H activation. Compared to the 16 structures from our previous study,30 the ten newly selected structures possess lower reaction energies as well as lower activation energies. |
Regarding the structures of the ligands of the ten selected complexes, as can be seen from Table 1, the linker unit connecting the two coordinating phosphorous atoms varies and includes flexible alkyl chains but also more rigid aromatic or 2,3-O-isopropylidene-2,3-dihydroxy-1,4-bisbutyl structures. All structures have in common that the phosphorus atoms bear at least one aromatic substituent. A more in-depth analysis of geometric and electronic parameters based on the DFT-optimized [Rh(PLP)(CO)(Cl)] equilibrium structures reveals that the newly suggested ligand structures have larger buried B5 Sterimol parameters,63 as well as an on average lower dipole moment. In total 20 geometric and electronic descriptors exhibit significant differences between the original and the newly proposed bisphosphine set (see ESI†). Importantly, steric factors describing the accessibility of the Rh center such as the % buried volume64 are included in this descriptor set. The complexes with a lower predicted reaction energy also have a lower % buried volume, pointing to the fact that the Rh center is more accessible for the substrate. Overall, the dependence of the C–H activation reaction energy on the complex structure cannot be explained with a single factor, instead, multiple geometric and electronic parameters influence the C–H activation process. This underlines the complexity when searching for new C–H activation catalysts.
When looking at the synthetic accessibility of the suggest structures it can be concluded that although the ligands look complex at first, the structures can be realistically synthezised as both the linker as well as the substituted phosphine part can be independently synthezised. E.g. the 2,3-O-isopropylidene-2,3-dihydroxy-1,4-bisbutyl linker found in L1, L5, L6 and L9 can be synthezised following the procedure from Kagan and Dang,65 while the unsymmetric phosphines can be synthesized following Singh and Nicholas.66 Thus, the structures suggested in Table 1 can be synthezised and subsequently their performance can be experimentally validated.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4dd00037d |
This journal is © The Royal Society of Chemistry 2024 |