Structure and Dynamics of Concentrated Hydrochloric Acid Solutions

Jul 1, 2010 - and NaCl) into the acid solution increases the stability of the hydrated proton pairs ... hydrochloric acid (HCl) in water has been expl...
1 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 9555–9562

9555

Structure and Dynamics of Concentrated Hydrochloric Acid Solutions Jianqing Xu,†,‡ Sergei Izvekov,‡ and Gregory A. Voth*,†,‡ Department of Chemistry, James Franck Institute, and Computation Institute, UniVersity of Chicago, 5735 S. Ellis AVe, Chicago, Illinois 60637; and Center for Biophysical Modeling and Simulation and Department of Chemistry, UniVersity of Utah, Salt Lake City, Utah 84112-0850 ReceiVed: March 19, 2010; ReVised Manuscript ReceiVed: June 9, 2010

Self-consistent iterative multistate empirical valence bond (SCI-MS-EVB) simulations for four different concentrated hydrochloric acid (HCl) concentrations (0.43, 0.85, 1.68, and 3.26 M) were conducted to study the structure and dynamic properties of the aqueous multiple excess proton systems. The present, more extensive simulations, confirm that the hydrated excess protons have a tendency to form metastable contact ion pairs by positioning the hydronium oxygen lone pair sides toward one another, in agreement with earlier work [Wang, F.; Izvekov, S.; Voth, G. A. J. Am. Chem. Soc. 2008, 130, 3120]. This unusual behavior results from the special “amphiphilic” nature of the hydrated excess proton. Increasing the acid concentration to 3.26 M was found to have minor effects on the solvation structures of the hydrated proton and chloride counterion, while the average lifetime of the hydrated proton contact ion pairs was reduced. The addition of salts (KCl and NaCl) into the acid solution increases the stability of the hydrated proton pairs and reduces the proton diffusion. 1. Introduction For over two centuries, significant scientific effort has been devoted to the investigation of proton solvation and transport in aqueous media due to its essential role in many areas of chemistry, biology, and materials science.1-10 As one of the most fundamental but important issues, the solvation structures of hydrochloric acid (HCl) in water has been explored both experimentally and theoretically.11-19 However, a definitive molecular level description, especially at intermediate HCl concentrations, has not yet been provided. The unusual amphiphilic property of the hydrated excess proton20 also increases the complexity of this issue, since the interaction between hydrated excess protons is predicted to be more complex than standard cation-cation electrostatic interactions.11 The amphiphilic nature of the hydrated excess proton is due to the directionally asymmetric ability of the core hydronium entity to form hydrogen bonds with neighboring water. On one side, the three hydrogen atoms of the hydronium can form very strong hydrogen bonds with the surrounding water molecules to make an Eigen cation, H9O4+, thus having hydrophilic character. However, the net charge of the hydronium is positive, and it is energetically unfavorable for the hydronium oxygen (OHy) to be a proton acceptor. Therefore, the lone pair region on the opposite side does not form a hydrogen bond with water, thus exhibiting “hydrophobic” character.20,21 This resulting “amphiphilic” nature of the hydrated proton also disrupts the hydrogen-bonding network in the bulk water. As a result, it was predicted that the hydronium concentration near the water-vapor interface can be enhanced in acidic aqueous solution, because the hydrophobic lone pair at the interface can be oriented toward the vapor region to minimize its disturbance in the bulk water.20 This prediction has subsequently received experimental support22-34 and attracted significant theoretical attention.11,21,35-44 † ‡

University of Chicago. University of Utah.

In order to investigate the solvation structure of multiple hydrated protons, Wang et al.11 carried out a computer simulation study employing the self-consistent iterative multistate empirical valence bond (SCI-MS-EVB)45 method for the HCl in bulk water solution with the HCl concentration between 0.43 and 3.26 M.11 It was observed that hydrated protons can form contact ion pairs in bulk water, with their hydrophobic OHy atoms pointing toward each other.11 It is reasoned that the hydrated protons have a tendency to transiently form ion pairs in order to minimize their hydrophobic influence in bulk water, but this result is surprising considering the very strong electrostatic repulsive force between the positively charged cations. However, the centers of excess charges on the hydrated proton Eigen cations (H9O4+) are each shifted ∼0.5 Å away from the OHy atom toward the hydrogen atoms, thereby increasing the effective charge separation while maintaining the molecular structure of the ion pair. This shift is due to the net positive charge of the electronic defect associated with the hydrated excess proton being delocalized over multiple water molecules, instead of residing on one particular hydronium core. This shifted distance between the charge centers (0.5 + 0.5 ) 1 Å) has the effect of reducing the electrostatic repulsive force between hydroniums, which increases the likelihood of forming transient hydrated proton contact ion pairs. These results were also independently supported by the results of a Car-Parrinello molecular dynamics (CPMD) simulation.11 The underlying MS-EVB model that was employed by Wang et al. was at an intermediate stage of development. It was a combination of MS-EVB246 and MS-EVB347 models. The recently refined MS-EVB3 model employed in the present work has been improved to yield a more accurate solvation structure of the hydrated excess proton in the aqueous solution. In particular, the solvent depletion on the lone pair side of hydronium oxygen, which was overestimated in the MS-EVB2 model, is more precisely described in the new model. The purpose of the current study is therefore to simulate concentrated HCl solutions with the MS-EVB3 model, and to provide

