Acid-Catalyzed Ring Opening of Furan in Aqueous Solution - Energy

Article pubs.acs.org/EF

Cite This: Energy Fuels 2018, 32, 4139−4148

Acid-Catalyzed Ring Opening of Furan in Aqueous Solution Xiao Liang, Brian S. Haynes, and Alejandro Montoya* School of Chemical and Biomolecular Engineering, The University of Sydney, Sydney, New South Wales 2006, Australia S Supporting Information *

ABSTRACT: The mechanism and energetics of furan decomposition via the opening of the five-membered ring structure in dilute acid solution were investigated using density functional Car−Parrinello molecular dynamics combined with metadynamics simulations. Diffusion of an acidic proton from the aqueous medium leading to the formation of a protonated furan is found to be the rate-limiting step of the ring-opening process, with protonation at the Cα position being 7 kcal mol−1 less activated than that at the Cβ position. Protonation at the Cα position results in the formation of 2,5-dihydro-2-furanol or 2,3-dihydro-3-furanol via nucleophilic attack by a solvent molecule. Subsequent protonation of the furanols at the ring oxygen initiates the opening of the furanic ring, leading to the formation of 4-hydroxy-2-butenal. in mildly acidic solution23 and 2-methylfuran in dioxane−D2O solution with dilute HClO4 acid22 supported the Cαprotonation pathway. Nevertheless, kinetic studies22,23 suggest that the pre-equilibrium of Cα protonation is faster than the subsequent cleavage of the furanic ring, conflicting with observations by Stamhuis et al. The activation energy for α-hydrogen exchange of 2methylfuran in a water−dimethyl sulfoxide mixture with dilute HClO4 between 70 and 100 °C was determined to be 21.7− 24.1 kcal mol−1.16 In contrast, no hydrogen exchange between solvent and Cβ was observed.16,23 It is proposed that the Cα positions are more readily protonated in both gas phase and aqueous solution,24,25 mainly as a result of a higher stability of the Cα-furanium ion, which has greater charge delocalization than its isomer Cβ-furanium ion.26 The ring oxygen has lower proton affinity than the Cα and Cβ positions.27 It has been shown that protonation at the ring oxygen only slightly alters the electron density distribution of the furanic ring, thereby reducing the possibility of a subsequent nucleophilic attack.28,29 While furan is detected only in a low concentration during the decomposition of monosaccharides in acid solutions, its ring configuration provides a framework to study the furanic ring-opening reaction without the interference of formyl and hydroxymethyl groups associated with the more prevalent furanic products. Therefore, furan has been selected as the model furanic compound in molecular modeling studies. In the gas phase, the relevant stabilities of furanium ions were computed at the G3 level of theory, showing that the most favored protonation position is Cα atoms, which is 10.5 and 27.7 kcal mol−1 more exothermic than that at Cβ and heterooxygen atoms, respectively.30 The kinetics of Cβ protonation of furan were separately studied in the gas phase at the B3LYP/6311G(d,p) level of theory,31 yielding the result that the proton

1. INTRODUCTION Furans represent a family of chemicals that can be derived from cellulose and hemicellulose. Furfural and 5-hydroxymethylfurfural (5-HMF) are valuable furanic compounds obtained from renewable resources because they are essential building blocks for other chemicals with wide applications in the automotive and chemical industry.1,2 These compounds are produced by dehydration of monosaccharides, such as glucose, fructose, and xylose. Furans are susceptible to decomposition via opening of their five-membered ring structures.3−6 The decomposition reaction pathways become predominant, especially in the presence of acid catalysts.7−10 Therefore, elucidation of the reaction mechanisms and kinetics of the furanic ring opening in acidic aqueous solution is fundamental to gaining insight into the process of sustainable production of chemicals from biomass. The stability of furan and its derivatives in acidic solutions has been studied under various acid concentrations and temperatures. The stability of the furanic ring depends significantly upon the nature of its substituents. Furanics with electron-attracting substituent groups are rather stable in acidic solution,11 while those with electron-donating groups are readily converted into their corresponding 1,4-diketones.12,13 In hot dilute aqueous acids, it has been accepted that the major reaction of furanics is hydrolytic ring opening that starts from protonation of one of the ring carbon atoms.14 The selectivity of the protonation site of a furanic ring has been debated in the literature. A rate-limiting protonation at the Cβ positions was proposed on the basis of several experimental studies of acid-catalyzed ring cleavage of furan and 2-methoxy-, 2-methyl-, and 2,5-dimethylfurans.15−17 Although the product distribution was subsequently revised,18 the Cβ protonation is used to explain the formation of levulinic acid and polymerization from 5-HMF.19,20 On the other hand, Stamhuis et al. have proposed that protonation at the Cα position is ratelimiting among the ring-opening reactions of furan and 2,5dimethylfuran in aqueous HClO4 solution at 25 °C.21 Additional kinetic isotope experiments confirmed that hydrogen exchange between solvent and Cα follows a two-step A− SE2 mechanism.16,22 Two separate experimental studies of furan © 2017 American Chemical Society

Special Issue: 6th Sino-Australian Symposium on Advanced Coal and Biomass Utilisation Technologies Received: October 22, 2017 Revised: December 20, 2017 Published: December 20, 2017 4139

DOI: 10.1021/acs.energyfuels.7b03239 Energy Fuels 2018, 32, 4139−4148

Article

Energy & Fuels Scheme 1. Possible Reaction Pathways of Acid-Catalyzed Ring Opening of Furan

furanol, as shown in Scheme 1. The second step comprises the acid-catalyzed ring opening via C−O bond dissociation of the furanol species to acyclic products.

