Modelling Reactions at the Active Sites of Chiral Ruthenium Catalysts Using Density Functional Theory

Journal Archive

Platinum Metals Rev., 2007, 51, (2), 54
doi: 10.1595/147106707X180891

Modelling Reactions at the Active Sites of Chiral Ruthenium Catalysts Using Density Functional Theory

NEW APPROACH TO UNDERSTANDING OF CATALYTIC REACTIONS

S. A. French
Johnson Matthey Technology Centre,
Blounts Court, Sonning Common, Reading RG4 9NH, U.K.
Email: frencsa@matthey.com

Article Synopsis

Selectivity is a key success factor in the chiral catalyst technologies market. Understanding the fundamental processes that occur when a reagent interacts with a homogeneous single site catalyst, both in its approach and at the active site, is therefore critical to the rational design of new catalysts. Ruthenium-based asymmetric hydrogenation catalysts have been considered as part of a collaborative research project. [(S)-XylBINAP-RuH₂-(S,S)-DPEN], first developed by Noyori (1–3), is studied as the parent or prototype model of a series of efficient hydrogenation catalysts, among them the catalysts based on the P-Phos, PhanePhos and ParaPhos ligand families (4).

This work addresses homogeneous asymmetric hydrogenation of C=O bonds, from ketones to alcohols, which has applications in the industrial production of pharmaceutical intermediates. The catalyst studied, [(S)-XylBINAP-RuH₂-(S,S)-DPEN], is shown in Figure 1.

Fig. 1

The structure of [(S)-XylBINAP-RuH₂-(S,S)-DPEN], asymmetric hydrogenation catalyst

Each reaction mechanism may be understood in terms of the energies of intermediates and the roles of ligands and additives, as determined by state-of-the-art computational techniques. The knowledge gained will then be exploited for the design and synthesis of ligands for improved catalyst systems. Advances in experimental techniques will allow rapid identification of lead ‘designed’ catalysts by automated parallel screening. With this scheme in mind, a consortium was assembled from leading experts (both industrial and academic) in all areas of the workflow. The partners are the Royal Institution of Great Britain, the University of Liverpool, the University of Southampton, Johnson Matthey, AstraZeneca, GlaxoSmithKline and Pfizer. The aims of the project were to implement an evolutionary improvement in ligand/catalyst design strategies:

COMPUTATION ⇔ SYNTHESIS ⇔ ACCELERATED TESTING

This computation-guided approach for catalyst discovery is expected to be more efficient, faster delivering and more revealing on the molecular aspects of a catalytic cycle than one-at-a-time synthesis or combinatory methodologies, which usually screen catalysts at random (5); see Scheme I.

Scheme I

Roles of research consortium members

Since the project's conception there has been a step change in the ability of industry to perform high-throughput screening. This acceleration has enormously reduced the time required to identify the right catalyst for any desired transformation from a library of existing catalysts or ligands. The preparation of the library of ligands and catalysts remains, however, the bottleneck in this process. For example, the selection of the correct substituents at phosphorus in any family of bidentate phosphine ligands is largely a matter of trial and error, with each component of the family requiring independent and often time-consuming synthesis. The design and synthesis of new ligand backbones is even more time-consuming, and there is no certainty that the new ligands will be effective in the desired transformation. Understanding the factors which govern the relationship between the structure of the ligand and its efficacy in catalysis will accelerate the process of ligand design by evaluation through computer simulation (‘in silico’) of a large number of structural variations, among which only the most promising structures will actually be synthesised.

The project, supported financially by the U.K. Department of Trade and Industry's ‘Manufacturing Molecules Initiative’ (6), is focused on two industrially important organic reactions: (a) the production of chiral alcohols via the asymmetric reduction of ketones; and (b) C–C bond forming reactions such as the Heck reaction. Molecular modelling has so far focused on Reaction (a). Computer simulations have been used, at the molecular-mechanical (7) and, as reported here, quantum-mechanical levels, to investigate the structures of proposed catalysts, and to probe the reaction mechanism.