10.1021/jp102516h  2010 American Chemical Society Published on Web 07/01/2010

9556

J. Phys. Chem. B, Vol. 114, No. 29, 2010

Xu et al.

additional analysis of the structure and dynamics of the hydrated proton contact ion pairs, as well as the solvation structure and diffusion properties of the concentrated acid systems. Additionally, due to the ability of salts to modify the hydrophobic effect of nonpolar molecules in water (e.g., the “salting out” effect), KCl and NaCl doped HCl acid solutions are also simulated. 2. Methodology 2.1. SCI-MS-EVB Simulation. The electronic charge defect from the hydrated excess proton is delocalized over several water molecules, which is called a “complex” in the MS-EVB model.8,10 The core structure of the complex exists in the limited forms of the Eigen (H9O4+)48 and Zundel (H5O2+)49 cations. Proton hopping often occurs by means of Grotthuss shuttling mechanism,1,3,6-8,10 which is a transition (through the water hydrogen-bonding network) from one Eigen complex, via the Zundel complex structure, to another Eigen complex centered on a neighboring oxygen atom, accompanied by chemical bonding rearrangement.50 The MS-EVB method8,10,45-47,51-54 treats the hydrated excess proton by a linear combination of “basis states”, corresponding to possible chemical bonding topologies around the excess proton. A quantum-like Hamiltonian matrix is then built with the diagonal terms representing the bonded and nonbonded interactions of these states with the off-diagonal terms describing the couplings between states. At each time step, the Hamiltonian matrix is diagonalized and the lowest eigenvalue is determined. The square of the corresponding eigenvector provides the weight of each EVB state, and the forces on each nuclei are calculated by the Hellmann-Feynman theorem. The self-consistent iterative (SCI-MS-EVB) method is an important generalization in order to include multiple hydrated proton in a simulation.45 Due to the delocalized nature of the charge defect associated with the excess proton, it is spread out over several water molecules in the MS-EVB complex. Therefore, to characterize the position of the excess charge, the concept of the center of excess charge (CEC) of the hydrated proton charge defect has been introduced, given by8,10,46,47 N

rCEC(t) )

i (t) ∑ ci2(t)rCOC

(1)

i

i where rCOC (t) is the instantaneous position of the center-ofcharge of the hydronium in the ith EVB state, and ci2(t) is the probability of each EVB state. In the subsequent discussion, the term “hydronium” refers to a single protonated water molecule in a particular EVB state, while the term “hydroniumlike state” refers to the EVB state in the MS-EVB representation time t having the largest EVB amplitude. 2.2. Simulation Details. All the simulated systems in this work follow the procedure described in the paper of Wang et al.11 Specifically, four systems including 2, 4, 8, and 16 HCl ion pairs in 256 H2O molecules (0.43, 0.85, 1.68, and 3.26 M, respectively) were studied. The simulations of each concentration consisted of 20 independent SCI-MS-EVB trajectories. Each trajectory was equilibrated at 298 K in the canonical ensemble (constant NVT) for 500 ps using a Nose´-Hoover thermostat with a relaxation time 0.5 ps. A 0.5 fs MD time step and the Ewald summation method were used in all the simulations. Subsequently, each trajectory was sampled in the microcanonical ensemble (constant NVE) for 2 ns (40 ns in total for 20 trajectories), 1 ns (20 ns), 500 ps (10 ns), and 250 ps (5 ns) for the 0.43, 0.85, 1.68, and 3.26 M systems, respectively. To study

the effects of added salts to these solutions, half of the HCl molecules in the 1.68 M (8 HCl) and 3.26 M (16 HCl) solutions were replaced by an equivalent amount of KCl and NaCl. The simulation lengths of each acid-salt system were the same as the acid only trajectories. The MS-EVB potential parameters were chosen from the MSEVB3 model.47 The salt parameters were taken from Dang’s work,55 which avoids an anomalous salt crystallization effect.56 The error bars in the figures reported later are the standard errors calculated from 20 independent data points of each HCl molarity. 2.3. State-Averaged Radial Distribution Function. In the MS-EVB simulations, the assignment of hydronium oxygen (OHy) depends on the EVB states. Given the distribution of EVB states, the simple radial distribution function (RDF) calculation in terms of hydronium-like OHy atoms is not obvious. In order to capture the delocalized nature of the hydrated excess proton charge defect, the “state-averaged” RDF was introduced.11 It can be described by using the single proton case, in which each EVB state represents a certain probability, ci2, that the particular water is the hydronium. Therefore, the RDF between the OHy and water oxygen (OW) is a linear combination of OHy-OW RDFs from each possible EVB state weighted by its coefficient: N

g′OHy-OW(r) )

i ∑ ci2gOH -OW(r) i)1

y

(2)