transfer from an adjacent hydronium ion to the Cβ position was 5.3 kcal mol−1 exothermic and nearly barrierless (0.2 kcal mol−1).31 The energy barrier for the proton shift from the α to β position, assisted by one water molecule, was 16.7 kcal mol−1.31 The reaction mechanism and energetics of Brønstedacid-catalyzed hydrolysis of dimethylfuran in aqueous solution have been recently explored using electronic structure calculations at the M062X/6-311+G(d,p) theory of level using an implicit water solvent model.32 Although the activation energy of α protonation was computed to be 3 kcal mol−1 lower than that of β protonation, it was proposed that the favored reaction pathway follows a general acid catalysis mechanism by rate-limiting Cβ protonation. This mechanism was proposed as a result of the absence of an energy barrier in the water binding reaction of β-furanium ion and the low energy barrier (9 kcal mol−1) of the subsequent ring-opening step compared to a much higher activation energy of 27 kcal mol−1 for the ring opening of α-furanium ion.32 A broad energy mapping of possible reaction pathways initiating at the Cβ and Cα positions of a furanic ring in solution using explicit water molecules is required to reveal the important dynamic effect of water molecules in the protonation of furanics and, therefore, to obtain additional information on the rate-limiting step of the ring-opening process. In this work, the mechanisms and energetics of the acid-catalyzed ring cleavage of furan are explored using density functional Car− Parrinello molecular dynamics (CPMD) combined with metadynamics. This computational approach has been successfully applied to some organic reaction systems and has been proven to provide accurate energetic information for aqueous-phase reaction systems.33−37 The ring opening of furan will be explored in two steps. The first step involves the protonation of furan and hydration of the furanium species. Protonation at the α position leads to the formation of 2,5dihydro-2-furanol and 2,3-dihydro-3-furanol, while protonation at the β position leads to the formation of 2,3-dihydro-2-

2. COMPUTATIONAL METHODOLOGY The molecular dynamics (MD) simulations were carried out using the CPMD program.38 The Goedecker pseudo-potentials39 and the Becke, Lee, Yang, and Parr (BLYP) functional40,41 were used to describe the potentials exerted by the core electronic shells of ions and the valence electrons, respectively. A plane wave cutoff of 70 Ry was used for the Kohn−Sham 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. A NVT ensemble was used at 300 K with a Nosé−Hoover chain thermostat.42 A periodic unit cell with cubic dimensions of 11 Å containing 35 H2O molecules, 1 furanic species, 1 hydronium ion (H3O+), and 1 chloride anion was used. The corresponding density of the solution is 939 g L−1, and the concentration of H3O+ is 1.25 mol L−1. The cell with the H2O molecules was pre-equilibrated using a Monte Carlo method.43 Then, an optimized furan or furanol molecule, together with one hydronium and one chloride species, was inserted into the box containing H2O molecules. The isolated furanic species were optimized using Gaussian 0944 at the B3LYP/6-311++G(d,p) level of theory; a 10 ps CPMD simulation was then conducted on the created supercell to equilibrate the system before performing the metadynamics simulations. The equilibrium of the systems was examined by plotting the root-mean-square deviation (RMSD) of each trajectory against simulation time, as presented in Figure S1 of the Supporting Information. The RMSD was observed to fluctuate in a small range (±0.15 Å) after 4 ps in all cases. In the first set of calculations, furan was allowed to evolve into dihydrofuranols via acid-catalyzed furan hydration reactions. To map the free energy landscape of this process, we have used two collective variables (CVs) involving the protonation C−H (CV1) and the reacting C−OH bonds (CV2), as presented in Figure 1. In the second set of calculations, the dihydrofuranols were allowed to evolve into the corresponding acyclic products via ring-opening reactions. The free energy landscape of the acid-catalyzed ring openings was evaluated using two CVs that describe the protonation of the heteroatom O1 (CV3) and the cleavage of the cyclic ether bond (CV4), as shown in 4140


Article

Energy & Fuels

Here, dij is the distance between the corresponding atoms, and d0 is the reference distance with the defined values of 2.0 Å for the distance between oxygen and carbon atoms and 1.4 Å for the distance between oxygen and hydrogen atoms. The high power constants p = q = 6 are used to distinguish between the coordinated and non-coordinated states. The CN value between two atoms indicates whether a covalent bond exists (1, bond; 0, no bond). The parameters of fictitious masses (m) and force constants (k), which control the dynamics of the fictitious CVs, are selected to be 100 amu and 3.0 au for CV1 and CV3 and 600 amu and 8.0 au for CV2 and CV4 in all metadynamics computations. These parameters have been shown to be able to adiabatically separate the ionic and electronic degrees of freedom and restrain the CVs close to their coupled fictitious variables (auxiliary particles) of organic reaction systems46,47 and have been previously used to study similar reaction systems.36,37 A normal Gaussian bias potential was used in our metadynamics simulations, with the width of the Gaussian potential set to 0.05. The heights of the Gaussian hills were determined by the shape of the underlying energy surface and set to be at most 0.6 kcal mol−1 in all simulations, except for the 2,5-dihydro-2-furanol ring-opening reaction. The average hill heights in these simulations were in the order of 0.4 kcal mol−1. The maximum height of Gaussian hills for the simulation of 2,5-dihydro-2-furanol ring opening was set to be 0.2 kcal mol−1 as a result of the relatively low energy barrier of this reaction. The addition of bias potential was allowed, with a minimum time

Figure 1. Relevant nomenclatures of atoms and definitions of CV1, CV2, CV3, and CV4 in the CPMD metadynamics computations of acid-catalyzed ring opening of furan. Bonds shown in red are used to define the CVs for the acid-catalyzed hydration of furan, while those in blue are for the CVs used to characterize ring opening of dihydrofuranols. Figure 1. All hydrogen atoms of the solvent system are considered in mapping the changes of CV1 and CV3. In each simulation, the CV is described as the sum of the coordination numbers (CNs) between relevant atoms. A CN between two atoms i and j is determined from eq 1.45

CN(i , j) =

1 − (dij/d0) p 1 − (dij/d0) p + q

