Towards predictive computational catalysis – a case study of olefin metathesis with Mo imido alkylidene N-heterocyclic carbene catalysts Free

Olefin metathesis is a cornerstone in chemistry for C–C bond formation under simultaneous redistribution of alkyl fragments.^1–5  It is used in the chemical industry, for example in the Shell Higher Olefin Process^6,7  and in the manufacturing of novel polymer-based materials.^8–10  More recent applications are in drug discovery,^11–14  total synthesis,^15–18  chemical biology,¹⁹  and biomass conversion.^20,21 

In olefin metathesis, one distinguishes between cross metathesis, that is the redistribution of alkylidene fragments, ring closing metathesis, where C–C bond formation results in a ring, and Ring Opening Metathesis Polymerization (ROMP), where a polymer is generated (see Scheme 1). Herrison and Chauvin were the first ones to propose a mechanism in 1971.²²  They proclaimed that the key species is a metal alkylidene compound that reacts with an olefin to form a metallacyclobutane intermediate. The latter undergoes ring opening to form the new olefin under regeneration of the metal alkylidene catalyst (compare Scheme 1).

Olefin metathesis would not have reached its importance without the development of highly efficient catalysts that allow for carbon–carbon bond formation under mild conditions. Despite the progress in recent years catalyst development is still ongoing. The ultimate goal is to develop a catalyst that is highly selective and highly active, but also very stable as characterized by high turnover numbers (TON) and turnover frequencies (TOF). In recent years, the experimental efforts have more and more been accompanied by theoretical investigations to understand the origin of activity and selectivity.^23–26 

Today, two main classes of homogeneous catalysts for olefin metathesis exist: late transition-metal catalysts in low oxidation states, typically with ruthenium as the active centre, also termed Grubbs type,²  and early transition-metal catalysts in high oxidation states, with molybdenum or tungsten as the active centre, also termed Schrock type.^5,27–29 

Examples of prototypical catalysts are shown in Scheme 2 (upper panel); Grubbs type catalysts are mostly 5-fold coordinated with tricyclohexylphosine (PCy₃) and/or N-heterocyclic carbene (NHC) ligands, whereas Schrock type catalysts have imido and alkoxy substituents. While Mo-based catalysts are more reactive, Ru-containing catalysts of the Grubbs or Grubbs–Hoyveda type are more stable, even when exposed to air and moisture³⁰  and have been more widely applied.³¹  However, a major breakthrough was achieved when Buchmeiser et al. reported the Mo imido alkylidene NHC class of catalysts with a 5-fold coordinated Mo-centre, an imido and a NHC ligand in addition to two alkoxy or triflate groups and the alkylidene (see Scheme 2 (lower panel)).³²  Variants can be cationic 4-fold coordinated or cationic 5-fold solvent coordinated species (Scheme 2).^5,33,34  These catalysts show an unprecedented functional group tolerance in the substrate, polymerizing monomers with hydroxyl, carboxyl, aldehyde, ether, and amine functionalities in addition to high activity.^32,35,36  The most active catalyst modifications reported by the Buchmeiser group have TONs up to 1 200 000,^35–38  while variations in the ligand sphere, substrate or solvent substantially affect the reactivity and selectivity. To decipher the origins of these different reactivities quantum chemical studies are of great help and also provide directions for catalyst development.^23,24,26,39 

Indeed, over the past decades, quantum chemical studies have advanced the understanding of homogeneous transition-metal catalysts by providing atomistic insight into elusive intermediates and by predicting the thermodynamics and the kinetics of the reaction.^40–42  The final goal in computational chemistry is to calculate the time evolution of the reaction under ambient conditions, that is, at finite temperature, pressure, concentration, and in explicit solvent with all reaction species – including counter ions – present. However, such ab initio molecular dynamics (MD) simulations are not feasible (yet) for large transition-metal catalysts in the condensed phase.