where N is the number of states in the EVB complex and giOHy-OW(r) is the OHy-OW RDF from a state |i〉 with probability ci2. This definition can also be extended to the multiple excess proton case. An example of the two proton summation is: N1

g′OHy-OHy(r) )

N2

ij ∑ ∑ ci2cj2gOH -OH (r) y

i)1 j)1

y

(3)

where the N1 and N2 are the number of states in each complex. For a system having n protons, a total of n(n - 1)/2 pairs in eq 3 must be included. 2.4. Lifetime of Hydrated Proton Contact Ion Pair. In order to quantify how long the hydrated proton contact ion pairs can persist in the simulations, the chemical equation

(4)

RhP

can be utilized to calculate the hydrated proton pair lifetime. The reactant R denotes the formation of hydrated proton contact ion pair, and the product P represents its dissociation products. Assuming a first-order kinetic rate equation, the macroscopic concentration difference of the ion pair decays exponentially57

∆cR(t) t ) exp ∆cR(0) τ

( )

(5)

where ∆cR(t) ) cR(t) - 〈cR〉 is the time-dependent difference between the instantaneous concentration and the equilibrium concentration of the reactant R. The inverse relaxation time τ-1 is the sum of the forward and back rate constants

Structure of Concentrated Hydrochloric Acid Solutions

τ-1 ) kRfP + kPfR

J. Phys. Chem. B, Vol. 114, No. 29, 2010 9557

(6)

A step function h(t) is defined in order to distinguish the R and P of a particular hydrated proton contact ion pair: at time t, if the distance between the OHy atoms of two hydroniumlike states is smaller than a critical distance r*, the system configuration is R and h(t) ) 1; otherwise it is P and h(t) ) 0. The critical distance r* was chosen based on the first minimum of the OHy-OHy RDF, which is around 4.1 Å (see details below), such that

hR(t) )

{

1 0

r(t) e r* r(t) > r*

(7)

Naturally, this definition can be extended to the state-averaged definition as well:

hRij (t) )

{

ci2(t)cj2(t) 0

rij(t) e r* rij(t) > r*

(8)

where ci2(t) means the probability of a tagged state |i〉 in complex A at time t, and cj2(t) is the probability of state |j〉 in complex B. The configuration of a hydrated proton pair can be described by N1

h′R(t) )

N2

∑ ∑ hRij (t)

(9)

i)1 j)1

The time autocorrelation function of h′R(t) can be defined in the usual way

CR(t) )

〈(h′R(t) - 〈h′R〉)(h′R(0) - 〈h′R〉)〉 〈(h′R(0) - 〈h′R〉)2〉

(10)

According to the fluctuation-dissipation theorem or Onsager’s regression hypothesis57

∆cR(t) t ) C(t) ) exp ∆cR(0) τ

( )

(11)

at sufficiently long times. By building the autocorrelation function CR(t) based on the h′R(t) function from the MS-EVB trajectories and fitting the linear part of the ln(CR(t)), the lifetime τ of the hydrated proton contact ion pair can be estimated. 3. Results and Discussion 3.1. Hydrated Proton Contact Ion Pair. A convenient way to quantify the hydrated proton contact ion pair is to calculate the radial distribution functions (RDF) as a function of the OHy-OHy separation. Figure 1 shows various RDF results (any RDFs related to OHy atom refer to the state-averaged RDF described in the methodology part) for the two excess proton (0.43 M) system. As emphasized in the Introduction, the lone pair side of the OHy has been proposed to have a hydrophobic character. As a result, compared to the water-water (OW-OW) RDF, a deeper well in OHy-OW RDF after the first peak is observed.47

Figure 1. State-averaged RDF results for a 0.43 M acid (2 HCl + 256 H2O) system. OHy, OW, CEC, and Cl- represent hydronium-like oxygen, water oxygen, center of excess charge, and chloride anion, respectively.

