J. Phys. Chem. B 2007, 111, 4453-4459
4453
Ab Initio Molecular Dynamics Study of Carbon Dioxide and Bicarbonate Hydration and the Nucleophilic Attack of Hydroxide on CO2 Kevin Leung,*,† Ida M. B. Nielsen,‡ and Ira Kurtz§ Sandia National Laboratories, MS 1415, Albuquerque, New Mexico 87185, Sandia National Laboratories, MS 9158, LiVermore, California 94551, and UniVersity of California at Los Angeles DiVision of Nephrology, Los Angeles, California 90095 ReceiVed: December 10, 2006; In Final Form: February 16, 2007
We apply ab initio molecular dynamics (AIMD) to study the hydration structures of the carbon dioxide molecule and the bicarbonate and carbonate anions in liquid water. We also compute the free energy change associated with the nucleophilic attack of the hydroxide ion on carbon dioxide. CO2 behaves like a hydrophobic species and exhibits weak interactions with water molecules. The bicarbonate and carbonate ions are strongly hydrated and coordinate to an average of 6.9 and 8.7 water molecules, respectively. The energetics for the reaction in the gas phase are investigated using density functional theory and second-order Møller-Plesset perturbation theory (MP2) in conjunction with high-quality basis sets. Using umbrella sampling techniques, we compute the standard state, aqueous phase free energy difference associated with the reaction CO2 + OH- f HCO3after correcting AIMD energies with MP2 results. Our predictions are in good agreement with experiments. The hydration structures along the reaction coordinate, which give rise to a predicted 9.7 kcal/mol standard state free energy barrier, are further analyzed.
I. Introduction When carbon dioxide (CO2) dissolves in water, it forms varying amounts of carbonic acid (H2CO3), carbonate (CO32-), and bicarbonate (HCO3-) anions:
CO2 + OH- T HCO3-
(1)
CO2 + H2O T H2CO3
(2)
HCO3- + OH- T CO32- + H2O
(3)
HCO3- + H+ T H2CO3
(4)
CO32- + H+ T HCO3-
(5)
These species are involved in pH homeostasis in the ocean1 and in atmospheric water,2 geological processes such as the dissolution of calcite,3 and biological processes such as pH regulation and transepithelial transport by membrane transporters.4,5 The equilibrium population of hydrated CO2, HCO3-, and CO32- in water, as well as the rates and mechanisms of their interconversion, is strongly pH dependent.6-11 Above pH 8-9, eq 1 is believed to supersede eq 2 as the predominant mechanism for nucleophilic attack of CO2.6,8 Infrared experiments on dissolved CO2 in liquid water suggest that the carbon atom in CO2 coordinates to the oxygen of H2O molecules.12 Density functional theory (DFT) calculations performed in the gas phase have confirmed that the most stable CO2-H2O complex involves a CCO2-Ow coordination (the subscript “w” refers to water), although the binding energy is * Corresponding author. E-mail:
[email protected]. † Sandia National Laboratories, Albuquerque. ‡ Sandia National Laboratories, Livermore. § University of California at Los Angeles Division of Nephrology.
small, on the order of 2.9 kcal/mol.13 On the other hand, there have been surprisingly few theoretical studies of the CO2 hydration structure in liquid water, although supercritical carbon dioxide/water interfaces, pertinent to “green” solvent industrial applications, have received more attention.14 Classical force field simulations have shown that the CO2 hydration structure exhibits hydrophobic ordering, but the results are sensitive to the size of the Lennard-Jones parameters.15 While eq 1 exhibits no reaction barrier in the gas phase,16 a considerable enthalpic barrier has been reported in liquid water.8 Thus, hydration effects on the reactant, product, and transition state complexes are of paramount importance to the reaction thermodynamics and kinetics. Increasing the rate of OH- attack on CO2 is, in fact, an important biological function provided by carbonic anhydrase enzymes.17-20 Previous theoretical work has applied ab initio and DFT methods to model the CO2OH- reaction complex,16,21 while the aqueous environment has been treated with various approximations. The dielectric continuum approximation of water has so far been unable to predict reasonable standard state free energy differences (∆G°) and barriers (∆G°†) associated with the reaction.22,23 Quantum mechanics/molecular mechanics (QM/MM) free energy perturbation treatment, using the nonpolarizable TIP3P and other models for water, proves to be more successful, but it has been suggested that water polarizability may resolve some discrepancy between modeling and experiments.23-25 In this work, we will apply the ab initio molecular dynamics (AIMD) technique to elucidate the CO2 hydration structure in water, as well as those of HCO3- and CO32-. All molecules, including water in the bulk region, are treated using DFT. DFT accounts for water polarizability and has been shown to yield hydration stuctures for organic anions that agree with experiments, even where widely used force fields fail.26 Understanding these hydration structures is crucial for modeling the chemical
10.1021/jp068475l CCC: $37.00 © 2007 American Chemical Society Published on Web 04/05/2007
4454 J. Phys. Chem. B, Vol. 111, No. 17, 2007 reactions eqs 1-5 as well as future simulations of the pKa’s of bicarbonate and carbonic acid. We also use AIMD-based potential of mean force (PMF) methods to compute ∆G° and ∆G°† associated with eq 1. Equation 1 is one of the simplest nucleophilic reactions in liquid water, and our study will help establish an important benchmark for AIMD simulation of chemical reactions. Studying this reaction is pertinent not only to biological and industrial processes,27 but it also yields predictive capability toward manipulating CO2 reaction rates in confined aqueous environments such as protein cavities4,5 and, potentially, in nanopores. Analogous AIMD simulations of eqs 2-5 will be performed in the future. To date, there have only been a handful of explicit AIMD PMF calculations in bulk liquid water environments. All such studies apply umbrella sampling and related techniques.28 Most involve proton transfer reactions and the making/breaking of N-H and O-H covalent bonds.29-33 In such cases, widely applied DFT exchange correlation functionals predict barriers and zero temperature reaction energies within a few kilocalories per mole of high-level quantum chemistry results (see the appendix of ref 32). When the breaking and making of other types of covalent bonds are involved,34 achieving quantitative agreement with experiments is more challenging. A recent AIMD study predicts a free energy barrier of nucleophilic attack of OH- on formamide that is in good agreement with experiments,35 provided that gas phase DFT results are corrected with second-order Møller-Plesset perturbation (MP2)36,37 predictions. The choice of the DFT exchange correlation functional will affect not only the eq 1 reaction energetics13,16,21,22 but also the structure and self-diffusion constant of liquid water38-40 and the hydration structure of one of the key reactants (OH-). Both the Hammer et al. revised Perdew-Burke-Ernzerhof (RPBE)41 and the Becke-Lee-Yang-Parr (BLYP)42,43 functionals predict four-coordinated OH- in liquid water (with one of the coordinating H2O molecules possibly being an outer shell contribution), while the Perdew-Burke-Ernzerhof (PBE)44 functional predicts 3 H2O in the first shell of the anion.45-49 The higher coordination number appears more consistent with the experimental diffusion constant. We will adopt the RPBE functional in this work and correct the energetics with quantum chemistry predictions.32,35 The paper is organized as follows. Section II describes the method used. The hydration structures of different species are depicted in section III. Section IV describes the potential of mean force calculations and the hydration structure along the reaction coordinate for eq 1, and section V concludes the paper with further discussion of the usage of AIMD in modeling chemical reactions in liquid water. II. Methods AIMD calculations are performed using the Vienna ab initio Simulation Package (VASP),50 the projector augmented wave method,51 and associated pseudopotentials.52 The RPBE calculations apply the PBE pseudopotentials;53 as discussed above, this functional appears to predict the correct coordination number for the hydroxide ion. As we will correct our RPBE energies using MP2 predictions, the effect of using PBE rather than RPBE pseudopotentials will be suitably adjusted in the end. The simulation cell is cubic with a linear dimension of 12.417 Å and contains 63 water molecules and a solute (either CO2, HCO3- with the C-OH bond in various stages of dissociation, CO32-, or OH-). The deuterium mass is substituted for all protons to reduce the time step size (0.5 fs for CO2 and HCO3-
Leung et al. hydration studies, 0.25 fs in all other cases) needed for energy conservation in this Born Oppenheimer dynamics AIMD simulation. This substitution does not affect the equilibrium properties. A 400 eV energy cutoff and a 10-6 eV accuracy convergence criterion are enforced at each time step. We have verified that increasing the energy cutoff in the VASP RPBE calculations changes the gas phase reaction energies by less than 0.2 kcal/mol. A thermostat keeps the systems at T ) 300 K; at this condition, RPBE predicts reasonable liquid water structure.38 The energy drifts in all trajectories are ∼1 K/ps. The effect of this drift on the potential of mean force will be discussed in section IV-B1. The trajectory lengths are 28, 13, and 9 ps for CO2, HCO3-, and CO32- after applying 2 ps equilibration runs on configurations obtained from classical force field molecular dynamics simulations. A longer, 5 ps AIMD equilibration run using a 0.5 fs time step is performed for OH-, after which statistics are collected for 8 ps using a 0.25 fs time step. For the nucleophilic reaction eq 1, the distance between CCO2 and OOH- is chosen as the reaction coordinate RCO. To compute the potential of mean force, ∆W(RCO), umbrella sampling constraints, Ui(RCO) ) Ai(RCO - Ri)2, are applied in 12 sampling windows. Ri ranges from 1.5 to 3.5 Å. Ai ) 2 or 3 eV/Å2, except that Ai ) 6 eV/Å2 in the window where ∆W(RCO) exhibits the sharpest rise. The trajectory used for each hydration study/ umbrella sampling window is at least 9 or 10 ps, after equilibration for 1-2 ps using a larger time step; the initial configuration is taken from a snapshot of an adjoining window. In windows where Ri > 2.2 Å, a penalty function, V(ROHw) ) B(Ro - ROHw)4 is applied whenever any water proton, Hw, is less than a distance of Ro ) 1.5 Å from the OH- oxygen. Here, B ) 50.0 eV. This prevents proton transfer between bulk water and the designated OH- and permits the construction of an unambiguous reaction coordinate.54 The effect of this bias is removed by reweighing the energy of each configuration with the factor exp[βV(ROHw)], where β ) 1/kBT. In other words, the umbrella sampling histogram reweighing is performed using exp[βV(ROHw) + βU(RCO)]. For the dissociation of HCO3- (eq 1), we introduce a Cartesian coordinate system with the x-axis going through the CCO2 and OOH- atoms without loss of generality. We further constrain the CCO2, OOH-, and the two OCO2 atoms to lie in the x-y plane to eliminate free rotation about the x-axis.56 Otherwise, unlike for instance ref 35, all atoms are free to move. Therefore, to reference ∆G obtained from ∆W(RCO) to standard states (i.e., 1.00 molar reactants and product), we apply55
∆G° ) -kBT ln{4πC0
∫0R
dr I(r)r2 exp[-β∆W(r)]}
TST
(6)
Here, all lengths are in units of angstroms, C0 is the standard molar concentration, 1/(55.56 × 29.90 Å3), and I(r) is 1 - Θ(r - RTST), where Θ is the Heaviside function, and the top of the barrier of ∆W(RCO) is located at RCO ) RTST. Equation 6 omits all symmetry factors and rotational entropies of dissociated species. It applies if ∆W(RCO) embodies adequate sampling of all degrees of freedom.57 For the product state of eq 1, namely, HCO3-, the in-plane O-C-O wag (θO-O, Figure 1) and C-O-H bend (θOH) about the CCO2-OOH- axis are severely restricted bending vibrations, which are adequately sampled using AIMD. In the reactants (Figure 1b), however, CO2 and OH- exhibit free rotation. The relatively short AIMD trajectories are not suited to sampling such diffusive motions. Custom-made solutions to this problem need to be devised;35 even cluster-based quantum chemistry methods encounter this problem.24,25 As will be discussed, we predict that OH- is
Carbon Dioxide and Bicarbonate Hydration
J. Phys. Chem. B, Vol. 111, No. 17, 2007 4455
Figure 1. In-plane O-C-O wag, defined as the angle between the bisector of the O-C-O bond and an axis perpendicular to the C-O bond axis, and C-O-H bend (θO-O and θOH, respectively). As the dissociation reaction proceeds, both turn from bending vibrations (panel a) to free rotations (panel b).
collinear with the CCO2-OOH- axis as the dissociation proceeds. To estimate the rotational entropy contributions in the dissociated state, we impose two more potentials only in the umbrella sampling windows where dissociation has taken place (i.e., when ∆W(RCO) has reached a constant, asymptotic value): Y(θO-O) ) C(sin θO-O - 1)2 and Z(θOH) ) D(cos θOH - 1)2. During histogram reweighing to compute ∆W(RCO), new Boltzman factors associated with Y(θO-O) and Z(θOH) are applied in addition to those for V(ROHw) and U(RCO). These terms restrict the in-plane O-C-O wag and the C-O-H bend to within a small region of angular space, which can be sampled using AIMD. The potential of mean force is then corrected by comparing the cos θO-O and cos θOH space sampled against all available rotational space, namely, 2.00 in both cases.56 For validation, we also consider a “transitional” window where ∆W(RCO) has not attained its asymptotic value. Then, ∆W(RCO) should be independent of the presence/absence of Y(θO-O) and Z(θOH) as long as C and D are chosen not to overly restrict the natural bending motions. We find that C ) 6 eV and D ) 3 eV satisfy this criterion. We propose that this may be a general method for AIMD calculations of the free energy of molecular dissociation reactions. Similarly, to reference ∆G† to standard states, we use
∆G° ) -kBT †
2 ln{4πC0R TST δ
exp[-β∆W(RTST)]} (7)
where δ is by convention 1 Å.58 Gas phase calculations are performed for the species involved in eq 1 using the Gaussian 03 suite of programs59 and the augmented correlation-consistent basis sets aug-cc-pVDZ, augcc-pVTZ, and aug-cc-pVQZ.60 Optimum structures are determined with the Hartree-Fock37 (HF) method, density functional theory (B3LYP43,64 and PBE44), and second-order MøllerPlesset perturbation theory36,37 (MP2) employing the aug-ccpVDZ basis. Single point energies are computed at the optimum aug-cc-pVDZ structures using the aug-cc-pVTZ and aug-ccpVQZ sets. An estimate for the infinite basis set limit for the MP2 reaction energy for eq 1 is obtained by extrapolation of the HF and MP2 correlation energies. HF energies are extrapolated from the DZ, TZ, and QZ energies using the exponential form61 HF E HF X ) E ∞ + A exp(-BX)
(8)
where X is the highest angular momentum represented in the basis set. MP2 correlation energies are extrapolated from the TZ and QZ energies using the two-point formula62
E (2) ∞
)
3 (2) 43E (2) QZ - 3 E TZ
43 - 33
(9)
Corrections for basis set superposition error are included in several cases (cf. section IV-A) and are computed using the
Figure 2. Pair correlation function between water atoms and carbon dioxide atomic sites: (a) C; (b) O. Solid and dashed lines refer to the water O and H sites, respectively.
Figure 3. Pair correlation function between water atoms and bicarbonate atomic sites: (a) O(1); (b) O(2); (c) O(3); (d) H. Solid and dashed lines refer to the water O and H sites, respectively. HCO3- is found to coordinate to a total of 6.9 water molecules. Inset: labeling of the atoms in HCO3-. Within AIMD simulation time scales, all three oxygens in HCO3- are distinct. The reaction coordinate RCO is the distance between C and O(3).
counterpoise correction of Boys and Bernardi.63 A potential energy profile for eq 1 is determined by computation of the energy ∆E(RCO) for the CO2-OH- system for values of the reaction coordinate, RCO, ranging from 1.4 to 4.0 Å in 0.1 Å increments. The energy profile is computed at the HF, B3LYP, and PBE levels of theory using the aug-cc-pVDZ basis. For each RCO, constrained geometry optimizations are performed, fixing the value of RCO and optimizing all other degrees of freedom. Note that the aqueous phase results (next section) are performed with the RPBE functional; the PBE results are quoted just for comparison purposes. III. Hydration Results Figure 2 depicts the site-site pair correlation functions between carbon dioxide and water. These g(r) indicate a predominantly hydrophobic arrangement of water around CO2. The pair correlation functions between the carbon and oxygen sites of CO2 and the oxygen atoms in water, respectively, gCOw(r) and gOOw(r), are similar to the classical force field results of Panhuis et al., when the scaling of Lennard-Jones parameters is set to unity.15 The OCO2-Hw g(r) exhibits small values at r < 2.5 Å, which is a typical cutoff distance for counting hydrogen bonds associated with electronegative groups (see Figures 3 and 4 below). Spatial integration of this g(r) to r ) 2.5 Å yields 0.3 coordinating water molecules per oxygen, confirming its weak hydrogen bonding with water protons. In addition to possible hydrogen bonding, we examine the coordination between CO2 carbon and water oxygen atoms. AIMD and Panhuis et al. both predict broad gCOw(r) first peaks
4456 J. Phys. Chem. B, Vol. 111, No. 17, 2007
Leung et al.
Figure 4. Pair correlation function between water atoms and (a) carbonate and (b) hydroxide oxygen sites. Solid and dashed lines refer to the water O and H sites, respectively. The two ions exhibit total hydration numbers of 8.7 and 5.0. A penalty function has been applied to preserve the identity of the OH- ion (see text).
centered at r ) 3.95 and 3.85 Å, respectively, and a gOOw(r) first peak at r ) 3.35 Å. One subtle, but important, difference is that the AIMD gCOw(r) indicates an appreciably higher probability that an Ow is in the vicinity of CCO2. Thus, AIMD predicts a shoulder at r ) 3.10 Å, while ref 15 does not. Integrating these gCOw(r) reveals that, on average, AIMD and force fields predict 0.60 and 0.32 water molecules residing within r ) 3.1 Å of the CO2 carbon atom. This coordinate feature seen in AIMD simulations is consistent with infrared experiment inferences.12 At the same time, since the predicted CO2-H2O binding energy13 is much smaller than that of a water-water hydrogen bond, it is consistent that such coordination is manifested as a small shoulder in a liquid water environment. The g(r) between water and HCO3- are shown in Figure 3 (see the inset for the oxygen site labels). Within AIMD time scales, rotation about the C-O(3) bond axis is not observed, HCO3- remains a planar molecule throughout, and the three bicarbonate oxygens are distinct. Integrating gOHw(r) and gHOw(r) to a standard hydrogen bond cutoff distance of 2.5 Å, which roughly corresonds to the first minimum in many of these g(r), O(1), O(2), O(3), and H are found to coordinate with 2.4, 2.8, 0.7, and 1.0 H2O molecules, respectively, yielding a total hydration number of 6.9. Unlike O(1), O(2) is not partially shielded from water by a nearby -OH group, so its hydration number is slightly higher. The hydration structures of CO32- and OH- are depicted in Figure 4. The three oxygen atoms of the divalent carbonate anion are equivalent, and they exhibit strongly structured hydration shells containing approximately 2.9 water molecules each. Thus, CO32- exhibits a total coordination number of 8.7. Note that there is much less short-ranged coordination between C and water oxygen atoms for CO32- and HCO3- (not shown) than for CO2 (Figure 2a) because the hydrogen bonding between the two anions and water sterically hinders the weak C-Ow interaction. The g(r) between OH- and water is obtained by applying the aforementioned penalty function, V(ROHw) ) B(Ro - ROHw),4 on any water proton, Hw, whenever ROHw < Ro ) 1.5 Å, and then unbiasing the results by reweighing. As can be seen in Figure 4, this procedure yields a sharp, discontinuous cutoff in gOHw(r) at r ) 1.3 Å. Thus, the effect of imposing V(ROHw) is similar to the 1.3 Å hard sphere cutoff employed in ref 35, achieved by reflecting momenta at the ROH ) 1.3 Å boundary. Comparing the OH- gOHw(r) with that computed in unconstrained RPBE AIMD simulations,45 we find that the first peak height and position are similar in the two cases. The unconstrained simulation correlation function extends to ROH ) 1.17 Å instead of exhibiting a sharp cutoff at ROH ) 1.3 Å. The hydration number is Nw ) 4.0 in both cases, as is that reported in AIMD simulations using the BLYP functional.46 In addition,
Figure 5. (a) Gas phase potential energy as a function of RCO for eq 1. Solid line, B3LYP predictions; dashed line, Hartree Fock; dotdashed line, PBE. (b) Aqueous phase potential of mean force in the dissociation of HCO3- into CO2 and OH- as a function of the reaction coordinate RCO. A proton penalty function is applied to maintain the chemical identity of the OH- ion (see text). (c) Mean hydration number (Nw) of OOH- in each of the 11 umbrella sampling windows.
the OH- proton hydrogen bonds to an average of 0.3 water molecules in the constrained case (not shown). As the nucleophilic reaction eq 1 proceeds, the hydration environments of the CO2 and OH- oxygen atoms exhibit large changes. O(1) and O(2) proceed from each forming 0.3 hydrogen bonds (Figure 2b) to Nw ) 2.4 and 2.8 (Figure 3a and b), while O(3) evolves from Nw ) 4 (Figure 4b) to Nw ) 0.7 (Figure 3c). This leads to significant differences between the gas phase and aqueous phase reaction mechanisms, as will be discussed in the next section. IV. OH- Nucleophilic Attack on CO2 A. Gas Phase Results. First, we consider the potential energy ∆E(RCO) ) EHCO3-(RCO) - ECO2 - EOH- as RCO varies for the CO2-OH- complex in vacuum, optimizing all other degrees of freedom. Figure 5a shows that our ∆E(RCO) results for Hartree Fock, PBE, and B3LYP are all monotonic and do not exhibit any barrier, in agreement with ref 16. At the optimal HCO3geometry, the C-O(3)-H angle is 102°. When RCO ) 4 Å, this angle becomes 82°; in other words, the nascent OH- and CO2 remain almost parallel to each other. This is in agreement with previous quantum chemistry calculations where the CO2OH- complex is optimized in the gas phase,23,25 but as will be shown, AIMD simulations with explicit water solvent predict a significant difference in the OH- orientation. The zero point energy differences between HCO3- and the dissociated products, not included in Figure 5a, are 3.66, 3.65, 3.40, and 4.38 kcal/ mol, respectively, for MP2, B3LYP, PBE, and the HartreeFock methods. The attractive energy decays slowly to zero with RCO because of the charge- and dipole-induced dipole interactions. Table 1 compares the reaction energy for eq 1 predicted using several DFT exchange correlation functionals, including RPBE, PBE, and B3LYP, as well as the MP2 method. Our highestlevel computation of the reaction energy is the infinite basis set limit MP2 result, ∆E ) -45.0 kcal/mol, obtained by extrapolation of HF and MP2 correlation energies (cf. section II for details). This value is somewhat lower than a previous ∆E value25 of -44.1 kcal/mol, computed at the MP2 level with
Carbon Dioxide and Bicarbonate Hydration
J. Phys. Chem. B, Vol. 111, No. 17, 2007 4457
TABLE 1: Energy Changes in Vacuum (kcal/mol) Associated with eq 1a method ∆E:
RPBE
PBE
B3LYP
MP2
-37.5
-43.9
-45.6
-45.0
a Zero point energy corrections are excluded. RPBE and PBE results were computed using the VASP code and PAW pseudopotentials. The B3LYP value was computed with the aug-cc-pVTZ basis set, including the counterpoise correction, and the MP2 result is an estimate for the infinite basis set limit (see text for details).
TABLE 2: Basis Set Convergence of the Reaction Energy (kcal/mol) for eq 1 in the Gas Phasea aug-cc-pVnZ B3LYP HF MP2
aug-cc-pVnZ CPC
n)D
n)T
n)Q
n)D
n)T
n)Q
∞
-47.9 -52.1 -43.1
-46.2 -51.8 -44.7
-45.8 -51.6 -44.9
-46.5 -51.0 -38.1
-45.6 -51.6 -42.2
-45.7 -51.5 -43.6
NA -51.4 -45.0
a
Figure 6. (a) Side and (b) top views of a snapshot of the OH- hydration shell, containing four water molecules in a quasi-square planar geometry coordinated to the OH- oxygen, when RCO ) 3.16 Å. Black circle, C; small white circle, H; large white circles, O(OH-); striped circles, O(CO2) and O(water) (striped in opposite directions). In this snapshot, the proton on the OH- is not strongly coordinated to any water molecules. In panel a, the water molecule whose O atom is closest to the H atom in OH- (3 Å away) is also shown, although it is not part of the hydration shell.
TABLE 3: Energy Corrections (kcal/mol) Computed in the Gas Phase Applied to Correct the Free Energy Differences and Barriers Associated with eq 1a ∆G° ∆G°†
s.s.
ZPE
MP2
C-O-H/O-C-O angle
+4.01 +2.04
+3.66 +0.70
-7.50 -2.65
+2.58 +2.58
Zero point vibrational energies are not included. Optimum B3LYP/ aug-cc-pVDZ geometries were employed for the B3LYP computations, and MP2/aug-cc-pVDZ geometries were used for the HF and MP2 methods. Superscript CPC denotes counterpoise-corrected values, and ∞ represents the infinite basis set limit (see text for details).
a The four columns are, respectively, conversion to standard states (1 M reactant and product); zero point energy correction; difference between MP2 (extrapolated to infinite basis set) and RPBE (VASP); and orientation entropy due to OH- and CO2 (see text).
a 6-311++G** basis set, which is similar in quality to the augcc-pVTZ set. We note that reaction energies for eq 1 computed at the MP2 and MP4 level with the 6-311++G** basis set have been found to be in close agreement,25 suggesting that the correlation treatment at the MP2 level is probably adequate. To further investigate the basis set requirements for accurate computation of the reaction energy for eq 1, we list in Table 2 computed gas phase reaction energies, including counterpoisecorrected values, for a number of methods and basis sets as well as the extrapolated infinite basis set limits. Basis set convergence is rapid at the HF and B3LYP levels, as evidenced by the small DZ f TZ f QZ changes in the reaction energy as well as the small size of the counterpoise correction. At the MP2 level, the uncorrected reaction energy also converges rapidly with basis set improvement. The MP2 counterpoise correction, however, is sizable, even at the QZ level, and moreover, application of the correction produces energies further away from the estimated basis set limit. For this system, we then conclude that the counterpoise correction should not be applied at the MP2 level of theory and that results obtained with the aug-cc-pVTZ basis, excluding the counterpoise correction, are within a few tenths of a kilocalorie per mole of the estimated basis set limit. B. Aqueous Phase Results. 1. AIMD Potential of Mean Force. The aqueous phase potential of mean force, ∆W(RCO), is depicted in Figure 5b. Judging from the well depth of ∆W(RCO) at RCO ) 1.38 Å, ∆G for OH- + CO2 f HCO3- is predicted to be -12.5 ( 1.4 kcal/mol in liquid water. The cumulative statistical uncertainty is estimated by dividing the trajectory in each window into four equal parts and reporting twice the standard deviation. ∆G between reactants and products in water is strongly reduced from the gas-phase energy difference. This is due to the larger hydration free energy of the reactant OH- compared with HCO3-, which is a much larger anion. A free energy barrier now appears at RTST ) 2.08 Å, and it arises from the disruption of the OH- hydration shell. This is underscored by the drastic reduction of the number of water molecules coordinated to OOH- along the reaction coordinate (Figure 5c). The long-range attraction between OHand CO2 appears to vanish beyond RCO ) 3 Å due to the waterinduced screening of electrostatic interactions. Furthermore, the
bulk water molecules provide additional polarizable centers and make the electrostatic interactions isotropic. Note that, when 2.75 Å < RCO < 3.2 Å, OH- orientation sampling (θOH, Figure 1b) is extremely restricted. The snapshot in Figure 6 is taken from one of these sampling windows. They show that, for OH- to assume the bulk coordination geometry, that is, four H2O molecules in a quasi-square pyramidal geometry,45-49 its proton points away from the CO2, and the four first hydration shell H2O molecules are layered between the anion and CO2. Substantial deviation from this collinear OHorientation will apparently cause one of the four coordinating H2O molecules to overlap with CO2. This is in contrast with the gas phase results, discussed above, where the OH- and the CO2 remain almost parallel to each other as RCO varies. In this snapshot, the OH- proton is not strongly coordinated to any water molecules. We have also examined the effect of the drift in AIMD total energy on the computed potential of mean force by repeating umbrella sampling simulations for the windows 1.47 Å < 1.60 Å and 1.60 Å < 1.88 Å using a 0.5 fs time step. The total energy drift is now significantly larger, ∼7 K/ps. However, this larger drift has little effect on the change in PMF ∆∆G between the end points of the windows. For the first window, the 0.25 fs (0.50 fs) time step yields ∆∆G ) 5.68 kcal/mol (5.35 kcal/ mol). For the second window, ∆∆G ) 9.75 kcal/mol (9.81 kcal/ mol). The difference of using a larger time step is within numerical uncertainty. Since the PMF calculations reported in this work utilize the smaller (0.25 fs) time step, leading to ∼O(1) K/ps total energy drift, our test indicates that such a drift in total energy has little effect on the computed PMF depicted in Figure 5b. 2. Corrections to Standard State ∆G. Several corrections to the potential of mean force in Figure 5b are required; they are summarized in Table 3. First, we apply eq 6 to reference ∆G to standard states. Using RTST ) 2.08 Å, we obtain a correction of +4.0 kcal/mol. Next, we correct for the rotational entropy. When RCO > 3.39 Å (data associated with the last umbrella sampling window in Figure 5b), significant in-plane O-C-O wag and C-O-H bend occur with respect to the CCO2-OOHaxis if these motions are unconstrained (see Figure 1 and section II). We thus apply the Y(θO-O) and Z(θOH) potentials and add
4458 J. Phys. Chem. B, Vol. 111, No. 17, 2007 -kBT log[(2/xOH-)(2/xCO2)] to ∆G°. Here, xOH- and xCO2 are the fractions of angular space actually sampled by the C-O-H bend and in-plane O-C-O wag with the constraining potentials in place. We define x to accommodate 95 percentile of the raw angular distribution data, such that xOH- ) 0.08 and xCO2 ) 0.66 at T ) 300 K. This correction amounts to 2.6 kcal/mol. We adopt the MP2 prediction of +3.66 kcal/mol as the difference in zero point energy corrections between reactants and product (see section IV-A). Finally, errors arising from the RPBE functional and the associated pseudopotentials are estimated using gas phase data. We use the MP2 limit value of -45.0 kcal/mol (cf. section IV-A) as the standard and add the MP2-RPBE difference, -7.5 kcal/mol, to the aqueous phase ∆G° as a postprocessing correction for our AIMD potential of mean force as has been done in the AIMD literature.32,35,65 Adding these four corrections, we obtain a standard state free energy change for eq 1, ∆G° ) -9.8 kcal/mol. This estimate is in excellent agreement with the widely accepted experimental value of ∆G° ) -9.4 kcal/mol at T ) 298 K in dilute salt solutions.7,10 3. Corrections to Standard State ∆G†. Next, we consider the free energy barrier ∆G°†. Reported experimental values for the activation enthalpy vary8,11 but are generally within (2 kcal/ mol of each other. We adopt the ∆G°† ) 12.1 kcal/mol estimate fitted by Peng and Merz25 from experimental data.8 (Pohorecki and Moniuk11 have compiled the eq 2 rate constant measured by different authors in slightly different environments; from these reported values, the inferred ∆G°† agree to within 0.4 kcal/ mol.) Our AIMD/RPBE PMF calculations predict that eq 1 exhibits ∆G† )7.1 kcal/mol (Figure 5b). Corrections similar to those for ∆G° are applied, as follows. (1) To reference to standard states, we apply eq 7 which yields +2.0 kcal/mol. We also add the rotational entropy penalty for restricted OH- and CO2 rotations of 2.6 kcal/mol discussed above. (2) The zero point correction at the transition state can be estimated using ref 25, which devises a reduced dimensional Hamiltonian to circumvent the lack of a barrier in the gas phase. The Hartree-Fock predicted gas phase “enthalpy” at room temperature listed in ref 25 is mostly due to zero point corrections. At RCO ≈ 2.08 Å,66 this correction is ∼0.7 kcal/ mol.25 (3) The effect of using different functionals on ∆G°† is estimated by freezing RCO ≈ 2.08 Å66 and allowing all other degrees of freedom to relax. The CO2/OH- pair in vacuum is used as a reference. The difference between RPBE (VASP, plane wave basis) and MP2 (Gaussian, aug-cc-pVTZ basis) is found to be -2.65 kcal/mol. Adding these corrections, the overall ∆G°† becomes 9.7 kcal/ mol, a little below the experimental value of 12.1 kcal/mol.7 Interestingly, Merz et al. strongly oVerestimate the free energy barrier, predicting a value of 19.2 kcal/mol when neglecting both solvent polarizability and internal relaxation of the CO2OH- complex in the presence of solvent molecules. 4. Validity of Our Reaction Coordinate. We note that, while the free energy change ∆G° is a state function, independent of the reaction coordinate RCO, the free energy barrier ∆G°† depends on the reaction pathway. To assess whether RCO is a good reaction coordinate, we have applied standard chemical kinetics theory and examined the “transmission coefficient”67 by taking random configurations and velocities from trajectories in the barrier top sampling window where RCO crosses 2.10 Å but removing the umbrella constraints.66 All five random configurations where eq 1 proceeds in the reverse direction, that is, RCO exhibits an initial negative velocity as it crosses the dividing surface at 2.10 Å, have ROH monotonically decreasing,
Leung et al. yielding the HCO3- product without recrossing the barrier (i.e., RCO reaching 2.10 Å again). All five random configurations where eq 1 proceeds in the forward direction lead to monotonically increasing RCO, yielding CO2 and OH- without recrossing the dividing surface. From chemical kinetics theory,67 this suggests that our choice of a reaction coordinate and our estimate of ∆G°† are reasonable. However, the above discussion ignores the fact that we have applied a penalty function to preserve the identity of OH-, thus preventing proton hopping.46 The OH- hydration is sterically frustrated when it is in close proximity of a hydrophobic species like CO2. Thus, as OH- approaches CO2 from afar, it experiences an unfavorable hydration free energy barrier. If we release the constraint in windows where Ri > 2.2 Å, the OH- would readily extract a proton from a nearby H2O, effectively switching identity with a water molecule and creating a new OH- farther away from the CO2. In other words, proton hopping introduces an additional pathway that allows OH- to “roll downhill” from our computed AIMD ∆W(RCO). While proton hopping configurations constitute a minority of the AIMD trajectory46-48 and will not affect the predicted ∆G° by more than 1 kBT, proton hopping may increase ∆G°† and may render RCO a less welldefined reaction coordinate. Such effects will be investigated in the future. V. Conclusion The energetics for the CO2 + OH- f HCO3- reaction in the gas phase have been studied using density functional theory (B3LYP and PBE) and the MP2 method, and the basis set dependence of the reaction energy has been investigated. Extrapolating to the infinite basis limit, we find a 45.0 kcal/ mol MP2 dissociation energy for HCO3-. The difference in zero point energies between the reactions and products of eq 2 is 3.66 kcal/mol. All gas phase calculations, using different functionals or quantum chemistry methods, predict that no barrier exists for this reaction in the gas phase. We have performed ab initio molecular dynamics (AIMD) simulations of CO2, HCO3-, and CO32- in liquid water. As expected, CO2 is predominantly a hydrophobic species and only weakly coordinates to water oxygens and protons through its C and O atoms, respectively. HCO3- and CO32- exhibit a far greater tendency to coordinate with water molecules, forming on average 6.9 and 8.7 hydrogen bonds, respectively. Using umbrella sampling, we also computed the potential of mean force associated with the CO2 + OH- f HCO3- reaction in liquid water. After correcting for zero point energies, density functional theory exchange correlation functional errors (referenced to MP2 results), and rotational entropies, the standard state free energy change and barrier, ∆G° and ∆G°†, are found to be -9.8 and 9.7 kcal/mol respectively. The statistical uncertainty is about 1.4 kcal/mol. Our ∆G° value agrees very well with the experimental value of -9.4 kcal/mol, and the computed ∆G°† is in reasonable agreement with the experimental value of 12.1 kcal/mol. This work also highlights the postprocessing protocols currently needed to correct AIMD potential of mean force calculations to yield reasonable aqueous phase chemical reaction free energies. The pertinent corrections include those due to the DFT exchange correlation functionals, zero point energies, and insufficient sampling of orientational degrees of freedom. The effects of proton hopping may be significant for the nucleophilic attack of OH- on CO2 in water, but they have not been explicitly considered; this will be examined in the future.
Carbon Dioxide and Bicarbonate Hydration Acknowledgment. This work was supported by the Department of Energy, NIH grants DK63125, DK58563, and DK07789, the Max Factor Family Foundation, the Richard and Hinda Rosenthal Foundation, and the Fredricka Taubitz Fund. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL8500. Supporting Information Available: A movie depicting the CO2 + OH- f HCO3- reaction, starting from a transition state C-O distance of ∼2.1 Å and proceeding to conclusion without recrossing the barrier (legend: green, carbon; red, oxygen; white, hydrogen). This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Caldeira, K.; Wickett, M. E. Nature 2003, 425, 365. (2) Plane, J. M. C. Ann. Geophys. 2000, 18, 807. (3) Orr, J. C.; Fabry, V. J.; Aumont, O.; et al. Nature 2005, 437, 681. (4) Kurtz, I.; Petrasek, D.; Tatischev, S. J. Membr. Biol. 2004, 197, 77. (5) Pushkin, A.; Kurtz, I. Am. J. Physiol. 2006, 290, F580. (6) Pinsent, B. R. W.; Roughton, F. J. W. Trans. Faraday Soc. 1951, 47, 263. (7) Harned, H. S.; Davis, R., Jr. J. Am. Chem. Soc. 1943, 65, 2030. (8) Himmelblau, D. M.; Babb, A. L. AIChE J. 1958, 4, 143. (9) Pinsent, B. R. W.; Pearson, L.; Roughton, F. J. W. Trans. Faraday Soc. 1956, 52, 1512. (10) Danckwerts, P. V.; Melkersson, K. A. Trans. Faraday Soc. 1962, 58, 1832. (11) Pohorecki, R.; Moniuk, W. Chem. Eng. Sci. 1988, 43, 1677. (12) Falk, M.; Miller, A. G. Vib. Spectrosc. 1992, 4, 105. (13) Abashkin, Y.; Mele, F.; Russo, N.; Toscano, M. Int. J. Quantum Chem. 1994, 52, 1011. (14) Kuznetsova, T.; Kvamme, B. Phys. Chem. Chem. Phys. 2002, 4, 937. da Rocha, S. R. P.; Johnston, K. P.; Westacott, R. E.; Rossky, P. J. J. Phys. Chem. B 2001, 105, 12092. (15) In Het Panhuis, M.; Patterson, C. H.; Lynden-Bell, R. M. Mol. Phys. 1998, 94, 963. (16) Jo¨nsson, B.; Karlstro¨m, G.; Wannerstro¨m, H. J. Am. Chem. Soc. 1978, 100, 1658. (17) Davis, R. P. J. J. Am. Chem. Soc. 1958, 80, 1958. (18) Bra¨r, M.; Pe´z-Lustres, J. L.; Weston, J.; Anders, E. Inorg. Chem. 2002, 41, 1454. (19) Merz, K. M., Jr.; Banci, L. J. Am. Chem. Soc. 1997, 119, 863. (20) Schroo¨der, D.; Schwarz, H.; Schenk, S.; Anders, E. Angew. Chem. 2003, 42, 5087. (21) Liang, J. Y.; Lipscomb, W. N. J. Am. Chem. Soc. 1986, 108, 5051. (22) Davidson, M. M.; Hillier, I. H.; Hall, R. J.; Burton, N. A. Mol. Phys. 1994, 83, 327. (23) Nemukhin, A. V.; Topol, I. A.; Grigorenko, B. L.; Burt, S. K. J. Phys. Chem. B 2002, 106, 1734. (24) Peng, Z.; Merz, K. M., Jr. J. Am. Chem. Soc. 1992, 114, 2733. (25) Peng, Z.; Merz, K. M., Jr. J. Am. Chem. Soc. 1993, 115, 9640. (26) See, for example: Leung, K.; Rempe, S. B. J. Am. Soc. Chem. 2004, 126, 344. (27) Kieke, M. L.; Schoppelrei, J. W.; Brill, T. B. J. Phys. Chem. 1996, 100, 7455. Takemura, F.; Matsumoto, Y. Chem. Eng. Sci. 2000, 55, 3907. (28) See, for example: Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (29) Sprik, M. Chem. Phys. 2000, 259, 139. (30) Sprik, M. J. Phys.: Condens. Matter 2000, 12, A161. Davis, J. E.; Doltsinis, N. L.; Kirby, A. J.; Roussev, C. D.; Sprik, M. J. Am. Chem. Soc. 2002, 124, 6594. (31) Ivanov, I.; Klein, M. L. J. Am. Chem. Soc. 2002, 124, 13380. (32) Leung, K.; Rempe, S. B. J. Chem. Phys. 2005, 122, 184506.
J. Phys. Chem. B, Vol. 111, No. 17, 2007 4459 (33) Lee, J. G.; Asciutto, E.; Babin, V.; Sagui, C.; Darden, T.; Roland, C. J. Phys. Chem. B 2006, 110, 2325. (34) For example, see: Meijer, E. J.; Sprik, M. J. Phys. Chem. A 1998, 102, 2893. (35) Blumberger, J.; Klein, M. L. Chem. Phys. Lett. 2006, 422, 210. (36) Møller, C.; Plesset, M. S. Phys. ReV. 1934, 46, 618. (37) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; McGrawHill: New York, 1989. (38) Pratt, L. R.; Asthagiri, D.; Kress, J. D. Phys. ReV. E 2003, 68, 041505. (39) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 120, 300. (40) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Chem. Phys. 2005, 122, 014515. (41) Hammer, B.; Hansen, L. B.; Norskov, J. K. Phys. ReV. B 1999, 59, 7413. (42) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (43) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1998, 37, 785. (44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (45) Asthagiri, D.; Pratt, L. R.; Kress, J. D.; Gomez, M. A. Proc. Natl. Acad. Sci. 2004, 101, 7229. (46) Tuckerman, M. E.; Marx, D.; Parrinello, M. Nature 2002, 417, 925. (47) Tuckerman, M. E; Chandra, A.; Marx, D. Acc. Chem. Res. 2006, 39, 151. (48) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. Science 1997, 275, 817. (49) Chen, B.; Ivanov, I.; Park, J. M.; Parrinello, M.; Klein, M. L. J. Phys. Chem. B 2002, 106, 12006. (50) Kresse, G.; Furthmu¨ller, J. Phys. ReV. B 1996, 54, 11169; Comput. Mater. Sci. 1996, 6, 15. (51) Blochl, P. E. Phys. ReV. B 1994, 50, 17954. (52) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758. (53) See http://cms.mpi.univie.ac.at/vasp/. (54) Unconstrained proton hopping between OH- and H2O in liquid water occurs at O(1) ps intervals, neither slow enough to be avoided on AIMD time scales nor fast enough to enable adequate sampling throughout the simulation cell. (55) Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. Biophys. J. 1997, 72, 1047. (56) The azimuthal angle associated with the free O-O rotation about the CCO2-OOH- axis in CO2, which becomes a vibration (hindered rotation) in HCO3-, is not sampled due to the three-atom coplanar constraint. Note that this constraint implies neglecting an additional out-of-plane bending mode of both HCO3- and CO2. However, the bending frequencies in these two species are sufficiently similar that the free energy correction to this constraint, ∼kBT log(ω1/ω2), is negligible. (57) With our system and constraints, the end product (CO2 and OH-) rotational entropy would yield a factor of 2/(2 × 4π), which would have to be included in the prefactor of C0 if we have frozen these degrees of freedom. For OH- in the fully dissociated state, the OH- azimuthal rotation about the CCO2-OOH- axis is a free rotation. However, since this motion describes a very small solid angle because θOH is close to zero, it should be well sampled within AIMD. (58) Hammes, G. G. Principles of Chemical Kinetics; Academic Press: New York, 1978; pp 52. (59) Frisch, M. J.; et al. Gaussian 03, revision C.02; Gaussian Inc.: Wallingford, CT, 2004. (60) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1994, 100, 2975. (61) Jensen, F. J. Chem. Phys. 1999, 110, 6601. Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Olsen, J. Chem. Phys. Lett. 1999, 302, 437. (62) Halkier, A.; Helgaker, T. Jørgensen, P.; Klopper, W.; Koch, H.; Olsen, J.;Wilson, A. K. Chem. Phys. Lett. 1998, 286, 243. (63) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (64) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (65) We have neglected the difference in hydration free energies predicted by different quantum mechanical methods because these involve subtle balancing between water-solute and water-water interactions. (66) For this computation, we have used an early estimate of RCO ) 2.1 Å for the barrier top position (obtained before the full set of statistics were accumulated) instead of RCO ) 2.08 Å, but we expect the resulting differences to be negligible. (67) Chandler, D. J. Chem. Phys. 1978, 68, 2959.