(1)

Figure 2. Free energy landscape of acid-catalyzed hydration of furan to 2,5-dihydro-2-furanol in aqueous phase at 300 K. The energy scale (kcal mol−1) is given in the inset to the left. Insets R, TS, and PI (with relative free energy in parentheses) in the energy landscape correspond to the reactant, transition state, and product basins, respectively. The average lengths for several important bonds of each state and their corresponding standard deviations are presented in angstroms. The redundant solvent molecules of every molecular configuration are deleted to emphasize the atomic configuration of the furanic ring. 4141


Article

Energy & Fuels

Figure 3. Free energy landscape of acid-catalyzed hydration of furan to 2,3-dihydro-3-furanol in the aqueous phase at 300 K. separation of 100 MD steps (∼10 fs) and with the displacements in CVs being larger than 1.5 times the width of the Gaussian function. The addition of consecutive potentials was forced after 300 MD steps (∼30 fs) if the displacement had not exceeded the threshold. An average of the free energies over several independent simulations of the metadynamics biasing potential yields very accurate free energies.48 However, this strategy is impractical for a large molecular system, such as the systems presented here, as a result of the computational resources required. In our study, one metadynamics simulation was computed for each reaction step. Several earlier studies have shown that the free energies computed from a single metadynamics profile are reasonably accurate if the control parameters are selected properly and the studied CV space is well-sampled.49−53 Before reconstruction of the free energy landscapes with the historydependent bias, the evolution of each CV was examined to ensure that the CV spaces of most interest are well-sampled. Graphs of the CV variation as a function of metadynamics steps and the sampling of CV trajectories for each metadynamics simulation can be found in Figures S2−S11 of the Supporting Information. The statistical error from a single metadynamics sampling can be estimated following Bussi et al.54 d S w ⎛ δs 2π ⎞ ⟨ε 2⟩ = ⎜ ⎟ βDτg ⎝ S ⎠

(

exp −

2

∑ k≠0

π 2k 2δs 2 2S 2

2 2

π k

integration of the velocity autocorrelation function of the CV. The smaller diffusion coefficient of the two CVs in each metadynamics simulation is used for the determination of error. d is the dimensionality of the CV space; k is a possible d dimensional vector of integers with non-zero norm; and β is equal to (kBT)−1, in which kB is the Boltzmann constant and T is the temperature. Using these relevant parameters determined from each metadynamics simulation, the statistical uncertainties in our computations are estimated to be in the range of 0.4−1.5 kcal mol−1. Some intermediate states in this study are identified on the basis of small free energy differences that are close to the corresponding metadynamics errors, raising the possibility that these intermediates may not be real stable states during reactions. In addition to the statistical errors, it is known that the BLYP exchange−correlation functional systematically underestimates the reaction energy barriers. In our simulated molecular systems, the energy barriers are associated with proton transfer and formation/ dissociation of carbon−oxygen bonds. The BLYP uncertainties in computations of the energy barriers of formation/dissociation of carbon−oxygen bonds were reported to be 2−5 kcal mol−1.55,56 However, this functional predicts the energy barriers of proton transfer comparable to high-level ab initio computations (MP2 and MP4 methods).57−59 Other errors involved in our CPMD computations may be related to the pseudo-potential for the interaction between core and valence electrons and the control parameters, such as plane wave cutoff, fictitious electronic mass, and time step. These uncertainties were estimated to be 2−3 kcal mol−1 in earlier studies with the same pseudo-potential and similar control parameters as used in our present work.60,61 In the present work, we are interested in comparing the relative free energy difference between α and β

) (2)

where S is the size of the free energy well, w is the Gaussian hill height, δs is the hill width, τg is the metadynamics time step, and D is the diffusion coefficient of the CV and can be calculated from the 4142


Article

Energy & Fuels

Figure 4. Free energy landscape of acid-catalyzed hydration of furan to 2,3-dihydro-2-furanol in the aqueous phase at 300 K. protonation of furan to determine the preferred reaction pathway of acid-catalyzed ring opening of furan rather than a direct comparison to experimental activation energies. The molecular graphics of selected snapshots were rendered with visual molecular dynamics (VMD).62 A hydrogen bond between a molecular pair is identified when the oxygen−oxygen (O···O) distance is less than 3.5 Å and, at the same time, the angle ∠O−H···O is greater than 120°. These geometric criteria have been used in previous studies of aqueous-phase systems.63,64 The average bond lengths and their corresponding standard deviations of the important states along each studied reaction pathway are determined from all of the configurations falling into the region where the free energy is within ∼1.5 kcal mol−1 of each minimum.

O1 and C5−O1 bond lengths are 1.39 and 1.40 Å, slightly longer than the reported gas phase experimental value of 1.364 Å.66 The protonation of furan proceeds via migration of the acid proton from the solvent toward the Cα position, leading to the formation of α-protonated furan. The energy basin R2 in Figure 2 represents a state of α-furanium, and TS1 represents the transition state configuration leading to this energy basin. The free energy change in moving from R1 to R2 (R1 → R2) is 18.1 kcal mol−1 endergonic, and the free energy barrier R1 → TS1 is 25.1 kcal mol−1. The bond length H1−C5 is 1.21 Å in TS1 and 1.10 Å in R2, significantly shorter than in R1. The H1−C5 bond length in R2 is in close agreement with the reported gas-phase value of 1.09 Å for α-furanium.25 The C5−O1 bond increases from 1.40 Å in R1 to 1.50 Å (7.1%) in TS1 and further elongates to 1.53 Å (9.3%) in R2, while the C2−O1 decreases to 1.32 Å (5.0%) in TS1 and R2. The small differences in geometrical parameters between TS1 and R2 indicate that α protonation of furan proceeds via a “late” transition state. There are considerable geometrical changes of the furanic ring upon protonation mainly as a result of an electron redistribution in the furanic ring, thereby activating the second Cα (C2) and Cβ (C4) positions for a nucleophilic attack of H2O. The nucleophilic attack on the second Cα position and regeneration of the acid proton produces 2,5-dihydro-2-furanol, corresponding to state PI1 in Figure 2. Inspection of the free energy landscape shows that the formation of furanol is 22.2 kcal mol−1 activated (R1 → TS2) and 4.3 kcal mol−1 endergonic (R1 → PI1). The average H···OH distance in TS2 of the reacting H2O is 1.17 Å, much longer than the corresponding length in pure liquid water.67 This indicates that

