Changjie Lu ,Weiqiang Tang ,Zijiang Dou ,Peng Xie ,Xiaofei Xu,* ,Shuangliang Zhao,,*

1 State Key Laboratory of Chemical Engineering and School of Chemical Engineering,East China University of Science and Technology,Shanghai 200237,China

2 Guangxi Key Laboratory of Petrochemical Resource Processing and Process Intensification Technology and School of Chemistry and Chemical Engineering,Guangxi University,Nanning 530004,China

Keywords:Solvent effect Reaction density functional theory Pyruvic acid Tautomerism Isomerization

ABSTRACT It is important to study the solvent effect on keto-enol tautomerism that has applications in many areas of chemical engineering.In this work,we use a multiscale reaction density functional theory (RxDFT) to study the keto-enol tautomerism and isomerization of pyruvic acid.The results show that both effects of solvation and water assistance could reduce the reaction barriers.The water molecule participates the reaction as a catalyst to accept/give the protons with forming a hexagonal ring in the transition state.As a result of this temporary and intermediate hexagonal ring,the solute configuration undergoes a small variation during the reaction,giving a diminished contribution to the intrinsic reaction free energy.The solvent distribution shows a local ordering behavior near the solute that also reduces the contribution of solvation effect to the reaction barrier.Water assistance plays a major role in both pre-reaction and postreaction process.In terms of the driving force for the reaction,the effects of both solvation and water assistance are important.

1.Introduction