Today, full-fledged catalysts can routinely be described by dispersion-corrected Density Functional Theory (DFT) in implicit solvent – providing a static picture of the reactivity. As transition-metal complexes are more challenging for DFT, the accuracy of the employed density functional needs to be evaluated on a case-by-case basis.^43–45  In addition, other sources of error may stem from (i) a neglect of the conformational diversity when the most stable conformer of each species is not identified,^46,47  (ii) a neglect of counter ions that form ion pairs with the catalysts,³³  and (iii) an approximate description of solvation, where explicit solute–solvent interaction is not accounted for. These sources of error are often overlooked in state-of-the-art investigations but can significantly impact the accuracy of the calculated results. In order to evolve from being descriptive to being predictive, quantum chemical investigations need to aim for the best possible description of a chemical problem. This means not only an electronic structure method has to be invoked – in practise this means choosing an adequate density functional and basis set combination – but also conformational diversity and interaction with the environment (counter ions, explicit solvent molecules) have to be accounted for. In this chapter, we shortly review the theoretical background of quantum chemistry including recent developments such as improved semi-empirical methods and local coupled cluster methodologies for Mo-based olefin metathesis catalysts. Next, we describe concepts and challenges to include conformational diversity, counter ions, and explicit solvent interactions in mechanistic investigations and to derive chemical descriptors to predict reactivities. Recent advances of these aspects are illustrated at the example of olefin metathesis with Mo imido alkylidene NHC catalysts.

To describe the making and breaking of bonds an explicit treatment of electrons is required. Consequently, the Schrödinger equation must be solved, which describes – in the Born–Oppenheimer form – the interaction between moving electrons and static nuclei. Finding an approximate solution to this equation is the ultimate goal in electronic structure theory. Especially transition-metal clusters are best described with wave function methods because of their intricate electronic structure. But unfavourable scaling of such methods limits their applications. A notable exception is the domain based local pair-natural orbital singles- and doubles coupled cluster (DLPNO-CCSD(T)) approach developed be Neese and co-workers, that allows single point calculations of large transition-metal complexes.^48,49 

In practice, we are often limited to DFT, which offers more favourable scaling. However, improvement of the accuracy is less straightforward because of the approximate nature of the exchange–correlation functional. While the quest for an accurate density functional is still an active field of research (for selected literature reviews see e.g. ref. 50–53), some milestones were achieved over the last years. For example, the development of hybrid functionals with local range separation that target a balance between a reduction of the self-interaction error and an incorporation of static correlation.^54–57  While first results are promising, their performance for transition-metal complexes has yet to be evaluated.

Most important for the modelling of transition-metal catalysts is the inclusion of dispersion interactions, either through empirical correction terms that are added a posteriori to the energy or through new density functionals that are (re)-parametrized to account for these types of interactions.^58,59  Furthermore, the class of double-hybrid density functionals has emerged such as the B2PLYP,⁶⁰  where in addition to a part of the exchange being described by Hartree–Fock, a part of the correlation is modelled by a virtual orbital-dependent term derived from Møller–Plesset perturbation theory. As an alternative approach density functionals with a large number of parameters, that can be adjusted to experimental or high-level computational data, have been developed that also account for dispersion interactions.^61–63 

Despite extensive benchmarking of these new density functionals, the transferability of results still poses a challenge because the accuracy is system dependent. Consequently, a case-by-case evaluation of their performance is still necessary. While geometries are fairly robust and show moderate dependency on the employed density functional (and basis set), final energies are more critical. To obtain reliable values not only a modern density functional is required but also a large enough basis set must be employed (see e.g. ref. 64).

With the rise of high throughput computational screening protocols and automated reaction mechanism exploration tools, the need for cheap yet accurate electronic structure methods became apparent. Semi empirical quantum mechanical methods make use of a number of integral approximations and minimal basis sets to speed up calculation time by orders of magnitude but at the cost of the accuracy, especially for transition-metal complexes.^65–68  Tight binding DFT (DFTB) offers an attractive compromise between accuracy and cost.^69–71  While the DFTB3 implementation uses pair specific parametrization schemes that ultimately limit the applicability,⁷²  GFN2-xTB circumvents this by using global and element-specific parameters.⁷³  It does not rely on “classical” force-field type correction terms and allows for “on the fly” calculation of parameter sets.⁷³  By including anisotropic second order density fluctuations and by providing parameters for the majority of the periodic table, GFN2-xTB has emerged as the semi empirical quantum mechanics method to be used in high throughput calculations.⁷⁴  Still, when employed to transition-metal catalysts results have to be carefully evaluated or require further treatment as GFN2-xTB and DFT energies may not coincide (vide infra).