The hydronium-hydronium (OHy-OHy) RDF in Figure 1 has a small peak around 3.2 Å. At this separation, not even a single water molecule can be accommodated between the two predominantly Eigen cation-like structures, which indicates that two hydrated protons are in a close proximity. The formation of the contact ion pair is facilitated due to the subtle balance between the OHy-OHy amphiphilic effect20 and the reduced cation-cation electrostatic repulsive interaction due to delocalization of the two hydrated proton charge defect CECs.11 Indeed, the CEC-CEC RDF around 4.2 Å suggests a preferred distance of an additional angstrom apart between the charge defect centers of the two proton complexes when a contact ion pair forms compared to the OHy-OHy distance. The typical structure of the ion pair can be seen in Figure 2, which shows a snapshot of two hydronium-like structure (blue) with their oxygen (OHy) atoms pointing toward each other in order to minimize the unfavorable interaction of their lone pair side with water, while the CEC (yellow balls) are further apart from each other in order to minimize the repulsive electrostatic interaction. This latter effect is a “nonclassical” feature of the hydrated proton ion pairs due to their ability to delocalize the net positive charge defects via the “special-pair dance” dynamical motions.50 Unlike Wang et al.’s results,11 however, the hydrated proton contact ion pair observed in the current study is not the absolute maximum in the OHy-OHy RDF for the two hydrated excess proton system. The height of the first peak in the present OHy-OHy RDF at 3.2 Å is slightly smaller than unity, indicating a more metastable structure than seen in ref 11. This result is presumably due to the somewhat reduced amphiphilic behavior of the final MS-EVB3 hydrated proton model.21,47 However, compared to the result of a classical H3O+ simulation that allows no proton CEC delocalization, the OHy-OHy RDF in the SCIMS-EVB model shows a significantly greater probability for the hydrated proton contact ion pair to form at small distances (Figure 3). As stated earlier, this is because the SCI-MS-EVB model can describe the delocalized nature of the excess proton charge defects, thus maximizing the separation of charge centers (CECs) while also maximizing the pairing of the hydronium lone pair sides to minimize their disruption of the water hydrogen-bonding network. 3.2. Ion Pair Structures at Different HCl Concentrations. Since the transient hydrated proton contact ion pair can be formed in the 2 excess proton (0.43 M) case, it is important to explore whether higher HCl concentrations will produce more or less stable pairs, and how the proton solvation structures will respond to the change of the acid concentration. Figure 4 displays the OHy-OHy, CEC-CEC, Cl--Cl-, and OHy-Cl- RDFs for the four studied acid concentrations (0.43,

9558

J. Phys. Chem. B, Vol. 114, No. 29, 2010

Xu et al.

Figure 2. A snapshot showing a typical hydrated proton contact ion pair configuration. Blue-colored balls are the hydronium. The yellow transparent balls are the CEC and the cyan balls are the Cl-.

Figure 3. A comparison between the OHy-OHy RDF from a classical hydronium model (or 1 state model), and the OHy-OHy RDF (stateaveraged) as well as the CEC-CEC RDF from the SCI-MS-EVB model. Both systems have the acid concentration 0.43 M.

0.85, 1.68, and 3.26 M). When the concentration of excess protons increases, the water structure is increasingly disrupted. On one hand, one might predict that more stable hydrated proton

pairs could form at higher HCl concentrations. On the other hand, the increased charge repulsion would suggest otherwise. The result is that adding more HCl into the solution has little effect on the hydrated proton contact ion pairs. Instead of seeing enhanced pairing, the introduction of more excess protons has a weak tendency to decrease the OHy-OHy pair correlation (Figure 4a), presumably due to the larger net charge repulsion. In addition, the increase in HCl molarity leads to an increased interaction between hydrated protons and chloride anions. This also appears to lower the probability of hydrated proton pairing, while increasing the probability of chloride anion pairing. The peaks around 5.2, 7.6, and 9.6 Å in the Cl--Cl- RDFs correspond to the chloride anions being separated by one, two, and three water molecules, respectively.19 The correlation between the hydrated proton structures and the chloride anions is also strong. 3.3. Lifetimes of Hydrated Proton Contact Ion Pairs. Although the metastable structures of the hydrated proton contact ion pairs at different concentrations were reported earlier via by the state-averaged RDFs, their dynamical stability is also

Structure of Concentrated Hydrochloric Acid Solutions

Figure 4. RDFs of (a) OHy-OHy, (b) CEC-CEC, (c) Cl--Cl-, and (d) OHy-Cl- for the four different HCl concentrations: 0.43, 0.85, 1.68, and 3.26 M. The dotted lines shown in panel (d) represent the coordination number (right Y axis) of chloride anions around each hydrated proton. The large dots on the coordination number curves denote the values at the first coordination shell.

J. Phys. Chem. B, Vol. 114, No. 29, 2010 9559

Figure 6. (a) Plot of hR(t) function verses time for one of the 0.43 M trajectories from SCI-MS-EVB simulation. (b) State-averaged treatment of the hR(t) function.

Figure 7. Semilog scale autocorrelation function ln(CR(t)) based on the function h′R(t) for the four acid concentrations studied. The linear fit at long time scale for each concentration gives an estimate of the lifetime τ. Figure 5. Ion pair “reactant” step function hR(t) over the time of the CPMD trajectory. hR ) 1 represents that two OHy atoms are within the critical distance r* (R configuration), while hR ) 0 indicates that the distance between two OHy atoms are larger than the critical distance r* (P configuration).

of considerable interest. To address this question, the lifetime of the hydrated proton contact ion pairs was calculated via the methodology described in eqs 4-11. The CPMD trajectory of 1 M HCl solution from Wang et al.’s work11 (2 HCl and 110 waters; the technical details of the CPMD simulations are found in ref 11) was also utilized to provide an independent estimate of ion pair lifetime. As described in eq 7, when the distance between two oxygen atoms of hydroniums is less than the critical distance r*, then hR(t) ) 1, which indicates the formation of a contact ion pair; otherwise hR(t) ) 0. The location of the first minimum in the OHy-OHy RDFs in Figure 4a is around 4.1 Å. This value was therefore chosen to be the critical distance r* to calculate the hR(t) function. Figure 5, for example, shows the hR(t) plot describing the dynamics configurations of the contact ion pair in the CPMD simulation. The black vertical lines represent the frequent exchange between configurations R and P (reactant and product). Due to the high computational cost, however, the CPMD simulation consisted of only a single 72 ps trajectory. Therefore, the statistics were not sufficient to build the correlation function given in eq 10. However, from the available data the lifetime of the contact pair in the CPMD simulation was estimated to be 14 ps for this HCl concentration (1 M). A more converged and independent lifetime estimate was obtained from fitting the h′R(t) autocorrelation function in the

