A computational study of Grubbs-type catalysts: Structure and application in the degenerate metathesis of ethylene

Journal of Molecular Structure: THEOCHEM 763 (2006) 37–42 www.elsevier.com/locate/theochem

A computational study of Grubbs-type catalysts: Structure and application in the degenerate metathesis of ethylene Ingrid T. Sabbagh, Perry T. Kaye * Department of Chemistry, Rhodes University, Artillery Road, Grahamstown, East Cape 6140, South Africa Received 13 December 2005; accepted 5 January 2006 Available online 10 March 2006

Abstract The ground-state geometries of the first-generation Grubbs pre-catalyst and its active monophosphine derivative have been explored at the molecular mechanics (MM), semi-empirical (SE) and density functional (DFT) levels; the theoretical data for the pre-catalyst have also been compared with the corresponding dimensions of known crystalline analogues. DFT studies of the ligand-dissociation energies for the first- and second-generation Grubbs pre-catalysts and the degenerate metathesis of ethylene using a truncated model of the first-generation catalyst have been undertaken and the results compared with available data. Transition states have been located for each step of the metathesis reaction and the calculated free energies of activation afford new insights into the rate-determining processes. q 2006 Elsevier B.V. All rights reserved. Keywords: Grubbs catalysts; Metathesis; Ethylene; Molecular mechanics; Semi-empirical PM3; Density functional theory; Dmol3

1. Introduction The discovery of efficient ruthenium alkylidene catalysts has been critical in the development of metathesis methods in modern organic synthesis. Notable amongst these are the first-generation tricyclohexylphosphine- (1) and secondgeneration imidazolidine-ruthenium catalysts (2) developed by Grubbs and co-workers [1]. The molecular complexity of the catalytic cycle and the difficulty of isolating putative intermediates have, however, hindered experimentally based mechanistic studies of cross-metathesis reactions. Consequently, several computational studies have been undertaken, using truncated models, to probe experimentally inaccessible mechanistic aspects. These include:- a CarParrinello ab initio molecular dynamics study of the behaviour of the model complex [(PH3)2Cl2RuaCH2] 3 and its reaction with ethylene at various temperatures [2]; DFT studies [3–5]; and, more recently, QM/MM studies by Adlhart and Chen [6,7]. Computer modelling methods constitute a powerful tool in the design and development * Corresponding author. Tel.: C27 46 6038254; fax: C27 46 6225109. E-mail address: [email protected] (P.T. Kaye).

0166-1280/$ - see front matter q 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.theochem.2006.01.014

of new catalytic systems, but the applicability of high-level methods is often limited by the size of the structures involved—a limitation typically addressed by the use of truncated models. With ongoing interests in the development of metal selective ligands [8], on one hand, and the elucidation of organic reaction mechanisms [9], on the other, we have embarked on a programme involving the rational design and synthesis of novel ruthenium alkylidene metathesis catalysts. In this communication, we report the results of studies aimed at: (i) exploring the applicability of rapid and readily accessible molecular mechanics (MM) and semiempirical (SE) methods in modelling ruthenium complexes; (ii) the use of DFT methods to evaluate the reliability of the MM and SE approaches; (iii) demonstrating application of the density functional package, ‘Dmol3’ [10], in determining ligand-dissociation energies for the first- and second-generation Grubbs pre-catalysts (1 and 2, respectively) and in elucidating the mechanism of the degenerate metathesis of ethylene using a truncated model 3 of the firstgeneration Grubbs pre-catalyst 1.

38

I.T. Sabbagh, P.T. Kaye / Journal of Molecular Structure: THEOCHEM 763 (2006) 37–42