Even less expensive are classical force fields, where the electronic structure is no longer explicitly described but only considered implicitly by definition of atom types and bonding patterns. As such, they cannot be used for describing chemical reactions but may be revisited when exploration of the potential energy surface (PES) is concerned.

If the dynamics of a catalyst and its time-evolution is of interest, (ab initio) molecular dynamics (MD) is the method of choice. In the Born–Oppenheimer approximation the nuclei are propagated in space and time.⁷⁵  In the absence of nuclear quantum effects this propagation is modulated by Newton's law of motion: at each step the electronic structure is converged in a self-consistent way and the forces are obtained from the negative gradient of the energy with respect to the nuclear coordinates. Even if nuclei are treated classically and orbitals are propagated without reaching self-consistency as done in the Car–Parrinello approach,^76,77  molecular dynamics are still very expensive. It is not yet feasible to simulate a condensed phase system, for example, a large transition-metal catalyst in solvent with all reactants present, up to the μ- or millisecond time scale. Yet, molecular dynamics approaches have gained interest for exploration of the PES and for generation of low energy conformers.⁷⁸  To that end, semi empirical quantum mechanical methods or classical force fields are used in implicit solvent or fully condensed phase calculations (vide infra).

However, even in classical MD, the timescales reached may not be sufficient to sample the regions of interest of the PES, and a number of enhanced sampling techniques is available.⁷⁹  They have in common that they drive the system away from its current energy minimum allowing exploration of larger regions of the PES. One such technique is accelerated MD⁸⁰  (aMD) with a speed-up up to a factor of 2000.⁸¹  In aMD, certain degrees of freedom (e.g., dihedral angles) are boosted, when being below a pre-defined energy threshold. As a consequence, the energy wells become shallower and transitions barriers are (artificially) lowered, thus, facilitating transitions between various conformer states. The advantage is that this approach requires little knowledge about the system itself but only a definition of the threshold energy, the degrees of freedom to be boosted (in most cases the dihedral angles), and the boosting strength. The disadvantage is that the original PES is modified and needs to be re-weighted to regain the original surface.⁸²  However, this is not such a problem if only energy minima are of interest. Another technique is the so-called metadynamics approach,⁸³  where a collective variable is defined before starting the simulation. Once a certain conformation has been visited a bias potential is added to the collective variable. Thereby, the potential well is successively filled with these Gaussian shaped bias potential functions – pushing the system away from its local minimum. The difficulty here is to define a collective variable suited to describe the system efficiently. However, regaining the original PES and calculating the resulting free energy is more straightforward.⁸⁴ 

Conformational diversity is often neglected in computational studies, where a mechanism is proposed based on one conformer only. This might be a sufficient approach for small and rigid catalysts or when the coordination environment of the transition metal undergoes little changes during the reaction and remains similar to its X-ray crystal structure. However, already in 2011, Besora et al. pointed out the importance of conformer searches when modelling transition-metal reaction mechanisms.⁴⁷  For a medium sized flexible palladium catalyst for Suzuki–Miyaura cross-coupling they found that neglecting the conformational diversity yielded errors up to 40 kJ mol⁻¹. In addition, a reaction may not necessarily progress from the lowest energy conformer. Instead, it may progress from a conformer that has a slightly higher energy than the overall minimum but that has a lower overall reaction barrier.⁴⁶  Therefore, multiple low lying conformers need to be tested for further reactions.⁸⁵  To this end, single conformer reaction studies provide an incomplete picture and potentially wrong energetics. Yet very few studies address this source of inaccuracies.

To implement conformational searches different strategies can be followed. Besora and co-workers used Monte Carlo Multiple Minimum searches as implemented in MacroModel, which uses a classical force field (MM3), to determine the low energy conformers of their transition-metal catalysts in the gas phase. Initial molecular mechanics conformers were then re-optimized with full DFT.⁴⁷ 

