The reaction mechanism of polyalcohol dehydration in hot pressurized water †

The use of high-temperature liquid water (HTW) as a reaction medium is a very promising technology in the field of green chemistry. In order to fully exploit this technology, it is crucial to unravel the reaction mechanisms of the processes carried out in HTW. In this work, the reaction mechanism of 2,5-hexanediol dehydration in HTW has been studied by means of three diﬀerent ab initio simulations: the string method, metadynamics and molecular dynamics in real time. It is found that the whole reaction involving protonation, bond exchange and deprotonation occurs in a single step without a stable intermediate. The hydrogen bonded network of the surrounding water has a vital role in assisting an eﬃcient proton relay at the beginning and at the end of the reaction. It is confirmed that the reaction is energetically most favorable in the S N 2 pathway with an estimated barrier of 36 kcal mol (cid:2) 1 , which explains the high stereoselectivity and the reaction rate observed in experiment. The mechanistic insights provided by our study are relevant for a prominent class of reactions in the context of sustainable biomass processing, namely dehydration reactions of polyalcohol molecules.


Introduction
2][3][4] Indeed, water offers many advantages over organic solvents because it is non-toxic, inexpensive, nonflammable and readily available. 2,5Moreover, the solubility of small organic compounds in high-temperature liquid water (HTW) is largely enhanced with respect to the solubility in room-temperature water. 2 Hence, the use of HTW holds great promise in the design of cleaner, safer and environmentally benign chemical processes.
The reaction rates and product selectivities of the processes carried out in HTW can be very sensitive to temperature and water density or pressure, 1,2 which means that these process variables can be tuned to fully harness the special properties of HTW.It is clear that a detailed understanding of the mechanism and the role played by water in reactions is of paramount importance in order to find the optimum conditions for a given chemical process.In broader terms, this detailed understanding can also have relevant implications in the context of the origin of life 6,7 and in the context of energy and fuels. 8However, elucidation of the mechanism and the corresponding role of water can be a very complex task due to the multifaceted nature of water-substrate interaction, e.g.hydrogen bonding, polarity, acidity and basicity, hydrophobicity, etc. Furthermore the reaction mechanisms in HTW can be very different from the mechanisms operating in organic solvents or in room-temperature water. 9,10It is well known that computational chemistry provides invaluable tools when it comes to unravelling reaction mechanisms.2][13][14][15] Therefore, there is an urgent need to carry out in silico studies in order to better understand how reactions take place in aqueous media under extreme thermodynamic conditions.Here, we present a computational study aimed at unravelling the mechanism of the stereoselective transformation of 2,5-hexanediol (HDO) into dimethyltetrahydrofuran (DMTHF) by means of an intramolecular dehydration process in hot pressurized water. 168][19][20][21] In this context, dehydration of polyalcohol compounds constitutes one of the most important families of reactions. 22,235][26][27][28][29] The addition of CO 2 to HTW leads to the fact that the pH of the aqueous medium lowers when this gas dissolves in water.It has been observed that an increase in the acidity of the aqueous medium can result in the acceleration of dehydration reactions since these transformations are acid-catalyzed processes.The use of CO 2 offers environmental advantages over the use of mineral acids.High-pressure CO 2 /HTW technology is increasingly being used in acid-catalyzed reactions 4,30 and will most likely play a prime role in novel green methodologies of biomass processing. 20,31ne of the simplest dehydrations of polyalcohol compounds in HTW reported to date is the dehydration of HDO to DMTHF (see Fig. 1). 16,32It is experimentally observed that this reaction is stereoselective in most aqueous solutions with strong acids such as HBr although exceptional cases exist such as sulfuric acid. 32The relevance of this reaction lies in the fact that it provided evidence that the CO 2 /HTW technology can be employed to carry out stereoselective dehydration reactions.The experiments performed by Shirai and coworkers 16 showed that the dehydration reaction of 2R,5R-HDO and 2S,5S-HDO in CO 2 -enriched HTW furnishes the cis isomer of DMTHF with a selectivity of more than 85%.Here we employ different computational techniques, i.e. ab initio molecular dynamics simulations, metadynamics, determination of minimum energy paths by means of the string method, applied to different model systems, including microsolvation of an HDO molecule and full solvation of an HDO molecule in the condensed phase, with the aim of elucidating the mechanism of the dehydration reaction of HDO and providing a rationale to its stereoselectivity.Our simulations demonstrate that the lowest free energy pathway of the reaction corresponds to an S N 2 mechanism, by virtue of which the cleavage of the C-O bond takes place simultaneously with the formation of a new C-O bond.Prior to this bond exchange, one of the alcohol groups, which is the leaving group in the S N 2 reaction, gets protonated by means of a Grotthuss relay mechanism. 33,34Yet the corresponding protonated HDO molecule is not a stable intermediate.Concerning the role of water, our simulations show that an increase in the number of water molecules in the first solvation shell of HDO results in an increase of the rate of the dehydration reaction.Overall, the details of the reaction mechanism at the atomistic level revealed by our simulations provide valuable insights not only for the specific reaction considered herein but also for the whole family of dehydration processes carried out in high-temperature liquid water.