3. RESULTS AND DISCUSSION 3.1. Acid-Catalyzed Hydration of Furan. The free energy landscape of protonation of furan at the Cα position leading to the formation of 2,5-dihydro-2-furanol, together with selected molecular configurations along the reaction pathway, is presented in Figure 2. Selected bond lengths and corresponding standard deviations (Å) are provided in each configuration. The process involves a proton transfer from the solvent to the Cα position, followed by a nucleophilic attack of H2O and regeneration of the acid proton. Furan is solvated by two layers of water molecules, with the acid proton in the second solvation shell in the form of a H5O2+ species. The molecular configuration R1 shows that the acid proton connects with O1 via five hydrogen bonds. However, a direct interaction between carbon atoms of furan and solvent molecules is not observed. The hydrophobicity of furan is consistent with its low solubility in water (1 wt % at 25 °C).65 The closest distance between hydrogen of H2O and the ring oxygen O1 is 2.62 Å, indicating a weak interaction. The C2− 4143


Article

Energy & Fuels

Figure 5. Free energy landscape of acid-catalyzed ring opening of 2,5-dihydro-2-furanol in the aqueous phase at 300 K.

activated (R1′ → TS3) and 6.1 kcal mol−1 endergonic (R1′ → PI2). The free energy barrier corresponding to the formation of 2,3-dihydro-3-furanol (R1′ → TS3) is only 1.3 kcal mol−1 greater than that predicted for the formation of 2,5-dihydro-2furanol (R1 → TS2), indicating that there is little kinetic preference for either reaction pathway. Thermodynamically, the formation of 2,3-dihydro-3-furanol (R1′ → PI2) is slightly less favored because it is 1.8 kcal mol−1 more endergonic than the formation of 2,5-dihydro-2-furanol (R1 → PI1). Considering that the energy differences between these two reaction pathways are close to the uncertainties of the computational method, the reaction pathways will depend upon the stabilities of the various furanols. Figure 4 presents the free energy landscape and selected molecular configurations corresponding to the protonation of furan at the Cβ position, followed by a nucleophilic attack of water at an adjacent Cα position, forming 2,3-dihydro-2-furanol. There is no clear basin state for the formation of the βfuranium species. The distance between C2 in TS4 and the closest water molecule is in the order of 2.53 Å, significantly shorter than the corresponding lengths between water and αfuranium (3.22 Å in R2; Figure 2). The formation of the βfuranium species (R1″ → TS4) is 30.5 kcal mol−1 activated, which is 7 kcal mol−1 more activated than the formation of αfuranium. The energy gain after passing TS4 is 24.7 kcal mol−1, which is larger than the energy gain from α-furanium to the corresponding furanols, as seen in Figures 2 and 3, further indicating that β-furanium is less stable than α-furanium in solution. A ring cleavage mechanism via β protonation and the

the nucleophilic attack occurs concurrently with the regeneration of the acid proton. The C2−O1, C5−O1, and H1−C2 bond lengths in TS2 are very close to those in R2, suggesting an “early” transition state during the hydration process. Figure 3 shows the free energy landscape and selected molecular configurations of protonation of furan at the Cα position, leading to the formation of 2,3-dihydro-3-furanol. The mechanism of the formation of α-furanium is similar to that in the formation of 2,5-dihydro-2-furanol discussed above, with the process being 16.9 kcal mol−1 endergonic (R1′ → R2′) and 23.8 kcal mol−1 activated (R1′ → TS1′). These two values are 1.2 and 1.3 kcal mol−1 lower than those determined during the formation of 2,5-dihydro-2-furanol, respectively, but these differences are within the estimated statistical errors (0.4−1.5 kcal mol−1), as discussed in section 2. Inspection of R1′, TS1′, and R2′ in Figure 3 shows that these molecular configurations are essentially the same as those observed during the formation of 2,5-dihydro-2-furanol; i.e., the furan configuration (R1′) connects via five hydrogen bonds with the acid proton, with a distance between the closest H2O and O1 in the order of 2.73 Å. The closely similar free energy profiles and molecular configurations of the Cα-protonation process in the two distinct simulations indicate a good consistency and energy convergence in these calculations. A nucleophilic attack of H2O on the adjacent Cβ position of α-furanium (R2′) regenerates the acid proton, forming 2,3dihydro-3-furanol (PI2). The nucleophilic attack is concerted with the regeneration of the acid proton; i.e., the average H··· OH distance in TS3 is 1.19 Å. This process is 23.5 kcal mol−1 4144


Article

Energy & Fuels

Figure 6. Free energy landscape of acid-catalyzed ring opening of 2,3-dihydro-3-furanol to produce 3,4-dihydroxybutenol in the aqueous phase at 300 K.