Another strategy was reported by Lepori et al. who used classical MD simulations to determine low energy conformations of a beta-diketiminatocobalt(ii) alkyl complex.⁸⁶  In their study, the authors used the Metal Center Parameter Builder MCPB.py tool⁸⁷  in combination with a method by Seminario⁸⁸  for intramolecular forces to generate force field parameters for the metal complex. They ran short MD simulations at 1000 K under solvent free conditions to generate an initial ensemble of structures.⁸⁶  To reduce the number of conformers to a few representative ones, they employed a K-means clustering algorithm and re-optimized the resulting structures with DFT. While this methodology was successful, there are several drawbacks: first of all, metal force field parameters need to be provided and even if they are available the description of the metal is very simplistic and does not account for changes in the coordination geometry. However, such changes occur in many 5-fold coordinated transition-metal complexes, which adapt their coordination depending on the ligands to be either square-pyramidal or to be trigonal bipyramidal as well as anything in between.⁸⁹  Second, the classical MD runs were very short which was compensated by letting them run at very high temperatures. However, most force fields are not designed to simulate at such elevated temperatures; their performance may not have been tested and they may not be reliable. Third, solvent effects were neglected, which again, is a simplification. Nevertheless, such a classical MD approach can be useful if solvent molecules should be considered and conformations of a complex in explicit solution are of interest. One such example are complexes, where strong solute–solvent interactions are required to retain their shape.^90,91 

A breakthrough for conformer searches of (transition-) metal complexes was the development of the semi empirical quantum chemical method GFN2-xTB. It outperforms other semi empirical approaches and does not require an a priori parametrization of the system.^73,74  This method was the basis for the Conformer Rotamer Ensemble Tool (CREST) developed by the Grimme group that automates the exploration of the PES and hence, the conformer search.⁷⁸  In this approach, short metadynamics type simulations are performed with GFN2-xTB, where the cartesian coordinates serves as collective variable. Each simulation is analysed for its unique minimum structures that are collected and saved. Multiple runs are necessary to avoid being trapped in a local minimum. The tool is fully automated and the user is provided with a list of unique structures that fall within a certain energy range and differ by a predefined root mean square deviation (RMSD) and energy criterion. If two conformers differ by more than this RMSD they are saved as different entries. Depending on the settings the user is supplied with several hundred conformers. As the performance of GFN2-xTB varies, some systems may require calculation of single point energies with DFT on the GFN2-xTB structures,⁹²  while others may require a full re-optimisation of the GFN2-xTB structures with DFT. Thereby, it is often necessary to reduce the number of structures to a few representatives, for example, in a hierarchical cluster algorithm, and only consider those for further processing.^93,94 

It is well-known in transition-metal chemistry that counter ions are not mere spectators but that they can impact catalysis. Prominent examples are gold catalysts^95,96  but also palladium catalysts,⁹⁷  ruthenium catalysts for alkyne to vinylidene isomerization,⁹⁸  and a number of catalysts for olefin metathesis⁹⁹  and olefin polymerization¹⁰⁰  were reported to be influenced by the counter ions. In the case of the Mo Imido Alkylidene NHC complexes, dissociated triflate anions were detected by ¹⁹F NMR spectroscopy to be in close vicinity to the active catalyst.³³  Despite these experimental evidences very few theoretical studies attempt to include counter ions in their chemical models.^46,101  Yet, counter ions were found in numerous instances to be crucial to explain reactivities and their neglect would lead to incorrect reaction profiles.⁹⁷ 

Of course, inclusion of counter ions prompts the questions of (i) where to place them (ii) how to model them – as naked counter ions or as solvent coordinated species – and (iii) what level of theory to use to describe them. The answer depends on the chemical question to be studied and whether it is the impact of the counter ions on the mechanism that is of interest or the ion dynamics. While most ions are treated as “naked” species, alkali cations such as Na⁺ or Li⁺ are modelled as solvent coordinated complexes.^102,103 

