Understanding Competition of Polyalcohol Dehydration Reactions in Hot Water

Dehydration of biomass-derived polyalcohols has recently drawn attention in green chemistry as a prototype of selective reactions controllable in hot water or hot carbonated water, without any use of organic solvents or metal catalysts. Here we report a free energy analysis based on first-principles metadynamics and blue-moon ensemble simulations to understand the mechanism of competing intramolecular dehydration reactions of 1,2,5-pentanetriol (PTO) in hot acidic water. The simulations consistently predict that the most dominant mechanism is the proton-assisted SN2 process where the protonation of the hydroxyl group by water and the C-O bond breaking and formation occurs in a single step. However the free energy barriers are different between the reaction paths: those leading to five membered ether products, tetrahydrofurfuryl alcohol (THFA), are few kcal/mol lower than those leading to six membered ether products, 3-hydroxytetrahydropyran (3-HTHP). A slight difference is seen in the timing of the protonation of the hydroxyl group of THFA and 3-HTHP on their reaction pathways. The detailed mechanism found from the simulations shows how the reaction paths are selective in hot water and why the reaction rates are accelerated in acidic environments, thus giving a clear explanation of experimental findings for a broad class of competing dehydration processes of polyalcohols.


Introduction
The development of green and sustainable processes for the production of chemicals and energy carriers is one of the most pressing challenges faced by our society. In this context, biomass conversion in hot pressurized water (be it either sub-or supercritical water) has recently emerged as a very promising technology. [1][2][3][4][5][6] On the one hand, the use of renewable biomass resources associated with this technology can give rise to an environmentally beneficial reduction in the carbon footprint of chemicals and liquid fuel and can contribute to the development of a sustainable industrial society. [7][8][9][10] In addition to this, the use of water as a reaction medium for organic processes is in line with green chemistry practices because water is non-toxic, non-flammable, inexpensive, environmentally benign and readily available. 11 In contrast with the low solubility for many small organic molecules in room-temperature liquid water, the same compounds are much more soluble in liquid water above ∼ 200 • C by virtue of a drastic decrease of the dielectric constant. 12 Consequently, an increase in temperature of liquid water widens the scope of reactions that can be carried out, thus rendering both subcritical high-temperature liquid water (HTW, i.e. liquid water above ∼ 200 • C) and supercritical water very interesting reaction media for biomass conversion and organic synthesis reactions. [13][14][15][16] The use of hot pressurized water as a reaction medium, besides being aligned with green chemistry guidelines, is very appealing from the perspective of optimization of chemical processes because water temperature and water density or pressure can exert a strong influence on product selectivities and reaction rates. 17-20 An in-depth knowledge of the reaction mechanism and the role played by water is clearly very beneficial when it comes to identifying the optimal conditions for a given reaction. Detailed knowledge on how the special properties of hot pressurized water affect reactivity is also very important in other research areas, such as geochemistry 21,22 and origin of life. 23,24 It can be very difficult to elucidate the role of hot pressurized water in reaction mechanisms exclusively from experiments because water can interact in many and complex ways with the reacting system (through hydrogen bonding, its polarity, its capacity of acting as an acid or basic catalyst, hydrophobic effects, etc.) and because reaction mechanisms operating in room-temperature water can drastically differ from those operating at high temperatures and pressures. 25,26 In view of the proven utility of computational chemistry in establishing reaction mechanisms in room-temperature aqueous media, [27][28][29][30][31] it is clear that computational studies of reactivity in hot pressurized water, which are still scarce 25,26,32,33 in comparison with those at ambient temperature, can provide valuable insights into reaction mechanisms, thereby facilitating the design of optimal chemical processes in hydrothermal conditions. Here we elucidate the mechanism and the origin of selectivity of an intramolecular dehydration process that serves as a model system for polyalcohol dehydration reactions in hot pressurized water.
Among all processes relevant for biomass conversion, the dehydration of polyalcohol compounds is one of the most intensively investigated class of reactions. [34][35][36] The use of CO 2 -enriched subcritical high temperature liquid water (HTW) has recently emerged as a very promising technology for polyalcohol dehydration 37-44 and other biomass conversion or organic synthesis processes. [45][46][47] The addition of CO 2 to HTW leads to the formation of carbonic acid, which results in an acidic aqueous medium that is able to provide acid-catalyzed dehydration reactions. The use of CO 2 instead of mineral acids as a means of providing an acid medium responds to environmental concerns associated with the latter type of acids.
The chemistry of intramolecular dehydration reactions of biomass-derived carbohydrates is very complex due to the large number of alcohol groups of these molecules and, thus, the large number of different products that can be formed. Taking this into consideration, studies on simpler model systems are very relevant because they can provide easier ways of establishing the underlying reaction mechanisms. One of the most prominent experimental studies on model systems for intramolecular dehydration reactions was reported by the research group led by Shirai a few years ago. 39 In this study, 1,2,5-pentanetriol (PTO) was successfully converted into tetrahydrofurfuryl alcohol (THFA) and 3-hydroxytetrahydropyran (3-HTHP) by intramolecular dehydration in HTW at 573 K in the presence of 20 MPa of CO 2 , as shown in Figure 1. This conversion was observed to be highly selective since the yield ratio of THFA to 3-HTHP was 7:1 in the total product mixture. This led to the conclusion that five-membered cyclic ethers were formed faster than six-membered cyclic ethers. 39 In the present study, we carried out first principles molecular dynamics simulations in the condensed phase (under the same conditions as described in the experiment by Shirai and coworkers) combined with several enhanced sampling techniques with the goal of establishing the mechanism of PTO dehydration in HTW and elucidating the origin of its selectivity. It should be mentioned that the mechanism of dehydration reactions of dialcohols or glycerol in HTW has been investigated in previous computational works using different simulation techniques. 32,48,49 However, this is the first computational work dealing with selectivity issues when both five-membered and six-membered cyclic ethers can be formed.

