Revisiting the free energy profile for the nucleophilic attack of hydroxide on formamide in aqueous solution

Chemical Physics Letters 422 (2006) 210–217 www.elsevier.com/locate/cplett

Revisiting the free energy profile for the nucleophilic attack of hydroxide on formamide in aqueous solution Jochen Blumberger *, Michael L. Klein Center for Molecular Modeling and Department of Chemistry, University of Pennsylvania, 231 S. 34th Street, Philadelphia, PA 19104-6323, United States Received 15 November 2005; in final form 10 February 2006 Available online 13 March 2006

Abstract The free energy profile for the first step of formamide hydrolysis in alkaline aqueous solution is computed using Car–Parrinello molecular dynamics simulation combined with umbrella sampling. An activation free energy of 22.3 kcal/mol and a reaction free energy of 17.3 kcal/mol are estimated from the potential of mean force in unprecedented good agreement with the corresponding free enthalpies obtained from experiment, 21.2 and 16.5 kcal/mol, respectively. The errors introduced by the BLYP density functional and pseudopotentials are investigated as well as the approximations made to relate the potential of mean force to the free enthalpies measured in experiment.  2006 Elsevier B.V. All rights reserved.

1. Introduction Alkaline hydrolysis of amides is a reaction of outstanding importance in organic chemistry, biochemistry as well as in the chemical industry. In recent years, amide hydrolysis in basic aqueous solution has also been the focus of many computational studies [1–9]. One reason for the growing interest of studying amide hydrolysis is its important role as reference reaction for enzymatic hydrolysis of peptide bonds. Moreover, because many amide hydrolysis reactions are experimentally well characterized, they can be used for validation or calibration of electronic structure calculations and solvation models. The first and usually rate determining step of alkaline amide hydrolysis is the nucleophilic attack of hydroxide on the carbon atom of the carbonyl group which results in the formation of the tetrahedral intermediate (TI) (see Scheme 1). Several theoretical investigations aimed at computing activation free energies for TI formation, but only a few could reproduce experimental results reasonably well (see *

Corresponding author. Fax: +1 215 573 6233. E-mail address: [email protected] (J. Blumberger).

0009-2614/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.02.035

Table 1). In the earliest attempt Weiner et al. [1], computed the energy profile for TI formation using ab initio gas phase calculation and a molecular mechanics model for solvation. The activation energy obtained deviated from the experimental activation free enthalpy energy by less than 1 kcal/ mol which is a surprisingly good result given the crude model for solvation that could be afforded at that time. However, if entropy effects had been included the deviation would have been 3–4 kcal/mol. To date, the best agreement with experimental data were obtained by Strajbl et al. [4] and more recently by Pliego [8]. Both authors performed ab initio gas phase calculations and used implicit solvation models. Interestingly, all computations carried out using explicit solvation models deviated significantly more from experiment. Bakowies et al. [3] who used force field based molecular dynamics simulation in combination with free energy perturbation methods overestimated the free energy barrier by 6.1 kcal/mol. Deviations of similar magnitude but opposite sign were also reported for recent Car–Parrinello MD (CPMD) simulations by Zahn [6] and Cascella et al. [7] Car–Parrinello MD has been successfully applied to model solvation and reactive events in the liquid phase. In particular the solvation structure of aqueous hydroxide was predicted [10–15] in excellent agreement with latest

J. Blumberger, M.L. Klein / Chemical Physics Letters 422 (2006) 210–217 H

H

HO- +

C O

1.

HO

H2N

C

O-

NH2

2.

H

NH3 +

C

O

-O

TI Scheme 1.

211

force to experimental activation and reaction free enthalpy. Anticipating our result, we obtain an estimate for these quantities in good agreement with experiment (deviation about 1 kcal/mol). Our computational approach is discussed and compared to previous MD simulations. Further avenues of investigation are given at the end of this letter. 2. Computational method

Table 1 Activation free energy DAà and reaction free energy DA for the formation of the tetrahedral intermediate (TI) in basic aqueous solution. The substrate is formamide unless stated otherwise

Experiment Weinerc Bakowiesd Strajble Pliegog Zahnh Cascellak This workl,m This workl,n This workl,o

DAà (kcal/mol)

DA (kcal/mol)

21.2a 21.7 27.3 19 23.4 15.8i 15 19.0 24.0 22.3

16.5b 16.0 13f 12j 9 16.6 22.9 17.3

a