Among the first to include counter ions in computational studies were Ziegler and co-workers. They reported a series of publications using a quantum mechanics/molecular mechanics (QM/MM) approach to study a reaction mechanism, where the counter ion is treated at the MM level.^104–106  Later studies resorted to a full quantum chemical description to evaluate the impact of counter ions on the reaction mechanism.⁹⁸  Matsumoto and co-workers combined MD simulations to sample various conformations before selecting representative ion pair structures for mechanistic quantum chemical investigations.^107,108  The sampling step was crucial to identify the most relevant states and to calculate a reliable reaction pathway. It must be accounted for in any static evaluation of counter ion effects. Other studies aimed at assessing the dynamic behaviour of the counter ion or ion pair. For example, in 2006, Corea and Cavallo employed classical MD simulations to study the dynamics of metallocenium ion pairs in explicit benzene,¹⁰⁹  whereas a 2013 study used QM/MM/MD to investigate the dynamics of the agostic interactions in Pt(ii) complexes.¹¹⁰ 

Homogeneous transition-metal catalysis is almost exclusively performed in solution and the catalytic activity is often linked to the solvent the reaction is carried out in. Hence, any realistic description of such species must include solvent effects.¹¹¹  In computational chemistry, several descriptions of solvents have emerged. The most simple ones are implicit solvation models, where the solute is placed in a molecular shaped cavity that is surrounded by a dielectric continuum that mimics the solvent.¹¹²  These models only account for bulk solvent effects but they are fast and avoid sampling and averaging over solvent conformations. Several implementations have been reported and they mainly differ by how the solute–continuum interactions are calculated and how the cavity is defined:¹¹²  There is the Polarizable Continuum Model (PCM) by Tomasi and co-workers,^113,114  the Conductor Like Screening Model (COSMO) and its extension for “real solvents” COSMO-RS by Klamt,^115,116  and the SMD method by Marenich, Truhlar and Cramer.¹¹⁷  The SMD and the COSMO-RS are the only two models that also contain “non-electrostatic” terms for calculating the interation between solute and solvent.

While implicit solvation models have sucessfully been applied in transition-metal catalysis including olefin metathesis,^33,93  they fall short if explicit solute–solvent interactions are expected to play a role.¹¹⁸  The alternative is a fully condensed phase description that includes many explicit solvent molecules in the calculations and that also captures the solvent dynamics. However, a fully ab inito MD treatment of solute and solvent is limited to small systems only. Even if only the catalyst is described by quantum mechanics and the solvent is treated by molecular mechanics, such treatment becomes computationally prohibitive rather quickly.¹¹⁹ 

A low-cost alternative is the so-called microsolvation approach, where only a “few” relevant explicit solvent molecules are considered in a quantum chemical calculation. In principle, two strategies exist to set-up microsolvated structures: in a bottom-up approach one starts from the solute alone and successively adds solvent molecules,^120–122  whereas in a top-down approach one starts from a fully condensed phase calculation, for example a classical MD simulation, and selects a few (strongly interacting) solvent molecules.¹²³  In any case, the questions remain, where to position the solvent molecules, how to orient them, and how many to include. Our group has recently developed a top-down multiscale approach to quantum chemical microsolvation, where the solvent molecules are automatically placed based on the free energy of solvation.¹²⁴  We performed classical MD simulations in explicit solvent and the obtained trajectories were analysed with Grid Inhomogeneous Solvation Theory (GIST). GIST allowed us to calculate the solute–solvent interaction on a grid and to identify favourable solute–solvent interaction sites. Subsequently, solvent molecules were placed at positions, where the averaged trajectory showed the highest solvent density. After placing the first explicit solvent molecule, the density of that solvent molecule is subtracted from the total (solvent) density until all density is assigned. As a last step, the solvent molecules are ranked according to their interaction strength so that the thermodynamically most favoured molecules can be identified and can be used to set-up the microsolvated cluster structure.

Despite our and other methods to set up microsolvated clusters, microsolvation has mostly been applied to organic complexes, while their performance for transition-metal catalysts has yet to be evaluated.^121,122,124 

To determine the reaction mechanism of a transition-metal catalyst, educts, products, intermediates, and transition states must be identified in their relevant conformation, while solvent effects as well as potential interactions with counter ions must also be accounted for. Thereby, it is essential to test different reaction pathways and to check for deactivation routes to identify the most probable mechanism. Proposed reaction mechanisms should also be verified by testing various substrates¹²⁵  but this can be challenging when reaction routes are substrate specific (vide infra).