Method
Metadynamics 50 (MTD) simulations were used to study the free energy landscape of the aqueous PTO molecule with respect to the reactive regions toward the THFA and 3-HTHP products. Blue moon ensemble 51 (BME) simulations were used to study free energy changes with respect to the non-reactive regions among the contracted and extended conformations of the aqueous PTO molecule. The simulations were conducted in the Born-Oppenheimer scheme, both for those based on first-principles density functional theory (DFT) and for those based on the semiempirical density functional tight binding (DFTB) method. We used the PIMD 52 code in which MTD and BME methods are implemented in hierarchical parallelization linking to electronic structure codes, in this case, 53 Vienna ab initio simulation package (VASP) 54-57 and DFTB+. [58][59][60][61] The MTD and BME simulations of the aqueous PTO molecule were carried out in the canonical (NVT) ensemble with constant volume corresponding to the pressure of 20 MPa and constant temperature of 573 K, i.e., the same conditions as the previous exper-iment. 39 The temperature was strongly controlled by massive Nosé-Hoover chain (MNHC) thermostats 62-64 for efficient statistical sampling. The simulations were executed in cubic boxes with periodic boundary conditions. The box sizes were determined from preliminary classical molecular dynamics simulations in the NPT ensemble at 20 MPa and 573 K using OPLS force field. 65 In order to mimic the experiment in acidic condition, we added one proton into the simulation box. This resulted in the aqueous solution being more acidic in the simulations (pH < 1) than in the experiment (pH ∼ 4), but this limitation cannot be avoided for molecular simulations taking account of water molecules explicitly. We note that the results of the current study on PTO dehydration with 30-70 water molecules should not be affected significantly by the system size since the reactions take place locally. In fact, in the case of hexanediol dehydration the free energy barrier in a microsolvated cluster with only six water molecules was found to have the same value (36 kcal/mol) as that in bulk water with 70 water molecules with periodic boundary conditions. 32 Now we describe the simulation setup in detail. A summary is also shown in Table 1.

DFT
The DFT simulations were carried out for a system containing one PTO molecule, 30 water molecules and a proton. Uniform background charge was added to maintain the charge neutrality of the system. The electronic structure calculations were based on density functional theory with the PBE exchange correlation functional 66 and Grimme's D2 van der Waals correction 67 using a cut off at 400 eV in the plane wave (PW) basis functions. The core electrons were taken into account using the projector augmented wave (PAW) method. 68,69 To check the accuracy and validity of the PBE functional, a set of static calculations were performed based on the polarizable continuum model using a dielectric constant in hot water ϵ = 25, see SI. With respect to the relative energies for the confined/extended conformations of PTO and the relative energies for the 3-HTHP and THFA products in the respective optimized geometries ( Figure S5), it was found that PBE functional exhibits the same qualitative trend as the calculation methods of higher accuracy such as B3LYP and MP2 (Table S1).