Ref. [32]. DGà at T = 25 C. Calculated from experimental data. We have taken the reaction enthalpy measured for N,N-dimethylformamide at T = 25 C and pH 14, DG = 17.9 kcal/mol, [35] and subtracted 1.4 kcal/mol. This free energy correction is the difference of the experimental activation free enthalpy between formamide, 21.2 kcal/mol, [32] and N,N-dimethylformamide, 22.6 kcal/mol [35]. c Ref. [1]. Estimates of energy barrier and energy difference. MP2, MM water. d Ref. [3]. MP2/6-31 G*, classical force field. e Ref. [4]. B3LYP/AUG-cc-pVDZ, Langevin dipoles. f Estimated from Fig. 3 of Ref. [4]. g Ref. [8]. MP2 and CCSD, polarizable continuum model. h Ref. [6]. CPMD/BLYP, N-methylacetamide + Na+ + OH in 33 water. i Converted from 66 kJ/mol to kcal/mol. j Estimated from Fig. 1 of Ref. [6]. k Ref. [7]. CPMD/BLYP, formamide + OH in 55 water, steered MD. l CPMD/BLYP, formamide + OH in 29 water, umbrella sampling. m Determined from the potential of mean force, Fig. 3b. n Calculated using Eqs. (8 and 6), respectively. o Calculated using Eqs. (8 and 6), respectively, BLYP density functional error corrected. b

2.1. Umbrella sampling The free energy profile [or potential of mean force (PMF)] for TI formation was computed using the distance rOh between the carbon atom of formamide and the oxygen atom of hydroxide as reaction coordinate. Sampling of this coordinate was enhanced using the umbrella sampling technique [17]. The bias potential U was chosen to have a similar form as the free energy profile along the minimum free energy path obtained from a 2-dimensional metadynamics simulation, [18] see Fig. 1. The minimum free energy path was replaced by the coordinate rOh and the profile was scaled by a factor 1.8897. We have sampled the reaction coordinate using 10 windows with harmonic restraining potentials Vi = ki/ 2(rOh  rOh,i)2, ki/105 = 3, 3, 3, 3, 3, 3, 3, 6, 24, 12 a.u., rOh,i = 4.19, 3.94, 4.21, 3.88, 3.56, 3.06, 2.73, 2.08, 1.89, ˚ for windows i = 1–10. Hence, the total bias poten1.26 A tial for each window i added to the unperturbed Hamiltonian is Ui = U + Vi. The corresponding biased distributions pi of the coordinate rOh are then given by pi ðrOh Þ ¼ hdðr0Oh  rOh ÞiEþU i

ð1Þ

where E denotes the unperturbed DFT potential energy of the system, E + Ui the effective potential energy from which the ionic forces are derived and Æ    æ denotes the canonical average. The biased distributions are then combined using the weighted histogram analysis method [19] (WHAM) to yield the optimal unbiased distribution p(rOh) and free energy profile A(rOh), respectively. AðrOh Þ ¼ k B T lnhdðr0Oh  rOh ÞiE ¼ k B T ln pðrOh Þ

ð2Þ

In Eq. (2) kB is the Boltzmann constant and T the temperature and b = 1/(kBT). neutron diffraction data [16]. The large discrepancy for the free energy barrier of amide hydrolysis is therefore rather surprising and calls for an analysis and a fresh investigation of this problem. First objective of this work is a validation of pseudopotentials and the BLYP density functional used in the condensed phase simulation. This is done by comparing the gas phase reaction energy with results of all electron calculations at different levels of theory. Using a similar simulation protocol as in the two previous CPMD studies of Refs. [6,7] but a different sampling technique, we compute the free energy profile for the first step of formamide hydrolysis in alkaline aqueous solution. We then discuss the approximations made to relate the computed potential of mean

2.2. Reaction free energy and activation free energy The standard reaction free enthalpy DG0 for the reaction of formamide (A) and hydroxide (B) to TI (C) in aqueous solution is related to the potential of mean force A(n) as follows [20].  0  Z C rA rB 0 DG ¼ k B T ln dnIðnÞJ n expðbAðnÞÞ 4p rC þ pDV

ð3Þ

In general, A(n), n = (r, /, h, U, H, X), AðnÞ ¼ k B T lnhdðn0  nÞiE þ const;

ð4Þ

212

J. Blumberger, M.L. Klein / Chemical Physics Letters 422 (2006) 210–217

U (kcal/mol)

40

20

0 1

2

3 rOh (Å)

4

5

Fig. 1. Bias potential U used for umbrella sampling of the reaction coordinate rOh.