2. Computational methods Molecular mechanics calculations were performed using the Accelrys Cerius2w software package [11] on an SGI O2 UNIXbased platform and using the general-purpose universal force field [12]—a force-field parameterized for the full periodic table and validated for numerous structure types including metal complexes. Charges were included via the chargeequilibration method [13], applying the Qeq_charged 1.0 parameter set supplied with the programme. All structures were minimized to convergence using the conjugate gradient and/or Newton–Raphson methods. Following initial minimization, isothermal and anneal simulation dynamics methods (the latter with 300 cycles and a mid-cycle temperature of 1000 K) were employed to explore conformational space and locate the global minimum for each compound. The 10 highest- and 10 lowest-energy conformers obtained from the anneal dynamics routine were minimised to afford a common lowest energy structure as the global minimum. Semi-empirical calculations were performed using the Wavefunction, Inc., Spartan Pro’02 programme [14]; equilibrium geometries were obtained using the semi-empirical PM3 model, which is parameterized for transition metals. Density functional calculations were conducted using the Accelrys Dmol3 DFT code in MaterialsStudio (version 2.2) [10] on LINUX-based Pentium IV PC’s. All calculations involved use of the generalized gradient approximation (GGA) functional by Perdew and Wang (PW91) [15] and the ‘double

numerical plus polarization’ (DNP) basis set: a polarized split valence basis set of numeric atomic functions which are exact solutions to the Kohn–Sham equations for the atoms [16]. Convergence criteria for geometry optimizations were, generally, ˚ , 0.005 A ˚ the threshold values: 2!10K5 Hartree, 0.004 Hartree/A K5 and 1!10 Hartree for energy, force, displacement and selfconsistent field (SCF) density, respectively; for structures 1 and 11–13 and the transition states I–III (Fig. 3), the corresponding ˚ , 0.005 A ˚ and thresholds were: 1!10K5 Hartree, 0.002 Hartree/A K6 1!10 Hartree. All calculations employed a method based on Pulay’s [17] direct inversion of iterative subspace (DIIS) technique to accelerate SCF convergence, using, where necessary, a small electron thermal smearing value of 0.005 Hartree. Preliminary transition state geometries were obtained using the integrated linear synchronous transit/quadratic synchronous transit (LST/QST) method [18], and then subjected to full TS optimization using an eigenvector following algorithm. These geometries were confirmed using intrinsic reaction path (IRP) calculations, based on the nudged elastic band (NEB) algorithm [19], to map the pathways connecting the relevant reactant, TS and product geometries. All structures identified as stationary points were subjected to full frequency analyses, to verify their classification as equilibrium geometries (zero imaginary frequencies) or transition states (one imaginary frequency). The reported energies refer to Gibbs free energy corrections to the total electronic energies at 298.15 K, including zero-point energy (ZPE) corrections, unless stated otherwise.

Table 1 Calculated bond lengths for complex 1 and the corresponding X-ray data for the structural analogues 4 [20], 5 [21] and 6 [22]