Method
This study has been done in three stages.At the first stage we calculated the minimum energy paths (MEPs) of the reaction for the microsolvated model which consists of an 2R,5R-HDO molecule and water molecules with an extra proton, i.e.HDO(H 2 O) x H + (with x = 2, 4 and 6) under the free boundary conditions.To obtain the MEPs for the reactant and product minima chosen appropriately, we employ the string method 35,36 which is able to optimize the position of discrete images arranged along the path.This provides a rough estimate of the energy profile and the structural changes during the reaction.For comparison, the MEP calculation is done in the same way for the case where the proton is absent, i.e.][39][40][41] At the second stage we calculated the free energy landscape with respect to three collective variables (CV), d, f and n, which are found to characterize the reaction in the first stage; d = r 28 À r 27 is the difference in bond distances between C 2 -O 8 and C 2 -O 7 atoms, f is the dihedral angle with respect to the O 7 -C 2 -C 3 -C 4 atoms, and n = n 7 is the hydrogen coordination number of the O 7 atom (see Fig. 2 and eqn (1) below for the definitions of atomic numbering and the coordination number).For this purpose, we used metadynamics simulations [42][43][44] of the microsolvated model mentioned above, as well as the aqueous solution of HDO where periodic boundary conditions are applied.In metadynamics the free energy basin is filled with a collection of history-dependent Gaussian hills accumulated where the system has sampled during the simulation.The free energy landscape obtained provides information on thermodynamic ensembles that belong to the reactant minimum and the transition states.
At the third stage we carried out molecular dynamics simulations of the aqueous solution running forwards and backwards starting from the structures near the transition state, which are taken from the second stage.This provides real time information of typical trajectories during the whole reaction process, albeit they were limited in number due to computational expense.
For the microsolvated model all the calculations were carried out using the PIMD code 45,46 combined with Turbomole. 47For the electronic structure calculations, we employed the resolutionof-identity density functional theory (RIDFT) 48 with the PBE exchange correlation functional 49 using the SV(P) basis set. 50We also employed the resolution-of-identity second-order Møller-Plesset perturbation theory (RIMP2) 51 using the SV(P) basis set 50 in the MEP calculations.The MEPs described in terms of 120  discrete images from the HDO reactant to the DMTHF product were obtained using the string method.The three-dimensional metadynamics was carried out based on the Born-Oppenheimer scheme in the NVT ensemble with Nose ´-Hoover chain thermostats [52][53][54] implemented in the PIMD code following the setup by Ensing et al. 43 The temperature was controlled at 523 K, with the time step of Dt = 0.25 fs.The size of Gaussian hills was gradually diminished as the reaction approached completion, and the final height (h) and the widths (w) were set to (h, w d , w f , w n ) = (1.0 kcal mol À1 , 0.09 Å, 5.01, 0.05).For a computational speedup, we used 12 walkers and we restricted the CV space to d o 1.85 Å and 132.51 o f o 227.51.The metadynamics simulations were run for 34 ps and 57 ps for the cases of x = 4 and x = 6, respectively, and on average a new hill was added every 15 fs for each walker.
For the aqueous solution, all the calculations were carried out using the CPMD code. 55The system contained an HDO molecule, 70 water molecules and a proton in a cubic box with a fixed side length of 14.93 Å, which was estimated from a preliminary classical MD simulation in the NPT ensemble under the experimental conditions of 20 MPa and 523 K. 16 A uniform background charge was added to maintain the charge neutrality of the system under the periodic boundary conditions.A similar setup was employed in previous ab initio MD simulations for systems with an excess proton, albeit the number of water molecules used therein (32, 56 36 9 ) were less than the present one (70).The main limitation arising from the system size is that the acidity (pH = 0.3) is set much lower than the experiment (pH = 4). 16The electronic structure calculations were done based on density functional theory with PBE exchange correlation functional 49 and Grimme's D2 van der Waals correction 57 using the cut off at 25 Ryd in the plane wave (PW) basis functions.The core electrons were taken into account using Vanderbilt's ultrasoft pseudopotentials. 58The fictitious mass of the orbitals was set to 400 a.u.The three-dimensional metadynamics was carried out based on the Car-Parrinello scheme 59 in the NVT ensemble using Nose ´-Hoover chain thermostats, [52][53][54] where the temperature is controlled at 523 K, with the time step of Dt = 0.097 fs.The size of Gaussian hills was gradually diminished as the reaction approached completion, and the final height and the widths were set to (h, w d , w f , w n ) = (0.3 kcal mol À1 , 0.005 Å, 0.31, 0.005).For a computational speedup, we used 8 walkers, and we restricted the CV space to d o 2.5 Å.The metadynamics simulations were run for 100 ps, and on average a new hill was added every 44 fs.The twodimensional metadynamics was also carried out in a similar manner with the collective variables d and f being unrestricted in space.To obtain real-time trajectories, we ran molecular dynamics simulations based on the Car-Parrinello scheme 59 in the NVE ensemble with the initial velocity corresponding to 523 K, but without any temperature control during the simulation.
The coordination number is defined in two ways as the rational function, and where the index j runs for all the hydrogen and carbon atoms, respectively, in eqn ( 1 3 Results and discussion

Reaction path analysis
We begin by searching for probable reaction paths in the microsolvated model.After many trial and error static calculations, we learned the following things.Firstly, the HDO molecule has several stable conformations, either in the extended form or in the contracted form, as shown in Fig. 3.However, the conformation must be contracted at some point for the reaction to take place.Secondly, for a small microsolvated reactant, HDO(H 2 O) x H + , to be stable the extra proton must be hydrated by a sufficient amount of water molecules.In particular, the stable structures were found when the extra proton is absorbed into a water molecule as the hydronium ion, which is hydrogen bonded to three other water molecules, i.e.
In fact as this structure is not available in the case of x = 2, its most stable form resulted in the protonated HDO molecule, failing to describe a stable state of proton in water.Thus at least x = 4 water molecules are required to properly describe the HDO molecule in a stable reactant state.This also implies that at least four water molecules are explicitly involved in the reaction process.
Based on these observations, we next looked for reaction paths of the x = 4 and x = 6 microsolvated models that start from a configuration with the H 3 O + (H 2 O) 3 cluster included using the string method.Here we show the PBE results for the x = 6 case, as they were qualitatively similar to the PBE results for the x = 4 case (Fig. S1 in ESI †) and the RIMP2 results of the x = 6 case (Fig. S2 in ESI †).In Fig. 4, we show two types of minimum energy paths, Channels 1 and 2, that start from different conformations of HDO.In Channel 1, Structure I has an internal hydrogen bond with the dihedral angle being f = 601 (cf.Fig. 3(c)), while in Channel 2 the hydrogen bond is broken with f = 1801 (cf.Fig. 3(d)).Investigating the process near Structure III in Fig. 4, the C 2 -O 7 bond breaking and the C 2 -O 8 bond forming are found to be stepwise and concerted, respectively, which is the signature of S N 1 and S N 2 mechanisms, respectively.Thus, let us refer to them hereafter as the S N 1-and S N 2-pathways.
There are some aspects shared among S N 1 and S N 2 pathways.A proton is received before (Structure II) and released after (IV) the O-C bond breaking/forming, and thus the reactions take place via the proton-assisted mechanism, in accordance with the experimental suggestion. 16In the calculation it was found that this process performs efficiently by means of Grotthuss proton relay, showing that the hydrogen bond network between water molecules plays a significant role in the reaction.The reaction not only starts in (Structure I) but also ends in (V) in the H 3 O + (H 2 O) 3 cluster to recover the stable form of the hydrated proton.
The energy profile and the structure variables along the MEPs of the x = 6 microsolvated model are shown in Fig. 5.In Fig. 5(a) we can see that the energy profiles of both the S N 1-and S N 2-pathways are single peaked at their respective transition states (Structure III), corresponding to the O-C bond exchange configuration.This indicates that they are single step processes starting from the relevant HDO conformations.There is no stable intermediate before and after the O-C bond exchange, which is the main rate-determining step.However, a crucial difference is that the energy barrier of the S N 2-pathway, 40 kcal mol À1 , is much lower than that of the S N 1-pathway, 57 kcal mol À1 .This indicates that the S N 2-pathway is dominant, and therefore the reaction is stereoselective.
Although all the tendencies mentioned above are unchanged for the x = 4 microsolvated model, the energy barriers of the S N 1and S N 2-pathways are found to be 58 and 41 kcal mol À1 , respectively, which are both slightly higher compared to the ones of the x = 6 microsolvated model.This implies that both reactions Fig. 4 Structures along the minimum energy paths, Channels 1 and 2, obtained from the PBE/SV(P) calculations.Color codes: oxygen atoms in red, hydrogen atoms in white, carbon atoms in grey.As an exception, the oxygen atom bonded to an excess proton is marked in blue.The numbers I-V correspond to the positions in Fig. 5.The change in O-C bond distances, r 28 and r 27 , their difference d, the coordination number n, and the dihedral angle f along the MEP are displayed in Fig. 5(b)-(e).This supports the proton assisted S N 1-and S N 2-mechanisms mentioned above.For instance, the sum of r 28 and r 27 values is found much larger in the S N 1pathway (6.6 Å) than that of the S N 2-pathway (4.3 Å), since a shortlived carbocation is formed in the S N 1 reaction.
We also studied the water-assisted pathway of this reaction in the absence of a proton by carrying out the MEP calculation of the microcluster model HDO(H 2 O) 6 (see Fig. S3 in ESI † for details).The energy barrier is found to be 57 kcal mol À1 , which is comparable with an earlier study of QM/MM calculation on the water-assisted pathway of 1,4-butanediol dehydration reaction in hot water, i.e. 53 and 62 kcal mol À1 , for different pathways in the absence of a proton. 38However, the energy barrier is much larger than that of the HDO(H 2 O) 6 H + , 40 kcal mol À1 , in the presence of an excess proton, suggesting that the water-assisted pathway is not the dominant one in this case.

Free energy analysis
As mentioned above, the variables d, n, and f are able to distinguish the reaction processes in terms of bond exchange, protonation, and S N 1/S N 2 characterization, respectively.Setting them to be the collective variables, three-dimensional metadynamics simulations of HDO aqueous solution were carried out at the temperature of 523 K (see Fig. S4 in ESI † for details).The free energy profiles were calculated from the basin of the HDO reactant to the transition state.
In Fig. 6, we display the free energy surface with respect to d, n and f.Fig. 6(a) shows that the free energy minimum is located at about (n = 1, d 4 1.5 Å) corresponding to the ensemble of structures in which the O 7 H group is neutral and the C 2 Á Á ÁO 8 distance is well separated, i.e. the original HDO molecule and a proton left in bulk water.
Meanwhile the transition state can be identified when we detect an event of the O-C bond exchange crossing the d = 0 line.Such a transition state is found around (n = 2, d = 0), where the O 7 atom is bonded to two hydrogen atoms, forming an OH 2 + group.
An important outcome is that the O-C bond exchange always occurs after protonation.The region of O-C bond exchange without protonation, (n = 1, d = 0), has never been visited because the free energy therein is too high.From the free energy profiles in Fig. 6(a), we also see that there is no minimum corresponding to the protonated HDO molecule, which is consistent with the finding of the previous MEP analysis.Therefore, protonation and O-C bond breaking are consecutive processes.Note that by construction metadynamics always finds the reaction channel with the smallest free energy barrier.
Analyzing the trajectory of all the walkers in metadynamics, it was found that the O 7 H group is protonated by an extra proton in water.We were unable to detect either an OH À intermediate in water or an O 8 À intermediate in HDO.Neither the unprotonated water molecules nor the OH group of the HDO molecule are the main source of protonation in this case.This finding is consistent with the result from the MEP calculations mentioned above that the water-assisted pathway is not a dominant one.Therefore we conclude that the proton is catalytic in this system, which must be relevant to the experimental fact that dehydration reactions in hot water can be accelerated by the acid environment, in particular under the subcritical conditions.In Fig. 6 The free energy profile of the x = 4 and x = 6 microsolvated models is shown in Fig. 7(a) and (b), respectively (also see Fig. S5 in the ESI † for another view of this free energy profile).In this case, the free energy has been obtained with a constraint of 132.51 o f o 227.51 focusing on the important region for the S N 2 reaction.It could be seen that the free energy profile is qualitatively similar to the result of the aqueous solution, Fig. 6.In Table 1, we summarize the reaction barriers in the energy of MEP and in the free energy at 523 K, with respect to the x = 4      and x = 6 microsolvated models and aqueous solutions.In the x = 4 and x = 6 microsolvated models, the energy and free energy barriers of the S N 2 reaction obtained are similar in value.This supports the reliability of metadynamics and that the set of collective variables (d, f, n) has been a reasonable choice.We can see that the barrier height is slightly decreased when the number of water molecules surrounding the HDO reactant is increased, which implies that the aqueous environment helps in accelerating the reaction.This is true not only for the S N 2-pathway but also for the S N 1 pathway.We can also see that the energy barrier is larger for the backward reaction, which is consistent with the fact that the backward reaction is not found in experiment.
The free energy landscape shown in Fig. 6 has covered only the reactive region within the contracted HDO configuration (d o 2.5 Å).In order to study the free energy landscape of the region covering the conformational change, another metadynamics simulation of HDO aqueous solution was carried out in two dimensions, d and f (with The result is shown in Fig. 8.The free energy minima found at (d, f) = (3.0Å, 3001), (3.3 Å, 601), and (3.5 Å, 1801) correspond to the structures shown in Fig. 3(a), (b) and (c), respectively.This indicates that the extended conformations of HDO are the most stable ones in aqueous solution.Compared with the extended conformation at (d, f) = (3.5 Å, 1801), the contracted conformation with the internal hydrogen bond, (d, f) = (1.9Å, 601) or (2.1 Å, 3001), is less stable by about 3 kcal mol À1 .Meanwhile, from Fig. 6 the free energy in the S N 2 transition state is 33 kcal mol À1 more than that of the HDO reactant in the contracted conformation with the internal hydrogen bond.The sum of these free energy differences, 36 kcal mol À1 , is the final estimate of the S N 2 reaction barrier in the aqueous solution.This is in reasonable agreement with the experimental estimation of the free energy barrier, DA = 39 kcal mol À1 , which is obtained by combining the experimental rate constant (with the half life t h E 30 min at T = 523 K) 16 and the transition state theory: where h and k B are Planck's constant and Boltzmann's constant, respectively.

Trajectory analysis
We Meanwhile, the behavior of the coordination numbers, N 7 and N 8 , is noticeably different from the one of the S N 2 MEPs, n 7 and n 8 .In the real-time trajectories, we find a large amount of thermal vibration that allows the proton to jump back and forth between the HDO/DMTHF molecule and its neighboring water molecules, as can be seen in the snapshots of Fig. 10 and 11.This process can last for ps time scale until the proton is finally absorbed from/released into the bulk water via the proton relay.This occurs when the hydrogen bond network of water is transformed such that H 3 O + is surrounded by three other water molecules.In light of these trajectories and others shown in the ESI, † the timing of the creation of such a hydrogen bond network seems to be a statistical event.

Conclusions
In this paper the dehydration reaction of HDO in hot acidic water has been studied theoretically in three different analyses, i.e.MEP analysis by the string method, free energy analysis by metaydnamics simulations, and trajectory analysis by molecular dynamics simulations.The former two have been performed in the protonated microsolvated model, while the latter two have been performed in aqueous solution.All these computational results consistently indicate that the S N 2 pathway is more favorable than the S N 1 pathway, and thus supports the recent experimental result that this reaction proceeds with a high stereoselectivity in acidic water. 16The qualitative mechanism of the S N 2 pathway is schematically shown in Fig. 13, with an emphasis on how the water molecules are involved in the reaction.The free energy barrier for the S N 2 reaction is estimated to be 36 kcal mol À1 in the aqueous solution, which is in reasonable agreement with a value converted from the experimental rate constant combined with the transition state theory, 39 kcal mol À1 .
The results obtained in these calculations do not contradict in any way to the reaction mechanism suggested in the experimental paper. 16However, our simulations have been able to add useful insights into the reaction mechanism with respect to the change in the energies and the structures step by step for each of the reaction processes.The reaction mechanism along with the respective time scales is summarized below.
1.The HDO molecule can fluctuate back and forth between different conformations in HTW.However, the extended conformation is more favorable than the contracted conformation from the viewpoint of free energy.
2. Few ps before of arriving at the transition state, the HDO molecule forms the contracted conformation.
3. At times less than a few ps before arriving at the transition state, one of the OH groups of HDO becomes protonated by the surrounding water via the Grotthuss relay mechanism.4. At the same time as step 2 or 3, the HDO molecule must form the f = 1801 conformation.
5. A bond is exchanged between O-C and OÁ Á ÁC at the same time in the S N 2 fashion, forming the ether ring.This occurs very quickly within a few tens of fs.This is the transition state of the whole reaction.
6.The bond angles are slightly shifted to adapt to the stable protonated DMTHF molecule.
7. The proton is picked up by water again in the Grotthuss relay, finalizing the reaction.It takes up to several ps for a proton to be released into bulk water after departing from the transition state.
The reaction proceeds in a single step.There is no stable intermediate throughout the reaction.
Finally it is important to comment on the role of water.On the quantitative basis, it is found that the energy barrier of the transition state is reduced as hydrogen bonded water molecules are increased, thus suggesting that water provides a field that accelerates the reaction.We provide a qualitative explanation on the role of water as an efficient solvent assisting the reaction.
The reaction proceeds only by way of the protonated HDO intermediate, capturing one proton at the early stage and releasing one proton at the last stage, wherein the proton must be stablized.Since the proton in the bulk can be stablized as a hydronium ion when it is hydrogen bonded to three water molecules, at least four water molecules are necessary to describe the whole reaction process.
As the proton is captured by the HDO reactant via the Grotthuss relay mechanism, the proper creation of a hydrogen bonded network is important.Once the HDO molecule is protonated it survives for a few ps although not being a stable intermediate.This is because for the protonated HDO molecule to lose the proton the hydrogen bond network must be rebuilt again.As the protonated HDO state lasts for more time, it provides more chance for an incoming energy to have the O-C bond elongated.If the system is in the contracted conformation, this immediately leads to the reaction.Once this happens, the extra proton is released into bulk water almost spontaneously without the need for energy, although it takes a few ps to form the hydrogen bond network.
The mechanistic insights offered by our work are not only relevant for the stereoselective dehydration reaction of HDO but also for the whole family of dehydration of polyalcohols in hot acidic water.We hope that our results will facilitate the design of more efficient and sustainable processes in HTW aimed at processing biomass derivatives.supercomputer in Center for Computational Science and E-systems at Tokai Research Establishment, JAEA.We are grateful to Dr A. Malins in JAEA for proofreading the manuscript of this paper.We also thank Prof. T. Sasaki at The University of Tokyo for useful discussions.
) and(2).The former was used in metadynamics simulations describing the number of O-H bonds, while the latter was used in the trajectory analysis of molecular dynamics simulations describing the sum of O-H and O-C bonds for a given oxygen atom i.The O-H and O-C bond length parameters were set to be d OH ij = 1.4 Å and d OC ij = 1.65 Å, respectively.

Fig. 5
Fig. 5 The change along the minimum energy paths, Channels 1 (blue) and 2 (red), obtained from the PBE/SV(P) calculations.(a) The relative energy, and the structural variables (b) f, (c) n 7 (solid line) and n 8 (dashed line), (d) r 27 (dashed line) and r 28 (solid line), and (e) d.The numbers I-V correspond to the ones in Fig. 4.

Fig. 6 Fig. 7
Fig. 6 Three-dimensional free energy surface of HDO aqueous solution with respect to d, f and n 7 , obtained from the PBE/PW metadynamics calculations.The views are in a direction (a) parallel to f and (b) parallel to n 7 .
(b), we show another side view of the free energy profile with respect to the variables f and d.In this case, three minima are found at f = 3001, 601 and 1801 with d E 2 Å.The minima at f = 601, 3001 and f = 1801 correspond to the thermodynamic ensemble of structures where the HDO molecule is in a contracted conformation with and without the internal hydrogen bond, respectively.They are comparable with the structures shown in Fig. 3(d) and (e), respectively.The transition state around (f = 1801, d = 0 Å) is characteristic of the S N 2 reaction.Combining the results of Fig. 6(a) and (b), it can be concluded that the reaction only occurs when n = 2 and f = 1801 is realized at the same time.The S N 2 reaction always undergoes via the protonated state of HDO in the conformation without the internal hydrogen bond.The reaction proceeds neither for the unprotonated state in this conformation nor the protonated state in any other conformation.

a
The sum of contributions from the conformational change and reaction.b Ref. 16. c Estimation combined with transition state theory.

Fig. 8
Fig. 8 Two dimensional free energy surface of HDO aqueous solution with respect to f and d, obtained from the PBE/PW metadynamics calculations.The focus is on the region of conformational changes of the HDO molecule.

Fig. 9
Fig. 9 The change of structural variables N 8 (left panel), r 27 and r 28 (upper right panel), and f (lower right panel) along the forward trajectory shown in Fig. 11.The color of N 8 indicates the position of an excess proton.It is red when O 8 is protonated, green when a water molecule hydrogen bonded to HDO is protonated, and blue otherwise.

Fig. 10
Fig. 10 Snapshots of a backward trajectory of aqueous solution from the transition state to the reactant HDO obtained from the PBE/PW molecular dynamics calculations of aqueous solution.The water molecules are shown only for those close to HDO.Color codes: oxygen atoms in red, hydrogen atoms in white, carbon atoms in grey.As exceptions, the oxygen atom bonded to the excess proton is marked in blue, and the hydrogen atoms that are hydrogen bonded to other water molecules are marked in yellow.

Fig. 11
Fig. 11 Snapshots of a forward trajectory of aqueous solution from the transition state to the product DMTHF obtained from the PBE/PW molecular dynamics calculations of aqueous solution.The water molecules are shown only for those close to DMTHF.Color codes: oxygen atoms in red, hydrogen atoms in white, carbon atoms in grey.As exceptions, the oxygen atom bonded to the excess proton is marked in blue, and the hydrogen atoms that are hydrogen bonded to other water molecules are marked in yellow.
studied the reaction process in real time by means of molecular dynamics simulations launched from the transition state region of the S N 2 reaction.The trajectories that led to either the reactant HDO or the product DMTHF correspond to the ones that have evolved backward and forward in time, due to time reversibility.The typical examples for such backward and forward trajectories are shown below.The results of other trajectories calculated (five backward and three forward) are shown in Fig. S6-S13 in the ESI.† Fig. 10 and 12 are the snapshots and the structural variables N x , r 27 , r 28 and f, respectively, along the backward trajectory.Fig. 11 and 9 show the respective results along a forward trajectory.Fig. 12 and 9 show the time evolution of the bond lengths and the dihedral angles, which start from the values of the transition state, r 27 E r 28 E 2.3 Å and f E 1801.The time evolutions near the transition state seem to resemble the changes along the S N 2 MEP shown in Fig. 5.The bondlengths quickly relax to the values of the reactant HDO (r 27 E 1.4 Å, r 28 4 3 Å) and the product DMTHF (r 27 4 3 Å, r 28 E 1.4 Å) in the backward and forward trajectories, respectively, within the sub-ps scale.This indicates that the bond exchange from CÁ Á ÁO-C to C-OÁ Á ÁC proceeds on this time scale.

Fig. 12
Fig. 12 The change of structural variables N 7 (top panel), f (middle panel), and r 27 and r 28 (bottom panel) along the backward trajectory shown in Fig. 10.The color of N 7 indicates the position of an excess proton.It is red when O 7 is protonated, green when a water molecule hydrogen bonded to HDO is protonated, and blue otherwise.

Table 1
The S N 2 activation energies and free energies in kcal mol À1