- Email: [email protected]
FEP method
Chemical Physics Letters 547 (2012) 97–102
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Theoretical study of the hydrolysis of ethyl benzoate in acidic aqueous solution using the QM/MC/FEP method Thanayuth Kaweetirawatt a,b, Yohei Kokita b, Shiho Iwai b, Michinori Sumimoto b, Kenji Hori b,⇑ a b
Ube Technical Center (Asia) Co., Ltd., 140/8 Moo4, Tambol Tapong, Muang Rayong, Rayong 21000, Thailand Graduate School of Science and Engineering, Yamaguchi University, 2-16-1 Tokiwadai, Ube, Yamaguchi 755 8611, Japan
a r t i c l e
i n f o
Article history: Received 25 April 2012 In final form 1 June 2012 Available online 15 June 2012
a b s t r a c t The hydrolysis of ethyl benzoate in acidic condition was theoretically studied for models with two (2W) or three water (3WA) molecules at the B3LYP/6-311++G(d,p) levels of theory. Activation free energy of solvation in aqueous solution (DGàcal) was calculated using the QM/MC/FEP method. The value of the 2W model in aqueous solution was calculated to be smaller by more than 5.0 kcal mol1 than the observed value (26.0 kcal mol1 at 298 K). The position of the third water molecule in the 3WA model plays an essential role in producing the DGàcal value (26.4 kcal mol1) consistent with the experimental value. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction The persistence of an organic compound in the aquatic environment is determined by its water solubility, vapor pressure and chemical and biological reactivity. The hydrolysis reactions of many organic pollutants such as carboxylic acid esters, amides, and phosphate esters, have been most widely studied. For instance, the hydrolysis of substituted methyl benzoates in aqueous solution was studied by Steinberg and Lena [1]. Furusjo and co-workers studied the kinetics and mechanism of acid and alkaline hydrolysis of phenyl acetate [2] and the hydrolysis of ethyl acetate in aqueous hydrochloric acid [3]. The reaction mechanism of alkaline hydrolysis of phenyl benzoate in ethanol–water media was investigated by Zhao et al. [4]. In addition, Zhu et al. [5] and Patnaik et al. [6] also performed research on the kinetics and mechanism of the hydrolysis of phthalate esters. The hydrolysis of ethyl benzoate has been a popular subject of investigation. The rate constants of the alkaline and acid-catalyzed hydrolyses are presented in Hilal’s study [7]. Usually, the hydrolysis proceeds in mixed solvents such as acetone/water or methanol/ water [8]. The use of pure water as a solvent is very rare. Generally kinetic parameters such as rates of reactions depend on their activation free energies. Density Functional Theory (DFT) calculations can compute such values. If solvent effects are included, theoretical methods are one of the feasible ways to avoid experimental difficulties [9].
⇑ Corresponding author. E-mail address: [email protected] (K. Hori). 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.06.001
Recently, interest in the hybrid explicit/implicit method for simulating solute–solvent systems has grown [10–14]. According to such method, part of the solvent effects, typically including a few solvent molecules, is treated explicitly along with the solute, while several methods are used to take the implicit effects of solvents into account. Typical methods for calculating such effects are the Self-Consistent Reaction Field (SCRF) method or a Molecular Dynamic (MD) simulation [15,16]. A Monte Carlo (MC) simulation is an alternative that can calculate properties concerned with solvents by generating random configurations of a particular ensemble [17–21]. Almost all MC simulations use classical force fields. In order to avoid the difficulties in obtaining the classical force-field parameters, quantum mechanical (QM) calculations have been adopted to compute the ensemble energies [22,23]. This method, called a ‘QM/MC simulation’, has been further developed to calculate the free-energy surfaces of chemical reactions in solution, i.e., the Free-Energy Perturbation (FEP) theory has been introduced to the QM/MC simulation, the ‘QM/MC/FEP method’, to calculate the difference in the free energy of solvation (DDGsol) between two perturbed structures [24,25]. In this contribution, we investigated the hydrolysis of ethyl benzoate to produce benzoic acid and ethanol in acidic aqueous solution. Moreover, our attention was focused on the activation free energy of the reaction. We attempted to find a suitable model with an appropriate level of theory to produce values consistent with
98
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
the experiment value. The implicit effects of solvent are calculated by both the QM/MC/FEP and the SCRF methods. 2. Experimental details 2.1. Materials Ethyl benzoate of analytical grade was purchased from Wako. For pH control, 1 M HCl was obtained from Sigma–Aldrich. Other reagents were analytical grade and without further purification. Laboratory grade water with 18 MX cm and was prepared using EASY pure RoDi water system. 2.2. Procedure Ethyl benzoate in aqueous solution was prepared by sonication. The pH of the solution was controlled to be 1.3 by adding HCl and the concentration of ethyl benzoate was adjusted to be 4.5 mM. The solution was observed to be homogeneous after putting in the sonicator for 5 min. The solution was divided into several screw cap vials. The vials were placed in a temperature-controlled water bath equipped with a magnetic stirrer. The temperature (T) range was varied between 303 and 333 K. Samples (1 ll) were withdrawn from a vial at appropriate reaction times. Even though the hydrolysis proceeded slowly in the low temperature such as 303 K, the measurement was completed within 3 days. The concentration of ethyl benzoate in the solution was determined using LC-MS (Shimadzu LCMS-2010EV). The thermodynamic parameters such as the activation enthalpy (DHà = 20.9 kcal mol1) and entropy (DSà = 0.017 kcal mol1 K1) were determined from the temperature dependencies of the rate constants (k) using an Eyring plot. R2 of the plot was 0.996. Then the DGàexp at a certain temperature can be computed according to the thermodynamic relation (Eq. (1)). As a result, a value of 26.0 kcal mol1 was obtained for DGàexp at 298 K.
DGzexp ¼ DHz T DSz
ð1Þ
3. Computational details 3.1. Gas-phase calculations All the geometry optimizations and activation free energies (DGàgas) in the gas phase at 298 K were calculated using the
B3LYP Functional [26] together with 6-311G(d,p), and 6311++G(d,p) basis sets in the GAUSSIAN 03 program [27]. We consider three interesting models, a two waters (2W) model [9,28,29], two types of three waters (3W) model, 3WA and 3WB, in this study. Their transition state (TS) structures are shown in Figure 1. A proton was also added to represent the specific acid catalysis. There are 200 intrinsic reaction coordinate (IRC) structures between s = 5.0 and 5.0 amu1/2 Bohr shown in Figure 2 [30,31]. 3.2. QM/MC/FEP calculations The solvent effects may cause change of the geometry in gas phase. However, our QM/MC/FEP program at present cannot optimize the structures in solution. Therefore, the DGàsol values from the QM/MC/FEP calculations were added to the DGàgas values to represent the activation free energies of reaction in aqueous solution. Other investigations adopted similar methods [32–35]. In order to calculate the activation free energies of solvation (DGàsol), the IRC structures (100 structures) were adopted as perturbed structures in the simulation process. An additional 50 interpolated structures that bridged the last IRC structure to the reactant were also included. The number of IRC and interpolated structures was decided according to the definition of the free energy perturbation theory (FEP). The calculations of FEP can converge properly when the difference between the two states is small enough. Therefore, there are 150 perturbed structures in the QM/MC/FEP calculations. MC simulations with the NPT ensemble used spherical droplets including a solute molecule in its center and many solvent molecules randomly placed within the droplet [22,23]. The 7.5 Å radius (r) of the droplet included 35 implicit water molecules for 2W model and 34 ones for 3W models. We have previously confirmed that this size of droplets is sufficient to obtain DGàsol values that are close to experimental ones [24]. The number of solvent molecules is dependent on the solvent density, which is 0.998 (g cm3) in this case. Additional parameters such as the maximum allowed displacement (a = 0.06 Å/step) and the maximum allowed angle of rotation (b = 0.12 rad/step) were also taken into account [24,25]. The 40,000 steps of the MC simulations in the Metropolis Sampling Algorithm were divided into two parts; the first half (20,000 steps) was used for equilibrating the droplet and the second half was used for calculating the energies of the ensemble. The energies in the latter part were calculated using the PM3 Hamiltonians of the MOPAC2000 program and were averaged to represent the energy of each droplet. This Hamiltonian gave reasonable results in previous studies [23–25].
Figure 1. TS structures of 2W, 3WA and 3WB models. The unit of length is Å.
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
99
Figure 2. The potential energy profile and geometry transformation of 2W model along the reaction coordinates calculated by B3LYP/6-311++G(d,p) level of theory. The unit of length is Å.
DGàsol values were obtained by integrating the implicit effects of solvents (DDGsol) from the reactant to the TS structure. The activation free energy in aqueous solution (DGzcal ) can be calculated by adding DGzsol to DGzgas [25].
DGzsol ¼
Z
TS
dsDDGsol
ð2Þ
reactant
DGzcal ¼ DGzgas þ DGzsol
ð3Þ
In the SCRF calculations, the conductor-like polarizable continuum model (CPCM) [36–38] was used to calculate the implicit effects of the solvent at 298 K. The optimizations at the B3LYP/6311G(d,p) and B3LYP/6-311++G(d,p) levels of theory were carried out using a dielectric constant of 78.39 to simulate an aqueous environment. The cavity of the United Atom Topological Model (UAHF) [39] was used in the calculations. Finally, the free energies were calculated by adding the zero-point vibrational energies. The TS structures for the 3W models were not optimized in the SCRF calculations. 4. Results and discussion 4.1. Differences between ethyl benzoate and methyl acetate First of all, we compared the hydrolysis of ethyl benzoate with that of methyl acetate under acidic conditions. The hydrolysis proceeds by the AAC2 mechanism [9,40]. The addition of a nucleophile to a carbonyl carbon accompanies a subsequent removal of alcohol. The carbonyl oxygen accepts a proton and this results an increase in the positive charge of the carbonyl carbon. For methyl acetate, the protonation facilitated the addition of water at that position to form a tetrahedral (TD) intermediate, as shown in Figure 3, i.e., the hydrolysis proceeds in two steps: the formation of the TD intermediate and its decomposition to produce a protonated acetic
Figure 3. Hydrolysis of ethyl benzoate and methyl acetate with two explicit-water molecules.
acid and a methanol [9]. Wolfe and co-workers confirmed that an acetic acid molecule can play a role in acting as both the proton donor and the acceptor. In this case, the hydrolysis of methyl acetate proceeds in only one step of reaction [41,42]. Figure 2 depicts the potential surface and the geometry transformation in the hydrolysis of ethyl benzoate calculated at the B3LYP/ 6-311++G(d,p) level of theory for the 2W model (Run 2). According to the results of Run 2, the C(1)–O(6) length increased along the reaction coordinates from 1.291 Å (s = 5.0 amu1/2 Bohr) 1/2 to 1.649 Å (s = 5.0 amu Bohr). Moreover, the length of the H(5)– O(6) bond decreased from 2.456 to 1.013 Å. These changes in the two lengths indicate that the cleavage of the C(1)–O(6) bond and the formation of the H(5)–O(6) bond proceed simultaneously. An optimization from the structure at s = 5.0 amu1/2 Bohr produced the protonated ethanol and benzoic acid. The results indicate that the TD intermediate does not form in the hydrolysis of ethyl benzoate. Therefore, the hydrolysis of the ethyl benzoate proceeds via
100
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
the one step mechanism. The other results of the 3W models also showed a similar transformation using the IRC calculations. 4.2. Two waters model There are six significant bonds between ethyl benzoate and two explicit-water molecules in the TS structures of the 2W model, TS(2W), as shown in Figure 1. In the gas phase, the inclusion of diffuse functions in the level of theory changes the geometries of the TS structure only slightly, so that they are all very similar to each other. For example, the geometrical parameters calculated at the B3LYP/6-311G(d,p) level of theory (Run 1) are different from those calculated in Run 2 within 0.014 Å. The largest difference is seen in the H(5)–O(6) length in Table 1. The roles of diffuse functions are the same in producing the TS structures calculated by the SCRF method. The geometry calculated at the CPCM-B3LYP/6-311G(d,p) level of theory (Run 4) are very similar to that of CPCM-B3LYP/6-311++G(d,p) level of theory (Run 5), different within 0.024 Å. For example the O(2)–H(3) length of Run 4 was calculated to be 1.747 Å while that of Run 5 turned out to be 1.771 Å. Therefore, it can be concluded that the TS structures are independent of the inclusion of diffuse functions. However, the TS structures produced by the SCRF method are obviously different from those of gas-phase calculations, especially the C(1)–O(2) and H(5)–O(6) lengths. For example, the lengths in Run 2 are 1.478 and 1.651 Å while those in Run 5 are 1.431 and 1.820 Å. The former differs from the later by 0.047 and 0.169 Å, respectively. Therefore, the implicit effects of solvent using the SCRF method do influence the TS structures of this reaction. As will be discussed later, the inclusion of a third explicit-water molecule changes this length markedly, i.e., the effects of the explicit solvent are more important than those of the implicit solvent included in the SCRF calculations. For the reactants in the 2W model in Figure 4, R(2W), the geometries depended on the diffuse functions very much. The bond lengths are also listed in Table 1. The diffuse functions cause the distances within each molecule to lengthen. For example, the C(1)–O(2) and H(3)–O(4) lengths of Run 2 are 2.786 and 1.846 Å, longer than those of Run 1 by 0.140 and 0.050 Å, respectively.
Similar elongations are also seen in the reactants in the SCRF method. The C(1)–O(2) length of Run 5 was determined to be 2.792 Å, longer than that of Run 4 by 0.164 Å. Moreover, the H(3)–O(4) length of the former (1.839 Å) is longer than that of the latter (1.805) by 0.034 Å. As the C(1)–O(2) and H(3)–O(4) lengths of Run 5 are very close to those of Run 2, meaning that the implicit effects of solvent using the SCRF method have no influence on the reactant geometry. The diffuse functions affected not only the geometry but also the energy. The activation free energy in the gas phase (DGàgas) of Run 1 without these functions is computed to be 26.0 kcal mol1, lower than that of Run 2 with the functions by 6.0 kcal mol1, as shown in Figure 5. Delchev and co-workers reported that the diffuse functions have strong effects on the energy calculations of the B3LYP method but are negligible for the BLYP method in the case of the complex ion [Cu(H2O)]+ [43]. The present results do seem to be consistent with those previously reported. In order to elucidate how large are the effects that the diffuse functions cause, the DGàgas value was calculated at the B3LYP/6311++G(d,p)//B3LYP/6-311G(d,p) level of theory (Run 3). That value of Run 3 was computed to be 29.6 kcal mol1, which is larger by 3.6 kcal mol1 than that of Run 1. The activation free energy of solvation (DGàsol) was then considered to obtain the activation free energy in aqueous solution, DGàcal. The QM/MC/FEP calculations produced large negative DGàsol values, depending on the perturbed structures, at each level of theory. For example, those of Run 1 and Run 2 were calculated to be 13.1 and 14.0 kcal mol1. Their difference is only 0.9 kcal mol1. Also, the DGàcal values of Run 1 and Run 2 in Figure 5 were calculated to be 12.9 and 18.0 kcal mol1, respectively. The value of Run 1 is lower by 5.1 kcal mol1 in comparison with that of Run 2. The DGàcal value of Run 3 was obtained as 16.5 kcal mol1. It is higher than that of Run 1 by 3.6 kcal mol1, but lower than that of Run 2 by only 1.5 kcal mol1. The DGàcal values of the 2W model calculated by the QM/MC/ FEP method are very consistent with those of the SCRF method. For example, the DGàcal values of Run 4 and Run 5 in Figure 5 are 9.2 and 17.4 kcal mol1, respectively. The DGàcal value of Run 5 differs from that of Run 2 by only 0.6 kcal mol1. As the observed
Table 1 Optimized geometries of 2W, 3WA and 3WB models. The unit of length is Å. Run
Model
Level of theory
1 2 4 5 6 7
2W 2W 2W 2W 3WA 3WB
B3LYP/6-311G(d,p) B3LYP/6-311++G(d,p) CPCM-B3LYP/6-311G(d,p) CPCM-B3LYP/6-311++G(d,p) B3LYP/6-311++G(d,p) B3LYP/6-311++G(d,p)
Reactant
TS structure
C(1)–O(2)
H(3)–O(4)
H(8)–O(10)
C(1)–O(2)
C(1)–O(6)
O(2)–H(3)
H(3)–O(4)
O(4)–H(5)
H(5)–O(6)
2.646 2.786 2.628 2.792 2.865 2.854
1.796 1.846 1.805 1.839 1.861 1.857
– – – – 1.632 1.653
1.478 1.478 1.432 1.431 1.432 1.491
1.419 1.420 1.421 1.419 1.451 1.433
1.458 1.466 1.747 1.771 1.488 1.459
1.059 1.057 1.005 1.004 1.053 1.06
1.016 1.020 0.999 0.999 1.046 1.031
1.666 1.651 1.804 1.820 1.529 1.586
Figure 4. Reactant structures of 2W, 3WA and 3WB models.
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
101
Figure 5. Activation free energies of hydrolysis of ethyl benzoate in gas phase and in aqueous at 298 K.
value (DGàexp) was 26.0 kcal mol1 at 298 K, the DGàcal values are not consistent with the experimental one. 4.3. Three waters models Due to the very large activation free energy of solvation in the 2W model, the reaction mechanisms were investigated more deeply. In TS(2W), the C(1) atom possesses two hydroxyl groups. One of the groups consists of the O(2) and H(9) atoms. These two atoms play a very important role in donating the proton H(3) to allow the reaction to proceed. It is easy to understand that the positions of the O(2) and H(9) atoms are changed significantly along the reaction coordinate by comparing the structures of R(2W) and TS(2W). The position of the other hydroxyl group consisting of O(7) and H(8) atoms in the TS(2W) is also very different from that in R(2W). After analyzing the population of Run 2 by the Merz-Singh– Kollman (MK) option, the charges of the O(2) and H(9) atoms were computed to be 0.69 and 0.50 in TS(2W) and they were found to be 0.77 and 0.37 in R(2W), respectively. The difference of the charge in the O(2) atom is 0.08. However, it is smaller than that of proton H(9), which is 0.13. Not only the positions, but also the charges of those atoms are dramatically changed, especially the proton H(9). In comparison, the charge of the proton H(8) was obtained to be 0.49 in TS(2W), but was calculated to be 0.45 in R(2W), different by only 0.04. Taking account of the positions and the charge differences of the hydroxyl groups along the reaction coordinate, it is possible that the implicit effects of solvent might not be included properly in the calculations of the 2W model. This resulted in overestimating the DGàsol value by as much as 14.0 kcal mol1 in Run 2. Therefore, one more explicit-water molecule has to be involved to interact with the hydroxyl groups in the TS structure in order to calculate the solvent effects more appropriately. In order to show that the position of a third explicit-water molecule plays a crucial role in calculating the solvent effects, two mechanisms including the third explicit-water molecule, 3WA (Run 6) and 3WB (Run 7), were proposed. An additional water molecule interacts with the proton H(9) in TS(2W) in the former mechanism whereas the water is added to interact with the proton H(8) rather than the proton H(9) in the latter. For both of the two types of three waters model, the calculations adopted only the B3LYP/6311++G(d,p) level of theory in the investigation. In the TS structure of the 3WA model, TS(3WA), the O(10) atom of the third explicit-water molecule interacts with the proton H(9), as depicted in Figure 1. The H(9)–O(10) length was calculated to be 1.664 Å. For TS(3WA), the geometrical parameters are significantly different from those of TS(2W) as shown in Table 1. For example,
the C(1)–O(2) length of Run 6 was computed to be 1.432 Å, shorter by 0.046 than that of Run 2. Moreover, the H(5)–O(6) length of the former turned out to be 1.529 Å, which is also shorter by 0.122 Å than that of the latter. On the other hand, the O(10) atom forms a hydrogen bond with the proton H(8), as shown by the H(8)–O(10) length (1.690 Å) in the TS structure of the 3WB model, TS(3WB). While the geometrical parameters of TS(3WA) are much different from those of TS(2W), the parameters of TS(3WB) are not much different to those of the latter. For example, the H(5)–O(6) length of Run 7 is computed to be 1.586 Å, shorter by 0.065 Å than that of Run 2. It is the largest difference between those two TS structures listed in Table 1. Moreover, there are many bond lengths in TS(3WB) which are longer than those in the TS(3WA). For example, the C(1)–O(2) and H(5)– O(6) lengths of the former were found to be 1.491 and 1.586 Å, while those of the latter turned out to be 1.432 and 1.529 Å. Even though the geometrical parameters of TS(3WA) and TS(3WB) are obviously different, the reactants of the three waters models, R(3WA) and R(3WB), are very similar within 0.021 Å, as shown in Figure 4. The biggest difference can be seen in the H(8)–O(10) length in Table 1. Furthermore, compared to the length of R(2W) in Run 2, the C(1)–O(2) lengths of R(3WA) (2.865 Å) and R(3WB) (2.854 Å) are longer by 0.079 and 0.068 Å, respectively. Consequently, the influence of the third water is to lengthen the distance between the ethyl benzoate and the first two waters. As plotted in Figure 5, the DGàgas values of Run 6 and Run 7 are calculated to be 33.6 and 35.0 kcal mol1. TS(3WA) is more stable than TS(3WB) by 0.8 kcal mol1. Figure 6 displays the DGàsol values of 3W models calculated using the QM/MC/FEP method. The DGàsol values of Runs 6 and 7 were computed to be 7.2 and 14.9 kcal mol1, respectively. Those DGàsol values are different by 7.7 kcal mol1. While the DGàsol value of 3WA model is far from that of the 2W model (Run 2), the value of the 3WB model is very close to that of the 2W model by less than 1.0 kcal mol1. Therefore, the position of the third water plays an essential role in producing the DGàsol values based on the QM/MC/FEP method. Corresponding to the DGàsol values, the DGàcal values of Runs 6 and 7 were 26.4 and 20.1 kcal mol1, respectively. We could not obtain the corresponding values of the models using the SCRF method since the optimizations for the TS structures were not successful. This seems to be a disadvantage of the SCRF method. 4.4. Comparison of DGàcal and DGàexp The DGàcal values of the 2W model were different from experiment (DGàexp) by 8.0 kcal mol1 at the B3LYP/6-311++G(d,p) level
102
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
Figure 6. Free energies of solvation (DGàsol) of 3WA and 3WB models at 298 K.
of theory. It means that two explicit-water molecules are adequate to make the reaction proceed via their TS structures but they are still not enough to represent the solvent effects of the reaction with regard to the DGàcal value. The position of the third explicit-water molecule is very important, as explained earlier. If the third water molecule was placed in the wrong position, the DGàcal value in the 3WB model was still not consistent with the experimental value. It differs from the observed free energy by more than 5.0 kcal mol1. Both the 2W and 3WB models seem to have the same problem, which is a large negative value of DGàsol. This propensity apparently comes from large changes in the charge and positions of the O(2)–H(9) group. On the other hand, when the third water molecule has been placed in the right position in the 3WA model, it assists the process of reaction and produces a value of DGàsol that is more reasonable. The outcome is that the DGàcal value for the 3WA model, 26.4 kcal mol1, is very consistent with the experimental free energy within less than 1.0 kcal mol1. 5. Conclusion Including diffuse functions in QM/MC/FEP calculations of the hydrolysis of ethyl benzoate directly affects the reactant geometries, especially the distance between the ethyl benzoate and two explicit-water molecules. The functions lengthen the distances between ethyl benzoate and the water molecules. However, these functions have only a small influence on the TS structures. The DGàgas value calculated by including the diffuse functions is higher than that when that are excluded. Two explicit-water molecules are not enough to correctly represent the mechanism of hydrolysis of ethyl benzoate since the DGàcal values calculated by both the QM/MC/FEP and the CPCM-SCRF method are not consistent with experiment. Therefore, a third explicit-water molecule is necessary. The third water molecule increases the distance between the ethyl benzoate and the first two water molecules in the reactant geometries. Moreover, the position of the third water in the TS structure plays a vital role in producing a DGàcal value from the QM/MC/FEP calculations that is highly consistent with the experimental value. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2012. 06.001.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
S.M. Steinberg, F. Lena, Water Res. 29 (1995) 965. E. Furusjo, L.G. Danielsson, Chemom. Intell. Lab. Syst. 50 (2000) 63. E. Furusjo, L.G. Danielsson, Anal. Chim. Acta 373 (1998) 83. Y. Zhao, G. Wang, W. Li, Z.L. Zhu, Chemom. Intell. Lab. Syst. 82 (2006) 193. Z.L. Zhu, J. Xia, J. Zhang, T.H. Li, Anal. Chim. Acta 454 (2002) 21. P. Patnaik, M. Yang, E. Powers, Water Res. 35 (2001) 1587. S.H. Hilal, United States Environmental Protection Agency, 2006. (EPA/600/R06/105). E.W. Timm, C.N. Hinshelwood, J. Chem. Soc. 3 (1938) 862. K. Hori et al., Tetrahedron 63 (2007) 1264. A. Okur, C. Simmerling, Annual Reports in Computational Chemistry, vol. 2, Elsevier, Amsterdam, 2006. B. Roux, T. Simonson, Biophys. Chem. 79 (1999) 1. M. Orozco, F.J. Luque, Chem. Rev. 100 (2000) 4187. G. Brancato, N. Rega, V. Barone, J. Chem. Phys. 128 (2008) 144501. E.F. Lima, J.W.M. Carneiro, C. Fenollar-Ferrer, S. Miertus, S. Zinoviev, N.C. Om Tapanes, D.A.G. Aranda, Fuel 89 (2010) 685. H. Freedman, H.N. Nguyen, T.N. Truong, J. Phys. Chem. B 108 (2004) 19043. H. Freedman, T.N. Truong, J. Phys. Chem. B 109 (2005) 4726. J.D. Madura, W.L. Jorgensen, J. Am. Chem. Soc. 108 (1986) 2517. J. Chandrasekhar, S.F. Smith, W.L. Jorgensen, J. Am. Chem. Soc. 106 (1984) 3049. J. Chandrasekhar, S.F. Smith, W.L. Jorgensen, J. Am. Chem. Soc. 107 (1985) 154. W.L. Jorgensen, J.F. Blake, D.C. Lim, D.L. Severance, J. Chem. Soc., Faraday Trans. 20 (1994) 1727. Y. Wu, Y. Xue, D.Q. Xie, C.K. Kim, F.S. Yan, J. Phys. Chem. B 111 (2007) 2357. T. Yamaguchi, M. Sumimoto, K. Hori, J. Chem. Phys. Lett. 460 (2008) 331. T. Yamaguchi, M. Sumimoto, K. Hori, J. Comput.-Aided Chem. 9 (2008) 62. K. Hori, T. Yamaguchi, K. Uezu, M. Sumimoto, J. Comput. Chem. 32 (2011) 778. T. Kaweetirawatt, T. Yamaguchi, H. Tsutomu, M. Sumimoto, K. Hori, J. Phys. Org. Chem. (2012), http://dx.doi.org/10.1002/poc.2944. D. Boese, J.M.L. Martin, J. Chem. Phys. 119 (2003) 3005. M.J. Frisch et al., GAUSSIAN 03 (Revision C.02), Gaussian Inc., Wallingford CT, 2004. L.C. Mether, D.V. Sagar, S.N. Naik, Renew. Sustain. Energy Rev. 10 (2006) 248. F. Ma, L.D. Clements, M.A. Hanna, Trans. ASAE 41 (1998) 1261. K. Fukui, Acc. Chem. Res. 14 (1981) 363. C. Gonzalez, H.B. Schlegel, J. Chem. Phys. 90 (1989) 2154. H. Zipse, G. Apaydin, K.N. Houk, J. Am. Chem. Soc. 117 (1995) 8608. G. Kaminski, W.L. Jorgensen, J. Phys. Chem. 102 (1998) 1787. O. Accevedo, W.L. Jorgensen, J. Am. Chem. Soc. 127 (2005) 8829. O. Accevedo, W.L. Jorgensen, J. Am. Chem. Soc. 128 (2006) 6141. O. Tapia, Theoret. Chim. Acta (Berl.) 47 (1978) 157. M. Chen, Z. Lin, J. Chem. Phys. 127 (2007) 154314. C.K. Kim, B. Park, H.W. Lee, C.K. Kim, Bull. Korean Chem. Soc. 32 (2011) 1985. Y. Takano, K.N. Houk, J. Chem. Theory Comput. 1 (2005) 70. M.B. Smith, J. March, March’s Advanced Organic Chemistry, fifth ed., 2001. Y. Hsieh, N. Weinberg, S. Wolfe, Can. J. Chem. 87 (2009) 539. Z. Shi, Y. Hsieh, N. Weinberg, S. Wolfe, Can. J. Chem. 87 (2009) 544. V.B. Delchev, H. Mikosch, J. Struct. Chem. 47 (2006) 979.
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Theoretical study of the hydrolysis of ethyl benzoate in acidic aqueous solution using the QM/MC/FEP method Thanayuth Kaweetirawatt a,b, Yohei Kokita b, Shiho Iwai b, Michinori Sumimoto b, Kenji Hori b,⇑ a b
Ube Technical Center (Asia) Co., Ltd., 140/8 Moo4, Tambol Tapong, Muang Rayong, Rayong 21000, Thailand Graduate School of Science and Engineering, Yamaguchi University, 2-16-1 Tokiwadai, Ube, Yamaguchi 755 8611, Japan
a r t i c l e
i n f o
Article history: Received 25 April 2012 In final form 1 June 2012 Available online 15 June 2012
a b s t r a c t The hydrolysis of ethyl benzoate in acidic condition was theoretically studied for models with two (2W) or three water (3WA) molecules at the B3LYP/6-311++G(d,p) levels of theory. Activation free energy of solvation in aqueous solution (DGàcal) was calculated using the QM/MC/FEP method. The value of the 2W model in aqueous solution was calculated to be smaller by more than 5.0 kcal mol1 than the observed value (26.0 kcal mol1 at 298 K). The position of the third water molecule in the 3WA model plays an essential role in producing the DGàcal value (26.4 kcal mol1) consistent with the experimental value. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction The persistence of an organic compound in the aquatic environment is determined by its water solubility, vapor pressure and chemical and biological reactivity. The hydrolysis reactions of many organic pollutants such as carboxylic acid esters, amides, and phosphate esters, have been most widely studied. For instance, the hydrolysis of substituted methyl benzoates in aqueous solution was studied by Steinberg and Lena [1]. Furusjo and co-workers studied the kinetics and mechanism of acid and alkaline hydrolysis of phenyl acetate [2] and the hydrolysis of ethyl acetate in aqueous hydrochloric acid [3]. The reaction mechanism of alkaline hydrolysis of phenyl benzoate in ethanol–water media was investigated by Zhao et al. [4]. In addition, Zhu et al. [5] and Patnaik et al. [6] also performed research on the kinetics and mechanism of the hydrolysis of phthalate esters. The hydrolysis of ethyl benzoate has been a popular subject of investigation. The rate constants of the alkaline and acid-catalyzed hydrolyses are presented in Hilal’s study [7]. Usually, the hydrolysis proceeds in mixed solvents such as acetone/water or methanol/ water [8]. The use of pure water as a solvent is very rare. Generally kinetic parameters such as rates of reactions depend on their activation free energies. Density Functional Theory (DFT) calculations can compute such values. If solvent effects are included, theoretical methods are one of the feasible ways to avoid experimental difficulties [9].
⇑ Corresponding author. E-mail address: [email protected] (K. Hori). 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.06.001
Recently, interest in the hybrid explicit/implicit method for simulating solute–solvent systems has grown [10–14]. According to such method, part of the solvent effects, typically including a few solvent molecules, is treated explicitly along with the solute, while several methods are used to take the implicit effects of solvents into account. Typical methods for calculating such effects are the Self-Consistent Reaction Field (SCRF) method or a Molecular Dynamic (MD) simulation [15,16]. A Monte Carlo (MC) simulation is an alternative that can calculate properties concerned with solvents by generating random configurations of a particular ensemble [17–21]. Almost all MC simulations use classical force fields. In order to avoid the difficulties in obtaining the classical force-field parameters, quantum mechanical (QM) calculations have been adopted to compute the ensemble energies [22,23]. This method, called a ‘QM/MC simulation’, has been further developed to calculate the free-energy surfaces of chemical reactions in solution, i.e., the Free-Energy Perturbation (FEP) theory has been introduced to the QM/MC simulation, the ‘QM/MC/FEP method’, to calculate the difference in the free energy of solvation (DDGsol) between two perturbed structures [24,25]. In this contribution, we investigated the hydrolysis of ethyl benzoate to produce benzoic acid and ethanol in acidic aqueous solution. Moreover, our attention was focused on the activation free energy of the reaction. We attempted to find a suitable model with an appropriate level of theory to produce values consistent with
98
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
the experiment value. The implicit effects of solvent are calculated by both the QM/MC/FEP and the SCRF methods. 2. Experimental details 2.1. Materials Ethyl benzoate of analytical grade was purchased from Wako. For pH control, 1 M HCl was obtained from Sigma–Aldrich. Other reagents were analytical grade and without further purification. Laboratory grade water with 18 MX cm and was prepared using EASY pure RoDi water system. 2.2. Procedure Ethyl benzoate in aqueous solution was prepared by sonication. The pH of the solution was controlled to be 1.3 by adding HCl and the concentration of ethyl benzoate was adjusted to be 4.5 mM. The solution was observed to be homogeneous after putting in the sonicator for 5 min. The solution was divided into several screw cap vials. The vials were placed in a temperature-controlled water bath equipped with a magnetic stirrer. The temperature (T) range was varied between 303 and 333 K. Samples (1 ll) were withdrawn from a vial at appropriate reaction times. Even though the hydrolysis proceeded slowly in the low temperature such as 303 K, the measurement was completed within 3 days. The concentration of ethyl benzoate in the solution was determined using LC-MS (Shimadzu LCMS-2010EV). The thermodynamic parameters such as the activation enthalpy (DHà = 20.9 kcal mol1) and entropy (DSà = 0.017 kcal mol1 K1) were determined from the temperature dependencies of the rate constants (k) using an Eyring plot. R2 of the plot was 0.996. Then the DGàexp at a certain temperature can be computed according to the thermodynamic relation (Eq. (1)). As a result, a value of 26.0 kcal mol1 was obtained for DGàexp at 298 K.
DGzexp ¼ DHz T DSz
ð1Þ
3. Computational details 3.1. Gas-phase calculations All the geometry optimizations and activation free energies (DGàgas) in the gas phase at 298 K were calculated using the
B3LYP Functional [26] together with 6-311G(d,p), and 6311++G(d,p) basis sets in the GAUSSIAN 03 program [27]. We consider three interesting models, a two waters (2W) model [9,28,29], two types of three waters (3W) model, 3WA and 3WB, in this study. Their transition state (TS) structures are shown in Figure 1. A proton was also added to represent the specific acid catalysis. There are 200 intrinsic reaction coordinate (IRC) structures between s = 5.0 and 5.0 amu1/2 Bohr shown in Figure 2 [30,31]. 3.2. QM/MC/FEP calculations The solvent effects may cause change of the geometry in gas phase. However, our QM/MC/FEP program at present cannot optimize the structures in solution. Therefore, the DGàsol values from the QM/MC/FEP calculations were added to the DGàgas values to represent the activation free energies of reaction in aqueous solution. Other investigations adopted similar methods [32–35]. In order to calculate the activation free energies of solvation (DGàsol), the IRC structures (100 structures) were adopted as perturbed structures in the simulation process. An additional 50 interpolated structures that bridged the last IRC structure to the reactant were also included. The number of IRC and interpolated structures was decided according to the definition of the free energy perturbation theory (FEP). The calculations of FEP can converge properly when the difference between the two states is small enough. Therefore, there are 150 perturbed structures in the QM/MC/FEP calculations. MC simulations with the NPT ensemble used spherical droplets including a solute molecule in its center and many solvent molecules randomly placed within the droplet [22,23]. The 7.5 Å radius (r) of the droplet included 35 implicit water molecules for 2W model and 34 ones for 3W models. We have previously confirmed that this size of droplets is sufficient to obtain DGàsol values that are close to experimental ones [24]. The number of solvent molecules is dependent on the solvent density, which is 0.998 (g cm3) in this case. Additional parameters such as the maximum allowed displacement (a = 0.06 Å/step) and the maximum allowed angle of rotation (b = 0.12 rad/step) were also taken into account [24,25]. The 40,000 steps of the MC simulations in the Metropolis Sampling Algorithm were divided into two parts; the first half (20,000 steps) was used for equilibrating the droplet and the second half was used for calculating the energies of the ensemble. The energies in the latter part were calculated using the PM3 Hamiltonians of the MOPAC2000 program and were averaged to represent the energy of each droplet. This Hamiltonian gave reasonable results in previous studies [23–25].
Figure 1. TS structures of 2W, 3WA and 3WB models. The unit of length is Å.
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
99
Figure 2. The potential energy profile and geometry transformation of 2W model along the reaction coordinates calculated by B3LYP/6-311++G(d,p) level of theory. The unit of length is Å.
DGàsol values were obtained by integrating the implicit effects of solvents (DDGsol) from the reactant to the TS structure. The activation free energy in aqueous solution (DGzcal ) can be calculated by adding DGzsol to DGzgas [25].
DGzsol ¼
Z
TS
dsDDGsol
ð2Þ
reactant
DGzcal ¼ DGzgas þ DGzsol
ð3Þ
In the SCRF calculations, the conductor-like polarizable continuum model (CPCM) [36–38] was used to calculate the implicit effects of the solvent at 298 K. The optimizations at the B3LYP/6311G(d,p) and B3LYP/6-311++G(d,p) levels of theory were carried out using a dielectric constant of 78.39 to simulate an aqueous environment. The cavity of the United Atom Topological Model (UAHF) [39] was used in the calculations. Finally, the free energies were calculated by adding the zero-point vibrational energies. The TS structures for the 3W models were not optimized in the SCRF calculations. 4. Results and discussion 4.1. Differences between ethyl benzoate and methyl acetate First of all, we compared the hydrolysis of ethyl benzoate with that of methyl acetate under acidic conditions. The hydrolysis proceeds by the AAC2 mechanism [9,40]. The addition of a nucleophile to a carbonyl carbon accompanies a subsequent removal of alcohol. The carbonyl oxygen accepts a proton and this results an increase in the positive charge of the carbonyl carbon. For methyl acetate, the protonation facilitated the addition of water at that position to form a tetrahedral (TD) intermediate, as shown in Figure 3, i.e., the hydrolysis proceeds in two steps: the formation of the TD intermediate and its decomposition to produce a protonated acetic
Figure 3. Hydrolysis of ethyl benzoate and methyl acetate with two explicit-water molecules.
acid and a methanol [9]. Wolfe and co-workers confirmed that an acetic acid molecule can play a role in acting as both the proton donor and the acceptor. In this case, the hydrolysis of methyl acetate proceeds in only one step of reaction [41,42]. Figure 2 depicts the potential surface and the geometry transformation in the hydrolysis of ethyl benzoate calculated at the B3LYP/ 6-311++G(d,p) level of theory for the 2W model (Run 2). According to the results of Run 2, the C(1)–O(6) length increased along the reaction coordinates from 1.291 Å (s = 5.0 amu1/2 Bohr) 1/2 to 1.649 Å (s = 5.0 amu Bohr). Moreover, the length of the H(5)– O(6) bond decreased from 2.456 to 1.013 Å. These changes in the two lengths indicate that the cleavage of the C(1)–O(6) bond and the formation of the H(5)–O(6) bond proceed simultaneously. An optimization from the structure at s = 5.0 amu1/2 Bohr produced the protonated ethanol and benzoic acid. The results indicate that the TD intermediate does not form in the hydrolysis of ethyl benzoate. Therefore, the hydrolysis of the ethyl benzoate proceeds via
100
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
the one step mechanism. The other results of the 3W models also showed a similar transformation using the IRC calculations. 4.2. Two waters model There are six significant bonds between ethyl benzoate and two explicit-water molecules in the TS structures of the 2W model, TS(2W), as shown in Figure 1. In the gas phase, the inclusion of diffuse functions in the level of theory changes the geometries of the TS structure only slightly, so that they are all very similar to each other. For example, the geometrical parameters calculated at the B3LYP/6-311G(d,p) level of theory (Run 1) are different from those calculated in Run 2 within 0.014 Å. The largest difference is seen in the H(5)–O(6) length in Table 1. The roles of diffuse functions are the same in producing the TS structures calculated by the SCRF method. The geometry calculated at the CPCM-B3LYP/6-311G(d,p) level of theory (Run 4) are very similar to that of CPCM-B3LYP/6-311++G(d,p) level of theory (Run 5), different within 0.024 Å. For example the O(2)–H(3) length of Run 4 was calculated to be 1.747 Å while that of Run 5 turned out to be 1.771 Å. Therefore, it can be concluded that the TS structures are independent of the inclusion of diffuse functions. However, the TS structures produced by the SCRF method are obviously different from those of gas-phase calculations, especially the C(1)–O(2) and H(5)–O(6) lengths. For example, the lengths in Run 2 are 1.478 and 1.651 Å while those in Run 5 are 1.431 and 1.820 Å. The former differs from the later by 0.047 and 0.169 Å, respectively. Therefore, the implicit effects of solvent using the SCRF method do influence the TS structures of this reaction. As will be discussed later, the inclusion of a third explicit-water molecule changes this length markedly, i.e., the effects of the explicit solvent are more important than those of the implicit solvent included in the SCRF calculations. For the reactants in the 2W model in Figure 4, R(2W), the geometries depended on the diffuse functions very much. The bond lengths are also listed in Table 1. The diffuse functions cause the distances within each molecule to lengthen. For example, the C(1)–O(2) and H(3)–O(4) lengths of Run 2 are 2.786 and 1.846 Å, longer than those of Run 1 by 0.140 and 0.050 Å, respectively.
Similar elongations are also seen in the reactants in the SCRF method. The C(1)–O(2) length of Run 5 was determined to be 2.792 Å, longer than that of Run 4 by 0.164 Å. Moreover, the H(3)–O(4) length of the former (1.839 Å) is longer than that of the latter (1.805) by 0.034 Å. As the C(1)–O(2) and H(3)–O(4) lengths of Run 5 are very close to those of Run 2, meaning that the implicit effects of solvent using the SCRF method have no influence on the reactant geometry. The diffuse functions affected not only the geometry but also the energy. The activation free energy in the gas phase (DGàgas) of Run 1 without these functions is computed to be 26.0 kcal mol1, lower than that of Run 2 with the functions by 6.0 kcal mol1, as shown in Figure 5. Delchev and co-workers reported that the diffuse functions have strong effects on the energy calculations of the B3LYP method but are negligible for the BLYP method in the case of the complex ion [Cu(H2O)]+ [43]. The present results do seem to be consistent with those previously reported. In order to elucidate how large are the effects that the diffuse functions cause, the DGàgas value was calculated at the B3LYP/6311++G(d,p)//B3LYP/6-311G(d,p) level of theory (Run 3). That value of Run 3 was computed to be 29.6 kcal mol1, which is larger by 3.6 kcal mol1 than that of Run 1. The activation free energy of solvation (DGàsol) was then considered to obtain the activation free energy in aqueous solution, DGàcal. The QM/MC/FEP calculations produced large negative DGàsol values, depending on the perturbed structures, at each level of theory. For example, those of Run 1 and Run 2 were calculated to be 13.1 and 14.0 kcal mol1. Their difference is only 0.9 kcal mol1. Also, the DGàcal values of Run 1 and Run 2 in Figure 5 were calculated to be 12.9 and 18.0 kcal mol1, respectively. The value of Run 1 is lower by 5.1 kcal mol1 in comparison with that of Run 2. The DGàcal value of Run 3 was obtained as 16.5 kcal mol1. It is higher than that of Run 1 by 3.6 kcal mol1, but lower than that of Run 2 by only 1.5 kcal mol1. The DGàcal values of the 2W model calculated by the QM/MC/ FEP method are very consistent with those of the SCRF method. For example, the DGàcal values of Run 4 and Run 5 in Figure 5 are 9.2 and 17.4 kcal mol1, respectively. The DGàcal value of Run 5 differs from that of Run 2 by only 0.6 kcal mol1. As the observed
Table 1 Optimized geometries of 2W, 3WA and 3WB models. The unit of length is Å. Run
Model
Level of theory
1 2 4 5 6 7
2W 2W 2W 2W 3WA 3WB
B3LYP/6-311G(d,p) B3LYP/6-311++G(d,p) CPCM-B3LYP/6-311G(d,p) CPCM-B3LYP/6-311++G(d,p) B3LYP/6-311++G(d,p) B3LYP/6-311++G(d,p)
Reactant
TS structure
C(1)–O(2)
H(3)–O(4)
H(8)–O(10)
C(1)–O(2)
C(1)–O(6)
O(2)–H(3)
H(3)–O(4)
O(4)–H(5)
H(5)–O(6)
2.646 2.786 2.628 2.792 2.865 2.854
1.796 1.846 1.805 1.839 1.861 1.857
– – – – 1.632 1.653
1.478 1.478 1.432 1.431 1.432 1.491
1.419 1.420 1.421 1.419 1.451 1.433
1.458 1.466 1.747 1.771 1.488 1.459
1.059 1.057 1.005 1.004 1.053 1.06
1.016 1.020 0.999 0.999 1.046 1.031
1.666 1.651 1.804 1.820 1.529 1.586
Figure 4. Reactant structures of 2W, 3WA and 3WB models.
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
101
Figure 5. Activation free energies of hydrolysis of ethyl benzoate in gas phase and in aqueous at 298 K.
value (DGàexp) was 26.0 kcal mol1 at 298 K, the DGàcal values are not consistent with the experimental one. 4.3. Three waters models Due to the very large activation free energy of solvation in the 2W model, the reaction mechanisms were investigated more deeply. In TS(2W), the C(1) atom possesses two hydroxyl groups. One of the groups consists of the O(2) and H(9) atoms. These two atoms play a very important role in donating the proton H(3) to allow the reaction to proceed. It is easy to understand that the positions of the O(2) and H(9) atoms are changed significantly along the reaction coordinate by comparing the structures of R(2W) and TS(2W). The position of the other hydroxyl group consisting of O(7) and H(8) atoms in the TS(2W) is also very different from that in R(2W). After analyzing the population of Run 2 by the Merz-Singh– Kollman (MK) option, the charges of the O(2) and H(9) atoms were computed to be 0.69 and 0.50 in TS(2W) and they were found to be 0.77 and 0.37 in R(2W), respectively. The difference of the charge in the O(2) atom is 0.08. However, it is smaller than that of proton H(9), which is 0.13. Not only the positions, but also the charges of those atoms are dramatically changed, especially the proton H(9). In comparison, the charge of the proton H(8) was obtained to be 0.49 in TS(2W), but was calculated to be 0.45 in R(2W), different by only 0.04. Taking account of the positions and the charge differences of the hydroxyl groups along the reaction coordinate, it is possible that the implicit effects of solvent might not be included properly in the calculations of the 2W model. This resulted in overestimating the DGàsol value by as much as 14.0 kcal mol1 in Run 2. Therefore, one more explicit-water molecule has to be involved to interact with the hydroxyl groups in the TS structure in order to calculate the solvent effects more appropriately. In order to show that the position of a third explicit-water molecule plays a crucial role in calculating the solvent effects, two mechanisms including the third explicit-water molecule, 3WA (Run 6) and 3WB (Run 7), were proposed. An additional water molecule interacts with the proton H(9) in TS(2W) in the former mechanism whereas the water is added to interact with the proton H(8) rather than the proton H(9) in the latter. For both of the two types of three waters model, the calculations adopted only the B3LYP/6311++G(d,p) level of theory in the investigation. In the TS structure of the 3WA model, TS(3WA), the O(10) atom of the third explicit-water molecule interacts with the proton H(9), as depicted in Figure 1. The H(9)–O(10) length was calculated to be 1.664 Å. For TS(3WA), the geometrical parameters are significantly different from those of TS(2W) as shown in Table 1. For example,
the C(1)–O(2) length of Run 6 was computed to be 1.432 Å, shorter by 0.046 than that of Run 2. Moreover, the H(5)–O(6) length of the former turned out to be 1.529 Å, which is also shorter by 0.122 Å than that of the latter. On the other hand, the O(10) atom forms a hydrogen bond with the proton H(8), as shown by the H(8)–O(10) length (1.690 Å) in the TS structure of the 3WB model, TS(3WB). While the geometrical parameters of TS(3WA) are much different from those of TS(2W), the parameters of TS(3WB) are not much different to those of the latter. For example, the H(5)–O(6) length of Run 7 is computed to be 1.586 Å, shorter by 0.065 Å than that of Run 2. It is the largest difference between those two TS structures listed in Table 1. Moreover, there are many bond lengths in TS(3WB) which are longer than those in the TS(3WA). For example, the C(1)–O(2) and H(5)– O(6) lengths of the former were found to be 1.491 and 1.586 Å, while those of the latter turned out to be 1.432 and 1.529 Å. Even though the geometrical parameters of TS(3WA) and TS(3WB) are obviously different, the reactants of the three waters models, R(3WA) and R(3WB), are very similar within 0.021 Å, as shown in Figure 4. The biggest difference can be seen in the H(8)–O(10) length in Table 1. Furthermore, compared to the length of R(2W) in Run 2, the C(1)–O(2) lengths of R(3WA) (2.865 Å) and R(3WB) (2.854 Å) are longer by 0.079 and 0.068 Å, respectively. Consequently, the influence of the third water is to lengthen the distance between the ethyl benzoate and the first two waters. As plotted in Figure 5, the DGàgas values of Run 6 and Run 7 are calculated to be 33.6 and 35.0 kcal mol1. TS(3WA) is more stable than TS(3WB) by 0.8 kcal mol1. Figure 6 displays the DGàsol values of 3W models calculated using the QM/MC/FEP method. The DGàsol values of Runs 6 and 7 were computed to be 7.2 and 14.9 kcal mol1, respectively. Those DGàsol values are different by 7.7 kcal mol1. While the DGàsol value of 3WA model is far from that of the 2W model (Run 2), the value of the 3WB model is very close to that of the 2W model by less than 1.0 kcal mol1. Therefore, the position of the third water plays an essential role in producing the DGàsol values based on the QM/MC/FEP method. Corresponding to the DGàsol values, the DGàcal values of Runs 6 and 7 were 26.4 and 20.1 kcal mol1, respectively. We could not obtain the corresponding values of the models using the SCRF method since the optimizations for the TS structures were not successful. This seems to be a disadvantage of the SCRF method. 4.4. Comparison of DGàcal and DGàexp The DGàcal values of the 2W model were different from experiment (DGàexp) by 8.0 kcal mol1 at the B3LYP/6-311++G(d,p) level
102
T. Kaweetirawatt et al. / Chemical Physics Letters 547 (2012) 97–102
Figure 6. Free energies of solvation (DGàsol) of 3WA and 3WB models at 298 K.
of theory. It means that two explicit-water molecules are adequate to make the reaction proceed via their TS structures but they are still not enough to represent the solvent effects of the reaction with regard to the DGàcal value. The position of the third explicit-water molecule is very important, as explained earlier. If the third water molecule was placed in the wrong position, the DGàcal value in the 3WB model was still not consistent with the experimental value. It differs from the observed free energy by more than 5.0 kcal mol1. Both the 2W and 3WB models seem to have the same problem, which is a large negative value of DGàsol. This propensity apparently comes from large changes in the charge and positions of the O(2)–H(9) group. On the other hand, when the third water molecule has been placed in the right position in the 3WA model, it assists the process of reaction and produces a value of DGàsol that is more reasonable. The outcome is that the DGàcal value for the 3WA model, 26.4 kcal mol1, is very consistent with the experimental free energy within less than 1.0 kcal mol1. 5. Conclusion Including diffuse functions in QM/MC/FEP calculations of the hydrolysis of ethyl benzoate directly affects the reactant geometries, especially the distance between the ethyl benzoate and two explicit-water molecules. The functions lengthen the distances between ethyl benzoate and the water molecules. However, these functions have only a small influence on the TS structures. The DGàgas value calculated by including the diffuse functions is higher than that when that are excluded. Two explicit-water molecules are not enough to correctly represent the mechanism of hydrolysis of ethyl benzoate since the DGàcal values calculated by both the QM/MC/FEP and the CPCM-SCRF method are not consistent with experiment. Therefore, a third explicit-water molecule is necessary. The third water molecule increases the distance between the ethyl benzoate and the first two water molecules in the reactant geometries. Moreover, the position of the third water in the TS structure plays a vital role in producing a DGàcal value from the QM/MC/FEP calculations that is highly consistent with the experimental value. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2012. 06.001.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
S.M. Steinberg, F. Lena, Water Res. 29 (1995) 965. E. Furusjo, L.G. Danielsson, Chemom. Intell. Lab. Syst. 50 (2000) 63. E. Furusjo, L.G. Danielsson, Anal. Chim. Acta 373 (1998) 83. Y. Zhao, G. Wang, W. Li, Z.L. Zhu, Chemom. Intell. Lab. Syst. 82 (2006) 193. Z.L. Zhu, J. Xia, J. Zhang, T.H. Li, Anal. Chim. Acta 454 (2002) 21. P. Patnaik, M. Yang, E. Powers, Water Res. 35 (2001) 1587. S.H. Hilal, United States Environmental Protection Agency, 2006. (EPA/600/R06/105). E.W. Timm, C.N. Hinshelwood, J. Chem. Soc. 3 (1938) 862. K. Hori et al., Tetrahedron 63 (2007) 1264. A. Okur, C. Simmerling, Annual Reports in Computational Chemistry, vol. 2, Elsevier, Amsterdam, 2006. B. Roux, T. Simonson, Biophys. Chem. 79 (1999) 1. M. Orozco, F.J. Luque, Chem. Rev. 100 (2000) 4187. G. Brancato, N. Rega, V. Barone, J. Chem. Phys. 128 (2008) 144501. E.F. Lima, J.W.M. Carneiro, C. Fenollar-Ferrer, S. Miertus, S. Zinoviev, N.C. Om Tapanes, D.A.G. Aranda, Fuel 89 (2010) 685. H. Freedman, H.N. Nguyen, T.N. Truong, J. Phys. Chem. B 108 (2004) 19043. H. Freedman, T.N. Truong, J. Phys. Chem. B 109 (2005) 4726. J.D. Madura, W.L. Jorgensen, J. Am. Chem. Soc. 108 (1986) 2517. J. Chandrasekhar, S.F. Smith, W.L. Jorgensen, J. Am. Chem. Soc. 106 (1984) 3049. J. Chandrasekhar, S.F. Smith, W.L. Jorgensen, J. Am. Chem. Soc. 107 (1985) 154. W.L. Jorgensen, J.F. Blake, D.C. Lim, D.L. Severance, J. Chem. Soc., Faraday Trans. 20 (1994) 1727. Y. Wu, Y. Xue, D.Q. Xie, C.K. Kim, F.S. Yan, J. Phys. Chem. B 111 (2007) 2357. T. Yamaguchi, M. Sumimoto, K. Hori, J. Chem. Phys. Lett. 460 (2008) 331. T. Yamaguchi, M. Sumimoto, K. Hori, J. Comput.-Aided Chem. 9 (2008) 62. K. Hori, T. Yamaguchi, K. Uezu, M. Sumimoto, J. Comput. Chem. 32 (2011) 778. T. Kaweetirawatt, T. Yamaguchi, H. Tsutomu, M. Sumimoto, K. Hori, J. Phys. Org. Chem. (2012), http://dx.doi.org/10.1002/poc.2944. D. Boese, J.M.L. Martin, J. Chem. Phys. 119 (2003) 3005. M.J. Frisch et al., GAUSSIAN 03 (Revision C.02), Gaussian Inc., Wallingford CT, 2004. L.C. Mether, D.V. Sagar, S.N. Naik, Renew. Sustain. Energy Rev. 10 (2006) 248. F. Ma, L.D. Clements, M.A. Hanna, Trans. ASAE 41 (1998) 1261. K. Fukui, Acc. Chem. Res. 14 (1981) 363. C. Gonzalez, H.B. Schlegel, J. Chem. Phys. 90 (1989) 2154. H. Zipse, G. Apaydin, K.N. Houk, J. Am. Chem. Soc. 117 (1995) 8608. G. Kaminski, W.L. Jorgensen, J. Phys. Chem. 102 (1998) 1787. O. Accevedo, W.L. Jorgensen, J. Am. Chem. Soc. 127 (2005) 8829. O. Accevedo, W.L. Jorgensen, J. Am. Chem. Soc. 128 (2006) 6141. O. Tapia, Theoret. Chim. Acta (Berl.) 47 (1978) 157. M. Chen, Z. Lin, J. Chem. Phys. 127 (2007) 154314. C.K. Kim, B. Park, H.W. Lee, C.K. Kim, Bull. Korean Chem. Soc. 32 (2011) 1985. Y. Takano, K.N. Houk, J. Chem. Theory Comput. 1 (2005) 70. M.B. Smith, J. March, March’s Advanced Organic Chemistry, fifth ed., 2001. Y. Hsieh, N. Weinberg, S. Wolfe, Can. J. Chem. 87 (2009) 539. Z. Shi, Y. Hsieh, N. Weinberg, S. Wolfe, Can. J. Chem. 87 (2009) 544. V.B. Delchev, H. Mikosch, J. Struct. Chem. 47 (2006) 979.