DFTB
The DFTB simulations were carried out for a system containing one PTO molecule, 70 water molecules, a proton and a chloride ion. In DFTB only the valence electrons are explicitly solved. The 3ob parameter set 70,71 was used to describe the DFTB Hamiltonian.

Metadynamics
A series of three-dimensional metadynamics was carried out following a similar protocol proposed by Ensing et al. 72 except that multiple walkers were used. 73 The pseudo-Hamiltonian Here, the first term is the Hamiltonian of the molecular system given by the sum of kinetic energy of atoms and the Born-Oppenheimer potential calculated with either the DFT or DFTB method. The second term is the kinetic energy of fictitious particles (walkers) where p i,k and µ i are the momentum and the mass of the walker, i and k being indices of the collective variable (CV) component and the walker, respectively. The third term is a harmonic potential that binds the position of walker, s i,k , to the actual CV value, is the force constant. The fourth term is a history-dependent potential given by the sum of small Gaussian functions (hills) deposited to the places where all the walkers have visited until time t, which can be expressed as where h and σ i are the hill height and the hill width, respectively, and N w is the number of walkers. The fifth term is the contribution from the MNHC thermostats. The free energy surface is estimated by where r b = 1.4 Å is a parameter to identify OH bonds. However, the atom numbers relevant to the CVs should be chosen differently for the respective reaction paths, see Figure 1. The atoms associated to the CVs for each reaction path are listed in Table 2.
We used N w = 12 walkers to speed up the calculation. The force constants were set is large enough to strongly bind the walkers to the actual CV positions. The hill size was diminished to (h, σ 1 , σ 2 , σ 3 ) = (0.2 kcal/mol, 0.08 Å, 0.033, 9 • ). A new hill was deposited every time (τ ) when a walker was displaced by 1.5|σ|. The fictitious mass was where τ (= 10 fs) is the characteristic time for fictitious particle to move 1.5 times the hill width and σ i is the hill width of CV. With this setup, the diffusibility of the walkers was controlled such that the average interval of hill deposition was 10.6 fs. The slow rate at which the hills were deposited ensures a sufficiently long simulation and, thus, a correct mapping of the free energy surfaces associated with the reactions, which gives access to the thermodynamically preferred reaction pathways and properly evaluated activation free energies.
To explore the free energy surface of the reactive region only (d < 2.38 Å), an artificial wall potential was added to V hills . With this prescription we were able to prevent the computation of MTD simulations from being too expensive and from being unstable in statistics.
However, the conformational change of PTO molecule was not fully covered within d < 2.38 Å. As shown below, the contracted conformations were in this range while the extended conformation was not. To remove this artifact, we used the blue moon ensemble simulation which is described next.

Blue moon ensemble
The blue moon ensemble simulations were used to compute the free energy difference with respect to the conformational change of PTO within the range 1.59 Å ≤ d ≤ 5.11 Å. The reactive region was not studied in the blue moon ensemble simulations because the reaction cannot be described by a single dimension with respect to d.
The pseudo-Hamiltonian of blue moon ensemble, in which the system is constrained to d =s(r), can be expressed as where the force constant was set to a large value, κ = 5600 kcal/mol Å −2 , so as to keep the actual CV,s(r), close to the target value, d. The free energy value at d = d α , was evaluated as a numerical integration by the trapezoidal rule where the thermal averages ⟨· · · ⟩ α are computed at constraint points of 1 ≤ α ≤ 12 with equal separations, d α = 1.59 Å +(α − 1) × 0.32 Å.

