ARTICLE pubs.acs.org/JPCA
Force Field Parametrization and Molecular Dynamics Simulation of Flexible POSS-Linked (NHC; Phosphine) Ru Catalytic Complexes Amirhossein Ahmadi,† Carl McBride,†,‡ Juan J. Freire,*,† Anna Kajetanowicz,§ Justyna Czaban,§ and Karol Grela§ †
Departamento de Ciencias y Tecnicas Fisicoquímicas, Facultad de Ciencias, Universidad Nacional de Educacion a Distancia, Senda del Rey 9, 28040 Madrid, Spain § Institute of Organic Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland ABSTRACT: In recent years, N-heterocyclic carbene (NHC) or phospine groups have been put forward as candidate catalysts ligands for olefin metathesis reactions to be performed using multistep methods. Some of these proposed ligands contain polyhedral oligomeric silsesquioxane (POSS) structures linked to NHC rings by means of alkyl chains. Some important properties for the prediction of catalytic activity, such as the theoretically defined buried volume, are related to the conformational characteristics of these complex ligands that can be studied through molecular dynamics simulations. However, the chemical structure of resulting catalytic complexes usually contains atoms or groups that are not included in the common forcefields used in simulations. In this work we focus on complexes formed by a catalytic metal center (Ru) with both phospine and POSS-linked NHC groups. The central part of the complexes contain atoms and groups that have bonds, bond angles, and torsional angles whose parameters have not been previously evaluated and included in existing force fields. We have performed basic ab initio quantum mechanical calculations based on the density functional theory to obtain energies for this central section. The force field parameters for bonds, bond angles, and torsional angles are then calculated from an analysis of energies calculated for the equilibrium and different locally deformed structures. Nonbonded interactions are also conveniently evaluated. From subsequent molecular dynamics simulations, we have obtained results that illustrate the conformational characteristics most closely connected with the catalytic activity.
’ INTRODUCTION Computational techniques offer the opportunity to explore many aspects of the behavior of molecular systems. In the present work, we have investigated the possibility of using molecular simulations in the study of particularly complex structures that act as catalysts in important organic chemistry reactions. Olefin metathesis is a remarkable topic in current chemistry because it is a very efficient and elegant method of formation of CdC double bonds.14 It is widely used not only in small-scale laboratory research chemistry, but also in large-scale industrial production.5 Metathesis owes its unique role in modern organic chemistry to well-defined Ru-based (pre)catalysts. Replacement of a phosphine ligand with a N-heterocyclic carbene (NHC) ligand in Grubbs first generation catalysts (Gru I)6 leads to NHC-based second generation catalysts (Gru II, Ind II, Ind II0 ) that typically display better thermal stability and activities compared to first generation catalysts (Figure 1).7 Although the metathesis (pre)catalysts have been successfully designed for stability,8 functional group tolerance,9 activity,10 and selectivity,11 enabling metathesis methodology to be broadly applied, this methodology does have some drawbacks. One of the most arduous tasks is the difficulty associated with the separation of the desired product(s) from the ruthenium catalysts, although r 2011 American Chemical Society
Figure 1. Examples of Grubbs (Gru) and indenylidene (Ind) type ruthenium catalysts. Generation two includes saturated and unsaturated mesityl-substituted (SImes and IMes) NHC rings.
they are typically used in small quantities, usually 15 mol %. This inconvenience prevents the wide employment of metathesis Received: July 1, 2011 Revised: September 16, 2011 Published: September 20, 2011 12017
dx.doi.org/10.1021/jp2062332 | J. Phys. Chem. A 2011, 115, 12017–12024
The Journal of Physical Chemistry A
ARTICLE
Figure 2. Complex with alkyl chain links of three carbons.
reactions in pharmaceutical synthesis, wherein permitted levels of such metallic contaminants are very low (for example, only 5 ppm or below a total of 20 ppm if two or more are present for oral administration of members of the platinum group [Pt, Pd, Ir, Rh, Ru, and Os]).12 Because even multistep methods, incorporating double purification on silica gel and treatment with activated carbon,13 fail to reach an acceptable level of ruthenium impurities, new techniques concerning separation are a subject of interest for many researchers.14 In previous work by some of the authors, the feasibility of connecting a second generation indenylidene catalyst with two POSS (polyhedral oligomeric silsesquioxane) moieties (for an example, see Figure 2) was explored. A catalyst of this type should be active and stable and, additionally, because of its bulky substituents, should be easily separated after reaction via filtration through polymeric membranes. However, although we were able to obtain the desired catalyst, the yield was very low (around 10%). Moreover, an activity test in the ring closing metathesis reaction of diethyl diallylmalonate gave rather poor results. These observations prompted us to study the optimal length of the linker between the aromatic ring and the POSS moiety. Different alkyl chain lengths imply different sizes of the substituents close to the catalytic center and a different degree of flexibility, which may have consequences for the catalytic activity. The main purpose of this paper is to establish the basis for performing adequate numerical calculations that can help in the understanding and prediction of the catalytic performance of these types of flexible complexes. A precise theoretical prediction of such features can only be carried out using the help of molecular simulations; in particular, using molecular dynamics (MD). This technique is able to describe the dynamic trajectory of the structure by solving classical mechanical equations.15,16 Given an array of atoms and a distribution of velocities (the latter associated to the system temperature), these equations are able to produce new positions and new velocities. These equations are numerically solved using a simulation time step sufficiently small so that forces can be considered to be constant and determined by the potential energy defined by the given set of atomic coordinates. To proceed with this iterative numerical scheme, a potential energy depending on the atomic coordinates should be considered. According to the molecular mechanics approach, different functions are defined to describe interactions between different atoms and atomic groups. Essentially, the parameters define equilibrium values and force constants for bonds and bond angles. Changes in torsional angles, which are directly related with the molecular flexibility, are also considered and they are described by means of adequate functions, having several maxima and minima to represent energy barriers. Moreover, van der
Waals and electrostatic interactions between atoms are described through relatively simple functions. All of this information is contained in the so-called force field files.17 In this work, we make use of the PCFF force field.18 PCFF is a member of a consistent family of second generation force fields. They have been parametrized considering a wide range of experimental properties for organic compounds containing hydrogen, carbon, nitrogen, oxygen, sulfur, phosphorus, halogen atoms and their ions, alkali metal cations, and several biochemically important divalent metal cations. However, PCFF (or any other similar force fields19) does not contain some of the important information needed for the definition of the central part of the complexes considered in this work. In particular, it does not contain information about the parameters corresponding to the carbene atom in the NHC ring or how to describe the somewhat loosely bonded structure around the metal atom. To rectify this paucity of force field data for the region of interest, our approach has been to perform basic quantum mechanics calculations on this central section. An energy minimization is first carried out to determine its equilibrium constants. The results are then compared with bibliographic experimental data corresponding to similar structures. Further calculations provide results for the energies corresponding to different locally deformed configurations. These results are then fitted to obtain the different new parameters. These parameters are then added to the standard PCFF force field, and subsequently, MD simulations of the catalytic complex structures immersed in an organic solvent are carried out. The analysis of the results allows us to investigate some structural differences in the central part of the complex due to its bulky outer parts and also provides some conclusions concerning conformational characteristics that may affect their catalytic performance.
’ METHODOLOGY We have made use of the Materials Studio (MS) software suit for some of the steps involved in the present study.20 The MS graphic module “Visualizer” has been employed to build the different complex structures. Quantum mechanics calculations are performed for the energy by means of the module “DMol3”. This module allows us to perform an ab initio calculation based on the density functional theory for the electronic structure and energy for organic and inorganic molecules. Because the feasibility of this calculation strongly depends on the number of atoms involved, we have restricted ourselves to the study of the complex central part (Figure 3) containing the atomic features that need to be parametrized. In DMol3, some different methodological options can be considered. We have undertaken a preliminary study to determine the ability of the different options available in DMol3 to 12018
dx.doi.org/10.1021/jp2062332 |J. Phys. Chem. A 2011, 115, 12017–12024
The Journal of Physical Chemistry A
ARTICLE
Figure 3. Optimized central part of the complex.
Figure 5. Fit of the DMol3 results for the energy variation with the carbeneRu distance to eq 1. Black squares, theoretical results; circle and spline line, eq 1 with the parameters shown in Table 2.
Figure 4. Variation of the potential energy, VT, with the torsional angle for n-butane. Solid line, ab initio calculations; symbols, PCFF calculations (total energy); dotted line, torsion barrier contribution; VΦ, also from PCFF.
predict the well-known parameters corresponding to the structure of simple alkyl molecules (bond lengths and angles and, particularly, the description of the potential energy corresponding to their torsional barriers). This verification has been carried out using the same procedures that have been subsequently used to study the complexes as described below. It reveals the need to use a gradient-corrected method (GC-PW91) in which all electrons are included,21,22 considering a second set of valence atomic orbitals and a polarization d-function on all nonhydrogen atoms (in the DMol3 module, this corresponds to the DND basis set). We should note that the functional density theory has been previously employed to obtain other force field parameters of some molecules.23 However, because the original PCFF schemes were largely based on the HF/6-31G* functional,18 we have performed a detailed comparison between the energy results obtained with DMol3 (with the indicated specifications) and the PCFF forcefield for different values of the main torsional angle of the n-butane molecule. We have verified that both sets of data are similar (see Figure 4). This seems to validate the use of the chosen functional and algorithm. It is important to estimate the contribution of the nonbonded interactions separately from the rest of the interactions,24 since repulsions between nonbonded atoms (hydrogen and carbon
atoms separated by four or five bonds) are not taken into account in the definition of the PCFF torsional barrier function (though they are reflected in the quantum calculations). However, these repulsive regions are not relevant from a simulation point of view, since they are implemented via nonbonded interactions. Thus, a direct fit of the DMol3 results to the PCFF torsional function data in the region of the minima should provide us with an adequate description of the torsional barrier. An initial geometry minimization yields an equilibrium structure for the central section of the molecule. This structure is shown in Figure 3. It should be noted that in most Grubbs structures the X-ray data25,26 show that hydrogens in the methine group are located in the plane formed by chlorine and Ru atoms. This type of discrepancy may be due to steric effects associated with the presence of bulky substituents. Further calculations are performed for different deformed configurations of the central structure. Because we are focusing on specific features within a many atom structure, we have not carried out a general study in terms of normal coordinates,23,27 but have restricted our energy calculations to a number of locally deformed configurations of bond distances or angles. (In the case of bond and torsional angle variations each particular deformation may involve several contributions). The results corresponding to different values of bond distances, bond angle values or torsional angles are fitted to expressions to describe these features included in the PCFF forcefield. In particular one has the following expressions Vb ¼ kb ðb beq Þ2
ð1Þ
for bond length, b, and Vθ ¼ kθ ðθ θeq Þ2
ð2Þ
for bond angle θ (in radians). The subscript eq denotes the equilibrium value. See Figures 5 and 6 for illustrations of the fits carried out to determine the parameters for bonds and bond angles. For the torsional barriers, we use the function VΦ ¼ kΦ ½1 þ cosðmΦ δÞ
ð3Þ
for the barrier corresponding to torsional angle Φ (m is the number of minima and δ is a phase angle). The different ki 12019
dx.doi.org/10.1021/jp2062332 |J. Phys. Chem. A 2011, 115, 12017–12024
The Journal of Physical Chemistry A
ARTICLE
Table 1. Comparison between our Theoretical Results for Bond Lengths (in Å) and Bond Angles (in Degrees) for the Central Structure and Experimental Data Obtained from Diffraction Experiments.a Bond lengths
a
b
c
d
Our results
Carbene-Ru
2.086
RuP
2.3975
2.090
2.069
2.11
2.082
2.407
2.4119
RuC
1.870
1.867
1.841
1.840
1.851
RuCl(1)
2.362
2.375
2.393
RuCl(2)
2.400
2.404
2.383
VVDW ðrÞ ¼ εij ½2ðrvij =rÞ9 3ðrvij =rÞ6
1.360
1.358
Carbene-N(1)
1.354
1.343
1.358
with the combination rules rvij ¼ ½ðrvii Þ6 þ ðrvjj Þ6 1=6 =21=6 εij ¼ 2ðεii εjj Þ1=2 ðrvii rvjj Þ3 =½ðrvii Þ6 þ ðrvjj Þ6
a
b
c
CarbeneRu-P
165.73
162.71
163.2
CarbeneRu-C
104.3
105.4
99.2
ð5Þ
d
Our results 176.1
95.7
91.2
168.6
170.8
146.9
Cl(1)Ru-C
87.1
104.0
106.5
Cl(2)Ru-C
104.3
106.5
PRuC PRuCl
97.1 89.86
92.6 90.4
N-Carbene-N
101.0
CarbeneRu-C(ring)
112.1
112.2
CRuP
97.1
92.6
ClRuP
89.9
86.9
ClRuCl
ð4Þ
2.39
1.366
CarbeneRu-Cl
parameters are constants associated to the simplest description of forces. The van der Waals interactions between atoms i and j at R distance are expressed as
2.39
Carbene-N(1)
Bond angles
Figure 6. Fit of the DMol3 results for the energy variation with the carbeneRuP bond to eq 2. Black squares, theoretical results; circle and spline line, eq 2 with the parameters shown in Table 2.
2.417
161.25
162.84
88.3
106.3
103.4
90.4
N(1)-carbene-Ru
125.4
128
N(1)-carbene-Ru
128.3
128
a
ð6Þ
It is worth stating that, from the point of view of performing molecular simulations, it is important to have precise information about the equilibrium values. Our equilibrium results can be compared with bibliographic data obtained from X-ray experiments for similar model structures.25,26,28 This comparison is shown in Table 1. It can be verified that our bond length results are remarkably close to the experimental data. Some larger differences are found for the angles centered on the metal atom, which, for the model molecules, try to adapt to the presence of more bulky substituents in the plane perpendicular to the phosphineRucarbene axis. Therefore, this axis shows a greater deviation from linearity than in the case of our central structure. Also, the angles formed by the Ru atom and perpendicular atoms clearly depend on the substituents (significant asymmetries can be noticed in the experimental data corresponding to the different model structures). Therefore, we conclude that our calculations for the simplified central structure of the complex yield reasonable values of bond lengths and bond angles. A very precise determination of the bond length and bond angle force constants, kb and kθ, is not quite so critical (actually, some simplified descriptions go as far as to use totally rigid bonds and angles). For the description of torsional barriers we have performed a full range sweep of the torsional angle values (in most cases, several torsional angles are coupled in a single atom displacement). As in the case of butane, we have tried to discern the cases where abnormal increases of energy along the barrier are due to the presence of strong repulsive interactions between nonbonded atoms, only performing fits to the energies
(a) Ref 26, Figure 1; (b) ref 26, Figure 2; (c) ref 25, Table 3; (d) ref 28, Figure 4.
corresponding to the unaffected range of angles. The nonbonded interactions should be computed separately from the rest of the contributions.24 Fortunately, the disposition of atoms in the central part of this simplified central structure (particularly the almost colinear carbeneRu and RuP bonds) avoids any significant presence of nonbonded interactions. Actually, the ab initio calculated variations of energy are always compatible with reasonable values of the torsional force constants. Figure 7 shows the fit corresponding to some coupled torsions on bonds involving the complex metal atom. The new torsional barrier parameters correspond to mainly rigid parts of the molecule. Therefore, the most relevant parameters to determine are the locations and the number of barrier minima (the flexible parts of the structure, particularly the alkyl chain links, contain groups whose interactions are fully described in the original PCFF forcefield.) These parameters are summarized in Table 2. As we explained above, we do not expect dramatic changes in the properties of these complexes due to nonbonded interactions in the region surrounding the Ru atom. However, PCFF does not include parameters for the Ru atom and some van der Waals potential-well values previously reported for Ru and similar metals,24,29 are an order of magnitude smaller than the values of εii shown in the PCFF file included in MS for Pd, Pt and other transition metals. For these reasons, we have performed an independent evaluation, using the GC-PW91 functional.21,22 To this end, we have calculated the interaction energies between these metal centers and an inert Ar atom, which is employed as a test particle.30 (The potential for ArAr interactions obtained 12020
dx.doi.org/10.1021/jp2062332 |J. Phys. Chem. A 2011, 115, 12017–12024
The Journal of Physical Chemistry A
ARTICLE
Table 2. Selection of New Parameters Incorporated to the Force Field that are Related to the Special Features Present in the Complex Central Structure
Rucarbene
2.084
140
RuP
2.42
85
RuCl
2.39
110
RuC
1.85
300
NC (ring) carbeneN
1.385 1.358
380 470
PC
1.839
220
θeq ()
bond angle
Figure 7. Fit of the DMol3 results for the energy variation with the torsional angle NcarbeneRuP to eq 3. Black squares: theoretical results. Circle and spline line: eq 3, with the combined contributions of the NcarbeneRuP, NcarbeneRuCl, and Ncarbene RuCl barrier parameters shown in Table 2 (the coupled variation of the torsional angle associated with a second N atom in the NHC ring is also included).
from experiments are closely reproduced by the aforementioned functional as provided in the software, whereas a similar reproduction for other noble gases would require a further nonimplemented parametrization.31) The metal parameters are then derived from fits of the results to eq 4. We have verified that the εij parameters for the PdAr and PtAr interactions obtained with this fitting procedure differ by less than 10% from those predicted applying the PCFF values and combination rules. Therefore, it seems apparent that the applied method can provide an adequate estimate of the van der Waals εij parameter for this type of atom. Performing a similar calculation for the interaction between Ru and Ar and comparing the results with those corresponding to the PdAr and PtAr cases, we have finally arrived at the van der Waals parameters of Ru that are shown in Table 2. These results are consistent with the range of parameter values exhibited by other transition metals in the PCFF file. The MD simulations are prepared following several steps. First, we build the whole molecular structure with the tools contained in the graphical MS platform. In particular, we consider the unsaturated (NHC; phospine) complex shown in Figure 2 with links formed by an alkyl chain of 3 carbons (n = 3) and similar structures with n = 5 and 7. At this stage, charges for the structure are calculated by using the charge equilibration method.32 We have verified that these results are close to those provided by our previous ab initio calculations for the central region. Charge values corresponding to the most significant central atoms are also incorporated in Table 2. Next, we build a simulation box having periodic boundary conditions using the MS module “Amorphous Cell”. Amorphous Cell uses different simulation techniques and an initial minimization to create an initial configuration having a reasonably low energy. Our systems are constituted by a single complex structure plus 170 organic solvent (dichloromethane) molecules contained in a box that has cubic periodic boundary conditions. The volume of this initial box is fixed so that density is about 70% of the
kb (kcal/mol/ Å2)
beq (Å)
bond
kθ (kcal/mol/rad2)
carbeneRuP
176
25
carbeneRuCl
88
75
carbeneRuC
91
60
ClRuCl
147
23
ClRuC RucarbeneN
106.5 128
27 35
carbeneNC(ring)
91
60
NC(ring)C(ring)
106
90
CRuP
92.6
40
ClRuP
90.4
55
NcarbeneRu
128
35
torsion angle
kΦ (kcal/mol/rad)
m
δ ()
0.7 0.7
2 2
0 180
NcarbeneRuP NcarbeneRuC NcarbeneRuCl
0.7
2
180
carbeneRuPC
0.13
3
180
CRuP-C
0.13
3
0
ClRuPC
0.13
3
0
carbeneRuCH
0.9
2
180
*C(ring)C(ring)*
4.07
2
0
*NC(ring)* *carbeneN*
2.25 2.25
2 2
0 0
Nonbonded Interactions van der Waals
εii (kcal/mol)
riiν (Å)
Ru
3
2.9 i
qi
Ru carbene in NHC ring
0.372 0.25
Cl
0.43
highest charges for central atoms
expected real density at room temperature. Subsequently we run a short MD simulation for each system using the MS module “Discover”. This simulation provides some output files needed for the final MD simulation runs, which are described below. A difference with respect to the DMol3 module is that both Amorphous Cell and Discover make use of a standard forcefield (such as PCFF). The introduction of our new set of molecular parameters in this procedure cannot be easily carried out by following the procedures implemented in the MS software. However, similar structures can be run by using conveniently modified molecules. Here these temporary structural modifications consist in 12021
dx.doi.org/10.1021/jp2062332 |J. Phys. Chem. A 2011, 115, 12017–12024
The Journal of Physical Chemistry A substituting the Ru atom by hexavalent S and bonding an additional H to the carbene atom in the NHC ring. For the MD simulation of the real systems, we use the open source molecular dynamics software DL-POLY.33 This software offers total flexibility when it comes to incorporating new forcefield parameters. Moreover, we have verified that its use in Linux platforms is consistently more efficient from a computational point of view (the simulations being about 23 times faster when using the same number of parallel processors). In DL-POLY, all the forcefield parameters to be used for the different molecular mechanics elements of a given system (bonds, bond angles, torsional angles, and nonbonded interactions) should be explicitly specified one by one in different input lines. This can be a problem for systems where the number of such elements is large (in our structure several thousand torsional angles have to be defined). To circumvent this problem, we have written a bespoke Linux shell script. This script detects all the elements system that are needed to be parametrized from the output connectivity file of MS Discover (the “mdf” file). Subsequently, the script reads the corresponding parameters from a forcefield file (in our case, an extension of the PCFF forcefield file in which we have incorporated our new parameters) and writes the corresponding input “FIELD” file for DLPOLY. We have also written another Linux script that transforms the MD output file containing the system coordinates from Discover (“car” file) into the input “CONFIG” input file in the format required to run the DL-POLY simulations. Before executing these Linux scripts, we have reversed the aforementioned temporary changes in the structures that were incorporated solely for the purpose of preparing and performing the short initial MD simulation with the MS modules. These changes are limited to a couple of atoms so any structural deformations they cause are easily corrected after running a few MD simulation steps with the real structures. The MD simulations with DL-POLY are performed in three steps. First, a preliminary equilibration run of 1 ns is performed at the initial lower density, maintaining constant volume and constant room (298 K) temperature (NVT). The lower density facilitates a fast equilibration. This is followed by a further 1 ns simulation at constant pressure (NPT) which yields a system with density close to the real density for this system. Then finally a 5 ns production simulation is performed for the real density system at constant volume. In all cases we use a 1 fs simulation time-step. The analysis of the results is carried out by calculating different properties included in the trajectory data form the production run, frames or “snapshots” being saved every 25 ps.
’ ANALYSIS OF RESULTS We have analyzed different conformational characteristics of the complexes. First of all, we have investigated the results for known bond lengths and bond angles along the trajectories. Most mean values for bonds and angles remain very close to the equilibrium values introduced in the force field (obtained, as previously mentioned, from a simplified central structure). However, we observe a noticeable increment in the Rucarbene and RuP distances, both growing by approximately 0.1 Å. This growth seems to alleviate tensions introduced by the presence of the bulky branches attached to the NHC rings. It should be noted that the constants for these particular bonds in the complex are particularly small compared with typical covalent bonds, so that these greater deformations can be allowed. Table 3 contains the
ARTICLE
Table 3. Structural Results for the Complex alkyl chain length, n
radius of gyration (Å2)
RuPOSS mean distance
3
11.23 ( 0.05
12.1 ( 0.1
5
12.98 ( 0.06
14.0 ( 0.2
7
12.58 ( 0.07
14.9 ( 0.3
Figure 8. Orientation (angle α) of the carbeneP axis in the complex relative to the axis formed by the two terminal C atoms of a probe propane molecule. The trajectory frames were saved every 25 ps.
mean quadratic radius of gyration of the complex and also the mean distance between the Ru atom and the center of masses of the POSS cages. We can observe the expected increase of these quantities with the link alkyl chain length (number of carbons, n) for the n = 3 and 5 cases. However, differences are smaller when comparing the n = 5 and the n = 7 (for which a slight decrease of the radius of gyration can actually be observed). This suggests that the n = 3 link is mainly rigid, while the new bonds incorporated in the n = 5, 7 links lead to a considerably higher degree of flexibility. We have also determined the orientation between an axis associated with the central part of the complex (carbeneP line) and the main chain axis of a potentially reactive molecule (propene) introduced into the system as a probe. The results, shown in Figure 8, demonstrate random orientations along the trajectory. Therefore, the probe molecule is able to rotate freely so that, in principle, the best orientation for a reaction can easily be accessed. We have also focused on the dynamic orientation of the two normalized vectors connecting the Ru atoms with both POSS cages 1 and 2, v1 and v2, characterized by the auto- and crosscorrelation functions Ci ðτÞ ¼ Ævi ðtÞ:vi ðτ þ tÞæt =Ævi ðtÞ:vi ðtÞæt
ð7Þ
and CP ðτÞ ¼ ÆPðtÞPðτ þ tÞæt =ÆP2 ðtÞæ
ð8Þ
with PðtÞ ¼ v1 ðtÞ:v2 ðtÞ
ð9Þ
The autocorrelation functions Ci(τ) exhibit a decay time of less than 1.5 ns (see Figure 9 for n = 3). Therefore, it is shown that the 12022
dx.doi.org/10.1021/jp2062332 |J. Phys. Chem. A 2011, 115, 12017–12024
The Journal of Physical Chemistry A
ARTICLE
Table 4. Results for the Buried Volume alkyl chain length, n
Figure 9. Self-correlation of the vectors (crosses and circles) joining the Ru atom with the center of the POSS cages.
Figure 10. Cross-correlation of the vectors joining the Ru atom with the center of the POSS cages.
molecule has completely reoriented over this period, which is consistent with our use of 5 ns of simulation trajectory. The crosscorrelation (Figure 10), however, exhibits a rather dynamically rigid structure, where both vectors maintain their mutual orientations, close to antiparallel, during the simulation trajectory. The buried volume is a useful property to predict the usefulness of a metal ligand complex for catalytic purposes.34 It measures the space occupied by the ligand in the first coordination sphere of the metal center. To calculate this property, we define a sphere of radius R centered on the metal atom. We construct a tridimensional mesh of cubic elements of volume, s3, whose centers are within the sphere. We assume that all these elements are inside the sphere (the number of elements should be high enough so that differences due to the incorrect assignment of the elements containing the surface as elements inside or outside the sphere can be neglected). Next we check whether the center of each element is within the Bondi radius (third column in Table 1 of ref 34) of any atom in the ligand, thus, defining if this element is occupied or free. The buried volume is calculated as the percentage of occupied elements. Because our structures are flexible, this calculation should be averaged over all the trajectory frames.
% buried volume
% buried volume
% buried volume by
by NHC ring ligand
by NHC ring and substituents
NHC ring, phosphine, and substituents
3
27 ( 3
50 ( 3
76 ( 1
5
32 ( 2
54 ( 2
76 ( 1
7
27 ( 3
50 ( 2
75 ( 1
For our calculation, we use the value R = 3.0 Å. The results are contained in Table 4. We consider three sets of results. In the first set we only take into account the atoms attached to the NHC. These numbers can be compared with results obtained for simpler rigid NHC ligands. They show a greater occupied volume for the intermediate case, correlated with the previously discussed changes in size and flexibility. However, the longer alkyl link tends to reduce the occupied volume. It seems that the extra flexibility in these links allow for a better accommodation of the bulky POSS cages. Our data are slightly higher to those previously reported for ligands with other relatively bulky substituents with simpler chemical structures that were reported to show34 percentages of buried volume ranging from 19 to 27%, see second column in Table 2 of ref 34. Therefore, the presence of the alkyl links and POSS cages do not seem to impose further steric restrictions with respect to accessibility. The second set of results corresponds to the buried volume simultaneously due to both the NHC ring ligand and the substituents attached to Ru (two chlorine atoms and a =CH2 group). Obviously, it has a higher value but it seems to have a smaller dependence of the alkyl chain length, being slightly greater for the intermediate (n = 5) alkyl chain link. We have also evaluated the simultaneous occupied volume of the whole complex structure around the Ru atom, before the reaction takes place. This reaches a significantly high value (about 75%) and does not show any significant variation with the link chain length. Because Ru, NHC-ring carbene, phosphorus, and the two chlorine atoms adopt an approximately planar configuration, the remaining free volume should be obviously located in the plane opposite to the =CH2 substituent, see Figure 3. In summary, our methodological implementation of different computational techniques has been able to provide some important quantitative data on the conformational characteristics related with the catalytic activity of a family of complexes formed by a metal and ligands with a relatively complex chemical. Naturally, this methodology can be applied to studies of other similar systems.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected]. Present Addresses ‡
Departamento de Quimica Fisica I, Facultad de Ciencias Quimicas, Universidad Complutense de Madrid, 28040 Madrid, Spain.
’ ACKNOWLEDGMENT This work has been supported by Project CP-FP 214095-2 “HiCat” of the 7th framework programme of the European Commission. J.J.F. also acknowledges Grant CTQ2010-16414 of the MICINN (Spain), and Grant P2009/ESP-1691 of the 12023
dx.doi.org/10.1021/jp2062332 |J. Phys. Chem. A 2011, 115, 12017–12024
The Journal of Physical Chemistry A CAM (Spain). We are grateful to Dr. D. Wolf and Dr. C. Azap, Evonik Industries, for helpful discussions on the calculation of properties within the framework of the HiCat project, and Prof. S. Pricl, University of Trieste, for her useful indications on parameter calculations.
’ REFERENCES
ARTICLE
(31) Kurita, N.; Hinoue, H.; Sekino, H. Chem. Phys. Lett. 2003, 370, 161. (32) Rappe, A. K.; Goddard, W. A. J. Chem. Phys. 1991, 95, 3358. (33) Smith, W.; Yong, C. W.; Rodger, P. M. Mol. Simul. 2002, 28, 385. (34) Poater, A.; Cosenza, B.; Correa, A.; Giudice, S.; Ragone, F.; Scarano, V.; Cavallo, L. Eur. J. Inorg. Chem. 2009, 13, 1759.
(1) Ivin, K. J. ; Mol, J. C. Olefin Metathesis and Metathesis Polymerization; Academic Press: San Diego, CA, 1997. (2) F€urstner, A. Angew. Chem., Int. Ed. 2000, 39, 3012. (3) Buchmeiser, M. R. Chem. Rev. 2000, 100, 1565. (4) Trnka, T. M.; Grubbs, R. H. Acc. Chem. Res. 2001, 34, 18. (5) Rouhi, M. A. Chem. Eng. News 2002, 80, 29. (6) Nguyen, S. T.; Grubbs, R. H.; Ziller, J. W. J. Am. Chem. Soc. 1993, 115, 9858. (7) Weskamp, T.; Schattenmann, W. C.; Spiegler, M.; Herrmann, W. A. Angew. Chem., Int. Ed. 1998, 37, 2490. Huang, J.; Schanz, H.-J.; Stevens, E. D.; Nolan, S. P. Organometallics 1999, 18, 5375. (8) Ritter, T.; Hejl, A.; Wenzel, A. G.; Funk, T. W.; Grubbs, R. H. Organometallics 2006, 25, 5740. (9) Fu, G. C.; Nguyen, S. T.; Grubbs, R. H. J. Am. Chem. Soc. 1993, 115, 9856. Chatterjee, A. K.; Morgan, J. P.; Scholl, M.; Grubbs, R. H. J. Am. Chem. Soc. 2000, 122, 3783. Love, J. A.; Morgan, J. P.; Trnka, T. M.; Grubbs, R. H. Angew. Chem., Int. Ed. 2002, 41, 4035. (10) Dias, E. L.; Nguyen, S. T.; Grubbs, R. H. J. Am. Chem. Soc. 1997, 119, 3887. Sanford, M. S.; Ulman, M.; Grubbs, R. H. J. Am. Chem. Soc. 2001, 123, 749. Sanford, M. J.; Love, J. A.; Grubbs, R. H. J. Am. Chem. Soc. 2001, 123, 6543. (11) Chatterjee, A. K.; Choi, T.-L.; Sanders, D. P; Grubbs, R. H. J. Am. Chem. Soc. 2003, 125, 11360. Ibrahem, I.; Yu, M.; Schrock, R. R.; Hoveyda, A. H. J. Am. Chem. Soc. 2009, 131, 3844. Malcolmson, S. J.; Meek, S. J.; Sattely, E. S.; Schrock, R. R.; Hoveyda, A. H. Nature 2008, 456, 933. Jiang, A. J.; Zhao, Y.; Schrock, R. R.; Hoveyda, A. H. J. Am. Chem. Soc. 2009, 131, 16630. (12) Welch, C. J.; Leonard, W. R.; Henderson, D. W.; Domer, B.; Childers, K. G.; Chung, J. Y. L.; Harmer, F. W.; Albaneze-Walker, J.; Sajonz, P. Org. Process Res. Dev. 2008, 12, 81. (13) Cho, J. H.; Kim, B. M. Org. Lett. 2003, 5, 531. (14) Clavier, H.; Grela, K.; Kirschning, A.; Mauduit, M.; Nolan, S. P. Angew. Chem., Int. Ed. 2007, 46, 6786. (15) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 1996. (16) Kotelyanskii, M.; Theodorou, D. Simulation Methods for Polymers; Marcel Dekker: New York, NY, 2004. (17) Burket, U.; Allinger, N. L. Molecular Mechanics; American Chemical Society: Washington, DC, 1992. (18) Sun, H.; Mumby, S. J.; Maple, J. R.; Hagler, A. T. J. Am. Chem. Soc. 1994, 116, 2978. (19) Sun, H. J. Phys. Chem. B 1998, 102, 7338. (20) Accelrys; Accelrys Software, Inc.: San Diego, CA. (21) Perdew, J. P.; Wang, Y. Phys. Rev. B 1992, 45, 13244. (22) Delley, B. J. Chem. Phys. 2000, 113, 7756. (23) Fermeglia, M.; Ferrone, M.; Pricl, S. Fluid Phase Equilib. 2003, 210, 105. (24) Allinger, N. L.; Zhou, X.; Bergsma, J. J. Mol. Struct. 1994, 312, 69. (25) Huang, J.; Stevens, E. D.; Nolan, S. P.; Petersen, J. L. J. Am. Chem. Soc. 1999, 121, 2674. (26) Broggi, M. J.; Urbina-Blanco, C. A.; Clavier, H.; Leitgeb, A.; Slugovc, C.; Slawin, A. M. Z.; Nolan, S. P. Chem.—Eur. J. 2010, 16, 9215. (27) Maple, J. P.; Hwang, M.-J.; Stockfish, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. J. Comput. Chem. 1994, 15, 162. (28) Rosen, E. L.; Sung, D. H.; Chen, Z.; Lynch, V. M.; Bielawski, C. W. Organometallics 2010, 29, 250. (29) Geremia, S.; Calligaris, M. J. Chem. Soc., Dalton Trans. 1997, 9, 1541. (30) Hill, J.-R. J. Comput. Chem. 1997, 18, 211. 12024
dx.doi.org/10.1021/jp2062332 |J. Phys. Chem. A 2011, 115, 12017–12024