It is important to explore how to control chemical reaction by regulating solvent conditions,which is an important topic in green chemistry.Solvents have a great impact on many aspects of reaction,such as reaction barrier,rates and even the mechanism [1–3].Keto-enol tautomerism (KET) is a chemical equilibrium between a keto form (＞HC―C(=O) and an enol form (＞C=C(OH)―) [4].The reaction involves the transfer of protons and the reorganization of bonding electrons.KET is sensitive to solvent conditions [1,5],such as catalyst,solvent quality,pH value,and ion concentration,etc.KET has important applications in many areas of chemical process [6,7].For example,phosphoenolpyruvate,an important intermediate in glycometabolism,has a high potential to transfer phosphate and energy in cell [8].The chemical process is just a KET of phosphorylated compound with the reactions of enzymes and ions[8].There are great significances if we could control the tautomerism by solvent conditions.

In a solvent,KET involves many degrees of freedom of the system.Description of KET in molecular scale is a challenging problem.There are three categories of theoretical methods.The first one is the zeroth-order approximation of solvent effect [2,3].The solvent is treated by continuum model.The solvent effect is accounted implicitly by electrostatic interaction with a dielectric constant.However,it may arise qualitative problems for describing the solvent effect.The second one is the hybrid quantum mechanism/molecular mechanics (QM/MM) [9]that couples classical mechanics with quantum mechanics.QM/MM takes the solvent effect directly into account and can quantitatively describe the solvent interaction by multiscale model during chemical reaction[10].But the applications of QM/MM method are subject to the restrictions of complex molecular model and time-consuming problem.The third one is the quantum mechanics coupled with classical statistical mechanics [11–15]which has derived a series of calculation methods.One of the most prominent and popular methods in this category is self-consistent field (SCF) method by the reference interaction site model (RISM) [16].Nevertheless,RISM-SCF as a method of using integral equation,which cannot to describe system dynamics.In those studies of the first two categories,simple models including a few water molecules were used to mimic the solvent condition.This kind of model cannot describe the complete physics of tautomerism.For instance,protons in KET could transfer over long distance accomplishing via hydrogenbonded networks [17].Today,it is still impossible for computer simulations to give a good description for this kind of large-scale tautomerism by the first two methods.

QM coupled classical density-functional theory (CDFT) could describe complex behaviors of KET.Recently,we proposed a multiscaled reaction density functional theory (RxDFT) [18,19],which combines quantum density functional theory (QDFT) and CDFT.RxDFT could predict the reaction barrier and molecular configurations by investigating several typical states (such as the initial,transition and final states)during the chemical reaction.The transport properties are not considered in RxDFT,but it could be described by dynamical DFT(DDFT)[20].The effect of intermolecular interactions and solvent density distribution can be considered in atomic details.While QDFT is a good method to study chemical reactions,it requires great computational resources.RxDFT provides a generic framework to consider effects of both reaction and solvents.It has the virtues of both efficiency and accuracy.The RxDFT method has been successfully applied to investigate the glycine tautomerization reaction pathways in aqueous solution[18]and the complete reaction pathways among different tautomers of 5-amino-1,3,4-oxadiazole-2(3H)-one (AOO) [21].We also predicted the solvent effects and the free energy profiles of the Type I and Type II SN2 reactions in aqueous solution by using the RxDFT method successfully [19].

In this work,we apply the RxDFT method to explore solvent effects on KET and isomerization reactions in pyruvic acid (PA).PA not only plays an important role in bioenergy metabolism,but also is a precursor for the synthesis of a variety of compounds.The additional reason to select PA is that the system of PA is wide-investigated by both experimental and theoretical methods [22,23].The data of molecular conformations and reaction barriers are available in literature.We could compare our results with those data to validate the effectiveness of RxDFT.The solvent effects on KET of PA will also be revisited by considering both direct reaction with a water molecule and indirect effect of solvent.

2.Model and Theory

2.1.Molecular model

Five types of enol isomers and two types of keto isomers are considered in this work;see both Table 1 and Fig.1.The stability of enol isomers increases from structure I to V.Structure B of keto isomer is more stable than structure A.Four cases of reaction are considered:I →II,III →IV,A →V,B →IV.The first two are isomerizations,and the last two are keto-enol tautomerisms.

We only consider the reaction of single solute molecule at T=298.15 K.Denote the solute molecule by M0for convenience.M0is fixed at the center of a cubic box with length of 50×10-10m.The box has a periodic boundary in all directions in contact with a bulk solvent reservoir.The solvent effect on every reaction is explored by comparing the results with solvent-free case.For solvent-on case,there is an environment of solvent molecules surrounding the solute molecule [24].The surrounding environment(solvent distribution) is related to a bulk water and the external potential of the solute.The data of bulk water are given in Section 2.4.

The solvent distribution is obtained by CDFT with regarding the solute’s potential as an external potential Vext(r,Ω) to solvent molecules.r is the position of beads,and Ω describes orientation of the solvent molecule.Vext(r,Ω)also serves as a boundary condition for the solute molecule.It is responsible for connecting the quantum properties with solvent effects in the study of chemical reaction.See Section 2.2 and 3 for details.Vext(r,Ω) comprises of a site-site Lennard-Jones (LJ) potential and Coulombic potential[18,19,25].It is given by

where rijis the distance between site i and site j.qiand qjare the partial charge on these sites.The LJ parameters of isomers and transition states in PA are obtained from the OPLS-AA force field database [26,27].The cross items of unlike pairs are evaluated by the Lorentz-Berthelot rule.The charge distribution of PA is evaluated by fitting the atomic charge data to molecular electrostatic potential,which can be obtained in quantum mechanical wave function of QDFT.

In order to investigate solvent effect directly,the assisting role of one explicit water molecule(OEWM)in the reactions is explored for both solvent-on (solution phase) and solvent-free (gas phase)cases;see Fig.3.Therefore,there are four cases for all reactions in total:M0in solution phase,M0+OEWM in solution phase,M0in gas phase and M0+OEWM in gas phase.For simplicity,we refer these four cases as ‘‘solution phase”,‘‘solution phase (1 water)”,‘‘gas phase”and ‘‘gas phase (1 water)”,successively.Fig.2 shows the initial (reactant),transition and final (product) states of OEWM-free case,and Fig.3 shows that of OEWM-on case.Note that for OEWM-on case,one hexagon ring can be formed between the OEWM and M0in the transition state;see Fig.3.

Table 1 Gibbs free energies (in Hartree) of pyruvic acid isomers

Fig.1.Schematic views of the enol (I to V) and keto (A and B) isomers of pyruvic acid.The white,gray,and red beads are H,C,and O atoms,respectively.

Fig.2.The optimal structures of initial (reactant),transition and final (product)states for the proton transfer reaction in gas phase.(a)I →II,(b)III →IV,(c)A →V,(d) B →IV.

There are several reasons to use only one single water molecule.(1)The assisting effect of one water molecule to the solute reaction is the major contribution to quantum effect.(2) The system complexity increases with the number of explicit water molecules.In reality,the water molecules form a hydrogen bond network with the solute.It is difficult to study the effect of hydrogen bond network on the chemical reaction.(3) The results of one water molecule could describe how water molecule plays the role in chemical reactions.It is also helpful to understand the real physics in the system of multi water molecules.(4) While most relevant studies are of one water molecule,it is convenient to compare our results with literature to validate the data.

2.2.Reaction density-functional theory

The reaction free energy FRcan be expressed as a summation of intrinsic reaction item and solvation item [18,19].

The reaction item ERis the intrinsic reaction free energy,which can be calculated by QDFT.For the reaction process of initial state(reactant)→transition state,For the reaction process of transition state →final state (product),The solvation free energy can be expressed as the grand potential difference between solute-on and solute-free cases.In grand canonical ensemble (μVT ensemble),it is given by[25]

Fig.3.The optimal structures of the initial,transition and final states for the proton transfer reactions in solution phase (1 water).(a) I →II,(b) III →IV,(c) A →V,(d)B →IV.

where x=(r,Ω) is the vector of position and orientation of solvent molecules.ρ(x) is the solvent density distribution surrounding the solute (M0).ρbis the bulk density.Vext(x) represents the external potential to the solvent system,which depends on the solute charge distribution and molecular configuration.The grand potential is given by [28]

where Λ denotes the thermal wavelength,and kBis the Boltzmann constant.Fexis excess free energy functional accounting for the solvent intermolecular interactions.The expression is given in Section 2.4.Under the condition of thermodynamic equilibrium,the solvent density distribution of the solute can be obtained by minimizing solvation free energy functional [25,29],

The equation is solved by Broyden-Fletcher-Goldfarb-Shanno method [30].The minimization could output the solvation free energy at equilibrium.The excess free energy Fextis a functional of ρ(x),which is evaluated in k-space by fast Fourier transformation (FFT) with FFTW3 library [31].

2.3.Evaluation of intrinsic reaction free energy

The intrinsic reaction free energies(ER)are calculated by QDFT.The second-order Moller-Plesset(MP2)perturbation theory is used with selecting a basis set of 6-311++G(2df,p) [23].The molecular configurations on the potential energy surface are optimized by using Berny’s algorithms [32].Single-point frequency calculations are then used to ensure that the optimized structures without imaginary frequency(reactants and products)or with only a single imaginary frequency (transition states).The authenticity of each transition state is further checked via the intrinsic reaction coordinate calculations [33].In addition,Tomasi’s polarized continuum model is applied at the MP2 level of theory to obtain structures of all isomers in aqueous solution [34].All calculations described above are carried by using the GAUSSIAN 09 package of programs[35].

2.4.Evaluation of excess free energy

The excess free energies of all reactions in pyruvic acid are calculated by CDFT.The molecular structures and charges of all states obtained by QDFT (see previous section)are used as an input data in CDFT.The excess free energy functional is given by [25,36,37]

with Δρ(x)is the density difference to the bulk state.ρbis the bulk density,and μbexis the excess item of bulk chemical potential.The two-body direct correlation function c(2)of bulk water is obtained by molecular dynamics simulation[36].FB[ρ(x)]is the bridge functional describing multi-interactions,given by [25]

3.Results and Discussion

3.1.Molecular configurations

Reaction is a barrier-crossing event.The state at the energy barrier (i.e.the transition state) corresponds to a saddle point in the landscape of free energy functional surface.The reaction proceeds forward only when reaction energy is greater than the reaction barrier.The value of reaction barrier quantifies the possibility and difficulty of the reaction.There are three typical states for all reactions:initial (reactant),transition and final (product) state.For pre-reaction process (initial state →transition state),the system tries to overcome the reaction barrier.For post-reaction process (transition state →final state),the system proceeds forward along the reaction and produces the new species at the final state.

It is important to optimize the molecular configurations at these three typical states.The optimal configurations of initial(reactant),transition and final(product)states for both gas and solution phase(1 water) are shown in Figs.2 and 3,respectively.Typical bond lengths and angles are listed in Table 2.In gas phase,the bond length of C6-O8 increases while that of C6-O7 decreases during the reaction.From initial to transition(TS1) state,the bond length of O7-H10 increases greatly.The angle of O7-C6-O8 decreases from 123.77° to 111.14°,while that of H10-O7-C6 decreases greatly from 106.31° to 70.70°.These behaviors of bond length and angle result in the formation of a quadrilateral;see Fig.2.O7 atom transmits a proton to O8 atom to form the quadrilateral.

Our data show that OEWM has a great effect on the behavior of bond angle during the reaction.In gas phase,the bond angle of CO-H changes by a value of only 6.35° during the reaction.In solution phase (1 water),the value is of 35.61°.In gas phase,the bond angle of O-C-O decreases about 10.20%.In solution phase(1 water),the value is of only 0.30%.The bond angle almost does not change for the reaction in solvent-on case.Moreover,the OEWM could give a stable transition state by forming a hexagonal ring,which gives a small tension value [21].Because of these two reasons,OEWM-on reactions have a much smaller barrier of reaction free energy;see Table 3.

The effect of OEWM is more important in the reaction of solution phase (1 water).A hydrogen bond can be formed between hydrogen (of M0) and oxygen (of the explicit water molecule)atoms in both initial and final states.There is a hexagonal ring in the transition state.We may regard the OEWM as a catalyst to the reaction.During this reaction,the hydrogen bond ruptures and forms successively for two times.A proton also exchanges between OEWM and M0for two times.Our results show that the assistance of water molecules plays an important role in both accepting and giving the protons.

3.2.Free energy profiles

Fig.4 and Table 3 show the data of relative free energy,which is defined as the difference of free energy between the consideredstate and the initial state (reactant):(ΔG ≡G-Ginitial).Note that the relative free energy at the transition state(ΔGTS≡GTS-Ginitial)is just the reaction barrier.The relative free energy at the final state(ΔGfinal≡Gfinal-Ginitial) quantifies the stability of product and the driving force for the reaction.

Table 2 The parameters of optimal configurations in the reaction I →II

Table 3 Relative free energy of all reactions at T=298.15 K

Fig.4.Relative free energy profiles for various reactions at T=298.15 K.(a) I →II and (b) III →IV are the isomerization reactions.(c) A →V and B →IV are the keto-enol tautomerism reactions.(1 kcal=4.1840 kJ).

The data listed in Table 3 shows that the prediction of RxDFT agrees well with the data of QDFT [23].The relative free energy for all reactions in gas phase are 1.83,0.53,12.81 and 10.81 kcal·mol-1.(List successively forI →II,III →IV,A →V and B →IV.The same in following.)Because of the eletronegativity difference between oxygen and carbon atom,oxygen atom has a stronger ability to adsorb electrons.The carbon–oxygen double bond in enol is more stable.So the PA of keto is more stable than that of enol in gas phase.The reaction barriers in gas phase are 32.40,30.91,69.03 and 66.91 kcal·mol-1,In comparison,the values in gas phase(1 water)are 11.69,11.14,39.35 and 40.49 kcal·mol-1.The OEWM significantly reduces the reaction barriers for all reactions by more than 20 kcal·mol-1.The OEWM also reduces the relative relative free energy of final state (ΔGfinal=Gfinal-Ginitial).The free energy barriers for KET are quite high while the ΔGfinalis positive.This means that the reaction is endothermic.So it cannot proceed forward if there is no solvent.

Results in solvent-on cases show that solvent is important to all reactions.The reaction barriers for all reactions in solution phase are 40.85,33.88,66.48 and 69.58 kcal·mol-1.These values deviate slightly to that in gas phase.The solvent environment surrounding the solute has small effect on proton transfer.In comparison,the reaction barriers in solution phase (1 water) are 6.94,8.74,29.37 and 32.53 kcal·mol-1.The OEWM reduces the reaction barrier by about 10 kcal·mol-1to that in gas phase.The assisting effect of water molecules plays a dominant role in the reaction.

We may decompose relative free energy (ΔG) of reaction in solution phase (1 water) into two items:(1) The solvation energy required to resist the mean potential of solvent environment.(2)The contribution of water assisting effect on the reaction.This part describes the configurational entropy of the solute molecule and the interaction effect to solvent molecules.Fig.5 compares these two contributions for all reactions.The data at the transition states(TS1 to TS4)show that the second items is more important for the reaction barrier.The data at the final states(II,IV,V,IV)show that the first and second items have comparable contributions to the stability of products.The assisting effect of water molecules plays a major role in both the pre-reaction process (initial state →transition state) and in post-reaction process (transition state →final state).In terms of product stability and driving force of reaction,the effects from both the assisting water molecules and solvent environment are important.

Fig.5.Comparison of solvation effect and OEWM’s contribution to relative free energy (ΔG=G-Ginitial ).The temperature is fixed at 298.15 K.Note that the direction of y-axis is reversed for the guide of view.

It is interesting to observe how water molecules distribute near the solute.Fig.6 shows the contour of radial distribution function for solvent in solution phase (1 water).The data of (I →TS1 and A →TS3) are shown here.The data of other reactions are similar and not shown.The water molecules are packed closely near the solute and the OEWM.Solvent distribution shows a local ordering behavior near the solute,which could reduce the contribution of solvation effect to the reaction barrier.The data in Fig.6 shows the clear evidence for the layering structure.

4.Conclusions

The keto-enol tautomerism and isomerization of pyruvic acid are studied by an RxDFT.The reaction of single solute molecule is explored for the effects of water assistance and solvent distribution.The intrinsic reaction free energy is evaluated by QDFT,and the contribution from solvent distribution is evaluated by CDFT.The RxDFT gives a good prediction for molecular conformations,reaction free energy profiles and barriers.Results show that both solvation effect and water assistance could reduce the reaction barriers.The water molecule participates the reaction as a catalyst to accept/give the protons with forming a hexagonal ring in the transition state.Because of the existence of the hexagonal ring,the solute configuration undergoes a small variation during the reaction.This leads to a diminished contribution to the intrinsic reaction free energy.As another result of water assistance,the solvent distribution shows a local ordering behavior near the solute that reduces the contribution of solvation effect to the reaction barrier.Water assistance plays a major role in the whole process of reaction.In terms of the driving force of reaction and product’s stability,the effects of both solvation and water assistance are important.

Fig.6.Contour of radial distribution function for solvent in solution phase(1 water).The data are a cross section of raw data on the plane of XOZ.The left panel shows the data of reactants,and the right panel shows the data of transition states.(1 ?=10-10 m).

With increasing the number of water molecules,the activation free energy barriers decrease first and then increase further [40–42].There is an optimal number for water molecules,at which the chemical reaction proceeds forward with the fastest speed.However,the role of hydrogen bond network in the chemical reaction still requires more studies.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Nos.21978079,and 21878078).