Results and discussion
Metadynamics Three-dimensional free energy profiles from DFT metadynamics are displayed in Figures   2 and 3 for the pathways leading to THFA (5PTO2 and 5PTO5) and 3-HTHP (6PTO1 and 6PTO5), respectively. (The corresponding results from DFTB metadynamics are also shown in Figure S1 and S2, respectively.) Each of these free energy profiles was obtained from the statistical average of four MTD runs. Gaining sufficient statistics helped us obtain qualitatively correct features of the free energy surfaces, although some roughness remains as an intrinsic error from the metadynamics sampling. (We will show that the errors of the reaction barrier later in Table 3 are about a few kcal/mol.) In Figures 2 and 3  At the transition state (when crossing the d = 0 line), we can see that n ≃ 2.0, which means that one of the hydroxyl groups is always protonated, and φ is in between 180 • and 270 • , which means that the backbone of PTO molecule is always in a specific conformation.
The protonation of the hydroxyl group in the process is in line with experimental evidence that the reaction rate is accelerated in acidic conditions by carbonated HTW. 39 Inspecting the trajectory, we find that the PTO conformation of 180 • < φ < 270 • allows the alignment of O· · · C-O atoms to enable the bond exchange via an S N 2 process. We did not find any free energy minima in the protonated region (n ∼ 2) before the bond exchange (d > 0). This means there are no stable protonated intermediates of PTO, and thus the protonation and bond exchange proceeds consecutively in a single step. All four paths found in this study, i.e., 5PTO2, 5PTO5, 6PTO1 and 6PTO5, share such a proton assisted S N 2 mechanism. This is also similar to what was found earlier for the dehydration reaction of hexanediol. 32 Therefore, this mechanism should be general for polyalcohol dehydration in HTW under acidic conditions.
In Figures 5 and 6, we show the molecular configurations along the 5PTO2 and 6PTO5 pathways. (The molecular configurations along the 5PTO5 and 6PTO1 pathways are also shown in Figures S6 and S7 MPa). In Figure 7, the origin of the free energy is set to that of the extended conformation.
The free energy changes due to the conformational change are of few kcal/mol, which is smaller than those of the bond breaking by an order of magnitude. However, we find that the free energy change from the extended to the confined conformations does depend on the reaction paths. Since the free energy change corresponds to the reversible work to fold the conformation, the energy change tends to increase as the folding stride, ∆d, increases.
Accordingly, the free energy in the confined conformations along the paths leading to 3-HTHP is larger than those of THFA. The dependence on the respective paths varies from 0.9 to 4.0 kcal/mol, which means that the free energy correction due to the conformational change cannot be neglected for correct estimation of reaction barriers.

Reaction Barrier
The free energy barrier of each reaction path, along with the contributions from metadynamics and blue moon ensemble simulations are summarized in Table 3 . 39 The measured rate constant was assumed to be equal to the rate constant from the transition state theory (TST), where k B is the Boltzmann constant, h is the Planck constant, and T is the temperature, 573 K. The free energy barriers thus estimated were ∆A exptl = 41.7 and 43.9 kcal/mol, respectively, for THFA and 3-HTHP.
The computed free energy barriers for the reaction path leading to THFA (33 kcal/mol for DFT and 33-35 kcal/mol for DFTB) are lower than those leading to 3-HTHP (35-36 kcal/mol for DFT and 37-43 kcal/mol for DFTB). Although the computed free energies are systematically lower than the experimental ones, both the results of DFT and DFTB support the experimental finding that THFA and 3-HTHP are the main and by-products, respectively. The difference between the computed free energy barriers between THFA and 3-HTHP products are 6 kcal/mol in DFT and 4 kcal/mol in DFTB, which are larger than those estimated from experiment, 2.2 kcal/mol. Possible origins for the discrepancy between the simulation and the experiment are multiple; it could be ascribed to the DFT and DFTB methods, the high pH condition in the computation, or the experimental estimation assuming the TST.
Still, the simulations provide some insights on a qualitative basis. It is estimated that the free energy barriers of two paths leading to THFA (5PTO2 and 5PTO5) are close to each other, 33 and 33-35 kcal/mol in DFT and DFTB, respectively. This suggests that reaction paths 5PTO2 and 5PTO5, which lead to S-THFA and R-THFA, are almost equally probable.
From this result it is predicted that the main product THFA would be a racemic mixture.
Meanwhile, the paths leading to 3-HTHP (6PTO1 and 6PTO5) are both stereoselective since they were S N 2 reactions. Therefore, for a 3-HTHP pathway, S-PTO leads to S-3-HTHP and R-PTO leads to R-3-HTHP. In fact, our previous studies on hexanediol (without the 3-OH) have shown that the reaction is stereoselective. 32