MS-EVB simulations from eq 10. The step function hR(t) and its state-averaged form h′R(t) are shown in Figure 6, utilizing one trajectory from the 0.43 M system as an example. The stateaveraged treatment (see methodology) of h′R(t) in Figure 6b is more correct to construct the autocorrelation function and extract the ion pair lifetime, because the proton charge defect delocalizes over several surrounding water molecules and each pure hydronium state only has a certain MS-EVB amplitude. Figure 7 illustrates the decay behavior of the log of the ion pair autocorrelation functions ln(CR(t)) at the four different HCl concentrations. Each concentration has several apparent exponential decay time scales, indicating multiple relaxation processes of the hydrated proton contact ion pairs. In the curves for the 0.43 and 0.85 M systems, it can be seen that there are at least two linear decay regions after one initial nonlinear region. The nonlinear part is intrinsically difficult to characterize,53 but the first linear region (not highlighted by a line) is likely related to the intrinsic proton hopping time, which can change the local R/P structures of the hydrated proton contact ion pairs on a short time scale. The second linear region should be most closely related to the correlation time (lifetime) describing how long the ion pair R can persist on average. The lifetimes by fitting the second linear areas of the curves of the 0.43 and 0.85 M systems are around 13.9 and 13.6 ps, which are consistent with the estimate from the CPMD simulation in this HCl concentration range. After increasing the molarity of HCl to 1.68 and 3.26 M, a fitting of the longer time scale linear regions of these two systems gives lifetimes of 8.5 and 5.2 ps, respectively. These

9560

J. Phys. Chem. B, Vol. 114, No. 29, 2010

Figure 8. Results for simulations in which half of the HCl molecules of 1.68 M (8 HCl) and 3.26 M (16 HCl) are replaced by the KCl and the NaCl. The OHy-OHy RDFs and log scale autocorrelation function ln(CR(t)) are evaluated and compared with the acid only solutions. All the results indicate an enhanced correlation of hydrated proton pairs after the introduction of salts.

results show that the average lifetime of the hydrated proton contact ion pair is decreased in more concentrated solutions. The reason for this decrease likely comes from the more crowded ionic environment at high acid concentration, which destabilizes the ion pairs through increased electrostatic interactions. Another reason might be due to the presence of more chloride counterions, because they have a tendency to penetrate into the first solvation shell of the hydrated protons. 3.4. Addition of Extra Salts. It is known that the solubility of a nonelectrolyte in an aqueous solution can be reduced by the addition of the electrolyte, called the “salting out” effect. For example, the solubility of hydrophobic solutes like methane58,59 and amphiphilic molecules like t-butanol60,61 can be decreased further by the addition of sodium halides to the solution. The solubility of proteins can also be changed by different salts, which is ordered in the Hofmeister series.62 A definitive understanding of this behavior at the molecular level is still elusive.63 However, rationalizing the salting-out mechanism is not the objective of current study. Instead, this section explores the effects of added salt on the ion pairing interactions between hydronium cations. Based on this motivation, half of the HCl molecules of the 1.68 M (8 HCl) and the 3.26 M (16 HCl) systems were replaced by NaCl and KCl. The hydrated proton pair stability was then calculated and compared with a pure acid solution at the same acid concentration. After addition of NaCl and KCl into the acid solutions, the stability of hydrated proton pairs in the 1.68 M HCl case (Figure 8a) seems to have a little effect. However, a larger stability enhancement of hydrated proton contact ion pair was observed in the 3.26 M HCl case (Figure 8b). Figure 8c,d shows the associated hydrated proton pair correlation functions (log scale of the autocorrelation function) corresponding to all of the cases studied. Both figures provide evidence that the lifetime of hydrated proton contact ion pair is increased if extra salts are present in the solution. This effect may be due to a change in the local water hydrogen bonding structure and dynamics because of the added Na+ or Ka+ cations, or more likely from changed electrostatic features. 3.5. Diffusion Properties. The difference in the solvation structures and hydrated proton ion pair dynamics at different concentrations will result in a change in proton diffusion. The

Xu et al.

Figure 9. Diffusion coefficients of the hydrated proton CEC, Cl- anion, and the hydronium-like oxygen at four HCl concentrations. The left inset shows the scaled simulation results (after the application of a constant correction factor)47 compared with experimental results.13 The right inset shows the separation of the CEC diffusion coefficient into its three contributions: hopping, vehicular, and their correlation. The blue and the green points are the CEC diffusion coefficients after the addition of the KCl and the NaCl into the acid solution. They are also separated into three components to investigate the origin of proton diffusion change upon addition of salt.

