794
J. Phys. Chem. B 2009, 113, 794–802
Hydrogen-Bonding Structure and Dynamics of Aqueous Carbonate Species from Car-Parrinello Molecular Dynamics Simulations P. Padma Kumar,*,† Andrey G. Kalinichev,*,‡ and R. James Kirkpatrick§ Department of Geology, UniVersity of Illinois at Urbana-Champaign, Urbana, Illinois 61801 ReceiVed: October 14, 2008; ReVised Manuscript ReceiVed: NoVember 17, 2008
A comprehensive Car-Parrinello molecular dynamics (CP-MD) study of aqueous solutions of carbonic acid (H2CO3), bicarbonate (HCO3-), carbonate (CO32-), and carbon dioxide (CO2) provides new quantitative insight into the structural and dynamic aspects of the hydrogen-bonding environments for these important aqueous species and their effects on the structure, H-bonding, and dynamical behavior of the surrounding water molecules. The hydration structures of the different carbonate species depend on their ability to accept and donate H-bonds with H2O. The H-bonds donated by the C-O-H sites of the carbonate species to water molecules are generally stronger and longer-lived than those accepted by these sites from water molecules. The structural relaxation among the water molecules is dominated by diffusional (translational) motion of H2O, whereas the H-bond reorganization is dominated by the librational motion of the water molecules and the carbonate species. The rates of structural relaxation of the H2O molecules and the rates of H-bond reorganization among them are slower in systems containing carbonate species, consistent with previous studies of simple salt solutions. The strengths and lifetimes of H-bonds involving the carbonate species positively correlate with the total negative charge on the species. H-bond donation from H2O to CO2 is weak, but the presence of CO2 noticeably affects the structure and structural relaxation of the surrounding H-bonding network leading to generally stronger H-bonds and slower relaxation rates, the behavior typical of a hydrophobic solute. 1. Introduction There is strong evidence that the local hydration structure around a solute molecule in aqueous solution and the response of the water molecules to the changes in the geometric and electronic structure of the solute have significant influence on the reactivity of solutes in the solution.1,2 Advances in understanding the thermodynamics and kinetics of aqueous solutions, however, require a dynamical treatment that goes beyond the conventional statistical approach and that also uses a molecular description of the solute rather than a dielectric continuum approach.1,2 The experimental microscopic temporal and spatial resolution of solution structure and dynamics sufficient to support such theoretical development has not been available in the past. Recent advances in time-resolved spectroscopic techniques and computer simulations are, however, now delivering such a molecular level description. Recent spectroscopic investigations of pure water3 and alkali halide and Na/Mg perchlorate (ClO4-) solutions4-6 provide valuable insight into the hydrogen-bonding and relaxation dynamics of water molecules in the bulk and in the solvation shell around solute molecules. A detailed microscopic picture of the solvation structure and dynamics of a variety of atomic and molecular species in aqueous solutions is also emerging from recent molecular dynamics simulations.7-17 * To whom correspondence should be addressed. E-mail: padmakumarp@ iitg.ernet.in;
[email protected]. † Current address: Department of Physics, Indian Institute of Technology Guwahati, Guwahati, Assam 781039, India. ‡ Current address: Department of Chemistry, Michigan State University, East Lansing, MI 48824. § Current address: College of Natural Science, Michigan State University, East Lansing, MI 48824.
Aqueous solutions containing dissolved CO2 and the various carbonate species are ubiquitous in the natural environment and are major players in a broad range of processes of critical importance in chemistry, biology, geochemistry, atmospheric sciences, astronomy, and environment sciences.18-22 Understanding the reactivity of aqueous carbonate species is especially critical to understanding the role of CO2 in the global carbon cycle and climate change and in developing viable carbon sequestration technology. Despite the importance of aqueous carbonate species, surprisingly little is known about their structural environments and dynamical behavior in aqueous solution. It has been known for many decades that when CO2 dissolves in water, about 1% of it converts to carbonic acid, H2CO3.23,24 It is also well-known that aqueous carbonic acid deprotonates to bicarbonate, HCO3-, and carbonate, CO32-, with pKa values of 6.4 and 10.3, respectively. Direct spectroscopic data for H2CO3 date from the early 1990s,25 and there is more recent spectroscopic evidence for its presence on calcium carbonate surfaces.26,27 Carbonic acid can be synthesized by irradiation of frozen films of CO2/H2O mixtures with highenergy He ions or protons25,28 and by protonation of HCO3-/ CO32- solutions at low temperatures.29 Crystalline H2CO3 can be sublimated and recondensed without decomposition to CO2 and H2O.30 There have been several recent quantum chemical studies of the structure, energetics, and reactivity of carbonate species using self-consistent field Hartree-Fock methods (SCFHF), coupled cluster methods (CCM), Møller-Plesset (MP) perturbation theory, and density functional theory (DFT).31-43 Because these computational techniques are expensive, they have been limited largely to gas phase systems consisting of a molecule of carbonic acid or carbon dioxide either in the bare state or microsolvated by a few water molecules.31-39 Quantum
10.1021/jp809069g CCC: $40.75 2009 American Chemical Society Published on Web 12/24/2008
Structure and Dynamics of Aqueous Carbonate Species chemical calculations for carbonate species in solution have either treated the solvent implicitly by mean-field or dielectric continuum approximations or have used empirical force fields.40-42 The hydrogen-bonded structures and vibrational frequencies of isolated H2CO3 dimers and larger clusters have also been studied by quantum chemical techniques.44,45 The stability of the bare carbonic acid molecule was a controversial topic43 until the mass spectroscopic studies by Terlouw et al. confirmed its existence in 1987.46 Theoretical calculations of the dissociation of gas-phase H2CO3 to CO2 and H2O have estimated the reaction to be exothermic by about 9-10 kcal/mol, but with an activation barrier of about 43 kcal/ mol, thus explaining the observed stability of gas phase H2CO3.33,37-39 More recent theoretical studies by Liedl and co-workers37,38 have shown that the activation barrier for H2CO3 dissociation decreases as the number of solvating water molecules increases, suggesting that hydrogen-bonding with H2O molecules plays a significant catalytic role in this reaction. Their transition state theoretical (TST) approach based on a coupled cluster method estimated an activation barrier of 43 kcal/mol for the dissociation of a bare H2CO3 molecule with this value decreasing monotonically to ∼21 kcal/mol when three H2O solvate the molecule.37,38 There have been several classical molecular dynamics (MD) and Monte Carlo (MC) computer simulation studies of CO2-H2O solutions,47-52 but this work focused largely on the thermodynamic properties and developing better equations of state for CO2-H2O mixtures. The molecular-scale structure and hydrogen-bonding in such solutions have attracted much less interest. Classical molecular simulations of this type employ predefined intermolecular potentials and can efficiently probe the structural and dynamic properties of bulk aqueous solutions, but their applicability is limited when chemical reactions occur. The development of “reactive” intermolecular potentials capable of modeling the formation and breaking of chemicals bonds in such situations has proven to be difficult.53,54 In addition, recent studies have observed that classical MD methods tend to overestimate the size of the solvation sphere.7,8,10,14 Explicit incorporation of ionic polarizability in the intermolecular potentials may improve the results but does not solve the difficulty in describing the chemical reactivity. Given the prohibitively large computational cost of traditional high-level quantum chemical methods in modeling bulk solvent environments and the limited applicability of classical molecular simulation techniques, several compromise ways to incorporate chemical reactivity into the molecular modeling of aqueous systems have been developed. The important advantage of these approaches over traditional quantum chemical techniques is in the explicit incorporation of finite temperature effects and the feasibility of performing calculations on large enough systems that the effects of the solution environment can be more realistically modeled. The QM/MM technique separates the solvent and solute system into two parts: (i) the QM (for quantum mechanical) part is treated using high-level quantum chemical methods but includes only the reaction center and the solvent molecules in its immediate vicinity and (ii) the MM (for molecular mechanics) part, which includes the rest of the system and is treated using a classical force field approach. Rode and co-workers16 have recently successfully applied this method in the form of quantum mechanical charge field molecular dynamics (QMCF-MD)55 to study the structure and dynamics of a number of oxo anions in water. A limitation of this approach is that the “classical” solvent molecules can frequently exchange positions with their “quantum” counterparts near the solute
J. Phys. Chem. B, Vol. 113, No. 3, 2009 795 during a typical simulation time of a few tens of picoseconds. This exchange partially eliminates the capability of some solvent molecules to react with the solute. The Car-Parrinello molecular dynamics technique (CPMD),55 sometimes also called ab initio molecular dynamics (AIMD), appears to be a more universal molecular modeling approach to the study of chemical reactions in many aqueous systems. In contrast to the QM/MM approach, all solvent molecules are treated equivalently in the CP-MD simulation, and all are capable of participating in chemical reactions in any region of the simulated system. Admitably, the gain is partially through the utility of density functionals developed within the local density approximation (LDA). CP-MD is superior to classical MD methods because the intermolecular interactions are computed “on the fly” at the DFT level, thus incorporating many-body interactions and ion polarizability. Leung et al.57 have recently used this approach to investigate the nucleophilic attack of hydroxide on CO2 in the aqueous environment, and Rustad et al.57 have used a combination of CP-MD with quantum chemical calculations of large supermolecular clusters to study carbon isotope fractionation in CO2(g), aqueous carbonate species, and carbonate minerals. Here we present a detailed CP-MD study of the structure and dynamics of hydrogen-bonding around molecules of carbonic acid (H2CO3), bicarbonate (HCO3-), carbonate (CO32-), and CO2, in aqueous solutions. The objective of these simulations is to provide new fundamental insight into the hydration structures, H-bond environments and lifetimes, and dynamical behavior of these important species, which are a required basis for future study of their molecular scale reactivity. In a recent study39 employing a metadynamics algorithm59 based on the Car-Parrinello technique,56,60 we reported the energetics of the conformational rearrangements of gas-phase carbonic acid and its dissociation to H2O and CO2. The results are in good agreement with previous theoretical studies, provide a detailed dissociation mechanism for gas phase H2CO3, and demonstrate the usefulness of the technique. These results, together with the present study, lay the foundation for ongoing efforts to understand the chemical reactivity of carbonate species in aqueous environment on a fundamental atomistic time and length scale. 2. Methods The CP-MD calculations here employ the Car-Parrinello molecular dynamics algorithm as implemented in version 3.11.1 of the CPMD software.61 The simulated systems are isoelectronic, with a single molecule of H2CO3, HCO3-, or CO32embedded in a box with 45 water molecules or a single molecule of CO2 embedded in 46 water molecules. Pure water is modeled using 48 water molecules. In all cases, periodic boundary conditions are used with the cubic cell edge of 11.3 Å. For pure water these parameters yield the experimental density of 1 g/cm3 (0.033 molecules/Å3). The starting configurations for the CPMD simulations were well equilibrated structures from 100 ps classical MD simulations using SPC62 and CVFF63 force field parameters. A further equilibration was carried out for ∼15 ps using CP-MD before 14 ps long equilibrium CP-MD trajectories were produced for detailed analysis of structure and dynamics. Note that with no counterions present the HCO3- and CO32systems are not electrostatically neutral, and the structures and dynamics obtained for them do not account for the effects of dissolved cations. For the sake of convenience, however, we use the terms solute and solution to refer to the carbonate species and aqueous systems containing all the carbonate species. These
796 J. Phys. Chem. B, Vol. 113, No. 3, 2009
Kumar et al.
two non-neutral systems are handled by adding a neutralizing charge to the simulation cell together with Ewald summation.61 The added charge is uniformly distributed in the simulation cell such that it has no spatial gradient and imparts no force on the ions and, hence, does not affect the dynamics of the ions. The valence electrons are treated explicitly within the DFT formalism employing the gradient-corrected Becke, Lee, Yang, and Parr (BLYP) functional.64,65 The interactions of the valence electrons with the nuclei and core electrons are represented by ultrasoft pseudopotentials originally proposed by Vanderbilt.66 The use of ultrasoft pseudopotentials (USPP) greatly reduces the number of plane waves necessary in DFT calculations relative to norm-conserving pseudopotentials.67 Therefore, in our simulations, the Kohn-Sham orbitals are expanded using plane wave basis sets,with a cutoff of 40 Ry. An electronic mass of 600 au and a time step of 4 au (≈0.1 fs) are used. A temperature of 310 K for the molecules and a fictitious kinetic energy of about 0.02 au for the electronic degrees of freedom are controlled using Nose´-Hoover thermostats. The physically uninteresting bulk translations of the system are removed every 100 MD steps. The H2CO3 molecule occurs in three different conformations, trans-trans (TT), cis-trans (CT), and cis-cis (CC), and we have undertaken calculations for the TT and CT conformers. The three conformers differ in the orientation of the hydroxyl groups with respect to the carbonyl group.33,39 Previous theoretical calculations suggest that in the gas phase TT is the most stable species but is only about 1-2 kcal/mol lower in energy than CT.33,39 This energy difference is small enough that their relative stability could be easily changed by H-bonding with water molecules, and it is not known which of these two species is the most stable in aqueous solution. The CC conformer is about 9 kcal/mol higher in energy than CT in the gas phase and is expected to be much less stable than the other two even in solution. This study focuses on the structural and dynamic aspects of hydrogen-bonding, and we use a uniform definition for the existence of an H-bond between two species that requires simultaneously satisfying three criteria: (i) the donor-acceptor O · · · O distance is less than 3.5 Å, (ii) the donor-acceptor H · · · O distance is less than 2.45 Å, and (iii) the hydrogendonor-acceptor angle is less than 30°.11 We examine the structural relaxation of hydrogen bonds in our simulated systems using two suitably defined time correlation functions. The function
CHB(t) ) 〈h(0) · h(t)〉/〈h〉
(1)
describes the structural relaxation of hydrogen bonds. Here h(t) is unity if an H-bond as defined above exists between a pair of donor-acceptor sites at time t and zero otherwise. A similarly defined function
CdHB(t) ) 〈h(0) · hd(t)〉/〈h〉
(2)
is based only on the distance criterion (i) and thus describes the structural relaxation of the H-bonds due to translational diffusion. Chandra11 has used these functions in a study of relaxation in NaCl and KCl solutions and found evidence for slower structural relaxation of H-bonds among water molecules in the presence of ions. We estimate the lifetimes of the H-bonds between specific donor-acceptor pairs using two similar time correlation functions:
SHB(t) ) 〈h(0) · H(t)〉/〈h〉
(3)
SdHB(t) ) 〈h(0) · Hd(t)〉/〈h〉
(4)
and
where H(t) is unity if a tagged pair of donor-acceptor sites remain continuously H-bonded for a period t and zero otherwise. SHB(t) probes the sum of all rotational and translational motions affecting H-bonding, whereas SdHB(t) probes only the translational contribution to the H-bond lifetimes, since it, like CdHB(t), is defined based only on the distance criterion (i) above. The power spectra of the carbonate species were calculated as Fourier transforms of their velocity autocorrelation functions over the 14 ps production trajectories of the CP-MD simulations. 3. Results and Discussion 3.1. Hydration Structure and Hydrogen-Bonding. The CPMD results show significant differences in the local hydration structure and hydrogen-bonding of the different carbonate species, as illustrated by the computed 3-dimensional spatial atomic density distributions for the hydrogens (gray) and oxygens (red) of water molecules around them (Figure 1a-e). In these diagrams, the first shells of hydrogen and oxygen include those directly H-bonded to the solute assuming a simplified H-bonding criterion of 2.4 Å between the donor and acceptor sites. The second shell includes the oxygens and hydrogens bonded to the atoms in the first shell. The enclosed volumes represent number densities of 0.109/Å3 or more around H2O (Figure 1f) and all the carbonate species except CO2. This threshold value corresponds to 3.5 times the average number density of water. For CO2 (Figure 1e) the number density visualization threshold is 0.0109/Å3 because this species makes too few H-bonds with H2O to be effectively visualized at the larger density, and there is only weak H-bond acceptance by the oxygens of CO2. The atomic density distributions of the hydrogen and oxygen of water in the hydration shell around an arbitrarily chosen H2O molecule from the pure water simulation (Figure 1f) are directly comparable with the CP-MD (BLYP/ pw) results of Mantz et al.68 for 64 water molecules, except that we use a 5% lower threshold density. The computed hydration structures involve three kinds of H-bonds: (i) OH-HO · · · Ow, donated by the hydroxyl groups (-C-OH) of the carbonate species to water, (ii) Ow-Hw · · · OH, donated by water to the hydroxyl oxygens of the carbonate species, and (iii) Ow-Hw · · · )O, donated by water to the carbonyl groups (-CdO). Symbols )O, OH, HO refer respectively to carbonyl oxygens, hydroxyl oxygens, and hydroxyl hydrogens of the carbonates. Ow and Hw refer to the oxygens and hydrogens of water molecules, respectively. The structure of the solute also affects the hydration structure, as demonstrated by the differences in the atomic density distributions for the TT and CT conformers of H2CO3 (Figure 1a,b). Thus, the hydration structures appear to result from the complex interplay of slightly different energetics and dynamics of the accepted and donated H-bonds between the different sites of the solute molecules and water molecules and among water molecules themselves. The atom-atom radial distribution functions (RDFs; g(r)) between different solute sites and H2O molecules provide a more detailed quantitative description of the hydration structures (Figure 2). For all the carbonate species except CO2, the C-Ow RDFs all contain a nearest-neighbor peak near 3.5 Å that
Structure and Dynamics of Aqueous Carbonate Species
J. Phys. Chem. B, Vol. 113, No. 3, 2009 797
Figure 1. 3-D atomic density distributions of hydrogens (gray) and oxygens (red) of water molecules H-bonded to carbonate species (or to an arbitrarily chosen water molecule in the case of pure H2O simulation) obtained from CP-MD simulations. The first shell of hydrogens (gray) or oxygens (red) of water molecules included in the atomic density calculation are within 2.4 Å of the oxygens or hydrogens of the solute molecule, which is a simplified H-bonding criterion. The second shell of density includes the water hydrogens or oxygens bonded to those in the first shell. Clockwise from top-left corner: (a) trans-trans (TT) H2CO3, (b) cis-trans (CT) H2CO3, (c) bicarbonate, HCO3-, (d) carbonate, CO32-, (e) CO2, (f) pure H2O. A uniform threshold number density of 0.109/Å3 is chosen for the isosurface in all plots, except for the case of CO2 where the threshold is set to 0.0109/Å3 (see text for details).
Figure 2. Atom-atom radial distribution functions between selected sites of the carbonate species and water molecules at 310 K. The symbols C, )O, OH, and HO refer respectively to carbon, carbonyl oxygen, hydroxyl oxygen, and hydroxyl hydrogen of the carbonate species. Ow and Hw refer to the oxygens and hydrogens of water molecules, respectively. Legends in (a) are common to all.
becomes progressively sharper as the charge of the carbonate species increases (Figure 2a). This result indicates a progressively stronger and more ordered solvation structure with increasing net solute charge. The C-Ow RDFs of TT- and CTH2CO3 conformers are similar. The C-Ow RDF of CO2 is quite
different, with a low-intensity feature (shoulder) at ∼2.75 Å and a broad peak between 3 and 5 Å. These structural results are in generally good agreement with calculations of Leung at al.57 and Rustad et al.58 The RDFs for the carbonyl oxygens and hydrogens of water ()O-Hw) show a first peak at ∼1.8 Å that shifts to shorter distances and becomes more intense with increasing negative charge of the species (Figure 2b). This trend clearly indicates a progressively increasing H-bond strength for this type of site with increasing net charge on the ion. Compaction of the 3-D hydration structure around the )O oxygens (Figure 1a-d) was also recently observed around the oxygens of perchlorite, sulfate, and phosphate in QMCF-MD simulations.16 The differences in the RDFs of the TT and CT conformers of H2CO3 are related to their structures and the observation that there are two groups of water molecules in the hydration shell: those donating H-bonds to the carbonate species and those accepting H-bonds from the carbonate species. Water molecules in these two groups do not make favorable H-bonds with water molecules in the other group. For the TT conformer, the water molecules donating H-bonds to )O are located between those accepting H-bonds from the -OH groups. For TT and CT conformers the H-bonds donated by the -OH groups to water molecules (HO · · · Ow) are stronger and presumably more “rigid” than those donated by H2O to the )O, as evidenced by their much shorter HO-Ow distances (Figure 2d). This configuration poses a significant structural constraint on the “basin” of water molecules around )O, resulting in fewer H-bonds accepted by the site. This effect is less pronounced for the CT conformer because the water molecules donating H-bonds to )O are located between those forming HO · · · Ow H-bonds (strong) on one side and those forming the relatively weaker OH · · · Hw bonds on the other side. The smaller steric hindrance could be partially responsible for the higher )O-Hw RDFs peaks for HCO3- and CO32-. As also seen in the C-Ow RDF discussed above, the
798 J. Phys. Chem. B, Vol. 113, No. 3, 2009
Figure 3. Oxygen-oxygen radial distribution functions for water, pure as well as containing different carbonate species from the CP-MD simulations at 310 K. Bins of size 0.02 Å are used in all cases. )
O-Hw RDF for CO2 is markedly different from the others. There is little intensity at distances less than 2 Å, suggesting very weak H-bonding between CO2 and H2O molecules. This conclusion is in qualitative agreement with the low solubility of CO2 in water and has been observed in previous MD simulations.57,58 The OH-Hw RDFs of HCO3- and the TT and CT conformers of H2CO3 (Figure 2c) have much smaller first peaks at a larger distance than the HO · · · Ow RDFs (Figure 2d). This result clearly shows that the H-bonds accepted by the -OH groups (OH · · · Hw) are much weaker than those donated by them. The donated HO · · · Ow H-bonds are shorter (stronger) for both conformers of carbonic acid than for HCO3-, which is indicative of the relative ease of deprotonation of carbonic acid, as shown by its lower pKa value. The RDFs between oxygens of water molecules (Ow), gOwOw(r), for the carbonate-containing solutions (Figure 3) are not greatly different from that of pure water. The Ow-Ow RDFs for the solutions, except that of CT-H2CO3, have slightly higher maxima and deeper minima, suggesting a somewhat more ordered solvent structure compared to pure water. It is not clear at this point whether the small structural differences of the different systems is due entirely to the differences in the first hydration shell around the solute molecules or to longer range effects. For pure H2O, the position (2.75 Å) and height (2.98) of the first peak of gOwOw(r) are in good agreement with the most recently published CP-MD simulations of 64 H2O molecules employing the BLYP/pw functional and at a comparable temperature.68,69 The values for the first minimum and second maximum for pure water are ∼3.35 Å, 0.6 and ∼4.45 Å, 1.35, respectively, and are also in good visual agreement with previous reports.68,69 The relatively flat first minimum and second maximum make it difficult to locate their positions accurately. The time-averaged H-bonding statistics of the various sites of the aqueous carbonate species computed using the three-point criterion discussed in section 2 show that the number of H-bonds accepted and donated by the different sites varies significantly from species to species (Figure 4). For the TT and CT conformers of H2CO3 and for HCO3-, the number of H-bonds donated by each -OH group to water (HO · · · Ow) is essentially unity, reflecting the large relative strength and stability of these H-bonds described above. The number of H-bonds accepted by the )O sites increases from CT-H2CO3 to HCO3-, to CO32-, reflecting the increasing negative charge on the )O sites. The
Kumar et al.
Figure 4. Statistics of the H-bonds between different sites of carbonate species with water molecules together with those between water molecules at 310 K averaged over the entire 14 ps long CP-MD trajectories. The numbers are per site of the carbonate species for the carbonate-water H-bonds. The number of H2O · · · H2O H-bonds is per water molecule.
OH sites of TT-H2CO3 accept a greater number of H-bonds than those of CT-H2CO3 and HCO3-. The number of H-bonds accepted by the two symmetrically inequivalent -OH groups of the CT conformer are also different: 0.27 per OH site for the one oriented toward the )O and 0.49 per OH for the one oriented away from it. Similarly, the )O of the bicarbonate closer to the HO accepts fewer H-bonds from water (2.19 per )O site) than the other )O (2.52 per site). These differences can be rationalized in terms of steric constraints arising from (i) the inability of the water molecules in the first hydration shell to form compact structures favorable for H-bonding among themselves and (ii) the relative strength of the H-bonds donated by the carbonate species, as discussed above. The presence of the carbonate species decreases the number of H-bonds between H2O molecules relative to pure water, except for CO2 which increases this number. Thus, CO2 behaves in this respect as a typical hydrophobic solute. 3.2. Structural Relaxation and H-Bonding Lifetimes. The hydrogen-bonding time correlation functions, CHB(t) and CdHB(t) (eqs 1 and 2), describe respectively the overall structural relaxation and the structural relaxation due only to translational diffusion. The decay of both functions is faster than exponential during the first 1-2 ps but approaches exponential at longer times (Figure 5). This behavior is in agreement with the results of previous classical MD simulations of simple salt solutions using the SPC/E water potential.11 For pure water, the CHB(t) and CdHB(t) relaxation times, τc and τdc , obtained from linear fits to the exponential portion of our results are 5.0 and 5.1 ps, respectively (Table 1). Given the technical differences between the methods, these values agree favorably with the values of 6.6 and 7.6 ps from the previous classical MD studies.11 For the carbonate solutions, the structural relaxation times (Table 1) are all larger than for pure water and increase from CT-H2CO3 to TT-H2CO3 to HCO3- to CO32-. The CT-H2CO3 solution exhibits a longer nonexponential regime than the others and relaxes significantly faster than that of the TT conformer. The behavior of the CT conformer is unexpected but consistent with the less ordered structure of its hydration shell (Figures 1 and 2). The logarithmic nature of the relaxation plots in Figure 5 suggests that longer simulation times are needed for precise estimates of the structural relaxation times for CT-H2CO3. The
Structure and Dynamics of Aqueous Carbonate Species
J. Phys. Chem. B, Vol. 113, No. 3, 2009 799
Figure 6. Time correlation functions, SHB(t) and SdHB(t) (see text for details), describing the lifetimes of H-bonds between water molecules for different solutions of carbonate species at 310 K.
Figure 5. Time correlation functions, CHB(t) and CdHB(t) (see text for details), showing the structural relaxation of H-bonds between water molecules for different solutions of carbonate species at 310 K. Legends in the lower panel are common to both.
longer relaxation times of the HCO3- and CO32- solutions correlate well with their increasing net charge and increasing H-bonding strength. This observation is again consistent with the results of classical MD simulations for simple salt solutions, in which increasing salt concentration causes slower relaxation.11 In both cases, the increasing ionic strength (or the charge/volume ratio) leads to strengthening of the ion atmosphere felt by the water molecules and slows the H-bonding relaxation. The CO2 system also exhibits increased structural relaxation times despite the relatively weak CO2 · · · H2O H-bonding. This behavior is consistent with the increased structural ordering of water molecules in the CO2 hydration shell (Figure 2), typical for a hydrophobic solute. The two relaxation times obtained from these results, τc and τdc (Table 1), are equal for all the systems studied here within typical statistical errors. This result demonstrates that the longer time structural relaxation (t > 2 ps) is due principally to translational diffusion of H2O molecules and not to their rotational or librational motion. The SHB(t) and SdHB(t) correlation functions (eqs 3 and 4) that probe H-bond lifetimes also exhibit nonexponential behavior at t < 0.5 ps and exponential decay at longer times (Figure 6, Table 2). These functions provide detailed insight into the H-bond dynamics. For all of our systems, the SdHB(t) lifetimes, τds , are significantly longer than the SHB(t) ones, τs, in contrast to the structural relaxation for which CdHB(t) decays at the same rate as CHB(t). This result demonstrates that librational modes play a more important role than translational diffusion modes in breaking of water H-bonds (H2O · · · H2O) in both pure water70 and the carbonate solutions. As for the structural relaxation, the τs and τds H-bond lifetimes between H2O molecules are generally longer in the solutions than in pure water (Table 2). Again, the only exception is the CT-H2CO3 system, which relaxes faster than pure water. We do not have a clear understanding of this result. A second CPMD run, starting from the final structure of the first run for CT-H2CO3, revealed that this feature is persistsent. This behavior
TABLE 1: Structural Relaxation Times τc and τcd of H-Bonds between Water Molecules Computed from the Correlation Functions CHB(t) and CdHB(t), Respectively system
τc (ps)
τdc (ps)
pure H2O TT-H2CO3 CT-H2CO3 HCO3CO32CO2
4.9 8.8 6.0 13.8 40.4 10.7
5.1 9.2 6.4 13.4 40.0 10.0
TABLE 2: Lifetimes of Different H-Bonds, τs, Computed from the Correlation Function SHB(t)a system
H2O · · · H2O (ps)
pure H2O TT-H2CO3 CT-H2CO3 HCO3CO32CO2
0.88 (2.57) 1.10 (3.32) 0.68 (2.28) 1.10 (3.61) 1.43 (6.82) 1.14 (4.04)
)
O · · · Hw (ps)
OH · · · Hw (ps)
0.31 0.51 0.99 2.05 ∼0.01
0.22 0.15 0.13
a For the H-bonds between water molecules, the relaxation times, τds , computed from SdHB(t) are given in parentheses. The SdHB(t) relaxation times, τds , of the other H-bonds and the τs values for HO · · · Ow H-bonds are not reported because of statistical uncertainty.
may be due to the system not being completely at equilibrium, but it may also be due to the natural behavior of the reactant (CT-H2CO3) being close to a reaction event (H2CO3 T HCO3+ H3O+). This reaction could have a low activation barrier but with the product unstable due to pH considerations. CO2 increases the H2O · · · H2O H-bond lifetimes despite the very weak CO2 · · · H2O H-bond interaction, again consistent with the structural relaxation described above and with the hydrophobic nature of CO2 in solution. For pure water, our τs and τds values from CP-MD (0.88 and 2.57 ps, respectively) are in reasonable agreement with those from the previous classical MD simulations by Chandra (0.54 and 1.75 ps, respectively).11 The lifetimes of H-bonds accepted from water molecules by the OH and )O sites of the carbonate species are shorter than those between water molecules, except for )O of CO32- (Table 2), and the lifetimes of OH · · · Hw H-bonds are shorter than those for )O · · · Hw bonds for all relevant species. This latter behavior is consistent with the relatively long OH · · · Hw bond lengths
800 J. Phys. Chem. B, Vol. 113, No. 3, 2009
Figure 7. Power spectra of the atomic motions for different carbonate species in aqueous environment at 310 K. Trajectories were stored every 0.5 fs, and the spectra are computed from the full 14 ps production run. Symbols: C, )O, OH, and HO refer respectively to carbon, carbonyl oxygens, hydroxyl oxygens, and hydroxyl hydrogens. The open and filled triangles are respectively the harmonic frequencies of the gasphase TT and CT conformers of H2CO3. The legends are common to all the panels.
(compared with the )O · · · Hw bonds) in the corresponding RDFs (Figure 2b,c). The OH · · · Hw H-bonds are, thus, clearly weaker than )O · · · Hw H-bonds for all the species studied. The increasing lifetimes of )O · · · Hw bonds from TT to CT to HCO3- to CO32- correlate positively with shorter )O · · · Hw bond distances (Figure 2b) and thus correspond to progressively stronger ) O · · · Hw H-bonds. These trends also parallel the increase in the total number of H-bonds accepted by )O (Figure 4). For all of the systems studied, the H-bonds donated by the C-OH groups to H2O (HO · · · Ow) are much longer lived than those accepted by the C-OH (OH · · · Hw) and )O sites ()O · · · Hw) from water molecules. In fact, very few events of HO · · · Ow H-bond breakage were observed during the 14 ps long simulations. This behavior is consistent with the relatively short HO · · · Ow bond lengths (Figure 2) and demonstrates that (HO · · · Ow) H-bonds are quite strong and may have average lifetimes of several picoseconds. Because the HO · · · Ow H-bonds in our simulations are few in number and have lifetimes comparable to the duration of the simulation runs, these lifetimes have large statistical uncertainties and are not reported in Table 2. In contrast, the H-bonds accepted by the C-OH groups of the carbonic acid and bicarbonate species (OH · · · Hw) are, in general, the weakest and the shortest lived in these systems. 3.3. Vibrational Density of States. The power spectra of atomic motions for the aqueous carbonate are calculated as Fourier transforms of the velocity autocorrelation functions for the respective atoms from the 14 ps production trajectories of the CP-MD simulations. They provide additional insight into the relationships among the structure, H-bonding, and spectroscopic characteristics of these species (Figure 7). However, a detailed comparison of the full vibrational spectra for all the species and discussion of the origins of the differences between computed and observed71,72 frequencies is beyond the scope of this paper. The large red shift of the O-H stretching frequencies of the TT and CT conformers of carbonic acid in solution relative to the computed gas-phase values39 of about 3660 cm-1 (Figure 7a) is a clear signature of the strong HO · · · Ow H-bonding that is evidenced in the RDFs and H-bond lifetimes discussed above.
Kumar et al. These O-H stretching bands are very broad, spreading over a frequency range from 2100 to 3100 cm-1. This result is in qualitative agreement with the pioneering IR study by Moore and Khanna,25 who reported vibrational frequencies of 2614, 2840, and 3030 cm-1 for the O-H stretching modes of H2CO3/ H2O(ice) at 250 K. Our results are also consistent with more recent FTIR spectroscopic data73-75 for solid amorphous H2CO3 below 230 K which show a broad O-H stretching band between 2600 and 3200 cm-1. Recent CBSB7 B3LYP quantum chemical calculations45 of the vibration frequencies for the hydrogenbonded H2CO3 dimers in implicit PCM continuum solventmimicking water also give 2728, 2820, and 3232 cm-1 for the IR-active stretching vibrations of these species. Over the 14 ps duration of the TT and CT H2CO3 solution simulations, there were a few instances of O-H breaking leading to the formation of contact ion pairs of HCO3- and H3O+. These transient deprotonation events contribute to the large computed red shifts and frequency ranges. In all such events the same proton recombined with the parent carbonic acid in a short span of time. At near-neutral pHs the HCO3concentration is much higher than that of H2CO3 and is thermodynamically favored. At the high concentrations in our simulations, however, formation of HCO3- from H2CO3 is unlikely to be favorable. The presence of a single H3O+ molecule corresponds to a very low macroscopic pH, which greatly inhibits the forward reaction. Microscopically, the back reaction may also be more common than the forward reaction because the latter requires a reorganization of the hydration structure while the former only requires a proton hop along a single H-bond. The O-H stretching modes of HCO3- occur near 3200 cm-1 (Figure 7a) and are less red-shifted and have a narrower frequency range than the H2CO3 conformers. These differences are consistent with the longer HO · · · Ow distances and hence with weaker HO · · · Ow H-bonds for HCO3- (Figure 2d). Our value of ∼3200 cm-1 is somewhat lower than that from recent quantum chemical calculations59 using the self-consistent reaction field method with the PCM continuum solvation model (3423.5 cm-1). However, both are significantly larger that the experimentally observed value for the solvated bicarbonate ion, 2620 cm-1.71,72 Accurate calculation of O-H vibrational frequencies for the hydrated bicarbonate ion remains a significant challenge.72 The CdO stretching of the various carbonate species have been the signature modes for identification of carbonic acid in several experimental IR studies of samples containing carbonic acid together with other carbonates species.25-27 Our calculations show that the CdO stretch is at about 1500 cm-1 for HCO3and at about 1600 cm-1 for CT- and TT-H2CO3 (Figure 7c,d). For CO32- our calculated CdO stretching band is at about 1300 cm-1 (blue line in Figure 7d), well separated from those of the carbonic acid species. Thus, the CdO stretching feature employed by Moore and Khanna25 (in their case, at 1705 cm-1) in identifying the presence of carbonic acid in an ensemble of other carbonates embedded in an ice matrix is clearly visible in our simulated solution environment. However, intensity at this frequency is also observed in the HO power spectrum (Figure 7a), suggesting that the mode at 1650 cm-1 is not a pure CdO stretch, but also involves motion of hydrogens of the -OH group. The )O and C power spectra of CO2 (Figure 7c,d) have distinct peaks near 2200 cm-1 for the asymmetric stretching of the CdO bond. The computed Hw vibrational density of states are nearly identical for all simulated systems, suggesting that the presence
Structure and Dynamics of Aqueous Carbonate Species of dissolved carbonate species largely affects the relatively low frequency translational and librational processes of H2O, for which 14 ps CP-MD simulations are not very sensitive. 4. Conclusions CP-MD calculations provide new quantitative insight into the structure and dynamics of hydrogen-bonding environments for hydrated carbonic acid (H2CO3), bicarbonate (HCO3-), carbonate (CO32-), and carbon dioxide (CO2) and the effects of these species on the structure, H-bonding, and dynamical behavior of the surrounding water molecules. The hydration structures of the different carbonate species depend in detail on their ability to accept and donate H-bonds with H2O. One important control on the hydration structure and dynamics is that H-bonds between water molecules that donate H-bonds to the carbonate species and those that accept H-bonds from the carbonate species are often difficult to form due to steric constraints in the first hydration shell around the carbonate species. H-bonds donated by the -OH groups of the carbonate species to water molecules are stronger and longer lived than those accepted by these sites from water molecules. The rates of structural relaxation of the water molecules and the rates of H-bond reorganization among them are slower in systems containing carbonate species, consistent with previous classical MD simulations of simple aqueous salt solutions.11 On time scales longer than 2 ps the rate of structural relaxation among the water molecules is dominated by diffusional motion of H2O, whereas the rate of H-bond reorganization of the water molecules and the carbonate species is dominated by librational motion. The strengths and lifetimes of H-bonds involving the carbonate species generally increase with increasing net charge in the system. Solvated CO2 experiences only very weak H-bond donation from surrounding water molecules, but its presence in solution significantly affects the structure and structural relaxation of water, resulting in stronger H-bonds among the water molecules and slower relaxation rates. In this sense, dissolved CO2 exhibits properties characteristic of a typical hydrophobic solute. Acknowledgment. This research was supported by the DOE BES Geoscience Program (Grants DE-FG02-00ER-15028 and DE-FG02-08ER-15929), and the computational resources were provided by the DOE National Energy Research Scientific Computing Center, National Center for Supercomputing Applications, and other NSF TeraGrid supercomputing facilities. P.P.K. thankfully acknowledges Axel Kohlmeyer and Nisanth N. Nair for useful technical discussions. We acknowledge the use of VMD software for visualizations of MD trajectories. References and Notes (1) Rossky, P. J.; Simon, J. D. Nature (London) 1994, 370, 263. (2) Hynes, J. T. Nature (London) 1994, 369, 439. (3) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698. Eaves, J. D.; Loparo, J. J.; Fecko, C. J.; Roberts, S. T.; Tokmakoff, A.; Geissler, P. L. Proc. Natl Acad. Sci. U.S.A. 2005, 102, 13019. (4) Kropman, M. F.; Bakker, H. J. Science 2001, 291, 2128. (5) Omta, A. W.; Kropman, M. F.; Woutersen, S.; Bakker, H. J. Science 2003, 301, 347. (6) Nickolov, Z. S.; Miller, J. D. J. Colloid Interface Sci. 2005, 287, 572. (7) Bruge´, F.; Bernasconi, M.; Parrinello, M. J. Am. Chem. Soc. 1999, 121, 10883. (8) Yarne, D. A.; Tuckerman, M. E.; Klein, M. L. Chem. Phys. 2000, 258, 163. (9) Rey, R.; Møller, K. B.; Hynes, J. T. J. Phys. Chem. A 2002, 106, 11993. (10) Raugei, S.; Klein, M. L. J. Chem. Phys. 2002, 116, 196. (11) Chandra, A. Phys. ReV. Lett. 2000, 85, 768.
J. Phys. Chem. B, Vol. 113, No. 3, 2009 801 (12) Silvestrelli, P. L.; Bernasconi, M.; Parrinello, M. Chem. Phys. Lett. 1997, 277, 478. (13) Sillanpa¨a¨, A. J.; Simon, C.; Klein, M. L.; Laasonen, K. J. Phys. Chem. B 2002, 106, 11315. (14) Leung, L.; Rempe, S. B. J. Am. Chem. Soc. 2004, 126, 344. (15) Ebner, C.; Onthong, U.; Probst, M. J. Mol. Liq. 2005, 118, 15. (16) Pribil, A. B.; Hofer, T. S.; Vchirawongkwin, V.; Randolf, B. R.; Rode, B. M. Chem. Phys. 2008, 346, 182. (17) Hammerich, A. D.; Buch, V.; Mohamed, F. Chem. Phys. Lett. 2008, 460, 423. (18) Haugan, P. M.; Dranke, H. Nature (London) 1992, 357, 318. (19) Siegenthaler, U.; Sarmiento, J. L. Nature (London) 1993, 365, 119. (20) Saito, T.; Kajishima, T.; Nagaosa, R. EnViron. Sci. Technol. 2000, 34, 4140. (21) Ridgwell, A.; Zeebe, R. E. Earth Planet. Sci. Lett. 2005, 234, 299. (22) Lackner, K. S. Annu. ReV. Energy EnViron. 2002, 27, 193. (23) Buytendyk, F. J. J.; Brinkman, R.; Mook, H. W. Biochem. J. 1927, 21, 576. (24) Soli, A. L.; Byrne, R. H. Mar. Chem. 2002, 78, 65. (25) Moore, M. H.; Khanna, R. K. Spectrochim. Acta 1991, 47A, 255. (26) Al-Hosney, H. A.; Grassian, V. H. J. Am. Chem. Soc. 2004, 126, 8068. (27) Baltrusaitis, J.; Schuttlefield, J. D.; Zeitler, E.; Jensen, J. H.; Grassian, V. H. J. Phys. Chem. C 2007, 111, 14870. (28) Brucato, J. R.; Palumbo, M. E.; Strazzula, G. Icarus 1997, 125, 135. (29) Hage, W.; Hallbrucker, A.; Mayer, E. J. Am. Chem. Soc. 1993, 115, 8427. (30) Hage, W.; Liedl, K. R.; Hallbrucker, A.; Mayer, E. Science 1998, 279, 1332. (31) George, P.; Bock, C. W.; Trachtman, M. J. Comput. Chem. 1982, 3, 283. (32) Nguyen, M. T.; Ha, T.-K. J. Am. Chem. Soc. 1984, 106, 599. (33) Wight, C. A.; Boldyrev, A. I. J. Chem. Phys. 1995, 99, 12125. (34) Sun, H.; Mundy, S. J.; Maple, J. R.; Hagler, A. T. J. Phys. Chem. 1995, 99, 5873. (35) Sadlej, J.; Makarewicz, J.; Chałasin´ski, G. J. Chem. Phys. 1998, 109, 3919. (36) Jena, N. R.; Mishra, P. C. Theor. Chem. Acc. 2005, 114, 189. (37) Loerting, T.; Tautermann, C.; Kroemer, R. T.; Kohl, I.; Hallbrucker, A.; Mayer, E.; Liedl, K. R. Angew. Chem., Int. Ed. 2000, 39, 892. (38) Tautermann, C.; Voegele, A. F.; Loerting, T.; Kohl, I.; Hallbrucker, A.; Mayer, E.; Liedl, K. R. Chem.sEur. J. 2002, 8, 66. (39) Kumar, P. P.; Kalinichev, A. G.; Kirkpatrick, R. J. J. Chem. Phys. 2007, 126, 204315. (40) Sato, H.; Matubayasi, N.; Nakahara, M.; Hirata, F. Chem. Phys. Lett. 2000, 323, 257. (41) Tossell, J. A. Geochim. Cosmochim. Acta 2005, 69, 5647. (42) Nguyen, M. T.; Raspoet, G.; Vanquickenborne, L. G.; Duijnen, P. Th. V. J. Phys. Chem. A 1997, 101, 7379. (43) Ludwig, R.; Kornath, A. Angew. Chem., Int. Ed. 2000, 39, 1421. (44) Ballone, P.; Montanari, B.; Jones, R. O. J. Chem. Phys. 2000, 112, 6571. (45) Tossell, J. A. Inorg. Chem. 2006, 45, 5961. (46) Terlouw, J. K.; Lebrilla, C. B.; Schwarz, H. Angew. Chem., Int. Ed. 1987, 26, 354. (47) Duan, Z. H.; Moller, N.; Wear, J. H. Geochim. Cosmochim. Acta 1992, 56, 3839. (48) Brodholt, J. P.; Wood, B. J. Am. Mineral. 1993, 78, 558. (49) Destrigneville, C. M.; Brodholt, J. P.; Wood, B. J. Chem. Geol. 1996, 133, 53. (50) Kuznetsova, T.; Kvamme, B. Phys. Chem. Chem. Phys. 2002, 4, 937. (51) Duan, Z. H.; Zhang, Z. G. Geochim. Cosmochim. Acta 2006, 70, 2311. (52) Tafazzoli, M.; Khanlarkhani, A. Fluid Phase Equilib. 2008, 267, 181. (53) Ma, Y.; Garofalini, S. H. J. Chem. Phys. 2006, 124, 234102. (54) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 1040. (55) Rode, B. M.; Hofer, T. S.; Randolf, B. R.; Schwenk, C. F.; Xenides, D.; Vchirawongkwin, V. Theor. Chim. Acta 2006, 115, 77. (56) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (57) Leung, K.; Nielsen, I. M. B.; Kurtz, I. J. Phys. Chem. B 2007, 111, 4453. (58) Rustad, J. R.; Nelmes, S. L.; Jackson, V. E.; Dixon, D. A. J. Phys. Chem. A 2008, 112, 542. (59) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562. Micheletti, C.; Laio, A.; Parrinello, M. Phys. ReV. Lett. 2004, 92, 170601. Bussi, G.; Laio, A.; Parrinello, M. Phys. ReV. Lett. 2006, 96, 090601. Ensing, B.; Laio, A.; Parrinello, M.; Klein, M. L. J. Phys. Chem. B 2004, 109, 6676.
802 J. Phys. Chem. B, Vol. 113, No. 3, 2009 (60) Marx, D.; Hutter, J. In Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; NIC: FZ Ju¨lich, 2000; pp 301-449. Remler, D. K.; Madden, P. A. Mol. Phys. 1990, 70, 921. (61) CPMD, Copyright IBM Corp 1990-2006. Copyright MPI fu¨r Festko¨rperforschung Stuttgart 1997-2001. (62) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction models for water in relation to protein hydration. In Intermolecular Forces; Pullman, B., Ed.; Riedel: Dordrecht, 1981; p 331. (63) Kitson, D. H.; Hagler, A. T. Biochemistry 1998, 27, 5246. (64) Becke, A. D. Phys. ReV. A 1988, 38, 3098. Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (65) Sprik, M.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1996, 105, 1142\ (66) Vanderbilt, D. Phys. ReV. B 1990, 41, 7892. (67) Laasonen, K.; Pasquarello, A.; Car, R.; Lee, C.; Vanderbilt, D. Phys. ReV. B 1993, 47, 10142.
Kumar et al. (68) Mantz, Y. A.; Chen, B.; Martyna, G. J. Chem. Phys. Lett. 2005, 405, 294. (69) Kuo, I.-F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I.; VandeVondele, J.; Sprik, M.; Hutter, J.; Chen, B.; Klein, M. L.; Mohamed, F.; Krack, M.; Parrinello, M. J. Phys. Chem. B 2004, 108, 12990. (70) Laage, D.; Hynes, J. T. Science 2006, 311, 832. (71) Davis, A. R.; Oliver, B. G. J. Solution Chem. 1972, 1, 329. (72) Rudolph, W. W.; Fischer, D.; Irmer, G. Appl. Spectrosc. 2006, 60, 130. (73) Hage, W.; Hallbrucker, A.; Mayer, E. J. Chem. Soc., Faraday Trans. 1996, 92, 3183. (74) Hage, W.; Hallbrucker, A.; Mayer, E. J. Chem. Soc., Faraday Trans. 1996, 92, 3197. (75) Winkel, K.; Hage, W.; Loerting, T.; Price, S. L.; Mayer, E. J. Am. Chem. Soc. 2007, 129, 13863.
JP809069G