Initially it was proposed to use activation energies calculated by considering transition states (TS) between reactant and products to compare the performance of catalysts. However, difficulties in simulating the reactants and products correctly caused the TS calculations to fail. Therefore an alternative strategy was implemented; a geometric constraint was applied and the reactant brought towards the reactive centre, exploring the pathway of the ketone molecule to the active site. Understanding of the correct relative positions of reactants and products, and further understanding the need to ‘lock’ conformations of ligands, has led to the capability of performing TS calculations on ‘cut down’ (i.e. simplified) model catalysts, as well as exploring the entry of reactants to a real catalyst system.

Methods

The processes of prime interest to us involve the breaking and creation of bonds, which means that the electronic structure as well as the molecular structure must be modelled. Traditional methods in electronic structure theory, in particular Hartree-Fock theory and its descendants, are based on the complicated many-electron wavefunction. The main objective of density functional theory (DFT) is to replace the many-body electronic wavefunction as the basic quantity by the electronic density. Within this study we have used the DFT code DMol³ (8, 9) for both model TS calculations and constrained optimisations.

All the DFT investigations were performed using the linear combination of atomic orbitals approximation, with a double numerical basis set augmented by polarisation functions (with a 5.5 Å cut-off). The calculations employed the gradient-corrected Perdew-Becker-Ernzerhof (PBE) exchange-correlation functional. The fine accuracy convergence criteria were used throughout for both electronic structure and atomic optimisation calculations. The criteria guarantee that the energy per bond, bond lengths and angles converge to approximately 0.1 eV, 0.01 Å and 1°, respectively.

Constrained Optimisation Calculations

To understand how the reactant molecule approaches the metal centre and what restrictions are placed on its passage, we have performed a large number of simulations to compile a trajectory of the path followed, and to compare the energy barriers that are encountered. At each stage, the geometry of the system is optimised with respect to one constraint, namely the distance between the hydride H on the ruthenium and the C of the carbonyl group in the ketone. (These species eventually become bonded to one another in the alcohol product.) The output from one simulation is used to generate the initial configuration for the next. Initially we used the results to understand the interaction between reactant and catalyst, to identify the most relevant reactants and products for TS calculations, which are described below. The application of the method was then extended to consider four ‘quadrants of attack’ of the ketone to the active site, to probe the potential energy surfaces of the reaction. This will provide vital information concerning the selectivity of the catalysts.

Transition State Calculations

The TS calculations have required a workflow applied to a cut down version of the catalyst (Figure 2), to arrive at a final model for the TS. The stages are as follows:

(a) Relaxation of the reactant;

(b) Construction of the product from the reactant using ‘chemical intuition’ to move atoms around;

(d) Full relaxation of the product;

(e) Linear synchronous transit (LST) method to approximate reaction path and provide input to full LST/quadratic synchronous transit (QST) calculation with conjugate gradient (CG) optimisation;

(f) Single-point calculation of TS with frequency analysis;

(g) Animation of negative mode to check which centres are involved;

(h) TS optimisation calculation where mode is followed;

(i) Single-point calculation of TS with frequency analysis;

(j) Animation of negative mode.

Fig. 2

Cut down model of the Noyori catalyst

Results and Discussion

Among the models used to rationalise the structure-activity relationship in asymmetric homogeneous catalysis, the so-called ‘quadrant approach’ is one of the simplest and most effective.

The space around the reactive centre is divided into four volumes, across which the substrate can bind to the metal centre in a number of different conformations. The ligand will prevent access to some quadrants by simple steric interaction, thereby forcing the substrate to bind to the metal in a preferred conformation that, upon transfer of the hydride from the metal to the substrate, becomes the precursor to the favoured enantiomer of the product. Such a model, although very simplistic, allows straightforward rationalisation of the sense of stereoinduction obtainable with a number of well-known hydrogenation catalysts such as Ru-BINAP, Rh-DuPhos and Rh-BisP* (10, 11).

Starting from this simplistic approach (see Figure 3) we have developed a more sophisticated model that takes into account the whole trajectory of the substrate into the ‘reactive pocket’ of the catalyst. It is well known that subtle modifications of the substituents at phosphorus can produce very significant changes in the activity and selectivity of the catalyst (one example of this being the so-called ‘meta-effect’). We suggest that the reason for these effects may reside not only in changes at the transition state, but also in the docking of the substrate into the reactive pocket, well before the bond breaking/bond forming interactions are established.

Fig. 3

Quadrants of interest showing the steric interactions between the reactant and ligands of the catalyst

Constrained Optimisations: Initial Quadrant