With the increase in computing power automated approaches have emerged, that autonomously explore entire reaction networks or elementary reaction steps.^126–130  For example, in the Chemoton program,¹²⁶  two species are aligned and let react with each other based on interaction sites in order to discover new reaction pathways with minimal human bias. While the applicability has been demonstrated by many examples from organic chemistry,^130–133  automated exploration of reactions mediated by transition-metal catalysts is still challenging but shows great promise.¹³⁴  However, the complexity, i.e., the high dimensionality of the often shallow potential energy surface, the multitude of potential interactions sites between catalyst and substrate, and the intricate electronic structures, render a fully automated exploration of the whole reaction network impracticable. A pragmatic protocol is presented by Wheeler and co-workers with their transition-state re-optimization toolkit AARON.¹³⁵  Given that a reaction pathway is already characterized, modified catalysts and substrates will be pre-aligned based on known transition-state structures and reoptimized.

As determination of the reaction mechanism is often a tedious process involving a lot of trial-and-error, any automation that minimises manual work is highly desirable.

Despite all difficulties, validation of results is essential and whenever applicable comparison with experiment should be aimed at. However, it is important to compare calculated and experimental data on equal footing. Quantum chemical electronic energies are typically calculated for one isolated molecule at 0 K in the gas phase, while experimentally obtained reaction or activation free energies are measured at T > 0 K for a finite concentration in solution. Thus, zero-point and thermal corrections have to be considered in the calculations. This can be achieved with a complex series of modifications of the standard rigid-rotator/harmonic oscillator model.^136,137  In addition, Boltzmann averaging over multiple low energy conformers may also be required to increase the accuracy. For the calculation of rate constants, concentration effects may have to be accounted for by microkinetic modelling.¹³⁸  For comparison of calculated free energies (of activation) to TOF or TON Kozuch and Shaik's energy span model can be invoked.¹³⁹  If experimental data is scarce, uncertainty quantification measures may aid to estimate the accuracy of the calculated electronic structure results.¹⁴⁰ 

To extract reactivity and selectivity trends in a series of catalysts or substrates, (quantum chemical) descriptors and chemical maps have been proven useful to formulate concepts and rationale for future catalyst design.^39,141  For example, quantum chemical ligand based descriptors such as the buried volume or the Tolman Electronic Parameters (TEP) have been used to characterize NHCs.^142,143  More involved is the calculation of chemical shift tensors that provide a direct link to the reactivity in olefin metathesis catalysts.^144–146  This approach has been used to characterize and distinguish metallacyclobutanes by their reactivity.¹⁴⁵  Chemical maps may reveal (underlying) features that are essential in catalyst optimization. Examples are the steric map of the surrounding of the metal centre to estimate the accessibility of the metal¹⁴¹  or the chemical space map that reveals how ligands influence the metal coordination geometry.¹⁴⁷  These efforts were complemented by the de novo design of Jensen and his group, who developed a number of tools for the in silico prediction of novel Ru based olefin metathesis catalysts.¹⁴⁸ 