Reaction Path Analysis
To understand the origin of the difference in free energy barriers, the minimum free energy paths (MFEPs) were analyzed. The MFEPs were obtained in the following manner from the free energy surfaces of DFT metadynamics simulations. The three-dimensional free energy surfaces, which are shown in Figures 2 and 3, were cut into slices with constant d values.
Within two-dimensional free energy surface of each slice, we found the point (n mfe (d), φ mfe (d)) where the free energy A mfe (d) becomes a minimum. The physical meaning of the free energy minimum is the most probable state, in this case, among the states from the canonical ensemble of molecular configurations. Therefore, given a d value, (n mfe (d), φ mfe (d)) are the most probable (n, φ) values that are found in the canonical ensemble.
In Figure 8(a), we display the free energy curves A mfe (d) along the MFEP, together with the free energy curves obtained from MTD and BME in the range 0 Å < d < 5.1 Å. The origin of these curves are aligned in the extended conformation of PTO. We can see that the free energy curves leading to THFA and 3-HTHP start to deviate below the range d ≃ 1.7 Å and above A ≃ 10 kcal/mol. In Figure 8(b), we display the plots of the coordination number along the MFEP n mfe (d) and A mfe (d). We find that the plots are very different between the MFEP leading to THFA and 3-HTHP above A ≃ 10 kcal/mol; the most probable of the protonation occurs at A ≃ 10 kcal/mol and A ≃ 30 kcal/mol in the cases of THFA and 3-HTHP reactions, respectively. Therefore, from the MFEP analysis we speculate that the timing of protonation before the C-O bond breaking and formation process may be causing the difference in the free energy barrier of the reaction leading to THFA and 3-HTHP. (In contrast, the MFEP of the DFTB calculation did not show much difference between the THFA and 3-HTHP pathways, see Figure S4; this may be ascribed to the failure of DFTB method to describe properly the change in proton affinities of PTO at different configurations.)

Conclusions
The competing reactions of PTO dehydration in HTW were studied by MTD and BME simulations. The results obtained from the DFT-and DFTB-based simulations indicate consistently that the free energy barriers leading to the five membered ether product, THFA, are lower than those leading to the six membered ether product, 3-HTHP. This is in accordance with experimental evidence that THFA and 3-HTHP are the main product and the byproduct respectively. It was found that the proton-assisted S N 2 process, in which the protonation of hydroxyl group and the C-O bond exchange occurs in a single step, is the most dominant mechanism. As this process was seen in all these four reaction paths studied, it is expected that this is common in polyalcohol dehydration in HTW under acidic conditions. The sim-ulations clarified the importance of hydrogen bond network of water enabling fast proton transfer such that it facilitates the protonation of hydroxyl group. This clearly explains why the reaction is accelerated in carbonated HTW. From the minimum free energy path analysis of DFT-based simulations, it was found that the reaction mechanism were slightly different between THFA and 3-HTHP reactions in the timing of the protonation of hydroxyl group, which could be a source of difference in their reaction rates.
The simulation predicts that the free energy barriers of the two S N 2 reaction paths from S-PTO leading to R-and S-TFHA are similar. Thus it is expected that the TFHA product should be racemized. On the other hand, S-PTO should mostly lead to S-3-HTHP since the stereoselectivity is kept in the two S N 2 reaction paths.
The results obtained in this work rely on the MTD simulations which detect a reaction path with lowest free energy barrier for a given range of CV space. Nothing is said so far about the possibility of the reactive process leading to a minor product which should, basically, experience a reaction path overcoming a higher free energy barrier. In fact, experiments found a deterioration of selectivity with minor products, such as polymers that arise presumably via inter-molecular dehydration. One of the next challenges in theoretical studies of polyalcohol dehydration in water would be to identify the minor product and understand the reaction mechanism leading to minor products.      : Three-dimensional free energy surface of 5PTO2 (top panels) and 5PTO5 (bottom panels) with respect to d, φ and n, obtained from DFT metadynamics simulations. The side views of free energy surfaces are displayed in a direction perpendicular to d, n (left panels) and parallel to d, φ (right panels). They are shown by different transparent colors with increasingly leveling contours, from black (confined conformation) to white (TS).