˚) Bond length (A

Molecular mechanics

Semi-empirical

DFT

X-ray data for 4

X-ray data for 5

X-ray data for 6

RuaC Ru–Cl1a Ru–Cl2a Ru–P1a Ru–P2a

2.022 2.414 2.412 2.608 2.603

1.851 2.357 2.357 2.313 2.315

1.879 2.455 2.453 2.472 2.490

1.839 2.401 2.395 2.397 2.435

1.811 2.388 2.375 2.368 2.390

1.852 2.400 2.393 2.415 2.422

a

The atom labels, 1 and 2, refer to the structures, permitting differentiation of the respective chlorine and phosphorus atoms.

I.T. Sabbagh, P.T. Kaye / Journal of Molecular Structure: THEOCHEM 763 (2006) 37–42

39

Table 2 Calculated bond angles for complex 1 and the corresponding X-ray data for the structural analogues 4 [20], 5 [21] and 6 [22] Angle (8)

Molecular mechanics

Semi-empirical

DFT

X-ray data for 4

X-ray data for 5

X-ray data for 6

Cl1–Ru–P1 P1–Ru–C P1–Ru–Cl2 Cl1–Ru–P2 C–Ru–P2 Cl1–Ru–C Cl1–Ru–Cl2 C–Ru–Cl2 P1–Ru–P2 Cl2–sRu–P2

89.9 95.4 88.8 91.1 90.6 96.6 175.1 88.2 173.8 89.7

88.0 95.6 90.9 91.4 97.7 97.9 174.7 87.4 166.7 88.4

91.2 94.5 89.7 83.6 103.0 93.6 159.6 106.6 162.0 83.6

87.2 97.5 91.5 90.8 101.2 88.7 167.6 103.7 161.1 86.5

91.65 98.07 87.59 88.16 97.26 89.16 173.50 97.34 142.29 90.86

89.51 98.59 89.09 91.07 98.92 97.80 173.44 88.74 162.23 88.33

3. Results and discussion The Grubbs first-generation catalyst 1 was modelled both de novo and using crystal structure data for the p-chloro analogue 4, the X-ray analysis of which had been reported by Grubbs and co-workers [20]; on minimisation, both models converged to the same structure. Isothermal and anneal dynamics routines at the MM level were then used to explore conformational space and locate the global minimum, which was subjected to geometry optimisation at the semi-empirical PM3 level using the SpartanPro package, followed by DFT calculations at the PW91/DNP level using Dmol3. The results of all three sets of calculations are compared with the X-ray structure data for the complexes 4, 5 and 6 in Tables 1 and 2. The X-ray structure of the p-chloro analogue 4 reveals a distorted square pyramidal coordination geometry about the metal centre, with the carbene moiety perpendicular to the P–Ru–P axis and an almost linear Cl–Ru–Cl angle [20]. It is, of course, possible that complex 1 and its p-chloro analogue 4 might favour somewhat different conformations in the crystal lattice. However, since no X-ray data is available for the Grubbs first-generation catalyst 1, it seems reasonable to use the reported [20] X-ray structure of the analogue 4 as a reference model. A search of the Cambridge Structural Database revealed that compound 4 was, in fact, the only ruthenium bis(tricyclohexylphosphino)dichlorobenzylidene complex for which X-ray crystal structure data had been reported. In this structure, and in the X-ray structures of the ethoxymethylene [21] and vinylogous [22] analogues 5 and 6, respectively, it is apparent that the metal coordination geometry is largely conserved (see Tables 1 and 2). While conformational differences between isolated (computed) and lattice-bound molecules cannot be precluded and while there are some obvious discrepancies, the tabulated data reveal that all three modelling approaches appear to afford reasonably credible and, from a design perspective, useful representations of the coordination geometry of the Grubbs first-generation catalyst 1. The SE and DFT methods performed equally well for bond length calculations (Table 1), with mean errors ˚ ; relative to structure 4) approximately, half of (0.057G0.002 A ˚ ). The DFT structure of complex 1 the MM value (0.118 A (Fig. 1) correlates most closely with previously computed bond angles [5] and with the X-ray geometry [20] of the structural

analogue 4. The mean absolute errors for the MM, semiempirical and DFT calculated bond angles were 6.5, 4.7 and 3.78, respectively (Table 2). In order to assess the capacity of the Dmol3 package to reproduce a co-ordinatively unsaturated ruthenium system, the active monophosphine derivative 7 of the pre-catalyst 1 was modelled at the DFT level (Fig. 1). The retention of the coordination geometry in the modelled monophosphine complex 7 is, in fact, consistent with reported computational data for several Grubbs-type catalysts [5], as well as the X-ray data for an isolated ruthenium monophosphine complex [20]. Attention was then turned to exploring the application of Dmol3 in calculating the energies of full Grubbs-type catalyst systems. The electronic ligand-dissociation energies (DE) were thus calculated for the dissociation of tricyclohexylphosphine (PCy3) from the Grubbs first- and second-generation pre-catalysts 1 and 2 (Fig. 2). The calculated value of DE for dissociation of the first-generation pre-catalyst 1 (34.8 kcal molK1) was found to be 79.6% of DE for the second-generation system 2 (43.7 kcal molK1). This not only compares very favourably with the corresponding result (80.8%) obtained by Adlhart and Chen (at the BP86 density functional level using the ADF programme) [7] and reasonably well with the 70.4% calculated by Cavallo (at the same level) [5], but also closely parallels experimental trends observed by Grubbs [23].

Fig. 1. Dmol3/GGA/PW91/DNP optimized structures of the pre-catalyst 1 and its active monophosphine derivative 7 (H atoms omitted for clarity).

40

I.T. Sabbagh, P.T. Kaye / Journal of Molecular Structure: THEOCHEM 763 (2006) 37–42

Fig. 2. Calculated ligand dissociation energies (DE) for precatalysts 1 and 2 (mass balanced for leaving group L).

Molecular mechanics (using Cerius2w) and semi-empirical (PM3; using Spartanw) studies were also undertaken for the monophosphine 7, the p-complex 9 and the metallacycle 10, permitting the co-ordination geometry, in each case, to be compared (Table 3) with the corresponding full structures modelled at the DFT (PW91/DNP) level using Dmol3. From the tabulated data, it is apparent that while the bond lengths calculated at each level are reasonably consistent, some of the bond angles vary quite widely. Moreover, geometry optimisation of the metallacycle 10 at the SE (PM3) level resulted in significant distortion of the tetrahedral geometry of the sp3 carbons in the 4-membered ring [24]. These observations indicate a need for caution in the use of both the MM and SE models. However, it should be noted that even the X-ray data reported for the ruthenium complexes 4–6 (Table 2) exhibit some significant bond-angle variations.

The computational time required for high-level calculations involving the full catalyst systems militates against their use in detailed mechanistic studies. Consequently, truncated analogues, such as complex 3 (Scheme 1), and model substrates, such as ethylene, are typically used to minimise computational time [2,3,7]. However, although facilitating computation, the use of truncated analogues may well fail to accommodate or reflect significant steric or electronic effects in the full systems. Moreover, theoretical data obtained for isolated molecules may

Table 3 Comparison of bond lengths and angles for the monophosphine 7, the p-complex 9 and the metallacycle 10, calculated at the MM and SE levels,a with the corresponding dimensions obtained at the DFT levelb p-complex 9

Monophosphine 7 MM

SE

DFT

Bond (A˚) RuaC Ru–Cl1 Ru–Cl2 Ru–P1 Ru–P2

–

–

–

Angle (8) Cl1–Ru–P1 P1–Ru–C P1–Ru–Cl2 Cl1–Ru–P2 C–Ru–P2 Cl1–Ru–C Cl1–Ru–Cl2 C–Ru–Cl2 P1–Ru–P2 Cl2–Ru–P2

90.8 96.7 90.1 – – 89.8 179.0 90.7 – –

92.47 93.98 90.82 – – 97.48 162.70 99.24 – –

92.44 96.78 93.08 – – 92.73 171.00 93.68 – –

a b c d

1.98 2.40 2.40 2.59

1.83 2.28 2.28 2.28

2.03 2.24 2.24 2.36

Using the Cerius2w and Spartanw programmes, respectively. For the full structures using Dmol3. To the ethylene p-bond. To CH2 in the metallacycle.

MM 2.02 2.42 2.42 2.59 2.16c 89.2 89.8 87.4 90.5c 95.8c 86.9 174.4 97.6 174.4c 92.4c

Metallacycle 10 SE

DFT

1.86 2.34 2.33 2.33 2.02c

1.89 2.42 2.42 2.45 2.29c

87.06 96.51 82.11 95.00c 93.48c 87.32 168.43 97.92 169.89c 94.98c

91.06 94.94 85.77 91.75c 98.78c 96.62 157.48 105.86 165.58c 86.23c

MM 2.24 2.40 2.40 2.58 2.21d 91.7 110.6 90.3 88.3d 62.9d 97.2 175.0 86.4 173.4d 90.3d

SE

DFT

2.10 2.30 2.35 2.36 1.94d

1.99 2.45 2.37 2.36 2.00d

83.41 99.26 82.32 121.53d 67.56d 93.69 164.73 93.94 151.28d 73.68d

80.39 136.07 95.10 85.81d 85.43d 86.97 175.47 96.67 134.61d 97.15d

I.T. Sabbagh, P.T. Kaye / Journal of Molecular Structure: THEOCHEM 763 (2006) 37–42

41

Scheme 1.

Fig. 3. DG298 profile for the degenerate metathesis of ethylene with the model catalyst 3.

not provide a completely satisfactory model for condensed phase reactions. The Dmol3 code was applied, at the PW91/DNP level, to explore the degenerate methathesis of ethylene with the truncated Grubbs catalyst model 3 via the widely accepted [4,7,25] dissociative pathway detailed in Scheme 1. The calculated Gibbs free energy (DG298) profile for this pathway is illustrated in Fig. 3, and the corresponding reaction energies are summarised in Table 4. The ground-state and transition-state structures were verified by vibrational analysis, and transitionstate structures were linked to the relevant intermediates by mapping the respective intrinsic reaction paths. The free energies (G), obtained by thermodynamic corrections to the electronic energies (E), are of particular interest in fundamental studies of this type since enthalpic and entropic effects are not reflected by the electronic energies (E). Such effects may be critical in steps involving a change in molecularity (e.g. 3/11 and 11/12; Table 3), where entropic effects render E values meaningless. The data in Table 4 correlates reasonably well with the results (Table 5) obtained by Adlhart and Chen [7] and by Cavallo [5] at the DFT/BP86 level, and no significant differences are apparent between the reported geometries and those obtained in the present study. In our study, however, we have also located the three relevant transition states, I–III (Figs. 3 and 4)—I and II, apparently for the first time [26]. In the truncated system examined here, the largest free-energy of activation (DG‡298 ; Table 4) is associated with formation of the monophosphine species 11, but transition state III is clearly the

highest point on the calculated free-energy profile and its formation is thus rate-determining. However, transition state III is only 1.0 kcal/mol higher in energy than transition state II in the preceding step, suggesting that both of these steps are likely to be kinetically significant. It should also be noted that standard states are assumed in the theoretical study whereas, in practice, the alkene substrate is generally present in large excess. Moreover, the results in Tables 4 and 5 are not solvent Table 4 Dmol3//PW91/DNP reaction energies for the degenerate metathesis of ethylene with the model catalyst 3 (Scheme 1, Fig. 4) Reaction

DE (kcal/mol)

Ea (kcal/mol)

DG298 (kcal/mol)

DG‡298 (kcal/mol)

3/11 11/12 12/13

22.07 K14.81 3.31

19.96 K2.71 10.39

9.44 K1.99 4.88

14.25 6.82 9.81

Table 5 Reaction energies for the degenerate metathesis of ethylene with the model catalyst 3, as calculated by Adlhart and Chen [7] and, in parentheses, by Cavallo [5] Reaction

DE (kcal/mol)

Ea (kcal/mol)

DG298 (kcal/mol)

DG‡298 (kcal/mol)

3/11 11/12

21.5 K14.8 (K11.5) 7.6 (5.6)

– –

8.2 K0.2

– –

12/13

14.8 (12.0)

9.1

13.7

42

I.T. Sabbagh, P.T. Kaye / Journal of Molecular Structure: THEOCHEM 763 (2006) 37–42

Acknowledgements The authors thank Sasol Technology Ltd, for a bursary (I.T.S) and Sasol Technology Ltd, the Technology and Human Resources for Industry Programme (THRIP; Project no. 2642) and Rhodes University for generous financial support, and the referee for helpful comments.

Fig. 4. Geometry-optimised structures of the transition states I–III.

corrected, but Cavallo [5] has shown that the inclusion of solvent effects has a negligible effect on the reaction profile. The geometry of the ruthenium complex in the late transition state I, corresponds closely to that of intermediate 11, the distance between the metal centre and the departing ˚ ) being greater than twice the Ru–P bondphosphine (5.597 A length in the precursor 3. The relative orientation of the interacting species in the transition state II is clearly supported by the vibrational analysis; there is little perturbation in the geometry of the ruthenium complex (compared to intermediate ˚ ) between the proximal ethylene 11) and the distance (4.507 A carbon and the metal centre is consistent with formation of an early transition state complex. Transformation of the p-complex 12 to the metallacycle 13, however, proceeds via a late transition state III, in which:- the ethylene moiety is linked as a s-complex to the metal centre; the Cl–Ru–Cl bond angle approaches 1808; the carbene methylene group is rotated through 908; and there is evidence of an incipient s-bond between the eclipsed, methylene termini of the acyclic ˚ ). precursor (C–C distance: 2.333 A The location of transition states I and II on the free-energy surface is significant as, in previous studies [7,25], formation of both the monophosphine 11 and the p-complex 12 have been found to be barrierless [26]. 4. Conclusions While application of the universal force field in the Cerius2 MM package has provided structures that approximate the X-ray crystal and DFT co-ordination geometries of the ruthenium complexes examined, caution is clearly needed in drawing detailed conclusions from the MM structures. Use of the Dmol3 code in exploring ligand–dissociation energies and reaction energetics has afforded results that compare favourably with those reported by other researchers. Furthermore, all three transition states involved in the widely accepted dissociative metathesis of ethylene using a truncated model of the first-generation Grubbs pre-catalyst 2 have been located. The calculated 1 kcal/mol free energy difference between transition states II and III suggest that both steps (11/12 and 12/13) are likely to contribute to the overall kinetics.

References [1] P. Schwab, R.H. Grubbs, J.W. Ziller, J. Am. Chem. Soc. 118 (1996) 100. [2] O.M. Aagaard, R.J. Meier, F. Buda, J. Am. Chem. Soc. 120 (1998) 7174. [3] F. Bernardi, A. Bottoni, G.P. Miscione, Organometallics 22 (2003) 940. [4] S. Fomine, S.M. Vargas, M.A. Tlenkopatchev, Organometallics 22 (2003) 93. [5] L. Cavallo, J. Am. Chem. Soc. 124 (2002) 8965. [6] C. Adlhart, P. Chen, Angew. Chem. Int. Ed. 41 (23) (2002) 4484. [7] C. Adlhart, P. Chen, J. Am. Chem. Soc. 126 (2004) 3496. [8] See, for example J.P. Hagemann, P.T. Kaye, J. Chem. Soc., Perkin Trans. 1 (1999) 341–347; A. Daubinet, P.T. Kaye, Synth. Commun. 32 (20) (2002) 3207–3217; P.T. Kaye, T.R. Tshikhudo, J. Chem. Res. (S) (2003) 179–181 (J. Chem. Res. (M) (2003) 0401–0418); P.T. Kaye, B.S.B. Gxoyiya, J.P. Hagemann, J. Chem. Res. 4 (2004) 252–256. [9] See, for example P.T. Kaye, M.J. Mphahlele, M.E. Brown, J. Chem. Soc., Perkin Trans. 2 (1995) 835–838; P.T. Kaye, R.S. Robinson, Tetrahedron 54 (1998) 13253–13256; P.T. Kaye, M.A. Musa, X.N. Nocanda, S. Robinson, Org. Biomol. Chem. 1 (2003) 1133–1138. [10] B. Delley, J. Chem. Phys. 92 (1990) 508; B. Delley, J. Phys. Chem. 100 (1996) 6107; B. Delley, J. Chem. Phys. 113 (2000) 7756. [11] Cerius2w, Version 4.0. [12] A.K. Rappe´, C.J. Casewit, K.S. Colwell, W.A. Goddard, W.M. Skiff, J. Am. Chem. Soc. 114 (1992) 10024. [13] A.K. Rappe´, W.A. Goddard, J. Phys. Chem. 95 (1991) 3358. [14] Wavefunction, Inc., 18401 Von Karman Avenue, Suite 370, Irvine, CA 92612. [15] J.P. Perdew, Y. Wang, Phys. Rev. B 45 (1992) 13244. [16] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [17] P. Pulay, J. Comput. Chem. 3 (1982) 556. [18] T.A. Halgren, W.N. Lipscomb, Chem. Phys. Lett. 49 (1977) 225. [19] G. Henkelman, H. Jonsson, J. Chem. Phys. 113 (2000) 9978. [20] P. Schwab, M.B. France, J.W. Ziller, R.H. Grubbs, Angew. Chem. Int. Ed. Engl. 34 (1995) 2039. [21] J. Louie, R.H. Grubbs, Organometallics 21 (2002) 2153. [22] S.T. Nguyen, R.H. Grubbs, J.W. Ziller, J. Am. Chem. Soc. 115 (1993) 9858. [23] M.S. Sanford, M. Ulman, R.H. Grubbs, J. Am. Chem. Soc. 123 (2001) 749. [24] Depending on the starting MM-minimised structure, SE (PM3) calculations afforded distorted metallacycles exhibiting bonding of ruthenium to the benzylic proton or (at higher energy) to one of the methylene protons. [25] S.F. Vyboishchikov, M. Bu¨hl, W. Thiel, Chem. Eur. J. 8 (2002) 3962. [26] Although transition states I and II each exhibit a single imaginary frequency, the presence of low-frequency (!50 cmK1) vibrations (three for TS I and one for TS II) may well affect the accuracy of the corresponding free-energy calculations.