Advances and challenges of computational modelling of homogeneous transition-metal catalysts are in the following being discussed at the example of two Mo imido alkylidene NHC catalysts, namely, Mo(N-2,6-Me₂-C₆H₃)(CHCMe₃)(I-^tBu)(OTf)₂ (1, I-^tBu = 1,3-di-t-butyl-1,3-dihydro-2H-imidazol-2-ylidene and Mo(N-2,6-Me₂-C₆H₃)(CHCMe₃)(IMes)(OTf)₂ (2, IMes = 1,3-dimesitylimidazol-2-ylidene). Both catalysts can occur in a syn and in an anti conformation with the syn conformer being more stable. For both compounds X-ray crystal structures (of the syn conformers) were reported.³² 

To address the structural diversity of Mo imido alkylidene NHC catalysts, we set out to generate conformer ensembles of 1 and 2 using CREST⁷⁸  in combination with GFN2-xTB.⁷³  To test the performance of GFN2-xTB, the X-ray crystal structures of 1 and 2 were optimized in implicit solvent and compared to full DFT optimisations (BP86/def2-SVP/dichloroethane). As both methods yielded similar structures for the two catalysts (see Fig. 1) with only small differences in the triflate orientations, the conformer generation with CREST was started. It was performed in implicit dichloromethane with an energy window of 30 kcal mol⁻¹, that is structures up to 30 kcal mol⁻¹ with respect to the most stable conformer are accepted, all others are discarded. To identify a new conformer, the energy threshold between two structures needed to be larger than 0.1 kcal mol⁻¹ and their RMSD must differ by more than 0.1 Å. These settings resulted in several hundred conformers.

The lowest energy GFN2-xTB conformer of each catalyst was reoptimized with DFT and compared to the optimized X-ray crystal structure. Surprisingly, results did not agree well. In the case of 1, GFN2-xTB found the optimized crystal structure to be 38.2 kJ mol⁻¹ higher in energy than the best CREST predicted conformer and the difference was even worse in the case of 2, where the X-ray crystal structure was 70.4 kJ mol⁻¹ higher in energy than the best CREST conformer. In stark contrast, re-optimization of these conformers with DFT (BP86/def2-SVP/dichloroethane) predicted the X-ray crystal structures to be the more favoured and the CREST conformer to be by 28.4 kJ mol⁻¹ (1) and by 16.8 kJ mol⁻¹ (2) higher in energy. As evident from the overlays displayed in Fig. 2, the major visible structural differences are the orientations of the two triflate groups.

The poor agreement between DFT and GFN2-xTB results prompted us to revisit the generated conformer ensemble. Rather than directly considering the predicted low energy conformers, we focused on obtaining conformers that are as structurally diverse as possible. Thereby, we did not only account for the fact that GFN2-xTB and DFT energies and structures differ but also take into consideration that the lowest energy conformer might not have the lowest energy barrier.

As it is not tractable to handle several hundred structures for mechanistic investigations, we reduced the total number of structures by grouping together those that are similar in their geometry using a so-called clustering algorithm. In the hierarchical average linkage clustering method, the two closest (most similar) structures are joined into one cluster, then the next closest clusters are grouped until all clusters are joined into one. The dendrogram of the hierarchical average linkage clustering is shown in Fig. 3 at the example of catalyst 1 in its syn conformation. On the y-axis the difference between two clusters is depicted calculated as an average distance between any pair of elements belonging to cluster A and B, respectively. By choosing a cut-off, here of 2.0 Å, 38 cluster structures were obtained, whereas the average distance between each two clusters is 2.0 Å or larger.

As visible in Fig. 4, where the cluster distribution and the resulting conformers of 1 in its syn and in its anti conformation are shown, the cut-off is critical for the number of resulting clusters. It is a trade-off between structural diversity and computational costs. But in any case, quite diverse structures were generated as evident from the overlay of the optimized cluster representatives in Fig. 4. Results for catalyst 2 were similar with 22 clusters for the syn conformer (cut-off 2.5 Å).

Resulting cluster representatives (centroids) were optimized with GFN2-xTB and various DFT methods to obtain energy minima and to compare structures and energies. While the results for both catalysts were very similar, the correlation between GFN2-xTB and DFT energies and structures is illustrated at the example of 2 in its syn conformation. The diagrams depicted in Fig. 5 impressively show that correlation between the GFN2-xTB optimized structures and the various DFT optimized structures is very poor (top panel Fig. 5). This finding indicates that GFN2-xTB energies (and structures) are not reliable for Mo imido alkylidene NHC catalysts. However, re-optimization even with a small basis set such as def2-SVP and without dispersion corrections resulted in good correlation with the fully optimized BP86/def2-TZVP/D3 structures (bottom left Fig. 5), and indeed, the correlation between the BP86/def2-TZVP/D3//BP86/def2-SVP single points and the fully optimized BP86/def2-TZVP/D3 structures is excellent (bottom right Fig. 5). A full optimization with def2-TZVP is therefore not necessary.

While the GFN2-xTB structures may not be accurate enough for Mo containing complexes, CREST can still be used to efficiently explore the PES and to generate an initial ensemble of structures. Clustering of the conformers ensures reduction of dimensionality and a structural diversity. While this protocol has successfully been applied by us (vide infra) and others,⁹⁴  further testing of various clustering algorithms is necessary.

In a recent study, our group investigated the reaction mechanism of 2 with methoxystyrene by DFT in implicit solvent (see Fig. 6).³³  In agreement with available ¹⁹F NMR data, our calculations predicted the reaction to proceed in an associative fashion, where the neutral pre-catalyst formed an unprecedented (neutral) olefin adduct 2₁ (Fig. 6) before dissociation of one of the triflates generated the catalytically active species 2₂ and the metathesis took place. The transition state towards the metallacyclobutane TS(2₂-2₃) was found to be rate-determining, with an overall barrier of ΔG = 115 kJ mol⁻¹. Ring opening and final product formation (2a) were found to have lower free energy barriers.

In order to obtain consistent free energies for the entire reaction pathway, it was necessary to explicitly include the dissociated counter ion. Calculations suggest that the negatively charged triflate and the cationic catalyst form an ion pair with a triflate to metal distance between 4 to 5 Å, which is in agreement with NMR data. The free energy of this ion pair was found to be more favourable than that of the separated ions indicating that the electrostatic attractions outperform the entropic penalty of the complex formation. It is noteworthy here that at the time of the study it was not possible to determine the most favourable position of the counter ion with respect to the catalyst by a comprehensive sampling of the phase space. Hence, we used chemical intuition to test several positions of the triflate and chose the one that yielded the lowest energy of the ion pair. Unexpectedly, the influence of the counter ion was rather consistent and the triflate lowered the relative electronic energy by around 50 kJ mol⁻¹.

To improve the sampling, one could use the previously described CREST in combination with GFN2-xTB and DFT re-optimization.

In another study, we investigated the E-selective ROMP of norbornene with 2 to yield the first insertion product.⁹³  Unlike methoxystyrene, the reaction did not proceed in an associative fashion but in a substrate induced dissociation of one triflate to yield the catalytically active species. As such, it is an example of a substrate specific mechanism. The reaction then proceeded via the known metallacyclobutane formation and reversion to form the first insertion product. Given that the catalyst can adopt a syn and an anti conformation and that the monomer can react with the catalyst in an ene_syn or ene_anti fashion, four different first insertion products can be formed.

The high flexibility of the alkylidene tail of the insertion product and the absence of an experimentally determined structure made it particularly challenging to find the lowest energy conformation based on chemical intuition alone. Indeed, in order to correctly predict the reaction to be exergonic, conformer searches of the first insertion products were essential. The multi-step computational protocol, using CREST to generate an ensemble as described previously, followed by clustering to reduce the number of structures, and re-optimization with full DFT, yielded conformers that were up to 30 kJ mol⁻¹ lower in energy than those predicted by chemical intuition.

While the standard reaction mechanism failed to explain the high E-selectivity, we detected a stereo-specific alternative reaction pathway. This route circumvented the formation of an olefin adduct but allowed direct addition of the olefin to yield the metallacyclobutane for one out of the four possible stereoisomers.⁹³  The direct addition was found to be by 30 kJ mol⁻¹ lower in free energy than the known mechanism – indicating once more the importance to test alternative pathways when evaluating the catalytic activity of transition-metal catalysts. Visualization of non-covalent interactions^149,150  further aided to rationalize the origin of E-selectivity of 2 by comparison of key intermediates to the non E-selective catalyst 1.⁹³ 

This chapter summarized recent advances and addressed current challenges in the computational modelling of homogeneous catalysts at the example of Mo imido alkylidene NHC complexes. An accurate description of the chemical system is required for calculations to become predictive. In this regard, significant advances have been made to include the conformational diversity and interactions with counter ions in quantum chemical protocols to improve the accuracy of the theoretical studies. However, the electronic structure of molybdenum still poses a challenge for fast semi empirical methods such as GFN2-xTB. Yet, semi-empirical methods are indispensable for screening of a large number of structures.

In addition to conformational flexibility and interaction with the environment, substrate or stereoisomer specific reaction pathways must be taken into consideration in automated high-throughput calculations. Only if they are accounted for, the reactivity of Mo imido alkylidene NHC catalysts can reliably be modelled and steric and/or electronic factors for rational design of superior catalyst can be deduced.

M. P. would like to thank the Austrian Science Fund (FWF) for financial support (P33528) and Sarah Schulz for performing several conformer studies described within this chapter.