Initially we have considered the catalyst [(S)-XylBINAP-RuH₂-(S,S)-DPEN] as an exemplar of the class of asymmetric materials that we are interested in understanding. As mentioned above, initially we focused on understanding the pathway of the reactant to the active site. With this in mind we chose to model the quadrant and orientation known to lead to the preferred product. Determining the conformational changes forced on the reactant or catalyst during approach would provide a better understanding of where reactants and products should be sited for TS calculations. The final structure from each of the constrained optimisation calculations, starting with a constraint (Ru–H - - - C=O) of 8 Å and reducing the separation between the reactant and catalyst in steps of 0.5 Å, is shown in an animation (12). It is clear from the animation and subsequent analysis of the potential energy surface for the pathway (Figure 4 and Figure 5) that we can begin to understand the complexity of these systems. There are two distinct energy barriers that the reactant must overcome before arriving at the active site of the catalyst. The reactant must first push into a pocket of the catalyst, before arriving at its final alignment. The C=O of the ketone and the Ru–N bonds lie parallel, thereby maximising orbital overlap with the hydrogen atoms that transfer to form the alcohol. The advantage of computational models is that the changes in geometrical structure as the ketone approaches can be observed ‘frame by frame’. It is then possible to follow the trajectory of approach, analyse the position of the barriers, and view the corresponding changes in atomic structure.

Fig. 4

The potential energy surface of the pathway for the reactant to approach the catalyst active site. The black circles show the positions on the trajectory that are illustrated in Figure 5

Fig. 5

Molecular models showing stages of the reactants’ pathway to the catalyst active site

When the reactant enters the pocket of the catalyst, which begins to occur when the constraint (Ru–H - - - C=O) is between 5 Å and 4 Å, the ketone-catalyst system stabilises. This is shown by the total energy of the system decreasing, before it has to overcome the largest barrier between 3.5 Å and 2.75 Å. The barrier is due to the interaction of the phenyl ring of the ketone with ligands of the catalyst; this interaction increases as the ketone is pulled down onto the active site. At a constraint of 2.75 Å there is a conformational change of the reactant, with the phenyl ring tilting so that all the carbons of the ring are no longer in the plane of the other atomic centres of the ketone. The driving force for this change of conformation is the formation of a hydrogen bond, which holds the ketone as it moves closer to the active site. The hydrogen bond distance reduces further as the reactant is pulled further towards the catalyst, but the position of the oxygen atom does not change greatly; the major movement is that of H(NH) of the catalyst.

For the reactant to leave the pocket there is then another barrier, which requires the alignment of the C=O bond of the ketone to the underlying Ru–N bond of the catalyst. To this end, the carbon moves further down, changing from sp² to sp³ hybridisation, until it is in the same plane as oxygen. This results in the C=O bond lying parallel to the Ru–N bond, maximising overlap with the two hydrogen atoms to be transferred from the catalyst to the reactant. Simultaneously, the Ru–H bond elongates; this would contribute to the barrier at a (Ru–H - - - C=O) constraint of 2.25 Å.

Constrained Optimisations: Other Quadrants

Having considered the approach of the reactant along the favoured quadrant, we then addressed the question of whether computational calculations contain sufficient detail to predict the selectivity of a specific catalyst. A catalyst that forms products with a high enantiomeric excess (ee) is highly desirable, as these reactions are important in a pharmaceutical context. Here the physiological reaction to one enantiomer may differ greatly from that to another. Using docking calculations, the aim is to investigate the various arrangements, and thereby provide insight into how catalysts may be optimised.

The difference in energy barriers between the different quadrants of approach has the potential to provide such discrimination. From Figure 2 and Figure 6 it is obvious that two of the possible orientations are sterically ‘favourable’ (quadrants Q1 and Q3) and two sterically ‘unfavourable’ (Q2 and Q4) (1). However, this is not sufficient to understand the selectivity, as Q1 and Q3 lead to (R) and (S) products, respectively. In fact what this simple analysis demonstrates is that the channels that would be followed by reactants approaching along Q2 and Q4 should be closed. However, of most importance is to understand why certain catalysts produce much higher ee than others. To this end we must understand the difference between Q1 and Q3.

Fig. 6

Starting configurations with a constraint of 5 Å showing the quadrants