not only depends on the separation distance r but on the position and orientation of B relative to A, determined by the spherical coordinates (r, /, h) and Euler angles (U, H, X), respectively. The constant in Eq. (4) has to be chosen so that the potential of mean force vanishes in the limit of infinite separation, limr!1A(n) = 0. In Eq. (3), C0 is the standard concentration 1 mol dm3 = 1/ ˚ 3), ri the symmetry number of solute 1660 (molecules A molecule i, b = 1/(kBT), Jn the Jacobian determinant, I(n) a step function which is unity if n is within configuration space assigned to C and 0 otherwise, and pDV is the reversible work due to change of molar volumes of solvation. Note that Eq. (3) slightly differs from the one given in Ref. [20]. The steric factor in the denominator obtained by integration over all 3 Euler angles, 8p2 in Ref. [20], was replaced by 4p because of the absence of a third rotation axis for hydroxide. The position and orientation of hydroxide relative to formamide are defined by introducing a body fixed frame centered on the C atom of formamide. The x-axis is defined by the vector connecting C and O atom of formamide and the y-axis by the component of the vector connecting C and N atom that is orthogonal to the x-axis. The z-axis is chosen to be orthogonal to x- and y-axes. The position of hydroxide is then defined by the distance rOh ” r, the angle / between x-axis and the projection of the vector along the C–Oh axis onto the x, y plane and the angle h between z-axis and the vector along the C–Oh axis. Similarly, the orientation of hydroxide is defined by the angle U between x-axis and the projection of the molecular axis of hydroxide onto the x, y plane, and the angle H between the z-axis and the molecular axis of hydroxide, X ” 0. Eq. (3) is a generalization of the simpler expression [21– 23] Z rA rB 0 0 DG ¼ k B T ln C 4p drIðrÞr2 expðbAðrÞÞ rC þ pDV ;

ð5Þ

which is valid only if A and B are spherically symmetric species (e.g., NaCl). In this special case, the reversible work

required to pull B apart from A does not depend on relative orientations. Hence, integration over angles is trivial and yields an explicit expression for the standard reaction free enthalpy as a function of the reversible work along r. By contrast, the nucleophilic attack studied here occurs at preferred orientations as dictated by anisotropic orbital interactions between hydroxide and formamide. In this case A(n) is a complicated 5-dimensional (5D) free energy function, which in principle could be computed from the corresponding probability distribution obtained from molecular dynamics simulation (MD). However, in practice the MD estimate for A(n) usually gives a poor approximation due to sampling problems and the high dimensionality of the distribution makes enhanced sampling methods impracticable. In order to relate DG0 to the radial PMF A(rOh) obtained from umbrella sampling we assume that in a range of angles [a1, a2], a = /, h, U, H A(n) is essentially independent on the orientation of hydroxide and equal to A(rOh). Outside this range A(n) is assumed to have significantly larger values giving zero contribution to the exponential on the right-hand side (RHS) of Eq. (3). The integral Eq. (3) can then be separated and estimated as follows DG0  DA0  k B T ln

Z

rTS

drOh r2Oh

 expðbAðrOh ÞÞ  k B T

0

 ln D  k B T ln C 0 ;

ð6Þ

where D¼

1 4p

Z

/2

Z

Z

h2

U2

Z

H2

d/ dh sin h /1

h1

dU dH sin H; U1

ð7Þ

H1

is a steric factor and limrOh !1 AðrOh Þ ¼ 0. The volume term is usually small at ambient pressure, i.e., DA0  DG0. All symmetry numbers were set equal to unity because the hydrogen atoms in A and C as well as the oxygen atoms in C are separated in configuration space by high barriers and are therefore distinguished in molecular dynamics simulations. Assignment of configuration space to product and reactant is unambiguous for the reaction studied here (as opposed to the pKa calculation of Ref. [23]) because the free energy curve exhibits a well defined maximum determining the transition state at a separation distance rTS (see Fig. 3). Hence, I(n) = 1 if rOh < rTS and zero otherwise. The range of angles used to compute the steric factor Eq. (7) and an estimation of the error introduced by the approximation Eq. (6) are discussed in Section 3. The activation free energy for formation of the tetrahedral intermediate, DAà, is estimated using the same approximation as for the computation of the reaction free energy. Adopting Eq. (6) the separation distance is fixed to rTS and the integration over rOh is omitted. DAz  AðrTS Þ  2k B T ln rTS  k B T ln D  k B T ln C 0

ð8Þ

J. Blumberger, M.L. Klein / Chemical Physics Letters 422 (2006) 210–217

2.3. Simulation details An alkaline aqueous solution of formamide is modeled by a periodically repeated cubic unit cell comprising one formamide molecule, one hydroxide ion and 29 water mol˚ which corecules. The box length was set equal to 10.032 A 3 responds to a density of 963 kg m . The net charge of the solution was compensated by a neutralizing homogeneous background charge. The biased distributions pi(rOh) are obtained using the CPMD code [24] which was modified to allow for the calculation of bias potential Ui and corresponding bias forces. During equilibration and production runs we have used a chain of Nose–Hoover thermostats [25] for the ions with a target temperature of 300 K, a MD time step of 5 a.u. (0.1209 fs) and a fictitious electron mass of 700 a.u. The initial configuration for each window i was taken from a snapshot of a reactive metadynamics trajectory. After discarding the first 2.5 ps of dynamics the distributions were averaged over trajectories of length 8– 10 ps. In order to prevent proton transfer between hydroxide and solvent molecules, the relative momentum vector between hydrogen atom H of a solvent molecule and the oxygen atom Oh of hydroxide was reflected if the H–Oh ˚ . This constraint is active distance was smaller than 1.3 A only in the rare event of an attempted hydrogen transfer and preserves total momentum and energy. An alternative coordination number constraint would have generated artificial forces at each MD step, which might have seriously affected the equilibrium structure of hydroxide. 2.4. Electronic structure method The pseudopotentials used in CPMD simulations and gas phase optimizations were constructed according to the Troullier–Martins scheme [26]. The 2s and 2p electrons of all second row atoms were included in the valence represented by l = s and l = p components with equal pseudization radii rc = 1.23, 1.12, 1.05 a.u. for C, N and O atoms. For hydrogen atoms we used a local s potential with rc = 0.5 a.u. The Kleinman–Bylander [27] transformation was used for O, N and C atoms. The orbitals were expanded in plane waves with a reciprocal space kinetic energy cutoff of 70 Ry unless stated otherwise. The exchange correlation energy was calculated using the generalized gradient approximation (GGA) for exchange according to Becke [28] and GGA for correlation according to Lee, Yang and Parr [29]. Gas phase geometry optimizations of formamide, hydroxide and TI were carried out with the CPMD code for two reciprocal space cutoffs, 70 and 130 Ry, until the maximum component of the nuclear energy gradient was less than 104 Hartree/bohr. The gas phase reaction energy was calculated also for the analytic Goedecker–Hutter pseudopotentials [30] (at 100 Ry) using the optimized structures obtained with Troullier–Martins pseudopotentials at 70 Ry.

213

The pseudopotentials and density functional were validated with all electron calculations in the gas phase using the GAUSSIAN program package [31]. The structures of formamide, hydroxide and TI were optimized at the BLYP, B3LYP and MP2 levels of theory using the 6-31 + G* basis set and applying the default convergence criteria. The final reaction energy was computed using the AUG-cc-pVTZ basis set. 3. Results and discussion 3.1. Gas phase calculations The results of the gas phase calculations are summarized in Table 2. We find that the pseudopotentials generated according to Troullier and Martins perform exceptionally well judging from the comparison with all electron (AE) calculations at the BLYP/AUG-cc-pVTZ level of theory. The dissociation energy DE(1) of TI into formamide and hydroxide is well converged at the relatively low cutoff 70 Ry and differs from the AE value by less than 0.1 kcal/mol. Similarly, the bond lengths of TI, formamide and of hydroxide are converged at 70 Ry and deviate by ˚ from the values of AE calculations at the less than 0.01 A BLYP/6-31 + G* level of theory. The only exception is ˚ larger the C–Oh bond length rOh which is at 70 Ry 0.03 A than the AE value. Increase of the cutoff to 130 Ry reduces ˚ . We have further tested the the overestimation to 0.01 A performance of Goedecker–Hutter pseudopotentials and find also in this case good agreement with AE calculation for the dissociation energy of TI. The performance of the BLYP functional in describing the TI dissociation in the gas phase is not as good. The dissociation energy obtained for CPMD/BLYP at 70 Ry is underestimated by 5.6 kcal/ mol compared to MP2/AUG-cc-pVTZ level of theory, ˚ and the C–Oh bond length is overestimated by 0.11 A (wrt MP2/6-31 + G*). Inclusion of exact exchange

Table 2 Gas phase calculations for the tetrahedral intermediate (TI) HCOOHNH 2 Method

rmin Oh

DE(1.86)

DE(1)

BLYP/TM/cut70 BLYP/TM/cut130 BLYP/GH/cut100 BLYP/6-31 + G* BLYP/AUG-cc-pVTZ B3LYP/6-31 + G* B3LYP/AUG-cc-pVTZ MP2/6-31 + G* MP2/AUG-cc-pVTZ

1.644 1.624

0.76

1.616

1.34 0.90 3.59 3.05 4.51 4.71

16.69 16.83 16.87 23.14 16.77 25.55 20.17 26.52 22.31

1.545 1.538

The equilibrium bond length rmin Oh between the C atom and the oxygen atom of OH was obtained from energy minimization of TI. DEðyÞ ¼ EðyÞ  Eðrmin Oh Þ denotes the optimized potential energy at a dis˚ relative to the energy minimum at rmin . The values for tance rOh = y A Oh DE(1) are counterpoise corrected. TM and GH stand for Troullier– Martins and Goedecker–Hutter pseudopotentials, respectively, and cut x denotes a plane wave kinetic energy cutoff x Rydberg. Bond lengths are ˚ and energies in kcal/mol. given in A

J. Blumberger, M.L. Klein / Chemical Physics Letters 422 (2006) 210–217

improves both, dissociation energy and bond lengths. The deviation at the B3LYP level of theory wrt MP2 is signifi˚. cantly smaller, 2.1 kcal/mol and 0.01 A In the gas phase, the dissociation is barrierless and the potential energy increases steadily with increasing separation distance rOh. This is in contrast to the aqueous phase where the barrier for dissociation of TI is about 5 kcal/ mol. In order to obtain an estimate for the error of the energy barrier in aqueous solution, we have calculated the energy difference between fully optimized TI in the ˚, gas phase and TI with rOh constrained to rTS = 1.86 A DE(1.86) in Table 2. The distance rTS corresponds to the free energy barrier of dissociation in solution (see Fig. 3b). The energy required for stretching of the C–Oh ˚ is bond from the respective minimum value to 1.86 A underestimated by 3.9 and 1.7 kcal/mol for CPMD/BLYP and B3LYP relative to MP2. The decrease of potential energy for the reverse process of bringing the isolated reac˚ is then underestimated tants to a C–Oh distance of 1.86 A by 5.6  3.9 = 1.7 and 2.1  1.7 = 0.4 kcal/mol for CPMD/BLYP and B3LYP, respectively. Note that in solution the reactant level (formamide + hydroxide) shifts below the one of TI. Therefore, on the basis of our DFT error analysis, the barrier for TI formation is expected to be 1.7 kcal/mol higher for CPMD/BLYP than for MP2, whereas the barrier for dissociation of TI is expected to be underestimated by 3.9 kcal/mol. However, these estimates only include the potential energy error in the gas phase due to the BLYP functional relative to MP2/AUGcc-pVTZ. Error contributions from entropy and hydrogen bonding with solvent molecules are not included. 3.2. Free energy calculation in solution 3.2.1. Potential of mean force The distributions of the order parameter rOh obtained from the 10 window potentials are shown in Fig. 2. They overlap sufficiently well giving a smooth free energy curve ˚, over the entire range of distances sampled, 1.35–4.25 A see Fig. 3b. Averaging over 8–10 ps per window was found to be sufficient because the distributions pi and the free energy curve A did not change significantly when calculated for trajectories of half length (see dotted lines in Figs. 2 and 3b). In fact, the change of the free energy barrier computed from half the statistical data were less than 0.2 kcal/mol. As one can see in Fig. 3 the free energy curve has a very small slope at the maximum separation ˚ confirming that rm was distance sampled, rm = 4.25 A chosen sufficiently large. The onset of the curve is at a dis˚ , the maximum corresponding to the tance of about 3.8 A ˚ and the minitransition state is located at rTS = 1.86 A mum corresponding to the tetrahedral intermediate (TI) ˚. at 1.58 A The part of the free energy curve describing the nucleophilic attack of hydroxide is well approximated by a parabola, see Fig. 3b. The correlation coefficient of the best quadratic fit function A = k/2(r  r0)2 for the interval

˚ is 0.998 and the force constant k is equal to 2.20–3.85 A ˚ 2. The part of the product well describ12.1 kcal mol1 A ing dissociation of TI is found to be linear rather than quadratic. The correlation coefficient of the fit function ˚ is 0.996 and A = a + (k/2)r for the interval 1.58–1.83 A ˚ 1. the force constant k is equal to 19.0 kcal mol1 A A comparison between Figs. 1 and 3b shows that the barrier height of minus the bias potential is twice as large as the one of the PMF. The bias potential therefore overcompensated the underlying free energy curve which explains the need of relatively many windows and the large force constants used for the harmonic restraining potentials Vi. This shows, however, that the PMF can be calculated accurately at affordable cost (see below) even if the bias potential is not chosen in an optimal way. 3.2.2. Reaction free energy and activation free energy The free energy difference obtained from the minima of the PMF corresponding to reactant and product, ˚ )  A(4.25 A ˚ ) = 16.6 kcal/mol, is virtually identiA(1.58 A cal with the reaction free enthalpy estimated from experiments, 16.5 kcal/mol (see Table 1). However, as explained in Section 2, for non-spherical reactants the experimental standard reaction free enthalpy cannot be simply related to the radial potential of mean force (Eq. (5)). The correct relation is given by Eq. (3) which we approximated by Eq. (6). Integration of the radial PMF over rOh, first term on the RHS of Eq. (6) yields a free energy difference of 17.2 kcal/mol, increasing the value obtained directly from the PMF by 0.6 kcal/mol. A slightly larger contribution comes from the steric factor, second term on the RHS of Eq. (6), which is estimated to be 1.3 kcal/mol (see below). The largest contribution is due to the concentration term, third term on the RHS of Eq. (6), which increases the free energy difference by 4.4 kcal/mol, to 22.9 kcal/mol.

10 9

8

0.1 probability

214

7

6

5

4 3 2 1

0 1

2

3 rOh (Å)

4

Fig. 2. Biased distributions of the distance rOh between the C atom of formamide and the oxygen atom of hydroxide. In the reactant state, the system comprises of 1 formamide molecule 1 hydroxide ion and 29 water molecules. The distributions were obtained from 10 biased Car–Parrinello MD trajectories of length 8–10 ps (solid and dashed lines for windows with odd and even index number) and 4–5 ps (dotted lines).

J. Blumberger, M.L. Klein / Chemical Physics Letters 422 (2006) 210–217

(a)

30

20

~ ~

E (kcal/mol)

40

A (kcal/mol)

(b)

A (kcal/mol)

20

10

0 2

10

3

4

rOh (Å)

0 1

2

3 rOh (Å)

4

5

Fig. 3. (a) Potential energy profile for the dissociation of the tetrahedral  intermediate (TI) in the gas phase, HCOOHNH 2 ! HCONH2 þ OH . The data are taken from Table 2, CPMD/TM/cut70 (+), BLYP/AUG-ccpVTZ//BLYP/6-31 + G* (s), B3LYP/AUG-cc-pVTZ//B3LYP/6-31 + G* (·), MP2/Aug-cc-pVTZ//MP2/6-31 + G* (h), rOh is the distance between carbon atom and the oxygen atom of OH. The potential energy of TI was set equal to the value of the free energy in aqueous solution. The energies ˚ are computed for infinite separation displayed at a distance rOh > 4 A using counterpoise correction to the basis set superposition error. (b) Free energy curve for the formation of the tetrahedral intermediate (TI) from formamide and hydroxide in aqueous solution. The free energy curve is obtained from the biased distributions shown in Fig. 2 using the weighted histogram analysis method (WHAM). The free energy curves depicted in solid (dotted) lines are obtained from trajectories of length 8–10 (4–5) ps. The inset shows the best quadratic fit (dashed line) of the free energy curve in the reactant well (solid line).

The additional terms involved increase the free energy difference by 6.3 kcal/mol relative to the value obtained directly from the PMF and seem to introduce a fairly large discrepancy with experiment. However, the value of 22.9 kcal/mol includes the error due to the BLYP density functional which overestimates the gas phase reaction energy for formation of TI by 5.6 kcal/mol compared to MP2 level of theory. Correcting for this artifact we obtain a free energy difference of 17.3 kcal/mol which is in good agreement with the experimental value. Thus, the error introduced by the density functional is largely canceled by the additional free energy terms of Eq. (6) explaining the agreement between the free energy difference obtained directly from the PMF and the experimental value. Similarly to reaction free energy, the free energy barrier for formation of TI compares very well with the experimental value. The barrier taken from the maximum of the PMF ˚ is 19.0 kcal/mol. Inclusion of the additional at rOh = 1.86 A

215

terms on the RHS of Eq. (8) increases the barrier by 5.0 kcal/mol while correction for the density functional error reduces the barrier by 1.7 kcal/mol giving a final value of 22.3 kcal/mol in good agreement with the experimental activation free enthalpy of 21.2 kcal/mol [32]. The fairly good agreement between density functional error corrected free energy difference and experiment might be due to a favorable approximation of the steric factor Eq. (7). However, as explained in the following, the uncertainty estimated for the corresponding steric free energy (1.3 kcal/ mol, second term on the RHS of Eq. (6)) is reasonably small and bound to ±(2–3) kcal/mol. The value 1.3 kcal/mol was obtained by taking the average of minimum and maximum steric free energy estimated to be 1.5 and 4.2 kcal/mol, respectively. The minimum value is obtained by integration over the full spheres, i.e., /1 = U1 = 0, /2 = U2 = 2p and h1 = H1 = 0, h2 = H2 = p. In this case A(n) of Eq. (3) is assumed to be equal to A(rOh) for all orientations of hydroxide, even for the least preferred ones, thereby providing a lower bound for activation and reaction free energies. The maximum value of 4.2 kcal/mol is estimated from the root-mean-square fluctuations ra of the angles a, a = /, h, U, H, as obtained from 10 ps of umbrella sampling at the transition state region, window 8 in Fig. 2, / = 232 ± 9.5, h = 28.8 ± 5.3, U = 288 ± 13.8, H = 87.2 ± 11.4 (all values in degree). The integration intervals in Eq. (7) are a1 = Æaæ  ra,a2 = Æaæ + ra and the steric factor obtained was multiplied by a factor of two to account for the possibility of nucleophilic attack from either side of the molecular plane of formamide. Assuming that A(n) is close to A(rOh) in the range Æaæ ± ra this choice provides an upper bound for activation and reaction free energies because all contributions from angles outside the small intervals ±ra are omitted. Thus, the uncertainty of free energies due to separation of radial and angular contributions, i.e., due to the approximation of Eq. (3) by Eq. (6), is expected to be ±1/2(4.2(1.5))  ± (2  3) kcal/mol. The uncertainty of the free energies due to the limited system size used is expected to be small. In the present simulation, the hydroxide ion has two solvation shells, which are assumed to account for most of the solvation free energy barrier of an infinitely diluted system. One might expect that barrier and free energy difference slightly increase with larger system size. An increase of the number of solvent molecules is likely to stabilize the hydroxide ion relative to the transition state because in the former the charge is more localized than in the latter. This effect should be rather small though, because the net charge does not change along the reaction path and the coordinated water molecules effectively screen the charge distribution. Indeed, as reported in Ref. [33], the PMF for dissociation of aqueous NaCl is virtually size independent in the range 106–862 water molecules per ion pair. Hence, while the PMF computed here for a 29 water molecule system might change slightly if a third solvation shell is included, we expect no major changes due to higher solvation shells and bulk solvation.

216

J. Blumberger, M.L. Klein / Chemical Physics Letters 422 (2006) 210–217

A further more subtle source for underestimation of the free energy barrier is due to the assumption of a 1-dimensional transition state [33,34]. Using transition path sampling the authors of Ref. [33] showed that the barrier for dissociation of the solvated ion pair NaCl was underestimated by about 1.7 kBT if the ion–ion separation distance was used as reaction coordinate. It was shown that in addition to separation distance preorganization of solvent molecules determines whether a configuration belongs to the transition state ensemble. The calculations of Ref. [33] suggest that the free energy required to reorganize the solvent is smaller for the ensemble of configurations that comprise the 1-dimensional transition state relative to the ones that comprise the transition path ensemble. One can expect that this shortcoming of a 1-dimensional solute coordinate is a general feature for chemical reactions in solution including formation of TI which could lead to an underestimation of activation free energy barrier by 1 kcal/mol or more. 3.2.3. Comparison to work of others The free energies obtained in this work are in better agreement with experiment than the ones reported in two previous CPMD simulation studies of alkaline hydrolysis of amides [6,7]. Zahn [6] reported, a free energy barrier of 66 kJ/mol = 15.8 kcal/mol for TI formation of Nmethylacetamide in basic aqueous solution. The experimental barrier for the similar N,N-dimethylacetamide is 24.1 kcal/mol [35] suggesting that the simulation of Ref. [6] underestimates the barrier by about 5 kcal/mol (taking into account the additional free energy terms of Eq. (8) and the error of the BLYP density functional which were not included in the estimate reported in Ref. [6]). This discrepancy is probably related to the small C–Oh distance of ˚ chosen as the ‘equilibrium value’ for the reactant 2.75 A state. At this distance the hydroxide ion is not fully solvated which leads to an underestimation for the free energy of desolvation, and consequently for the formation of the transition state. Indeed, the free energy required to decrease the C–Oh distance from the equilibrium value ˚ , to a distance 2.75 A ˚ is accordchosen in this work, 4.25 A ing to our simulation about 6.5 kcal/mol (see Fig. 3b). Another source of uncertainty is introduced by the presence of a counter ion Na+ which significantly perturbs the solvation structure of hydroxide in a small simulation cell as used in Ref. [6] (similar size as used here). We further point out that the system dimensions chosen in Ref. [6] give a density of 1.5 kg dm3, significantly higher than the experimental density of aqueous N-methylacetamide of 1.0 kg dm3. The reason for deviations between our study and the CPMD study of Cascella et al. [7] is less clear. The authors report a free energy barrier of 15 kcal/mol as obtained from the PMF for alkaline hydrolysis of formamide which is 4 kcal/mol lower than our estimate, 19.0 kcal/mol. Taking into account the additional free energy terms of Eq. (8) and the error of the BLYP density functional the free energy barrier obtained in Ref. [7] increases to 18.3 kcal/

mol. Thus, the deviation with the experimental value, 21.2 kcal/mol, is 2 kcal/mol larger than in our simulation and of opposite sign. Different too, is the location of the ˚ , a distance transition state reported to be at 2.20 A ˚ 0.34 A larger than our value. The system specification in Ref. [7] is identical with the one used here except for the number of water molecules which was 55/cell (29/cell in this work but similar density). This difference is not the reason for the discrepancy with our result because one would expect an increase rather than a decrease of the barrier with increased number of solvent molecules. Also the different pseudopotential used for hydrogen and the slightly larger fictitious electron mass used in Ref. [7] are not expected to have a significant effect. The deviation is most likely related to different sampling techniques used to obtain the free energy profile. Indeed, using steered molecular dynamics the authors of Ref. [7] could only enforce the reverse reaction, dissociation of TI, but not the nucleophilic attack of hydroxide. The lowest ˚ /ps) was still too high for pulling speed applied (0.2 A enforcing desolvation of hydroxide. The pulling speed was possibly also too high for reformation of the solvation shell of hydroxide in the reverse reaction, as indicated by the too small decrease of free energy after dissociation of TI. 4. Conclusion Computing a total of about 100 ps of biased Car–Parrinello MD trajectories, we obtained an estimate for reaction free enthalpy and activation free enthalpy of the first step of formamide hydrolysis in good agreement with experiment. The deviation with the experimental free energies is about 1 kcal/mol. However, due to non-spherical symmetry of the reactants, the free energies could be calculated only approximately from the radial potential of mean force. The uncertainty of our results is estimated to be 2– 3 kcal/mol. Judging from calculations in the gas phase the pseudopotential approach used in the plane wave implementation of CPMD and the moderate plane wave cutoff 70 Ry do not introduce any significant errors. However, all electron calculations showed that the BLYP density functional underestimates the dissociation energy by 5.6 kcal/mol relative to MP2 level of theory. The free energies in solution were corrected for this artifact. The satisfactory agreement with the experimental values suggests that solvation of hydroxide and in particular the energetics of desolvation of hydroxide is well described by the BLYP functional and the relatively small aqueous model system used. In the gas phase, the nucleophilic attack of hydroxide on formamide is steady down hill in energy. Therefore, the free energy barrier in aqueous solution has often been referred to as ‘completely solvent induced’ [1,4,7]. In a forthcoming publication [18], we analyze solvent dynamics, desolvation and solute reorganization along the entire reaction path. Moreover, we address the question whether hydrolysis occurs via a direct nucleophilic attack – as

J. Blumberger, M.L. Klein / Chemical Physics Letters 422 (2006) 210–217

assumed in this work and all previous computational studies – or via a general base mechanism. Acknowledgements We want to thank Dr. Bernd Ensing for sharing his umbrella sampling subroutine, for many fruitful discussions and for proofreading the manuscript. The referee of this manuscript is acknowledged for valuable comments. NIH is acknowledged for financial support. References [1] S.J. Weiner, U. Chandra Singh, P.A. Kollman, J. Am. Chem. Soc. 107 (1985) 2219. [2] Y.-J. Zheng, R.L. Ornstein, J. Mol. Struct. (Theochem) 429 (1998) 41. [3] D. Bakowies, P.A. Kollman, J. Am. Chem. Soc. 121 (1999) 5712. [4] M. Sˇtrajbl, J. Floria´n, A. Warshel, J. Am. Chem. Soc. 122 (2000) 5354. [5] X. Lopez, J. Ina˜ki Mujika, G.M. Blackburn, M. Karplus, J. Phys. Chem. A 107 (2003) 2304. [6] D. Zahn, Chem. Phys. Lett. 383 (2004) 134. [7] M. Cascella, S. Raugei, P. Carloni, J. Phys. Chem. B 108 (2004) 369. [8] J.R. Pliego Jr., Chem. Phys. 306 (2004) 273. [9] D. Cheshmedzhieva, S. Ilieva, B. Galabov, J. Mol. Struct. (Theochem) 681 (2004) 105. [10] M.E. Tuckerman, K. Laasonen, M. Sprik, M. Parrinello, J. Phys. Chem. 99 (1995) 5749. [11] M.E. Tuckerman, K. Laasonen, M. Sprik, M. Parrinello, J. Chem. Phys. 103 (1995) 150. [12] D. Marx, M.E. Tuckerman, J. Hutter, M. Parrinello, Nature 397 (1999) 601. [13] M.E. Tuckerman, D. Marx, M. Parrinello, Nature 417 (2002) 925.

217

[14] Z. Zhu, M.E. Tuckerman, J. Phys. Chem. B 106 (2002) 8009. [15] B. Chen, J.M. Park, I. Ivanov, G. Tabacchi, M.L. Klein, M. Parrinello, J. Am. Chem. Soc. 124 (2002) 8534. [16] A. Botti, F. Bruni, S. Imberti, M.A. Ricci, A.K. Soper, J. Mol. Liq. 117 (2005) 81. [17] G.M. Torrie, J.P. Valleau, Chem. Phys. Lett. 28 (1974) 578. [18] J. Blumberger, B. Ensing, M.L. Klein, Angew. Chem. Int. Ed., accepted. [19] S. Kumar, D. Bouzida, R.H. Swendsen, P.A. Kollman, J.M. Rosenberg, J. Comput. Chem. 13 (1992) 1011. [20] M.K. Gilson, J.A. Given, B.L. Bush, J.A. McCammon, Biophys. J. 72 (1997) 1047. [21] D. Chandler (Ed.), Introduction to Modern Statistical Mechanics, Oxford University Press, Oxford, 1987. [22] J. Pranata, W.L. Jorgensen, Tetrahedron 47 (1991) 2491. [23] J.E. Davies, N.L. Doltsinis, A.J. Kirby, C.D. Roussev, M. Sprik, J. Am. Chem. Soc. 124 (2002) 6594. [24] CPMD Version 3.10.0, The CPMD consortium, MPI fu¨r Festko¨rperforschung and the IBM Zurich Research Laboratory 2005. Available from: . [25] G.J. Martyna, M.L. Klein, M. Tuckerman, J. Chem. Phys. 97 (1992) 2635. [26] N. Troullier, J. Martins, Phys. Rev. B 43 (1991) 1993. [27] L. Kleinman, D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. [28] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [29] C. Lee, W. Yang, R. Parr, Phys. Rev. B 37 (1988) 785. [30] C. Hartwigsen, S. Goedecker, J. Hutter, Phys. Rev. B 58 (1998) 3641. [31] GAUSSIAN 03, Revision C.01, Frisch, M.J. et al., Gaussian, Inc., Wallingford CT, 2004. [32] H. Slebocka-Tilk, F. Sauriol, M. Monette, R.S. Brown, Can. J. Chem. 80 (2002) 1343. [33] P.L. Geissler, C. Dellago, D. Chandler, J. Phys. Chem. B 103 (1999) 3706. [34] R.P. Muller, A. Warshel, J. Phys. Chem. 99 (1995) 17516. [35] J.P. Guthrie, J. Am. Chem. Soc. 96 (1974) 3608.