ARTICLE pubs.acs.org/JPCB
Molecular Dynamics Study of Acid-Catalyzed Hydrolysis of Dimethyl Ether in Aqueous Solution Xiao Liang, Alejandro Montoya,* and Brian S. Haynes School of Chemical and Biomolecular Engineering, The University of Sydney, Sydney, NSW 2006, Australia ABSTRACT: The acid-catalyzed hydrolysis of dimethyl ether (DME) to methanol was examined using ab initio density functional metadynamics simulations. Diffusion of the acid proton from the aqueous medium, leading to the formation of a protonated DME, is 12.3 kcal mol1 activated and 9.3 kcal mol1 endothermic, indicating a greater affinity of the acid proton to water than to the ether group. Subsequent scission of the protonated ether bond is found to be 30.7 kcal mol1 activated, leading to the formation of a solvated methyl-carbocation, which is thermodynamically unstable. The methyl-carbocation reacts readily to form methanol and regenerate the acid proton.
1. INTRODUCTION The ether bond is a basic linkage unit, abundantly present in different biomass components and their derivatives. Hydrolyses of ethers assisted by an aqueous acid catalyst are now attracting a rapidly increasing interest because they provide insight into the decomposition pathways of biomass components into fuels.1,2 Several chemical factors seem to influence the rate of ether hydrolysis, including catalyst, solvent, the carbon chain size, configuration, and conformation.36 Dimethyl ether (DME) is the simplest alkyl ether because only one methyl group (CH3) on each side of the ether group (O) is present. The hydrolysis of DME provides useful information on the interaction of the acid proton with the ether group because the influence of the carbon chain is minimized. DME is a store and carrier of hydrogen,7 soluble in water (6.7 mol L1 at 20 °C8) and can be sustainably produced.912 Hence, it has been proposed as an alternative fuel in direct fuel cells.1317 A detailed understanding of DME hydrolysis is needed because this process constitutes an undesired side reaction that occurs under DME fuel cell operation conditions.16,18 In addition, protonation of primary alkyl ethers represents a model chemistry system that provides insight into the hydrolysis of alkoxysilanes in aqueous solution for the effective preparation of silicon particles.1921 DME hydrolysis is known to be an exceedingly slow reaction process in aqueous acid solutions even under sub- and supercritical water conditions.22 As a result, there are only limited data that permit elucidation of kinetic parameters, including some molecular modeling analyses. Static DFT calculations of DME interaction with clusters of up to five water molecules indicate that hydrophilic interaction with the ether oxygen and hydrophobic interaction with the methyl group coexist.23 Similarly, experimental and DFT computational studies on the conformations and energetics of clatharate hydrates of DME revealed that the ether oxygen of DME forms weak hydrogen bonds with solvent molecules.24,25 The interaction between ether oxygen and solvent is limited by the steric effect of methyl groups.25 A r 2011 American Chemical Society
static reaction pathway of acid-catalyzed DME hydrolysis has been reported using clusters of up to 23 water molecules.26 Inspection of the static hydrolysis configurations seems to indicate a migration of hydrogen bonds as the system moves from reactants, to transition states, to products.26 Available literature data do not provide dynamic information of the interplay between water, proton, and DME during the different stages of hydrolysis. The objective of this work is to determine the proton affinity to water relative to the ether group, the dynamics of the hydrogen bond network during hydrolysis, and, subsequently, the free energy of DME hydrolysis. In this study, we applied ab initio density functional theory molecular dynamics to study DME hydrolysis in order to obtain a robust statistical treatment of the system. Because the rate of DME hydrolysis is slow and exceedingly lengthy molecular dynamics computations would be required, we have used the CarParrinello metadynamics technique27 to explore the reaction mechanism and access the corresponding free energy landscape in a reasonable computational time. In the metadynamics approach, a repulsive biasing potential is established where the system would otherwise spend most of its time, thereby allowing the system to escape local energy minima and cross energy barriers.28 The metadynamics approach has successfully been applied to a number of reaction systems and has been proved to provide accurate thermodynamic information for liquid-phase reaction systems.2931
2. COMPUTATIONAL METHODOLOGY Density functional CarParrinello molecular dynamics simulations were carried out using the CPMD program.32 Goedecker pseudopotentials33 were used to describe the inner electronic shells of ions and the Becke, Lee, Yang, and Parr (BLYP) functional was used to describe the valence electrons.34,35 A plane-wave cutoff of 70 Ry was used for the KohnSham Received: March 1, 2011 Revised: May 6, 2011 Published: June 09, 2011 8199
dx.doi.org/10.1021/jp201951a | J. Phys. Chem. B 2011, 115, 8199–8206
The Journal of Physical Chemistry B
ARTICLE
orbitals. The equations of motion were integrated with a time step of 4 au (0.1 fs), and a fictitious mass of 800 amu for electrons was used. The CPMD simulations were conducted under the NVT ensemble at 300 K with a NoseHoover chain thermostat.36 A periodic unit cell with cubic dimensions of 10 Å was used. Similar unit cell sizes have been used in previous CPMD studies of acidcatalyzed hydrolysis of methylacetamide37 and methyl formate.38 Our simulated unit cell contains 30 water molecules, 1 DME molecule, 1 hydronium ion (H3Oþ), and 1 chloride ion. The chloride ion, acting as a counterion, was used to maintain electric neutrality of the system. During the simulations, the chloride anion is observed to remain always in the second solvation shell of DME without interacting with H3Oþ, nor participating in the reaction. The involvement of a chloride counterion to neutralize the simulated system has been used in several studies of acidcatalyzed reactions.37,39,40 The corresponding density of the solution is 1030 g L1, and the concentrations of DME and H3Oþ are both 1.66 mol L1. The initial configuration was obtained by inserting one optimized isolated DME molecule, one hydrogen, and one chloride atom into a box containing water molecules in equilibrium and deleting the water molecules that are overlapped by the DME molecule. The DME molecule was optimized by using Gaussian 0941 at the level of B3LYP/6-311þþG(d,p), while the equilibrium configuration of water molecules was achieved from a Monte Carlo simulation using the Dice program.42 A 2.5 ps Car Parrinello molecular dynamics calculation was then conducted on the created supercell to equilibrate the system (water þ DME þ HCl) before performing the metadynamics simulations. We studied the acid-catalyzed hydrolysis of DME in two steps. Step 1 involves protonation of DME and subsequent formation of a methyl-carbocation, methanol, and water, and step 2 involves formation of a second methanol molecule and regeneration of the acid proton species: CH3 OCH3 þ H3 Oþ f CH3 OH þ CH3 þ þ H2 O ðstep1Þ CH3 þ þ 2H2 O f CH3 OH þ H3 Oþ
ðstep2Þ
Potential free energy surfaces for steps 1 and 2 were generated from metadynamics simulations43 by biasing the dynamics with a history-dependent potential. This biasing potential can act as a coarse-grained representation of the system, providing free energies that depend on a set of predefined collective variables (CVs) that characterize the reaction coordinate. In this study, two CVs were selected for each step, as presented in Figure 1. CV1 was defined as the difference in the coordination number (CN) between C1 and O1 and between C1 and O2, whereas CV2 was described by the CN changes of H1 with respect to O1 (Figure 1). In step 2, the CN changes of C1 with respect to O2 and CN changes of H2 with respect to O2 were used. In all cases, the CN values were defined by the following equation44 CN ¼
1 ðdij =d0 Þp 1 ðdij =d0 Þp þ q
ðE1Þ
where dij is the distance between atoms i and j, d0 is the reference distance, and p = q = 6 are constants to distinguish between the coordinated and noncoordinated states. The CN value between two atoms indicates whether a covalent bond exits (1 = bond, 0 = no bond). The selected values of cutoff distance d0 for all selected
Figure 1. Selected collective variables in the acid-catalyzed hydrolysis of DME. Step 1: scission of the O1C1 ether bond and formation of the O2C1 bond (CV1) and protonation of DME (CV2). Step 2: formation of methanol (CV3) and regeneration of hydronium ion (CV4).
Table 1. Control Parameters of the Selected Collective Variables d0 (Å)
m (amu)
k (au)
CV1
2.0
600
8.0
CV2
1.5
100
3.0
CV3
2.0
600
8.0
CV4
1.5
100
3.0
CVs are shown in Table 1 coupled with the parameters of fictitious masses (m) and force constants (k), which controls the dynamics of the fictitious CVs. A normal Gaussian bias potential was used in our metadynamics simulations. The width of the Gaussian potential was set to 0.05 au in all simulations. The heights of the Gaussian hills in the protonation and dissociation of DME were determined by the shape of the underlying energy surface and set to be, at most, 6.2 kcal mol1 for the first 200 metadynamics steps and reduced to 1.8 kcal mol1 in the subsequent metadynamics steps. The average hill height in the simulation of step 1 was 1.5 kcal mol1. The heights of the Gaussian hills in the proton regeneration were set at 0.3 kcal mol1. The addition of bias potential was allowed with a minimum time separation of 100 MD steps and with the displacements in CVs being larger than 1.5 times the width of Gaussian function. The addition of consecutive potentials was forced after 300 MD steps if the displacement had not exceeded the threshold. Similar settings in previous studies were shown to be able to accurately predict the energetics and dynamics of hydrolysis of esters and sugars in aqueous solution.38,45 The molecular graphics were rendered with VMD.46 Inside our unit cell, DME is surrounded by two water shells, which hinders interaction of the DME with DME molecules located in adjacent periodic unit cells. Because DME hydrolysis produces two methanol species located inside our finite unit cell, the metadynamics free energy maps inevitably contain the methanolmethanol interaction with the proton species. We have estimated the required energy for infinite separation of the reactants and products by performing static electronic calculations on the reactant and products. Static calculations were carried out in Gaussian 0941 for optimizing the structures and estimating the solvation free energies of reactant species. The geometry optimizations were computed at the B3LYP/ 6-311þþG(d,p) level of theory. The nature of the structures at equilibrium was determined by performing a normal-mode 8200
dx.doi.org/10.1021/jp201951a |J. Phys. Chem. B 2011, 115, 8199–8206
The Journal of Physical Chemistry B
ARTICLE
Figure 2. Two-dimensional free energy landscape of the reaction, CH3OCH3 þ H3Oþ f CH3OH þ H3Cþ þ H2O, in the aqueous phase at 300 K (energy in kcal mol1). The energy scale is given in the inset to the right. Positions R, TS, and PI in the energy landscape correspond to selected key conformations of the reactant, transition state, and product basins, which are presented in Figure 3.
analysis. Each minimum is characterized by possessing all real vibrational frequencies with its Hessian matrix possessing all positive eigenvalues. The unscaled frequencies were used to obtain corrections for Gibbs free energies. Calculated gas-phase Gibbs free energies in a 1 atm ideal gas-phase standard state were converted to the molar Gibbs free energies in the gas phase via the relationship, ΔGg,1 M = ΔGg,1 atm þ RT ln(Q1 M/Q1 atm), where the term RT ln(Q1 M/Q1 atm) is 1.9 kcal mol1 at 298.15 K. Gibbs free energies of solvation in water solution (ΔsG) were calculated by means of the SMD model47 on the basis of gasphase geometry optimized at the same level of B3LYP/ 6-311þþG(d,p).
3. RESULTS AND DISCUSSION 3.1. Protonation and Dissociation of DME. Figure 2 shows the free energy surface for protonation and dissociation of DME as a function of CV1 (CO scission) and CV2 (proton transfer). Figure 3 shows selected snapshots collected from the metadynamics simulations corresponding to the energy surface in Figure 2. The average bond lengths and their corresponding standard deviations in every energy basin are presented for selected bonds in Table 2. The data are computed from all the configurations falling into the region where the free energy is within ∼1.5 kcal mol1 of each minimum. As is apparent from Figure 2, there are five states (R1R5) in the reactant basin in which the value of CV1 is between 0.8 and 1.0. The minimal energy reaction path begins from the global minimum R1, which is located at CV1 = 0.85 and CV2 = 0.01. The small value of CV2 in state R1 indicates a very weak interaction between the solvent and the DME molecule. Inspection of the R1 configuration in Figure 3 shows that the proton tends to form an H5O2þ species that is positioned in the second solvation shell of DME and solvated by water molecules exclusively. This indicates that the proton affinity to water is stronger
than to the ether oxygen. The average bond length of C1O1 in R1 is 1.49 Å with a small variability of 0.10 Å. The calculated bond length is slightly longer than the experimental CO bond length of 1.42 Å for DME in the gas phase.48 The average length H1 3 3 3 O1 at R1 is 3.13 Å, suggesting a weak hydrogen bond interaction. Meanwhile, H1 also donates another hydrogen bond toward the neighborhood water molecule with the length of 1.84 Å, consistent with observations of clatharate hydrates of DME.25 The hydrolysis reaction proceeds by diffusion of the proton via the hydrogen bond network toward DME, which leads to the formation of intermediates R2R5. First, water molecules start to interact with DME by crossing a small free energy barrier (8.4 kcal mol1) at TS1 (CV1 = 0.85 and CV2 = 0.09) and reaching a new reactant intermediate R2 (CV1 = 0.85 and CV2 = 0.19). The free energy at R2 is 4.5 kcal mol1 higher than that of R1. The protonated water species in configuration R2 remains in the second solvation shell of DME and interacts with DME via four hydrogen bonds (as seen in Figure 3). The average C1O1 bond length in R2 is 1.50 Å, very close to the length of this bond in R1. The value of CV2 increases as the proton further approaches DME. This process results in a reactant basin, R3 at CV1 = 0.86 and CV2 = 0.52. The free energy at R3 is 7.7 kcal mol1 higher than R1, and 3.2 kcal mol1 higher than R2. The corresponding energy barrier, at TS2, is 8.8 kcal mol1 with respect to the energy of R1. Inspection of the R3 configuration shows that the proton has moved to the first solvation shell, and the number of hydrogen bonds between protonated water species and DME has decreased to two. The relative changes of the ether bond length remain very small with respect to R1 and R2. Therefore, the activation energy for this step is interpreted as the diffusion energy barrier of the proton from the second to the first water shell due to a rather stronger attraction by water than by DME. Following the minimum free energy pathway in Figure 2, it is observed that the proton is next transferred to a point near O1, as shown in Figure 3, intermediate R4. State R4 is located at 8201
dx.doi.org/10.1021/jp201951a |J. Phys. Chem. B 2011, 115, 8199–8206
The Journal of Physical Chemistry B
ARTICLE
Figure 3. Snapshots from metadynamics simulations of the CH3OCH3 þ H3Oþ f CH3OH þ H3Cþ þ H2O reaction.
Table 2. Average Bond Lengths and Corresponding Standard Deviations (in Å) of Selected Configurations along the Reaction Pathway CH3OCH3 þ H3Oþ f CH3OH þ CH3þ þ H2 Oa C1O1
a
C1O2
H1O1
R1
1.49 ( 0.10
4.65 ( 0.47
3.13 ( 0.42
R2
1.50 ( 0.10
4.90 ( 0.83
1.78 ( 0.19
R3
1.50 ( 0.08
4.64 ( 0.88
1.52 ( 0.14
R4
1.54 ( 0.09
4.58 ( 0.67
1.34 ( 0.11
R5
1.53 ( 0.09
5.16 ( 0.89
1.19 ( 0.10
TS5 PI1
1.93 ( 0.14 2.21 ( 0.16
2.21 ( 0.11 1.86 ( 0.21
1.06 ( 0.03 1.01 ( 0.02
Atom nomenclature is presented in Figure 1.
CV1 = 0.80 and CV2 = 0.65, 9.4 kcal mol1 less stable than R1. The energy barrier for proton diffusion toward O1 in the first solvation shell is 10.5 kcal mol1 with respect to the energy of R1. There is just one hydrogen bond between the protonated water species and DME in R4. The H1 3 3 3 O1 distance is 1.34 Å and the C1O1 bond length increases to 1.54 Å due to the proton interaction with O1. The hydrolysis continues as CV2 increases, forming protonated DME, as observed in Figure 3, intermediate R5. Here, CV1 is 0.84 and CV2 is 0.79. The free energy of R5 is 7.5 kcal mol1 higher than that of R1, but lower than for R4 by a few kcal mol1. It is noteworthy that the energy gradient of R5 is significant along
the hydrolysis reaction coordinate, meaning that any small progress toward the dissociation of the ether bond requires a considerable energy input. For example, the free energy change is 6.1 kcal mol1 as the value of CV1 decreases from 0.84 to 0.76 and the CV2 increases from 0.79 to 0.88. The C1O1 bond length in R5 is 1.53 Å, which is 0.04 Å longer than that in R1. The distance between H1 and O1 decreases to 1.19 Å in R5, indicating that a proton has been transferred to DME; the transition state TS4 connecting R4 and R5 is at CV1 = 0.82 and CV2 = 0.69, with a free energy 10.4 kcal mol1 above R1. Dissociation of protonated DME occurs between the potential wells R5 and PI1 (Figure 2). The maximum value of the free energy along the minimum reaction pathway appears at CV1 = 0.13 and CV2 = 0.89, corresponding to TS5, 30.7 kcal mol1 above R5 or 38.2 kcal mol1 above R1. The snapshot of TS5 in Figure 3 shows a nascent methanol molecule and a methylcarbocation species. The methyl-carbocation configuration is nearly flat. The ether bond C1O1 has been stretched to 1.93 Å in TS5, which is 29.5% longer than the equilibrium bond length in DME. Coupled with the scission of the ether bond and the formation of the methyl-carbocation, a water molecule is observed to diffuse from the solvent to participate in the hydrolysis reaction. The reacting water molecule adsorbs to the methylcarbocation with a C1 3 3 3 O2 distance of 2.21 Å, thereby assisting in distributing the positive charge. The configuration of TS5 clearly indicates that acid hydrolysis of DME in aqueous solution follows an A2 bimolecular nucleophilic substitution mechanism, consistent with previous indications.5,49 8202
dx.doi.org/10.1021/jp201951a |J. Phys. Chem. B 2011, 115, 8199–8206
The Journal of Physical Chemistry B After passing through TS5, the system gains 1.1 kcal mol1 of energy and reaches a product intermediate potential well, PI1. The distance between C1 and O2 decreases to 1.86 Å, and a second adjacent H2O molecule is found to be tilted toward the methyl-carbocation, preparing to regenerate the acid proton and form the second methanol species. The stabilization provided by the two H2O molecules facilitates further separation of the methyl-carbocation from methanol; that is, the distance between O1 and C1 increases to 2.21 Å. 3.2. Regeneration of Proton Species. PI1 was selected as the initial configuration for a second metadynamics computation to
ARTICLE
simulate the formation of methanol and regeneration of the acid proton as a function of the collective variables (CV3 and CV4) presented in Figure 1, step 2. The free energy surface is presented in Figure 4, and selected snapshots along the reaction coordinate are presented in Figure 5. The average lengths and corresponding standard deviations of selected bonds in the system at the bottom of the various energy basins are presented in Table 3. The data are computed from all the configurations falling into the region where the free energy is within ∼1.5 kcal mol1 of each minimum. The lowest energy of the reactant basin is centered at CV3 = 0.61 and CV4 = 0.89, corresponding to the PI2 state. The nomenclatures PI1 and PI2 are used to distinguish the same intermediate state along the hydrolysis reaction in two respective simulations that correspond to different collective variables. By manually converting the average bond lengths of C1O2 and H2O2 of PI1 into CV3 and CV4, the energy of PI1 can be located at the free energy surface generated from the second metadynamics simulation, as the red square presented in Figure 4. The energy of PI1 is observed to be 0.3 kcal mol1 less stable than that of PI2. This slight energy difference is probably attributable to the different heights of the bias Gaussian function and collective variables selected in two calculations. The smaller average height of Gaussian hills and more targeted collective Table 3. Average Bond Lengths and Corresponding Standard Deviations (in Å) of Selected Configurations along the Reaction Pathway CH3þ þ 2H2O f CH3OH þ H3Oþa C1O2
Figure 4. Free energy surface of H3Cþ þ 2H2O f CH3OH þ H3Oþ in the aqueous phase at 300 K (energy in kcal mol1). Insets PI and TS correspond to the energy of a selected configuration of the product and transition state basins, which are presented in Figure 5.
a
H2O2
PI2
1.84 ( 0.06
1.05 ( 0.04
TS6 PI3
1.70 ( 0.07 1.53 ( 0.06
1.11 ( 0.02 1.12 ( 0.08
PI4
1.51 ( 0.08
1.88 ( 0.15
PI5
1.48 ( 0.11
4.26 ( 0.69
Atom nomenclature is presented in Figure 1.
Figure 5. Snapshots from metadynamics simulations of the H3Cþ þ 2H2O f CH3OH þ H3Oþ reaction. 8203
dx.doi.org/10.1021/jp201951a |J. Phys. Chem. B 2011, 115, 8199–8206
The Journal of Physical Chemistry B variables provide a more complete sampling of this intermediate state in the second simulation; however, the comparable energy and geometrical parameters of PI1 and PI2 indicates good consistency and energy convergence between the two separate simulations. The reaction continues as the bond length of C1O2 decreases, passing through a saddle point TS6 and reaching a product basin PI3 at CV3 = 0.84 and CV4 = 0.84. Inspection of PI3 in Figure 5 shows clearly that a protonated methanol is generated. The bond length of C1O2 decreases from 1.84 Å in PI2 to 1.53 Å in PI3. The protonated methanol is not thermodynamically stable as the proton can diffuse almost without a barrier (0.5 kcal mol1) back to the solvent medium to which it has a greater affinity. At product basin PI4, two methanol molecules and an acid proton are generated. PI4 is located at CV3 = 0.85 and CV4 = 0.20, which is 5.9 kcal mol1 more stable than PI3. Inspection of PI4 in Figure 5 shows that the proton is stabilized in the form of H5O2þ in the original DME solvation shell. The bond lengths of H2O2 and C1O2 are 1.88 and 1.51 Å, respectively. The final product, PI5, further separates the two methanol molecules and proton. The equilibrium bond length of C1O2
Figure 6. Optimized structure of DME 3 (H9O4þ) in the gas phase at B3LYP/6-311þþG(d,p). The selected bond length between DME and H2O is presented.
ARTICLE
in methanol is 1.48 Å with a standard deviation of 0.11 Å, consistent with reported equilibrium geometrical values of hydrated methanol using CPMD molecular dynamics.26,50 The snapshot of PI5 in Figure 5 indicates that the proton is connected to methanol via four hydrogen bonds, similar to that in the case of the DME configuration in R1. The direct interaction between methanol molecules is minimized in the periodic supercell. The global free energy of Rxn-1 calculated from available Gibbs free energies of formation in aqueous phase at 298.15 K and 1 M standard state of DME (63.0 kcal mol1),51,52 water (75.6 kcal mol1),52,53 and methanol (68.3 kcal mol1)51,52 is ΔGinf = 2.0 kcal mol1. The calculated free energy from our metadynamics calculation is 18.3 kcal mol1. CH3 OCH3 þ H2 O f CH3 OH þ CH3 OH
ðRxn-1Þ
As an estimate of the energy required bringing DME from R1 in Figure 3 to infinite dilution, we have computed from static electronic structure calculations at the B3LYP/6-311þþG(d,p) level the Gibbs free energy in the aqueous phase of Rxn-2 DME þ 3H2 O þ H3 Oþ f DME 3 ðH9 O4 þ Þ
ðRxn-2Þ
The optimized geometry in the gas phase of DME 3 (H9O4þ) is presented in Figure 6. The equilibrium geometry of DME 3 (H9O4þ) agrees well with the configuration of DME in the energy basin R3. The solvation free energy (ΔGsol) of DME 3 (H9O4þ) using the continuum solvation model SMD45 was computed to be 54.2 kcal mol1. Using the experimental values of ΔGsol of DME (1.9 kcal mol1),52 H2O (6.3 kcal mol1),52 and H3Oþ (110.2 kcal mol1),54 we computed the aqueous phase reaction Gibbs free energy of Rxn-2 to be 9.5 kcal mol1. Because the energy difference between R1 and R3 is 7.7 kcal mol1, the infinite separation of R1 is deduced to be 1.8 kcal mol1. This small energy for infinite separation seems to be consistent with the weak interaction of the proton with DME in R1. Therefore, the infinite separation of the final product PI5 is deduced to be 18.3 þ 1.8 2.0 = 18.1 kcal mol1.
Figure 7. Free energy at 300 K as a function of acid-catalyzed hydrolysis of DME. The nomenclature in the x axis corresponds to those used in Figures 25. R and P correspond to the reactants and products at infinite separation. Broken lines (-\\-) are used to connect two energy minima for which no TS is provided. 8204
dx.doi.org/10.1021/jp201951a |J. Phys. Chem. B 2011, 115, 8199–8206
The Journal of Physical Chemistry B Figure 7 summarizes the minimum free energy pathway of DME hydrolysis. The pathway includes the infinite separation of R1 and PI5. It is observed that diffusion of the acid proton from infinity to the ether oxygen is ∼12.3 kcal mol1 activated (R f TS3) and 9.3 kcal mol1 endothermic (R f R5). The scission of the protonated ether bond (R5 to PI1 step in Figure 7) is the step with the highest activation requirement (30.7 kcal mol1) and the proton regeneration (PI1 to PI5) is nearly barrierless, consistent with mechanistic descriptions of hydrolysis of alkyl ethers where the formation of a carbocation species is ratelimiting and regeneration of the proton species is a fast reaction step.3 The global free energy barrier is 40.0 kcal mol1. The large energy barrier is consistent with the slow conversion rates of DME to methanol under strong acid-catalyzed conditions.3,49 The high energy barrier for DME hydrolysis is mainly associated with the formation of a methyl-carbocation that, despite being readily solvated, is highly unstable, possibly due to a large positive charge concentration on the carbon atom. Our free activation energy compares relatively well with a reported 37.8 kcal mol1 (electronic) energy barrier for the hydrolysis of DME, which is inside a protonated water cluster of 23 water molecules.26 The free energy is comparable to the energy barrier due to a small entropy contribution because there are no drastic geometrical changes of the cluster and the number of water molecules remains constant. An electronic energy of 7.6 kcal mol1 is reported to be required to separate the solvated reactant cluster to a solvated H3Oþ and solvated DME to infinite separation.26 Because the process of infinite separation changes the number of molecules from reactant to products, this process is better represented by the release of 1.8 kcal mol1 free energy, as presented here.
4. CONCLUSION Metadynamics simulations were carried out to determine the pathways and free energy surface for the acid-catalyzed hydrolysis of DME in the aqueous phase. The energy required to bring a proton from infinity to the vicinity of the DME molecule is 11.2 kcal mol1. The endothermicity is associated with the greater proton attraction by the bulk solvent than by the DME molecule. Subsequent formation of the protonated DME is 1.0 kcal mol1 activated and 1.9 kcal mol1 exothermic. The rate of the protonation of DME is controlled by the slow step of proton transportation toward DME. The scission of the CO ether bond is the overall rate-limiting step, which follows an A2 type of reaction mechanism. The free energy of activation for the acid hydrolysis of DME is 40.0 kcal mol1. ’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT This work has been funded by the Australian Research Council through Discovery Grant DP1096802. The authors acknowledge the Australian National Computational Infrastructure (NCI) for the computational time allocated to this project.
ARTICLE
’ REFERENCES (1) Polydore, C.; Roundhill, D. M.; Liu, H.-Q. J. Mol. Catal. A: Chem. 2002, 186, 65. (2) Yu, Y.; Lou, X.; Wu, H. Energy Fuels 2008, 22, 46. (3) Burwell, R. L., Jr. Chem. Rev. (Washington, DC, U.S.) 1954, 54, 615. (4) Whalley, E. Trans. Faraday Soc. 1959, 55, 798. (5) Hughes, E. D.; Ingold, C. K. J. Chem. Soc. 1935, 244. (6) O’Reilly, K. T.; Moir, M. E.; Taylor, C. D.; Smith, C. A.; Hyman, M. R. Environ. Sci. Technol. 2001, 35, 3954. (7) Takeishi, K. Recent Advances in Energy & Environment, Proceedings of the 4th IASME/WSEAS International Conference on ENERGY and ENVIRONMENT, Cambridge, U.K., Feb 2426, 2009; WSEAS Press, 2009; p 449. (8) Yoo, J.-H.; Choi, H.-G.; Chung, C.-H.; Cho, S. M. J. Power Sources 2006, 163, 103. (9) Xu, M.; Lunsford, J. H.; Goodman, D. W.; Bhattacharyya, A. Appl. Catal., A 1997, 149, 289. (10) Khandan, N.; Kazemeini, M.; Aghaziarati, M. Appl. Catal., A 2008, 349, 6. (11) Lee, S.; Sardesai, A. Top. Catal. 2005, 32, 197. (12) Taniewski, M. Clean: Soil, Air, Water 2008, 36, 393. (13) Liu, Y.; Mitsushima, S.; Ota, K.-I.; Kamiya, N. Electrochim. Acta 2006, 51, 6503. (14) Mench, M. M.; Chance, H. M.; Wang, C. Y. J. Electrochem. Soc. 2004, 151, A144. (15) Yu, J.-H.; Choi, H.-G.; Cho, S. M. Electrochem. Commun. 2005, 7, 1385. (16) Serov, A.; Kwak, C. Appl. Catal., B 2009, 91, 1. (17) Muller, J. T.; Urban, P. M.; Holderich, W. F.; Colbow, K. M.; Zhang, J.; Wilkinson, D. P. J. Electrochem. Soc. 2000, 147, 4058. (18) Mizutani, I.; Liu, Y.; Mitsushima, S.; Ota, K.-i.; Kamiya, N. J. Power Sources 2006, 156, 183. (19) Stoeber, W.; Fink, A.; Bohn, E. J. Colloid Interface Sci. 1968, 26, 62. (20) El-Nahhal, I. M.; El-Ashgar, N. M. J. Organomet. Chem. 2007, 692, 2861. (21) Brochier Salon, M.-C.; Belgacem, M. N. Colloids Surf., A 2010, 366, 147. (22) Nagai, Y.; Matubayasi, N.; Nakahara, M. J. Phys. Chem. A 2005, 109, 3558. (23) Zeng, X.; Yang, X. J. Phys. Chem. B 2004, 108, 17384. (24) Monreal, I. A.; Cwiklik, L.; Jagoda-Cwiklik, B.; Devlin, J. P. J. Phys. Chem. Lett. 2010, 1, 290. (25) Kulig, W.; Kubisiak, P.; Cwiklik, L. J. Phys. Chem. A [Online early access]. DOI: 10.1021/jp111245z. Published Online: February 7, 2011. (26) Terleczky, P.; Nyulaszi, L. J. Phys. Chem. A 2009, 113, 1096. (27) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Chem. Phys. 2002, 116, 6941. (28) Laio, A.; Gervasio, F. L. Rep. Prog. Phys. 2008, 71, 126601. (29) Park, J. M.; Laio, A.; Iannuzzi, M.; Parrinello, M. J. Am. Chem. Soc. 2006, 128, 11318. (30) Cucinotta, C. S.; Ruini, A.; Catellani, A.; Stirling, A. ChemPhysChem 2006, 7, 1229. (31) Schreiner, E.; Nair, N. N.; Marx, D. J. Am. Chem. Soc. 2008, 130, 2768. (32) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471. (33) Goedecker, S.; Teter, M.; Hutter, J. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703. (34) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (35) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter. Mater. Phys. 1988, 37, 785. (36) Nose, S. J. Chem. Phys. 1984, 81, 511. (37) Zahn, D. J. Phys. Chem. B 2003, 107, 12303. (38) Gunaydin, H.; Houk, K. N. J. Am. Chem. Soc. 2008, 130, 15232. (39) Dong, H.; Nimlos, M. R.; Himmel, M. E.; Johnson, D. K.; Qian, X. J. Phys. Chem. A 2009, 113, 8577. 8205
dx.doi.org/10.1021/jp201951a |J. Phys. Chem. B 2011, 115, 8199–8206
The Journal of Physical Chemistry B
ARTICLE
(40) Liu, D.; Nimlos, M. R.; Johnson, D. K.; Himmel, M. E.; Qian, X. J. Phys. Chem. A 2010, 114, 12936. (41) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, € Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; S.; Daniels, A. D.; Farkas, O.; Fox, D. J. Gaussian 09, revision A1; Gaussian, Inc.: Wallingford, CT, 2009. (42) Coutinho, K.; Canuto, S. Dice, version 2.9; University of S~ao Paulo: Brazil, 2003. (43) Ensing, B.; De Vivo, M.; Liu, Z.; Moore, P.; Klein, M. L. Acc. Chem. Res. 2006, 39, 73. (44) Iannuzzi, M.; Laio, A.; Parrinello, M. Phys. Rev. Lett. 2003, 90, 238302. (45) Petersen, L.; Ardevol, A.; Rovira, C.; Reilly, P. J. J. Phys. Chem. B 2009, 113, 7331. (46) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33. (47) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2009, 113, 6378. (48) Tamagawa, K.; Takemura, M.; Konaka, S.; Kimura, M. J. Mol. Struct. 1984, 125, 131. (49) Long, F. A.; Paul, M. A. Chem. Rev. (Washington, DC, U.S.) 1957, 57, 935. (50) van Erp, T. S.; Meijer, E. J. Chem. Phys. Lett. 2001, 333, 290. (51) Goos, E.; Burcat, A.; Ruscic, B. Third Millennium Ideal Gas and Condensed Phase Thermochemical Database for Combustion with Updates from Active Thermochemical Tables 9; Private communication, Dec 2010.
[email protected]. (52) Cramer, C. J.; Truhlar, D. G. J. Comput.-Aided Mol. Des. 1992, 6, 629. (53) Linstrom, P. J., Mallard, W. G., Eds.; NIST Chemistry WebBook, NIST Standard Reference Database Number 69; National Institute of Standards and Technology: Gaithersburg, MD. http://webbook.nist. gov. (accessed Dec 9, 2010). (54) Pliego, J. R., Jr.; Riveros, J. M. Phys. Chem. Chem. Phys. 2002, 4, 1622.
8206
dx.doi.org/10.1021/jp201951a |J. Phys. Chem. B 2011, 115, 8199–8206