diffusion coefficient of the hydrated proton CEC can be obtained by fitting the linear part of the CEC mean square displacement (MSD) at long times:

DCEC ) lim tf∞

〈|rCEC(t) - rCEC(0)| 2〉 6t

(12)

Figure 9 depicts the diffusion coefficients of the hydrated proton CEC and the Cl- counteranions at the four concentrations. The diffusions of CEC and Cl- both decrease as the HCl concentration increases. The present SCI-MS-EVB simulation does not include nuclear quantum effects,47,52 so the diffusion coefficient of the proton is lower than the experimental values.13 However, when the acid concentration increases, the reduction in excess proton diffusion is consistent with the experimental results13 after the application of a constant correction factor (left inset in Figure 9).47 In order to better understand the reduction in the proton diffusion as a function of concentration, one can separate the total MSD of the proton into the contribution due to the proton hopping (Grotthuss mechanism), and the standard (vehicular) diffusion of hydronium-like cation, as well as the correlation between them.64 Although the CEC can denote the position of the excess charge, it is a continuous and smooth trajectory with time, and thus it cannot be decomposed. The most hydroniumlike entity is the hydronium structure in the MS-EVB state having the largest amplitude. For a time it follows a continuous vehicular diffusion like a regular ion. If a proton hops, the identity of the hydronium-like oxygen transfers from one water molecule to another, which is a discrete hopping process. At long time scales, the total MSD of the hydronium-like oxygen is very close to the MSD of the CEC,64 as shown in Figure 9. Therefore, it can be utilized to understand how the hopping and vehicular motions contribute to the overall hydrated proton diffusion. The displacement vector of the hydronium-like oxygen position can be written as a summation of the hopping and vehicular displacement:

Structure of Concentrated Hydrochloric Acid Solutions

∆b r ) b(t) r - b(0) r ) ∆b r hopping + ∆b r vehicular

J. Phys. Chem. B, Vol. 114, No. 29, 2010 9561

(13)

The MSD can be written as

〈∆b r · ∆b〉 r ) 〈∆b r hopping · ∆b r hopping〉 + 〈∆b r vehicular · ∆b r vehicular〉 + 2〈∆b r hopping · ∆b r vehicular〉 (14) Therefore, the diffusion coefficient of the hydronium-like oxygen can be decomposed into three components by fitting the linear part of each MSD

Dtotal ) Dhopping + Dvehicular + Dcorrelation

(15)