Figure 7 shows clearly that for [(S)-XylBINAP-RuH₂-(S,S)-DPEN], Q1 possesses a much lower barrier to approach and would therefore be expected to show high ee. This is confirmed experimentally where, depending upon experimental conditions, an ee of around 97% is achieved.

Fig. 7

The energy profiles for approach of the reactant to the catalyst along different quadrants. (Numbers are normalised to 8 Å, as here the reactant and catalyst do not interact)

With an understanding of the pathway of the reactant to the catalyst active site and of the stable position of the ketone along this path, we could return to considering the calculations of TS and therefore the activation energies for the hydrogenation processes. We are currently considering other catalysts and reactants to ascertain whether we can correlate the barrier for entry into the catalyst pocket with the selectivity of known systems.

Transition States

To make the best use of available computer resources, initial TS calculations were performed on a cut down model of a commercial catalyst (Figure 2). Acetone was initially considered as the reactant molecule. From our understanding of the approach of the ketone to the catalyst, we have been able to determine a stable configuration for the reactant. The constrained optimisation calculations showed that with the ketone hydrogen bonded to the catalyst, (R₂C=O - - - HNH) was at a minimum energy position when the (Ru–H - - - C=O) was held at 2.75 Å, with the hydrogen bond holding the reactant in place. We further optimised this structure to provide the starting point for our TS calculations. The product state has the alcohol physisorbed above the catalyst, and is stabilised by a hydrogen bond between COH - - - NH.

The TS forms a six membered ring (Figure 8). The hydride bond is elongated from 1.7 Å to 1.75 Å. It can be seen from Figure 9 that there is also a change in the hybridisation of carbon, which moves from sp² towards sp³. The TS has a (O)C - - - H(Ru) distance of 1.96 Å and (C)O - - - H(NH) of 1.86 Å. It is possible to visualise the imaginary mode that is associated with the TS, and it is found that the two hydrogen atoms and carbon are the atomic centres that move the most. Figure 10 shows the energy profiles for TS searches using the LST and QST calculation methods.

Fig. 8

Valence bond representation of transition state in hydrogenation of acetone

Fig. 9

Molecular configuration for transition state in hydrogenation of acetone; E_A = 3.07 kcal mol⁻¹

Fig. 10

Energy profiles of TS searches in hydrogenation of acetone

The activation energy that we have calculated is 3 kcal mol⁻¹, which is in agreement with previous calculations on similar-sized models of the commercial catalyst, while the calculated reaction energy is exothermic by 6 kcal mol⁻¹. The hydrogenation of acetone is therefore extremely facile, proceeding as follows:

Incipient bond formation is signalled by the shortening of the O - - - H distance (from 2.11 Å to 1.86 Å) in the N–H - - - O hydrogen bond, and of the Ru–H - - - C distance (from 3.02 Å to 1.95 Å);
Small changes in the same direction are observed in the other bond lengths (< 2%).

The structure of the TS, therefore, resembles much more that of the reactant complex [RuH₂–acetone] than that of the product complex [RuH₂–iPrOH]. This process is therefore a good example of the Hammond principle. This states that the structure of the transition state will resemble that of the product more closely than that of the reactant for endothermic processes, whereas the opposite is true for exothermic reactions. A previous computational study on a trans-dihydro(diamine)ruthenium(II) Noyori-type model catalyst has evaluated a reaction barrier for the hydrogenation of acetone of 3.6 kcal mol⁻¹ at B3LYP/6-31G** level (13). Our results are in apparent agreement with previous calculations on the formaldehyde/methanol transformation by the RuH(NH₂CH₂CH₂NH)(η⁶-benzene) complex performed at B3LYP (14) and generalised gradient approximation (GGA) (15). This shows that classical reaction barriers computed with GGA functionals are smaller than those obtained with B3LYP by about 2 kcal mol⁻¹.

We certainly anticipated that the methodology used would impact on the activation energy (E_A), and we are currently evaluating the effect of changing the density functional. Initial results show that increasing electron localisation by moving from GGA via hybrid to meta functionals leads to a slight increase in E_A. We are also considering other ketones, and attempting to build up a larger model of the catalyst system, so that we can make direct comparison with experimental data for industrially relevant systems.

Conclusions

The two complementary DFT simulation methodologies of transition state searches and constrained geometry optimisations are now yielding results that are of considerable importance to understanding catalyst behaviour, potentially leading to the prediction and design of new catalysts for the ketone hydrogenation reaction.