2-butenal is nearly thermo-neutral, being only 2.6 kcal mol−1 endergonic. Figure 6 shows the free energy landscape of the ring opening of 2,3-dihydro-3-furanol (PI2′ in Figure 6). The distance between the first water shell and furanol is 1.93 Å, which is only 0.07 Å shorter than in 2,5-dihydro-2-furanol. However, the acid proton is equilibrated at five hydrogen bonds from 2,3-dihydro3-furanol, one more hydrogen bond apart than in the case of 2,5-dihydro-2-furanol. The ring opening proceeds via a single step PI2′ → TS6 → P2. The energy barrier is 34.2 kcal mol−1, and the process is 14.2 kcal mol−1 endergonic. The ring opening occurs via the C2−O1 bond dissociation. The bond length of the new HP−O1 in TS6 is 1.12 Å, and the length of C2−O1 is 2.06 Å, indicating a concerted protonation and ringopening process. A H2O molecule at a point near C2 with an average distance of 2.35 Å is required to stabilize the TS6 configuration. The direct reaction product is 3,4-dihydroxybutenol, with the proton equilibrated in the second solvation shell. During our metadynamics simulations, it was observed that 3,4-dihydroxybutenol is rapidly converted to 4-hydroxy-2-butenal after it was formed. The conversion proceeded by losing the −OH group connected to C3 and the HP connected to O1. The conversion process from 3,4-dihydroxybutenol to 4-hydroxy-2butenal is not shown in the free energy landscape in Figure 6 because this was not clearly distinguished on the basis of the currently selected CV3 and CV4. This reaction was not studied further because the formation of 3,4-dihydroxybutenol is less thermodynamically favored than the formation of 4-hydroxy-2butenal from 2,5-dihydro-2-furanol.

possibility of an internal proton shift from the α to β position are unfavorable as a result of the instability of the β-furanium configuration. Therefore, we predict that 2,5-dihydro-2-furanol and 2,3-dihydro-3-furanol are the main products of furan reactions with acid protons and water. 3.2. Acid-Catalyzed Ring Opening of Furanol Intermediates. The free energy landscapes and some selected geometrical parameters of the acid-catalyzed ring openings of 2,5-dihydro-2-furanol and 2,3-dihydro-3-furanol as a function of proton transfer from H2O to the furanic ring are presented in Figures 5 and 6, respectively. The interaction of H2O with 2,5-dihydro-2-furanol (PI1′ in Figure 5) is significantly stronger than with furan; i.e., the closest distance is on average 0.72 Å shorter than the corresponding values in furan. PI1′ is connected with the acid proton via four hydrogen bonds, equilibrated in the second solvation shell. The ring opening initiates with an increased degree of interaction between O1 and the proton. This interaction leads to the protonation of the ring oxygen O1 and the dissociation of the C2−O1 bond (TS5), which is 13.4 kcal mol−1 activated. Inspection of TS5 shows that a new HP−O1 bond is almost formed (1.07 Å), indicating that the proton has been transferred from bulk H2O to the ring. The C2−O1 bond in TS5 increased from 1.48 Å in PI1′ to 2.05 Å, suggesting that TS5 is formed via a concerted process involving dissociation of C2−O1 and protonation of O1. The system gains 10.8 kcal mol−1 after passing TS5 and reaches the energy basin of 4hydroxy-2-butenal. The average distance of C2−O in 4hydroxy-2-butenal is 1.28 Å, indicating that an aldehyde group has been completely formed. Formation of 4-hydroxy4145


Article

Energy & Fuels Figures 2−6 indicate that furan decomposition in aqueous solution is driven by the formation of the α-protonated furan. The computed free energy barrier of ∼24 kcal mol−1 is consistent with reported activation energies of acid-catalyzed furanic ring reactions, e.g., 22.8 kcal mol −1 for 2,5dimethylfuran and 26.0 kcal mol−1 for furan in 0.2 mol L−1 HClO4 solution,15 17 kcal mol−1 for furan in 4.61 mol L−1 HClO4 solution,23 29.9 kcal mol−1 for furfural in dilute H2SO4 solution,8 and 22.8−26.4 kcal mol−1 for 5-HMF in various dilute acid solutions.9,68−70 These different activation energies suggest an important effect of the acid concentration or changes in the reaction mechanisms as a result of the side chains of the furanic ring. It should be noted that furan and the acid proton are both located inside a finite unit cell in our simulations. The free energy landscapes thus inevitably include the interactions between furan and proton species solvated in a second or third solvation shell. Previous metadynamics studies of acid-catalyzed sugar reactions have shown that the proton transportation from bulk solvent to the proximity of sugar molecules is endergonic, therefore contributing to the overall reaction energy barriers.36,37,60 Therefore, additional static calculations were carried out to estimate the energy difference between equilibrated furan in the presence of solvated H3O+ with a system in which furan and proton are individually solvated. The gas-phase equilibrium conformations were computed using the B3LYP/6-311++G(d,p) level of theory, and the free energies of solvation were calculated at the same level of theory using the solvation model based on density (SMD).71 The detailed description of the method for computing the free energies of infinite separation can be found in earlier publications.36,37 The energy of bringing together a solvated proton with furan is endothermic by 7 kcal mol−1, attributed to a greater proton attraction ability of bulk water compared to furan. The global free energy barrier in our simulations is therefore 7 + 23.8 (Figure 2) = ∼30 kcal mol−1, referenced from an infinite separation of solvated furan from the acid proton.

4. CONCLUSION CPMD simulations were carried out to determine the reaction pathways and corresponding free energy landscapes for the acid-catalyzed ring opening of furan in the aqueous phase. We observed that furan interacts with water only via a weak hydrogen bond linked to its ring oxygen. Protonation of furan is endergonic, indicating a greater proton attraction by bulk water molecules than furan. Protonation at the Cα position is energetically preferred to that at the Cβ position. The α protonation leads to the formation of 2,5-dihydro-2-furanol and 2,3-dihydro-3-furanol. Resultant 2,5-dihydro-2-furanol is subject to the protonation on the ring oxygen, triggering a ring opening to produce 4-hydroxy-2-butenal. In contrast, the ring opening of 2,3-dihydro-3-furanol is energetically less favored as a result of its high free energy barrier. We predict that a primary product of decomposition of furan in acid solution is 4hydroxy-2-butenal.

■

■