The right inset of Figure 9 depicts the decomposed diffusion coefficient of the hydronium-like oxygen for the various systems. The Dhopping drops most extensively even for the two low concentrations 0.43 and 0.85 M, while the vehicular parts are almost identical at these two molarities, with a slight change in the correlation. This result indicates that the decrease of proton diffusion at relatively low concentrations is mainly due to less proton hopping. When raising the HCl concentration (1.68 and 3.26 M), all the components go down, even for the regular vehicular diffusion. At the lowest concentrations, the relative pure water environment allows the proton to transport more easily. By contrast, the higher concentration results in a denser ionic environment, in which the effective solvation shell of each hydrated proton shrinks because of less available water, and the hopping events become hard to take place due to the disrupted solvation structures, since a new solvation shell needs to be formed for a successful proton hopping event.3,50,65 This result points to two possible conclusions: (1) the proton reduces its hopping rate, or (2) the hopping rate does not change much, but the proton hops back and forth more around the same average location. By counting the hopping rate at each concentration, the averaged number of hopping events for each proton per picosecond is 8.10 ( 0.03, 7.97 ( 0.03, 7.91 ( 0.03, and 7.82 ( 0.06 for 0.43, 0.85, 1.68, and 3.26 M HCl, respectively. These results do exhibit a small tendency for the hopping rate to decrease. However, compared to the larger decrease of the diffusion constant, the difference seems very small among the investigated four concentrations. This in turn implies that the hydrated protons are more likely to hop back and forth around the same location at higher HCl concentrations, which decreases the effective hopping diffusion at long time scales. The vehicular component of the hydrated proton diffusion is not as sensitive as the hopping component to the change in HCl concentration. It mostly shows a tendency of diffusion suppression for the 1.68 and 3.26 M HCl systems. When raising the HCl molarity, the ions can maintain the same vehicular diffusion for the lowest concentrations, but this will eventually decrease because of the relatively more crowded ionic environment at the higher concentrations. All of the studied concentrations also show a negative (anti) correlation between the hopping and vehicular diffusion, meaning that statistically the two diffusions can move in an opposite direction.64 This effect might be due to the influence from the chloride anions, since they are heavy and have a slower motion, but can significantly attract the positively charged hydrated protons. This behavior has been seen previously in the proton exchange membrane Nafion due to the sulfonate anion groups on the polymer side chains.64

The proton diffusion changes affected by NaCl and KCl were also evaluated. It is not surprising that the existence of extra salts can decrease the proton diffusion, since these additives can disrupt more water hydrogen bonds in the solution.66,67 In Figure 9, when comparing the 4 HCl + 4 KCl and 4 HCl + 4 NaCl systems with the 4 HCl only system, or comparing the 8 HCl + 8 KCl and 8 HCl + 8 NaCl systems with 8 HCl only system, the introduction of NaCl and KCl was indeed found to significantly diminish the hydrated proton diffusion. The right inset of Figure 9 displays the three components of the CEC diffusion coefficient when the extra salts are present. As expected, the hopping component is reduced the most significantly. 4. Conclusions In summary, various concentrated hydrochloric acid aqueous solutions (0.43, 0.85, 1.68, and 3.26 M) have been studied by the SCI-MS-EVB computer simulation methodology. The MSEVB3 model predicts a metastable formation of hydrated proton contact ion pairs, which is due to the competition between cation-cation electrostatic repulsion and OHy-OHy “amphiphilic” attraction. Increasing the acid concentration decreases the average lifetime of the hydrated proton contact ion pairs from 13.9 ps (0.43 M) to 5.2 ps (3.26 M). Moreover, the introduction of extra salts enhances the stability of the hydrated proton pairs. After increasing the molarity of HCl, the water hydrogenbonding network is disrupted by the introduction of more ions. The increase in the concentration induces a more crowded ionic environment with more complex interactions, which suppresses the overall hydrated proton diffusion by reducing both the proton hopping and vehicular diffusion. The hydrated proton diffusion is further decreased by the addition of the KCl and NaCl in the solutions due to the further disruption of the water hydrogenbonding network. Acknowledgment. This research was supported by the National Science Foundation (CHE-0719522). The computational resources for this project have been provided in part by a grant of computer time from the DOD High Performance Computing Modernization Program at Maui High Performance Computing Department of Defense Supercomputing Center. The authors thank Drs. C. Mark Maupin, Takefumi Yamashita, Luca Larini, and Brian Hopkins for many helpful discussions. References and Notes (1) de Grotthuss, C. J. T. Ann. Chim. 1806, LVIII, 54–74. (2) Bell, R. P. The proton in chemistry; Cornell University Press: Ithaca, NY, 1959. (3) Agmon, N. Chem. Phys. Lett. 1995, 244, 456–462. (4) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601–604. (5) Kreuer, K. D.; Paddison, S. J.; Spohr, E.; Schuster, M. Chem. ReV. 2004, 104, 4637–4678. (6) Cukierman, S. Biochim. Biophys. Acta 2006, 1757, 876–885. (7) Marx, D. Chemphyschem 2006, 7, 1848–1870. (8) Voth, G. A. Acc. Chem. Res. 2006, 39, 143–150. (9) Wraight, C. A. Biochim. Biophys. Acta 2006, 1757, 886–912. (10) Swanson, J. M. J.; Maupin, C. M.; Chen, H. N.; Petersen, M. K.; Xu, J. C.; Wu, Y. J.; Voth, G. A. J. Phys. Chem. B 2007, 111, 4300–4314. (11) Wang, F.; Izvekov, S.; Voth, G. A. J. Am. Chem. Soc. 2008, 130, 3120–3126. (12) Triolo, R.; Narten, A. H. J. Chem. Phys. 1975, 63, 3624–3631. (13) Dippel, T.; Kreuer, K. D. Solid State Ionics 1991, 46, 3–9. (14) Agmon, N. J. Phys. Chem. A 1998, 102, 192–199. (15) Cukierman, S. Biophys. J. 2000, 78, 1825–1834. (16) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M. A.; Soper, A. K. J. Chem. Phys. 2004, 121, 7840–7848. (17) Harsa´nyi, I.; Pusztai, L. J. Phys.: Condens. Matt. 2005, 17, S59– S65.

9562

J. Phys. Chem. B, Vol. 114, No. 29, 2010

(18) Botti, A.; Bruni, F.; Ricci, M. A.; Soper, A. K. J. Chem. Phys. 2006, 125, 014508-014509. (19) Heuft, J. M.; Meijer, E. J. Phys. Chem. Chem. Phys. 2006, 8, 3116– 3123. (20) Petersen, M. K.; Iyengar, S. S.; Day, T. J. F.; Voth, G. A. J. Phys. Chem. B 2004, 108, 14804–14806. (21) Iuchi, S.; Chen, H.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2009, 113, 4017–4030. (22) Petersen, P. B.; Saykally, R. J. J. Phys. Chem. B 2005, 109, 7976– 7980. (23) Petersen, P. B.; Saykally, R. J. Annu. ReV. Phys. Chem. 2006, 57, 333–364. (24) Petersen, P. B.; Saykally, R. J. Chem. Phys. Lett. 2008, 458, 255– 261. (25) Pegram, L. M.; Record, M. T. J. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 14278–14281. (26) Pegram, L. M.; Record, M. T. J. Chem. Phys. Lett. 2008, 467, 1– 8. (27) Radu¨ge, C.; Pflumio, V.; Shen, Y. R. Chem. Phys. Lett. 1997, 274, 140–144. (28) Shen, Y. R.; Ostroverkhov, V. Chem. ReV. 2006, 106, 1140–1154. (29) Tian, C.; Ji, N.; Waychunas, G. A.; Shen, Y. R. J. Am. Chem. Soc. 2008, 130, 13033–13039. (30) Baldelli, S.; Schnitzer, C.; Shultz, M. J. Chem. Phys. Lett. 1999, 302, 157–163. (31) Gopalakrishnan, S.; Liu, D.; Allen, H. C.; Kuo, M.; Shultz, M. J. Chem. ReV. 2006, 106, 1155–1175. (32) Buch, V.; Tarbuck, T.; Richmond, G. L.; Groenzin, H.; Li, I.; Shultz, M. J. J. Chem. Phys. 2007, 127, 204710–204715. (33) Tarbuck, T. L.; Ota, S. T.; Richmond, G. L. J. Am. Chem. Soc. 2006, 128, 14519–14527. (34) Levering, L. M.; Sierra-Hernandez, M. R.; Allen, H. C. J. Phys. Chem. C 2007, 111, 8814–8826. (35) Iyengar, S. S.; Day, T. J. F.; Voth, G. A. Int. J. Mass Spectrom. 2005, 241, 197–204. (36) Iyengar, S. S.; Petersen, M. K.; Day, T. J. F.; Burnham, C. J.; Teige, V. E.; Voth, G. A. J. Chem. Phys. 2005, 123, 84309. (37) Mucha, M.; Frigato, T.; Levering, L. M.; Allen, H. C.; Tobias, D. J.; Dang, L. X.; Jungwirth, P. J. Phys. Chem. B 2005, 109, 7617–7623. (38) Burnham, C. J.; Petersen, M. K.; Day, T. J. F.; Iyengar, S. S.; Voth, G. A. J. Chem. Phys. 2006, 124, 24327. (39) Petersen, M. K.; Voth, G. A. J. Phys. Chem. B 2006, 110, 7085– 7089. (40) Ishiyama, T.; Morita, A. J. Phys. Chem. A 2007, 111, 9277–9285. (41) Kofinger, J.; Dellago, C. J. Phys. Chem. B 2008, 112, 2349–2356. (42) Kuldin, K. N.; Car, R. J. Am. Chem. Soc. 2008, 130, 3915–3919.

Xu et al. (43) Vacha, R.; Horinek, D.; Berkowitz, M. L.; Jungwirth, P. Phys. Chem. Chem. Phys. 2008, 10, 4975–4980. (44) Buch, V.; Milet, A.; Vacha, R.; Jungwirth, P.; Devlin, J. P. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 7342–7347. (45) Wang, F.; Voth, G. A. J. Chem. Phys. 2005, 122, 144105. (46) Day, T. J. F.; Soudackov, A. V.; Cuma, M.; Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 2002, 117, 5839–5849. (47) Wu, Y. J.; Chen, H. N.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 467–482. Additions and Corrections: J. Phys. Chem. B 2008, 112, 7146. (48) Eigen, M. Angew. Chem., Int. Ed. 1964, 3, 1–19. (49) Zundel, G. The Hydrogen Bond-Recent DeVelopments in Theory and Experiments. II Structure and Spectroscopy; North-Holland: Amsterdam, 1976. (50) Markovitch, O.; Chen, H.; Izvekov, S.; Paesani, F.; Voth, G. A.; Agmon, N. J. Phys. Chem. B 2008, 112, 9456–9466. (51) Schmitt, U. W.; Voth, G. A. J. Phys. Chem. B 1998, 102, 5547– 5551. (52) Schmitt, U. W.; Voth, G. A. Isr. J. Chem. 1999, 39, 483–492. (53) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361– 9381. (54) Day, T. J. F.; Schmitt, U. W.; Voth, G. A. J. Am. Chem. Soc. 2000, 122, 12027–12028. (55) Dang, L. X. J. Am. Chem. Soc. 1995, 117, 6954–6960. (56) Joung, I. S.; Cheatham, T. E. J. Phys. Chem. B 2008, 112, 9020– 9041. (57) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (58) Docherty, H.; Galindo, A.; Sanz, E.; Vega, C. J. Phys. Chem. B 2007, 111, 8993–9000. (59) Graziano, G. J. Chem. Phys. 2008, 129, 084506-084509. (60) Paschek, D.; Geiger, A.; Herve, M. J.; Suter, D. J. Chem. Phys. 2006, 124, 154508–154513. (61) Jeufack, H. M.; Suter, D. J. Chem. Phys. 2007, 126, 144501– 144505. (62) Hofmeister, F. Arch. Exp. Pathol. Pharmakol. 1888, 24, 247. (63) Grover, P. K.; Ryall, R. L. Chem. ReV. 2004, 105, 1–10. (64) Petersen, M. K.; Voth, G. A. J. Phys. Chem. B 2006, 110, 18594– 18600. (65) Chandra, A.; Tuckerman, M. E.; Marx, D. Phys. ReV. Lett. 2007, 99, 145901. (66) Agmon, N.; Goldberg, S. Y.; Huppert, D. J. Mol. Liq. 1995, 64, 161–195. (67) Hertz, H. G.; Versmold, H.; Yoon, C. Ber. Bunsenges. Phy. Chem. 1983, 87, 577–582.

JP102516H