References

R. Noyori, Angew. Chem. Int. Ed., 2002, 41, (12), 2008 – Nobel Lecture LINK http://dx.doi.org/10.1002/1521-3773(20020617)41:12<2008::AID-ANIE2008>3.0.CO;2-4
R. Noyori and T. Ohkuma, Angew. Chem. Int. Ed., 2001, 40, (1), 40 LINK http://dx.doi.org/10.1002/1521-3773(20010105)40:1<40::AID-ANIE40>3.0.CO;2-5
T. Ohkuma, M. Koizumi, K. Muñiz,, G. Hilt, C. Kabuto and R. Noyori, J. Am. Chem. Soc., 2002, 124, (23), 6508 LINK http://dx.doi.org/10.1021/ja026136+
A. Zanotti-Gerosa,, W. Hems, M. Groarke and F. Hancock, Platinum Metals Rev., 2005, 49, (4), 158 LINK https://www.technology.matthey.com/article/49/4/158-165
A. Hagemeyer, B. Jandeleit, Y. Liu, D. M. Poojary, H. W. Turner, A. F. Volpe and W. H. Weinberg, Appl. Catal. A: Gen., 2001, 221, (1–2), 23 LINK http://dx.doi.org/10.1016/S0926-860X(01)00886-9
‘Manufacturing Molecules Initiative (MMI): A Source of Funding’, URN 02/527, U.K.Department of Trade and Industry, London, 2002
E. J. Palin, G. A. Grasa and C. R. A. Catlow, Mol. Simulat., 2006, 32, (10–11), 901s LINK http://dx.doi.org/10.1080/08927020600982491
B. Delley, J. Chem. Phys., 1990, 92, (1), 508 LINK http://dx.doi.org/10.1063/1.458452
B. Delley, J. Chem. Phys., 2000, 113, (18), 7756 LINK http://dx.doi.org/10.1063/1.1316015
I. D. Gridnev and T. Imamoto, Acc. Chem. Res., 2004, 37, (9), 633 LINK http://dx.doi.org/10.1021/ar030156e
J. M. Brown, ‘Hydrogenation of Functionalized Carbon-Carbon Double Bonds’, in “Comprehensive Asymmetric Catalysis”, eds.E. N. Jacobsen, A. Pfaltz and H. Yamamoto, Springer, 1999, Vol. 1, p. 121
Animation showing result from constrained optimisation calculations for approach of ketone reactant to asymmetric hydrogenation catalyst https://www.technology.matthey.com/images/Q1_new.gif
K. Abdur-Rashid, S. E. Clapham, A. Hadzovic, J. N. Harvey, A. J. Lough and R. H. Morris, J. Am. Chem. Soc., 2002, 124, (50), 15104 LINK http://dx.doi.org/10.1021/ja016817p
M. Yamakawa, H. Ito and R. Noyori, J. Am. Chem. Soc., 2000, 122, (7), 1466 LINK http://dx.doi.org/10.1021/ja991638h
J.-W. Handgraaf, J. N. H. Reek and E. J. Meijer, Organometallics, 2003, 22, (15), 3150 LINK http://dx.doi.org/10.1021/om030104t

Acknowledgements

The author would like to thank R. Catlow, E. Palin and D. Di Tommaso (Royal Institution); J. Xiao, Z. Chen, X. Wu and J. Ruan (University of Liverpool); A. Danopoulos and N. Stylianides (University of Southampton); F. King, F. Hancock and A. Zanotti-Gerosa (Johnson Matthey); P. Hogan and M. Purdie (AstraZeneca); P. Ravenscroft (GlaxoSmithKline); and P. Levett and A. Pettman (Pfizer).

The Author

Sam French holds a degree in Chemistry with Medicinal Chemistry from the University of Warwick, U.K., as well as a Ph.D. from the Royal Institution of Great Britain (RI). At the RI, he worked under the supervision of Prof. Richard Catlow FRS. He held three subsequent postdoctoral research assistant positions at the RI. All involved large collaborations with industrial input, specialising in QM/MM techniques and grid technology in chemistry and catalysis. He joined Johnson Matthey as a Senior Scientist in November 2004, to initiate and lead the Computational Chemistry Group.