RMSD in coordinates as a function of the simulation time for the equilibrated furan, 2,3-dihydro-3-furanol, and 2,5-dihydro-2-furanol chemistry systems described in section 2 (Figure S1), evolution of the CVs as a function of metadynamics steps in the simulation of acid-catalyzed hydration of furan to 2,5-dihydro-2-furanol (Figure S2), metadynamics sampling of the studied CV trajectories for the proton transfer (CV1) and formation/breakage of the >C−OH bond (CV2) in the simulation of acidcatalyzed hydration of furan to 2,5-dihydro-2-furanol (Figure S3), evolution of the CVs as a function of metadynamics steps in the simulation of acid-catalyzed hydration of furan to 2,3-dihydro-3-furanol (Figure S4), metadynamics sampling of the studied CV trajectories for the proton transfer (CV1) and formation/breakage of the >C−OH bond (CV2) in the simulation of acidcatalyzed hydration of furan to 2,3-dihydro-3-furanol (Figure S5), evolution of the CVs as a function of metadynamics steps in the simulation of acid-catalyzed hydration of furan to 2,3-dihydro-2-furanol (Figure S6), metadynamics sampling of the studied CV trajectories for the proton transfer (CV1) and formation/breakage of the >C−OH bond (CV2) in the simulation of acidcatalyzed hydration of furan to 2,3-dihydro-2-furanol (Figure S7), evolution of the CVs as a function of metadynamics steps in the simulation of acid-catalyzed ring opening of 2,5-dihydro-2-furanol to produce 4hydroxy-2-butenal (Figure S8), metadynamics sampling of the studied CV trajectories for the proton transfer (CV3) and formation/breakage of the C−O bond (CV4) in the simulation of acid-catalyzed ring opening of 2,5dihydro-2-furanol to produce 4-hydroxy-2-butenal (Figure S9), evolution of the CVs as a function of metadynamics steps in the simulation of acid-catalyzed ring opening of 2,3-dihydro-3-furanol to produce 3,4dihydroxybutenol (Figure S10), and metadynamics sampling of the studied CV trajectories for the proton transfer (CV3) and formation/breakage of the C−O bond (CV4) in the simulation of acid-catalyzed ring opening of 2,3-dihydro-3-furanol to produce 3,4dihydroxybutenol (Figure S11) (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Xiao Liang: 0000-0003-1201-9166 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

This work has been funded by the Australian Research Council through Discovery Grant DP1096802. The authors acknowledge The University of Sydney high performance computing (Artemis) resources and the National Computational Infrastructure (NCI) facility supported by the Australian Government for providing resources that have contributed to the results reported in this paper.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.energyfuels.7b03239. 4146


Article

Energy & Fuels

■

(22) Kankaanperä, A.; Kleemola, S.; Nielsen, P. H. Equilibrium protonation of a carbon base prior to its hydrolysis. Acid-catalyzed cleavage of the furan ring. Acta Chem. Scand. 1969, 23, 3607−3608. (23) Unverferth, K.; Schwetlick, K. Mechanism of furan hydrolysis. J. Prakt. Chem. 1970, 312, 882−886. (24) Houriet, R.; Schwarz, H.; Zummack, W.; Andrade, J. G.; Schleyer, P. V. R. The α- vs β-protonation of pyrrole, furan, thiophene and cyclopentadiene. Gas phase proton and hydrogen affinities. The bishomocyclopropenyl cation. Nouv. J. Chim. 1981, 5, 505−509. (25) Lorenz, U. J.; et al. Protonation of heterocyclic aromatic molecules: IR signature of the protonation site of furan and pyrrole. Int. J. Mass Spectrom. 2007, 267, 43−53. (26) Belen’kii, L. I. Positional selectivity in electrophilic substitution in π-excessive heteroaromatics. Adv. Heterocycl. Chem. 2010, 99, 143− 183. (27) Hiraoka, K.; Takimoto, H.; Yamabe, S. Stabilities and structures in cluster ions of five-membered heterocyclic compounds containing oxygen, nitrogen, and sulfur atoms. J. Am. Chem. Soc. 1987, 109, 7346−7352. (28) Gubina, T.; Pankratov, A.; Labunskaya, V.; Voronin, S.; Kharchenko, V. Study of the mechanism of recyclization of furans into thiophenes and selenophenes in conditions of acid catalysis. 6. Experiments with labeled atoms. Quantum chemical calculations of intermediates of recyclization and hydrolysis. Chem. Heterocycl. Compd. 1997, 33, 903−909. (29) Gubina, T.; Labunskaya, V.; Pankratov, A.; Voronin, S.; Kharchenko, V. Study of the mechanism of recyclization of furans into thiophenes and selenophenes in conditions of acid catalysis. 5. Direction of protonation of furans. Chem. Heterocycl. Compd. 1997, 33, 898−902. (30) van Beelen, E. S. E.; Koblenz, T. A.; Ingemann, S.; Hammerum, S. Experimental and Theoretical Evaluation of Proton Affinities of Furan, the Methylphenols, and the Related Anisoles. J. Phys. Chem. A 2004, 108, 2787−2793. (31) Zeng, K.; Cao, Z.-X. Protonation of pyrrole and furan by H3O+ and NH4+ in the gas phase: A density functional theory study. Chin. J. Chem. 2006, 24, 293−298. (32) Nikbin, N.; Caratzoulas, S.; Vlachos, D. G. On the Brønsted Acid-Catalyzed Homogeneous Hydrolysis of Furans. ChemSusChem 2013, 6, 2066−2068. (33) Park, J. M.; Laio, A.; Iannuzzi, M.; Parrinello, M. Dissociation Mechanism of Acetic Acid in Water. J. Am. Chem. Soc. 2006, 128, 11318−11319. (34) Cucinotta, C. S.; Ruini, A.; Catellani, A.; Stirling, A. Ab initio molecular dynamics study of the keto-enol tautomerism of acetone in solution. ChemPhysChem 2006, 7, 1229−1234. (35) Schreiner, E.; Nair, N. N.; Marx, D. Influence of Extreme Thermodynamic Conditions and Pyrite Surfaces on Peptide Synthesis in Aqueous Media. J. Am. Chem. Soc. 2008, 130, 2768−2770. (36) Liang, X.; Montoya, A.; Haynes, B. S. Local Site Selectivity and Conformational Structures in the Glycosidic Bond Scission of Cellobiose. J. Phys. Chem. B 2011, 115, 10682−10691. (37) Liang, X.; Montoya, A.; Haynes, B. S. Molecular Dynamics Study of Acid-Catalyzed Hydrolysis of Dimethyl Ether in Aqueous Solution. J. Phys. Chem. B 2011, 115, 8199−8206. (38) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471. (39) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (40) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (41) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (42) Nose, S. A unified formulation of the constant-temperature molecular-dynamics methods. J. Chem. Phys. 1984, 81, 511−519.

REFERENCES

(1) Lange, J.-P.; van der Heide, E.; van Buijtenen, J.; Price, R. FurfuralA Promising Platform for Lignocellulosic Biofuels. ChemSusChem 2012, 5, 150−166. (2) Rosatella, A. A.; Simeonov, S. P.; Frade, R. F. M.; Afonso, C. A. M. 5-Hydroxymethylfurfural (HMF) as a building block platform: Biological properties, synthesis and synthetic applications. Green Chem. 2011, 13, 754−793. (3) Kabyemela, B. M.; Adschiri, T.; Malaluan, R. M.; Arai, K. Glucose and fructose decomposition in subcritical and supercritical water: Detailed reaction pathway, mechanisms, and kinetics. Ind. Eng. Chem. Res. 1999, 38, 2888−2895. (4) Chuntanapum, A.; Yong, T. L.-K.; Miyake, S.; Matsumura, Y. Behavior of 5-HMF in Subcritical and Supercritical Water. Ind. Eng. Chem. Res. 2008, 47, 2956−2962. (5) Aida, T. M.; et al. Reactions of D-fructose in water at temperatures up to 400 °C and pressures up to 100 MPa. J. Supercrit. Fluids 2007, 42, 110−119. (6) Aida, T. M.; et al. Dehydration of D-glucose in high temperature water at pressures up to 80 MPa. J. Supercrit. Fluids 2007, 40, 381− 388. (7) Girisuta, B.; Janssen, L. P. B. M.; Heeres, H. J. Green chemicals a kinetic study on the conversion of glucose to levulinic acid. Chem. Eng. Res. Des. 2006, 84, 339−349. (8) Marcotullio, G.; Tavares Cardoso, M. A.; De Jong, W.; Verkooijen, A. H. M. Bioenergy II: Furfural destruction kinetics during sulphuric acid-catalyzed production from biomass. Int. J. Chem. React. Eng. 2009, 7, A67. (9) Asghari, F. S.; Yoshida, H. Kinetics of the Decomposition of Fructose Catalyzed by Hydrochloric Acid in Subcritical Water: Formation of 5-Hydroxymethylfurfural, Levulinic, and Formic Acids. Ind. Eng. Chem. Res. 2007, 46, 7703−7710. (10) Rose, I. C.; Epstein, N.; Watkinson, A. P. Acid-Catalyzed 2Furaldehyde (Furfural) Decomposition Kinetics. Ind. Eng. Chem. Res. 2000, 39, 843−845. (11) Dunlop, A. P. Furfural formation and behavior. Ind. Eng. Chem. 1948, 40, 204−209. (12) Skrabal, A.; Skrabal, R. The velocity of hydrolysis of enol ethers. Z. Phys. Chem. 1937, 181A, 449−468. (13) Gilman, H. Organic Chemistry; John Wiley and Sons, Inc.: New York, 1953; Vol. 4, pp 740. (14) Joule, J. A.; Keith, M. Heterocyclic Chemistry; Blackwell Publishing, Ltd.: Chichester, U.K., 2010; pp 347. (15) Kankaanperä, A.; Salomaa, P.; Mellingen, K.; Åkeson, Å.; Theorell, H.; Blinc, R.; Paušak, S.; Ehrenberg, L.; Dumanović, J. Kinetics of the hydrolytic cleavage of the furan ring. Acta Chem. Scand. 1967, 21, 575−577. (16) Salomaa, P.; Kankaanperä, A.; Nikander, E.; Kaipainen, K.; Aaltonen, R.; Swahn, C.-G. Electrophilic attack on the furan ring. Acidcatalyzed protodedeuteriation, protodetritiation, deuteriodeprotonation, and deuteriodetritiation at the α-position. Hydrolytic cleavage of the ring. Acta Chem. Scand. 1973, 27, 153−163. (17) Kankaanperä, A.; Aaltonen, R.; Liaaen-Jensen, S.; Tricker, M. J.; Svensson, S. General acid catalysis in the hydrolysis of a furan derivative. Acta Chem. Scand. 1972, 26, 2537−2540. (18) Garst, J. E.; Schmir, G. L. Hydrolysis of 2-methoxyfuran. J. Org. Chem. 1974, 39, 2920−2923. (19) Horvat, J.; Klaic, B.; Metelko, B.; Sunjic, V. Mechanism of levulinic acid formation. Tetrahedron Lett. 1985, 26, 2111−2114. (20) Horvat, J.; Klaic, B.; Metelko, B.; Sunjic, V. Mechanism of levulinic acid formation in acid catalyzed hydrolysis of 2(hydroxymethyl)furan and 5-(hydroxymethyl)-2-furancarboxaldehyde. Croat. Chem. Acta 1986, 59, 429−438. (21) Stamhuis, E. J.; Drenth, W.; van den Berg, H. Mechanism of reactions of furans. I. A kinetic study of the acid-catalyzed hydrolysis of furan and 2,5-dimethylfuran. Recl. Trav. Chim. Pays-Bas 1964, 83, 167−176. 4147


Article

Energy & Fuels (43) Coutinho, K.; Canuto, S. DICE: A Monte Carlo Program for Molecular Liquid Simulation, Version 2.9; University of São Paulo: São Paulo, Brazil, 2003. (44) 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, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (45) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302. (46) Qian, X. Mechanisms and Energetics for Acid Catalyzed β-DGlucose Conversion to 5-Hydroxymethylfurfurl. J. Phys. Chem. A 2011, 115, 11740−11748. (47) Liu, D.; Nimlos, M. R.; Johnson, D. K.; Himmel, M. E.; Qian, X. Free Energy Landscape for Glucose Condensation Reactions. J. Phys. Chem. A 2010, 114, 12936−12944. (48) Laio, A.; Gervasio, F. L. Metadynamics: A method to stimulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601. (49) Petersen, L.; Ardevol, A.; Rovira, C.; Reilly, P. J. Mechanism of Cellulose Hydrolysis by Inverting GH8 Endoglucanases: A QM/MM Metadynamics Study. J. Phys. Chem. B 2009, 113, 7331−7339. (50) Fu, C.-W.; Lin, T.-H. Theoretical Study on the Alkaline Hydrolysis of Methyl Thioacetate in Aqueous Solution. J. Phys. Chem. A 2011, 115, 13523−13533. (51) Barker, I. J.; Petersen, L.; Reilly, P. J. Mechanism of Xylobiose Hydrolysis by GH43 β-Xylosidase. J. Phys. Chem. B 2010, 114, 15389− 15393. (52) Alfonso-Prieto, M.; Biarnes, X.; Vidossich, P.; Rovira, C. The Molecular Mechanism of the Catalase Reaction. J. Am. Chem. Soc. 2009, 131, 11751−11761. (53) Comas-Vives, A.; Stirling, A.; Lledos, A.; Ujaque, G. The Wacker Process: Inner- or Outer-Sphere Nucleophilic Addition? New Insights from Ab Initio Molecular Dynamics. Chem. - Eur. J. 2010, 16, 8738− 8747. (54) Bussi, G.; Laio, A.; Parrinello, M. Equilibrium Free Energies from Nonequilibrium Metadynamics. Phys. Rev. Lett. 2006, 96, 090601. (55) Stirling, A.; Papai, I. H2CO3 Forms via HCO3- in Water. J. Phys. Chem. B 2010, 114, 16854−16859. (56) Blumberger, J.; Ensing, B.; Klein, M. L. Formamide hydrolysis in alkaline aqueous solution: Insight from ab initio metadynamics calculations. Angew. Chem., Int. Ed. 2006, 45, 2893−2897. (57) Geissler, P. L.; Van Voorhis, T.; Dellago, C. Potential energy landscape for proton transfer in (H2O)3H+: Comparison of density functional theory and wave function-based methods. Chem. Phys. Lett. 2000, 324, 149−155. (58) Pan, Y.; McAllister, M. A. Characterization of low-barrier hydrogen bonds. 4. Basis set and correlation effects: An ab initio and DFT investigation. J. Mol. Struct.: THEOCHEM 1998, 427, 221−227. (59) Xu, X.; Goddard, W. A. III. Bonding Properties of the Water Dimer: A Comparative Study of Density Functional Theories. J. Phys. Chem. A 2004, 108, 2305−2313. (60) Dong, H.; Nimlos, M. R.; Himmel, M. E.; Johnson, D. K.; Qian, X. The Effects of Water on β-D-Xylose Condensation Reactions. J. Phys. Chem. A 2009, 113, 8577−8585.

(61) Isayev, O.; Gorb, L.; Leszczynski, J. Theoretical calculations: Can Gibbs free energy for intermolecular complexes be predicted efficiently and accurately? J. Comput. Chem. 2007, 28, 1598−1609. (62) Humphrey, W.; Dalke, A.; Schulten, K. VDM: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (63) Lerbret, A.; et al. Molecular dynamics studies of the conformation of sorbitol. Carbohydr. Res. 2009, 344, 2229−2235. (64) Morando, M. A.; et al. NMR and molecular modeling reveal key structural features of synthetic nodulation factors. Glycobiology 2011, 21, 824−833. (65) McKillip, W. J.; Collin, G.; Höke, H.; Zeitsch, K. J. Furan and derivatives. Ullmann’s Encyclopedia of Industrial Chemistry; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 1981. (66) Liescheski, P. B.; Rankin, D. W. H. Molecular structure of furan, determined by combined analyses of data obtained by electron diffraction, rotational spectroscopy, and liquid crystal NMR spectroscopy. J. Mol. Struct. 1989, 196, 1−19. (67) Ichikawa, K.; Kameda, Y.; Yamaguchi, T.; Wakita, H.; Misawa, M. Neutron-diffraction investigation of the intramolecular structure of a water molecule in the liquid phase at high temperatures. Mol. Phys. 1991, 73, 79−86. (68) Heimlich, K. R.; Martin, A. N. A kinetic study of glucose degradation in acid solution. J. Am. Pharm. Assoc., Sci. Ed. 1960, 49, 592−597. (69) McKibbins, S. W.; Harris, J. F.; Saeman, J. F.; Neill, W. K. Kinetics of the acid catalyzed conversion of glucose to 5hydroxymethyl 2-furaldehyde and levulinic acid. For. Prod. J. 1962, 12, 17−23. (70) Girisuta, B.; Janssen, L. P. B. M.; Heeres, H. J. A kinetic study on the decomposition of 5-hydroxymethylfurfural into levulinic acid. Green Chem. 2006, 8, 701−709. (71) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 2009, 113, 6378−6396.

4148


Acid-Catalyzed Ring Opening of Furan in Aqueous Solution - Energy

Recommend Documents